source: sasview/src/sas/sasgui/perspectives/pr/pr.py @ 0391dae

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 0391dae was c10d9d6c, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 8 years ago

updated references to old paths.

  • Property mode set to 100644
File size: 47.1 KB
Line 
1"""
2    P(r) perspective for SasView
3"""
4################################################################################
5#This software was developed by the University of Tennessee as part of the
6#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
7#project funded by the US National Science Foundation.
8#
9#See the license text in license.txt
10#
11#copyright 2009, University of Tennessee
12################################################################################
13
14
15# Make sure the option of saving each curve is available
16# Use the I(q) curve as input and compare the output to P(r)
17
18import sys
19import wx
20import logging
21import time
22import math
23import numpy
24import pylab
25from sas.sasgui.guiframe.gui_manager import MDIFrame
26from sas.sasgui.guiframe.dataFitting import Data1D
27from sas.sasgui.guiframe.events import NewPlotEvent
28from sas.sasgui.guiframe.events import StatusEvent
29from sas.sasgui.guiframe.gui_style import GUIFRAME_ID
30from sas.sascalc.pr.invertor import Invertor
31from sas.sascalc.dataloader.loader import Loader
32import sas.sascalc.dataloader
33
34from pr_widgets import load_error
35from sas.sasgui.guiframe.plugin_base import PluginBase
36
37
38PR_FIT_LABEL = r"$P_{fit}(r)$"
39PR_LOADED_LABEL = r"$P_{loaded}(r)$"
40IQ_DATA_LABEL = r"$I_{obs}(q)$"
41IQ_FIT_LABEL = r"$I_{fit}(q)$"
42IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
43GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
44GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
45
46
47
48class Plugin(PluginBase):
49    """
50        P(r) inversion perspective
51    """
52    DEFAULT_ALPHA = 0.0001
53    DEFAULT_NFUNC = 10
54    DEFAULT_DMAX = 140.0
55
56    def __init__(self):
57        PluginBase.__init__(self, name="Pr Inversion")
58        ## Simulation window manager
59        self.simview = None
60
61        ## State data
62        self.alpha = self.DEFAULT_ALPHA
63        self.nfunc = self.DEFAULT_NFUNC
64        self.max_length = self.DEFAULT_DMAX
65        self.q_min = None
66        self.q_max = None
67        self.has_bck = False
68        self.slit_height = 0
69        self.slit_width = 0
70        ## Remember last plottable processed
71        self.last_data = ""
72        self._current_file_data = None
73        ## Time elapsed for last computation [sec]
74        # Start with a good default
75        self.elapsed = 0.022
76        self.iq_data_shown = False
77
78        ## Current invertor
79        self.invertor = None
80        self.pr = None
81        self.data_id = IQ_DATA_LABEL
82        # Copy of the last result in case we need to display it.
83        self._last_pr = None
84        self._last_out = None
85        self._last_cov = None
86        ## Calculation thread
87        self.calc_thread = None
88        ## Estimation thread
89        self.estimation_thread = None
90        ## Result panel
91        self.control_panel = None
92        ## Currently views plottable
93        self.current_plottable = None
94        ## Number of P(r) points to display on the output plot
95        self._pr_npts = 51
96        self._normalize_output = False
97        self._scale_output_unity = False
98
99        ## List of added P(r) plots
100        self._added_plots = {}
101        self._default_Iq = {}
102        self.list_plot_id = []
103
104        # Associate the inversion state reader with .prv files
105        from inversion_state import Reader
106
107        # Create a CanSAS/Pr reader
108        self.state_reader = Reader(self.set_state)
109        self._extensions = '.prv'
110        l = Loader()
111        l.associate_file_reader('.prv', self.state_reader)
112        #l.associate_file_reader(".svs", self.state_reader)
113
114        # Log startup
115        logging.info("Pr(r) plug-in started")
116
117    def delete_data(self, data_id):
118        """
119        delete the data association with prview
120        """
121        self.control_panel.clear_panel()
122
123    def get_data(self):
124        """
125            Returns the current data
126        """
127        return self.current_plottable
128
129    def set_state(self, state=None, datainfo=None):
130        """
131        Call-back method for the inversion state reader.
132        This method is called when a .prv file is loaded.
133
134        :param state: InversionState object
135        :param datainfo: Data1D object [optional]
136
137        """
138        try:
139            if datainfo.__class__.__name__ == 'list':
140                if len(datainfo) >= 1:
141                    data = datainfo[0]
142                else:
143                    data = None
144            else:
145                data = datainfo
146            if data is None:
147                msg = "Pr.set_state: datainfo parameter cannot "
148                msg += "be None in standalone mode"
149                raise RuntimeError, msg
150
151            # Ensuring that plots are coordinated correctly
152            t = time.localtime(data.meta_data['prstate'].timestamp)
153            time_str = time.strftime("%b %d %H:%M", t)
154
155            # Check that no time stamp is already appended
156            max_char = data.meta_data['prstate'].file.find("[")
157            if max_char < 0:
158                max_char = len(data.meta_data['prstate'].file)
159
160            datainfo.meta_data['prstate'].file = \
161                data.meta_data['prstate'].file[0:max_char]\
162                + ' [' + time_str + ']'
163
164            data.filename = data.meta_data['prstate'].file
165            # TODO:
166            #remove this call when state save all information about the gui data
167            # such as ID , Group_ID, etc...
168            #make self.current_plottable = datainfo directly
169            self.current_plottable = self.parent.create_gui_data(data, None)
170            self.current_plottable.group_id = data.meta_data['prstate'].file
171
172            # Make sure the user sees the P(r) panel after loading
173            #self.parent.set_perspective(self.perspective)
174            self.on_perspective(event=None)
175            # Load the P(r) results
176            #state = self.state_reader.get_state()
177            data_dict = {self.current_plottable.id:self.current_plottable}
178            self.parent.add_data(data_list=data_dict)
179            wx.PostEvent(self.parent, NewPlotEvent(plot=self.current_plottable,
180                                                   title=self.current_plottable.title))
181            self.control_panel.set_state(state)
182        except:
183            logging.error("prview.set_state: %s" % sys.exc_value)
184
185
186    def help(self, evt):
187        """
188        Show a general help dialog.
189
190        :TODO: replace the text with a nice image
191
192        """
193        from inversion_panel import HelpDialog
194        dialog = HelpDialog(None, -1)
195        if dialog.ShowModal() == wx.ID_OK:
196            dialog.Destroy()
197        else:
198            dialog.Destroy()
199
200    def _fit_pr(self, evt):
201        """
202        """
203        # Generate P(r) for sphere
204        radius = 60.0
205        d_max = 2 * radius
206
207        r = pylab.arange(0.01, d_max, d_max / 51.0)
208        M = len(r)
209        y = numpy.zeros(M)
210        pr_err = numpy.zeros(M)
211
212        total = 0.0
213        for j in range(M):
214            value = self.pr_theory(r[j], radius)
215            total += value
216            y[j] = value
217            pr_err[j] = math.sqrt(y[j])
218
219        y = y / total * d_max / len(r)
220
221        # Perform fit
222        pr = Invertor()
223        pr.d_max = d_max
224        pr.alpha = 0
225        pr.x = r
226        pr.y = y
227        pr.err = pr_err
228        out, cov = pr.pr_fit()
229        for i in range(len(out)):
230            print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
231
232        # Show input P(r)
233        title = "Pr"
234        new_plot = Data1D(pr.x, pr.y, dy=pr.err)
235        new_plot.name = "P_{obs}(r)"
236        new_plot.xaxis("\\rm{r}", 'A')
237        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
238        new_plot.group_id = "P_{obs}(r)"
239        new_plot.id = "P_{obs}(r)"
240        new_plot.title = title
241        self.parent.update_theory(data_id=self.data_id, theory=new_plot)
242        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
243
244        # Show P(r) fit
245        self.show_pr(out, pr)
246
247        # Show I(q) fit
248        q = pylab.arange(0.001, 0.1, 0.01 / 51.0)
249        self.show_iq(out, pr, q)
250
251    def show_shpere(self, x, radius=70.0, x_range=70.0):
252        """
253        """
254        # Show P(r)
255        y_true = numpy.zeros(len(x))
256
257        sum_true = 0.0
258        for i in range(len(x)):
259            y_true[i] = self.pr_theory(x[i], radius)
260            sum_true += y_true[i]
261
262        y_true = y_true / sum_true * x_range / len(x)
263
264        # Show the theory P(r)
265        new_plot = Data1D(x, y_true)
266        new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
267        new_plot.name = "P_{true}(r)"
268        new_plot.xaxis("\\rm{r}", 'A')
269        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
270        new_plot.id = "P_{true}(r)"
271        new_plot.group_id = "P_{true}(r)"
272        self.parent.update_theory(data_id=self.data_id, theory=new_plot)
273        #Put this call in plottables/guitools
274        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot,
275                                               title="Sphere P(r)"))
276
277    def get_npts(self):
278        """
279        Returns the number of points in the I(q) data
280        """
281        try:
282            return len(self.pr.x)
283        except:
284            return 0
285
286    def show_iq(self, out, pr, q=None):
287        """
288            Display computed I(q)
289        """
290        qtemp = pr.x
291        if not q == None:
292            qtemp = q
293
294        # Make a plot
295        maxq = -1
296        for q_i in qtemp:
297            if q_i > maxq:
298                maxq = q_i
299
300        minq = 0.001
301
302        # Check for user min/max
303        if not pr.q_min == None:
304            minq = pr.q_min
305        if not pr.q_max == None:
306            maxq = pr.q_max
307
308        x = pylab.arange(minq, maxq, maxq / 301.0)
309        y = numpy.zeros(len(x))
310        err = numpy.zeros(len(x))
311        for i in range(len(x)):
312            value = pr.iq(out, x[i])
313            y[i] = value
314            try:
315                err[i] = math.sqrt(math.fabs(value))
316            except:
317                err[i] = 1.0
318                print "Error getting error", value, x[i]
319
320        new_plot = Data1D(x, y)
321        new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
322        new_plot.name = IQ_FIT_LABEL
323        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
324        new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
325        title = "I(q)"
326        new_plot.title = title
327
328        # If we have a group ID, use it
329        if pr.info.has_key("plot_group_id"):
330            new_plot.group_id = pr.info["plot_group_id"]
331        new_plot.id = IQ_FIT_LABEL
332        self.parent.update_theory(data_id=self.data_id, theory=new_plot)
333
334        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
335
336        # If we have used slit smearing, plot the smeared I(q) too
337        if pr.slit_width > 0 or pr.slit_height > 0:
338            x = pylab.arange(minq, maxq, maxq / 301.0)
339            y = numpy.zeros(len(x))
340            err = numpy.zeros(len(x))
341            for i in range(len(x)):
342                value = pr.iq_smeared(out, x[i])
343                y[i] = value
344                try:
345                    err[i] = math.sqrt(math.fabs(value))
346                except:
347                    err[i] = 1.0
348                    print "Error getting error", value, x[i]
349
350            new_plot = Data1D(x, y)
351            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
352            new_plot.name = IQ_SMEARED_LABEL
353            new_plot.xaxis("\\rm{Q}", 'A^{-1}')
354            new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
355            # If we have a group ID, use it
356            if pr.info.has_key("plot_group_id"):
357                new_plot.group_id = pr.info["plot_group_id"]
358            new_plot.id = IQ_SMEARED_LABEL
359            new_plot.title = title
360            self.parent.update_theory(data_id=self.data_id, theory=new_plot)
361            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
362
363    def _on_pr_npts(self, evt):
364        """
365        Redisplay P(r) with a different number of points
366        """
367        from inversion_panel import PrDistDialog
368        dialog = PrDistDialog(None, -1)
369        dialog.set_content(self._pr_npts)
370        if dialog.ShowModal() == wx.ID_OK:
371            self._pr_npts = dialog.get_content()
372            dialog.Destroy()
373            self.show_pr(self._last_out, self._last_pr, self._last_cov)
374        else:
375            dialog.Destroy()
376
377
378    def show_pr(self, out, pr, cov=None):
379        """
380        """
381        # Show P(r)
382        x = pylab.arange(0.0, pr.d_max, pr.d_max / self._pr_npts)
383
384        y = numpy.zeros(len(x))
385        dy = numpy.zeros(len(x))
386        y_true = numpy.zeros(len(x))
387
388        total = 0.0
389        pmax = 0.0
390        cov2 = numpy.ascontiguousarray(cov)
391
392        for i in range(len(x)):
393            if cov2 == None:
394                value = pr.pr(out, x[i])
395            else:
396                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
397            total += value * pr.d_max / len(x)
398
399            # keep track of the maximum P(r) value
400            if value > pmax:
401                pmax = value
402
403            y[i] = value
404
405        if self._normalize_output == True:
406            y = y / total
407            dy = dy / total
408        elif self._scale_output_unity == True:
409            y = y / pmax
410            dy = dy / pmax
411
412        if cov2 == None:
413            new_plot = Data1D(x, y)
414            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
415        else:
416            new_plot = Data1D(x, y, dy=dy)
417        new_plot.name = PR_FIT_LABEL
418        new_plot.xaxis("\\rm{r}", 'A')
419        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
420        new_plot.title = "P(r) fit"
421        new_plot.id = PR_FIT_LABEL
422        # Make sure that the plot is linear
423        new_plot.xtransform = "x"
424        new_plot.ytransform = "y"
425        new_plot.group_id = GROUP_ID_PR_FIT
426        self.parent.update_theory(data_id=self.data_id, theory=new_plot)
427        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
428        return x, pr.d_max
429
430    def load(self, data):
431        """
432        Load data. This will eventually be replaced
433        by our standard DataLoader class.
434        """
435        class FileData(object):
436            x = None
437            y = None
438            err = None
439            path = None
440
441            def __init__(self, path):
442                self.path = path
443
444        self._current_file_data = FileData(data.path)
445
446        # Use data loader to load file
447        dataread = data
448        # Notify the user if we could not read the file
449        if dataread is None:
450            raise RuntimeError, "Invalid data"
451
452        x = None
453        y = None
454        err = None
455        if dataread.__class__.__name__ == 'Data1D':
456            x = dataread.x
457            y = dataread.y
458            err = dataread.dy
459        else:
460            if isinstance(dataread, list) and len(dataread) > 0:
461                x = dataread[0].x
462                y = dataread[0].y
463                err = dataread[0].dy
464                msg = "PrView only allows a single data set at a time. "
465                msg += "Only the first data set was loaded."
466                wx.PostEvent(self.parent, StatusEvent(status=msg))
467            else:
468                if dataread is None:
469                    return x, y, err
470                raise RuntimeError, "This tool can only read 1D data"
471
472        self._current_file_data.x = x
473        self._current_file_data.y = y
474        self._current_file_data.err = err
475        return x, y, err
476
477    def load_columns(self, path="sphere_60_q0_2.txt"):
478        """
479        Load 2- or 3- column ascii
480        """
481        # Read the data from the data file
482        data_x = numpy.zeros(0)
483        data_y = numpy.zeros(0)
484        data_err = numpy.zeros(0)
485        scale = None
486        min_err = 0.0
487        if not path == None:
488            input_f = open(path, 'r')
489            buff = input_f.read()
490            lines = buff.split('\n')
491            for line in lines:
492                try:
493                    toks = line.split()
494                    x = float(toks[0])
495                    y = float(toks[1])
496                    if len(toks) > 2:
497                        err = float(toks[2])
498                    else:
499                        if scale == None:
500                            scale = 0.05 * math.sqrt(y)
501                            #scale = 0.05/math.sqrt(y)
502                            min_err = 0.01 * y
503                        err = scale * math.sqrt(y) + min_err
504                        #err = 0
505
506                    data_x = numpy.append(data_x, x)
507                    data_y = numpy.append(data_y, y)
508                    data_err = numpy.append(data_err, err)
509                except:
510                    logging.error(sys.exc_value)
511
512        if not scale == None:
513            message = "The loaded file had no error bars, statistical errors are assumed."
514            wx.PostEvent(self.parent, StatusEvent(status=message))
515        else:
516            wx.PostEvent(self.parent, StatusEvent(status=''))
517
518        return data_x, data_y, data_err
519
520    def load_abs(self, path):
521        """
522        Load an IGOR .ABS reduced file
523
524        :param path: file path
525
526        :return: x, y, err vectors
527
528        """
529        # Read the data from the data file
530        data_x = numpy.zeros(0)
531        data_y = numpy.zeros(0)
532        data_err = numpy.zeros(0)
533        scale = None
534        min_err = 0.0
535
536        data_started = False
537        if not path == None:
538            input_f = open(path, 'r')
539            buff = input_f.read()
540            lines = buff.split('\n')
541            for line in lines:
542                if data_started == True:
543                    try:
544                        toks = line.split()
545                        x = float(toks[0])
546                        y = float(toks[1])
547                        if len(toks) > 2:
548                            err = float(toks[2])
549                        else:
550                            if scale == None:
551                                scale = 0.05 * math.sqrt(y)
552                                #scale = 0.05/math.sqrt(y)
553                                min_err = 0.01 * y
554                            err = scale * math.sqrt(y) + min_err
555                            #err = 0
556
557                        data_x = numpy.append(data_x, x)
558                        data_y = numpy.append(data_y, y)
559                        data_err = numpy.append(data_err, err)
560                    except:
561                        logging.error(sys.exc_value)
562                elif line.find("The 6 columns") >= 0:
563                    data_started = True
564
565        if not scale == None:
566            message = "The loaded file had no error bars, statistical errors are assumed."
567            wx.PostEvent(self.parent, StatusEvent(status=message))
568        else:
569            wx.PostEvent(self.parent, StatusEvent(status=''))
570
571        return data_x, data_y, data_err
572
573    def pr_theory(self, r, R):
574        """
575            Return P(r) of a sphere for a given R
576            For test purposes
577        """
578        if r <= 2 * R:
579            return 12.0 * ((0.5 * r / R) ** 2) * ((1.0 - 0.5 * r / R) ** 2) * (2.0 + 0.5 * r / R)
580        else:
581            return 0.0
582
583    def get_context_menu(self, plotpanel=None):
584        """
585        Get the context menu items available for P(r)
586
587        :param graph: the Graph object to which we attach the context menu
588
589        :return: a list of menu items with call-back function
590
591        """
592        graph = plotpanel.graph
593        # Look whether this Graph contains P(r) data
594        if graph.selected_plottable not in plotpanel.plots:
595            return []
596        item = plotpanel.plots[graph.selected_plottable]
597        if item.id == PR_FIT_LABEL:
598            #add_data_hint = "Load a data file and display it on this plot"
599            #["Add P(r) data",add_data_hint , self._on_add_data],
600            change_n_hint = "Change the number of"
601            change_n_hint += " points on the P(r) output"
602            change_n_label = "Change number of P(r) points"
603            m_list = [[change_n_label, change_n_hint, self._on_pr_npts]]
604
605            if self._scale_output_unity or self._normalize_output:
606                hint = "Let the output P(r) keep the scale of the data"
607                m_list.append(["Disable P(r) scaling", hint,
608                               self._on_disable_scaling])
609            if not self._scale_output_unity:
610                m_list.append(["Scale P_max(r) to unity",
611                               "Scale P(r) so that its maximum is 1",
612                               self._on_scale_unity])
613            if not self._normalize_output:
614                m_list.append(["Normalize P(r) to unity",
615                               "Normalize the integral of P(r) to 1",
616                               self._on_normalize])
617
618            return m_list
619
620        elif item.id in [PR_LOADED_LABEL, IQ_DATA_LABEL, IQ_FIT_LABEL, IQ_SMEARED_LABEL]:
621            return []
622        elif item.id == graph.selected_plottable:
623            if issubclass(item.__class__, Data1D):
624                return [["Compute P(r)",
625                         "Compute P(r) from distribution",
626                         self._on_context_inversion]]
627        return []
628
629    def _on_disable_scaling(self, evt):
630        """
631        Disable P(r) scaling
632
633        :param evt: Menu event
634
635        """
636        self._normalize_output = False
637        self._scale_output_unity = False
638        self.show_pr(self._last_out, self._last_pr, self._last_cov)
639
640        # Now replot the original added data
641        for plot in self._added_plots:
642            self._added_plots[plot].y = numpy.copy(self._default_Iq[plot])
643            wx.PostEvent(self.parent,
644                         NewPlotEvent(plot=self._added_plots[plot],
645                                      title=self._added_plots[plot].name,
646                                      update=True))
647
648        # Need the update flag in the NewPlotEvent to protect against
649        # the plot no longer being there...
650
651    def _on_normalize(self, evt):
652        """
653        Normalize the area under the P(r) curve to 1.
654        This operation is done for all displayed plots.
655
656        :param evt: Menu event
657
658        """
659        self._normalize_output = True
660        self._scale_output_unity = False
661
662        self.show_pr(self._last_out, self._last_pr, self._last_cov)
663
664        # Now scale the added plots too
665        for plot in self._added_plots:
666            total = numpy.sum(self._added_plots[plot].y)
667            npts = len(self._added_plots[plot].x)
668            total *= self._added_plots[plot].x[npts - 1] / npts
669            y = self._added_plots[plot].y / total
670
671            new_plot = Data1D(self._added_plots[plot].x, y)
672            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
673            new_plot.group_id = self._added_plots[plot].group_id
674            new_plot.id = self._added_plots[plot].id
675            new_plot.title = self._added_plots[plot].title
676            new_plot.name = self._added_plots[plot].name
677            new_plot.xaxis("\\rm{r}", 'A')
678            new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
679            self.parent.update_theory(data_id=self.data_id, theory=new_plot)
680            wx.PostEvent(self.parent,
681                         NewPlotEvent(plot=new_plot, update=True,
682                                      title=self._added_plots[plot].name))
683
684    def _on_scale_unity(self, evt):
685        """
686        Scale the maximum P(r) value on each displayed plot to 1.
687
688        :param evt: Menu event
689
690        """
691        self._scale_output_unity = True
692        self._normalize_output = False
693
694        self.show_pr(self._last_out, self._last_pr, self._last_cov)
695
696        # Now scale the added plots too
697        for plot in self._added_plots:
698            _max = 0
699            for y in self._added_plots[plot].y:
700                if y > _max:
701                    _max = y
702            y = self._added_plots[plot].y / _max
703
704            new_plot = Data1D(self._added_plots[plot].x, y)
705            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
706            new_plot.name = self._added_plots[plot].name
707            new_plot.xaxis("\\rm{r}", 'A')
708            new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
709            self.parent.update_theory(data_id=self.data_id, theory=new_plot)
710            wx.PostEvent(self.parent,
711                         NewPlotEvent(plot=new_plot, update=True,
712                                      title=self._added_plots[plot].name))
713
714    def start_thread(self):
715        """
716            Start a calculation thread
717        """
718        from pr_thread import CalcPr
719
720        # If a thread is already started, stop it
721        if self.calc_thread != None and self.calc_thread.isrunning():
722            self.calc_thread.stop()
723            ## stop just raises the flag -- the thread is supposed to
724            ## then kill itself. In August 2014 it was shown that this is
725            ## incorrectly handled by fitting.py and a fix implemented.
726            ## It is not clear whether it is properly used here, but the
727            ## "fix" of waiting for the previous thread to end breaks the
728            ## pr perspective completely as it causes an infinite loop.
729            ## Thus it is likely the threading is bing properly handled.
730            ## While the "fix" is no longer implemented the comment is
731            ## left here till somebody ascertains that in fact the threads
732            ## are being properly handled.
733            ##
734            ##    -PDB January 25, 2015
735
736        pr = self.pr.clone()
737        self.calc_thread = CalcPr(pr, self.nfunc,
738                                  error_func=self._thread_error,
739                                  completefn=self._completed, updatefn=None)
740        self.calc_thread.queue()
741        self.calc_thread.ready(2.5)
742
743    def _thread_error(self, error):
744        """
745            Call-back method for calculation errors
746        """
747        wx.PostEvent(self.parent, StatusEvent(status=error))
748
749    def _estimate_completed(self, alpha, message, elapsed):
750        """
751        Parameter estimation completed,
752        display the results to the user
753
754        :param alpha: estimated best alpha
755        :param elapsed: computation time
756
757        """
758        # Save useful info
759        self.elapsed = elapsed
760        self.control_panel.alpha_estimate = alpha
761        if not message == None:
762            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
763        self.perform_estimateNT()
764
765    def _estimateNT_completed(self, nterms, alpha, message, elapsed):
766        """
767        Parameter estimation completed,
768        display the results to the user
769
770        :param alpha: estimated best alpha
771        :param nterms: estimated number of terms
772        :param elapsed: computation time
773
774        """
775        # Save useful info
776        self.elapsed = elapsed
777        self.control_panel.nterms_estimate = nterms
778        self.control_panel.alpha_estimate = alpha
779        if not message == None:
780            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
781
782    def _completed(self, out, cov, pr, elapsed):
783        """
784        wxCallAfter Method called with the results when the inversion
785        is done
786
787        :param out: output coefficient for the base functions
788        :param cov: covariance matrix
789        :param pr: Invertor instance
790        :param elapsed: time spent computing
791        """
792        # Ensure hat you have all inputs are ready at the time call happens:
793        # Without CallAfter, it will freeze with wx >= 2.9.
794        wx.CallAfter(self._completed_call, out, cov, pr, elapsed)
795
796    def _completed_call(self, out, cov, pr, elapsed):
797        """
798        Method called with the results when the inversion
799        is done
800
801        :param out: output coefficient for the base functions
802        :param cov: covariance matrix
803        :param pr: Invertor instance
804        :param elapsed: time spent computing
805
806        """
807        # Save useful info
808        self.elapsed = elapsed
809        # Keep a copy of the last result
810        self._last_pr = pr.clone()
811        self._last_out = out
812        self._last_cov = cov
813
814        # Save Pr invertor
815        self.pr = pr
816        cov = numpy.ascontiguousarray(cov)
817
818        # Show result on control panel
819        self.control_panel.chi2 = pr.chi2
820        self.control_panel.elapsed = elapsed
821        self.control_panel.oscillation = pr.oscillations(out)
822        self.control_panel.positive = pr.get_positive(out)
823        self.control_panel.pos_err = pr.get_pos_err(out, cov)
824        self.control_panel.rg = pr.rg(out)
825        self.control_panel.iq0 = pr.iq0(out)
826        self.control_panel.bck = pr.background
827
828        # Show I(q) fit
829        self.show_iq(out, self.pr)
830
831        # Show P(r) fit
832        self.show_pr(out, self.pr, cov)
833
834    def show_data(self, path=None, data=None, reset=False):
835        """
836        Show data read from a file
837
838        :param path: file path
839        :param reset: if True all other plottables will be cleared
840
841        """
842        #if path is not None:
843        if data is not None:
844            try:
845                pr = self._create_file_pr(data)
846            except:
847                status = "Problem reading data: %s" % sys.exc_value
848                wx.PostEvent(self.parent, StatusEvent(status=status))
849                raise RuntimeError, status
850
851            # If the file contains nothing, just return
852            if pr is None:
853                raise RuntimeError, "Loaded data is invalid"
854
855            self.pr = pr
856
857        # Make a plot of I(q) data
858        if self.pr.err == None:
859            new_plot = Data1D(self.pr.x, self.pr.y)
860            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
861        else:
862            new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
863        new_plot.name = IQ_DATA_LABEL
864        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
865        new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
866        new_plot.interactive = True
867        new_plot.group_id = GROUP_ID_IQ_DATA
868        new_plot.id = self.data_id
869        new_plot.title = "I(q)"
870        wx.PostEvent(self.parent,
871                     NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
872
873        self.current_plottable = new_plot
874        # Get Q range
875        self.control_panel.q_min = min(self.pr.x)
876        self.control_panel.q_max = max(self.pr.x)
877
878    def save_data(self, filepath, prstate=None):
879        """
880        Save data in provided state object.
881
882        :TODO: move the state code away from inversion_panel and move it here.
883                Then remove the "prstate" input and make this method private.
884
885        :param filepath: path of file to write to
886        :param prstate: P(r) inversion state
887
888        """
889        #TODO: do we need this or can we use DataLoader.loader.save directly?
890
891        # Add output data and coefficients to state
892        prstate.coefficients = self._last_out
893        prstate.covariance = self._last_cov
894
895        # Write the output to file
896        # First, check that the data is of the right type
897        if issubclass(self.current_plottable.__class__,
898                      sas.sascalc.dataloader.data_info.Data1D):
899            self.state_reader.write(filepath, self.current_plottable, prstate)
900        else:
901            msg = "pr.save_data: the data being saved is not a"
902            msg += " sas.data_info.Data1D object"
903            raise RuntimeError, msg
904
905    def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
906                             bck=False, height=0, width=0):
907        """
908            Set up inversion from plotted data
909        """
910        self.alpha = alpha
911        self.nfunc = nfunc
912        self.max_length = d_max
913        self.q_min = q_min
914        self.q_max = q_max
915        self.has_bck = bck
916        self.slit_height = height
917        self.slit_width = width
918
919        try:
920            pr = self._create_plot_pr()
921            if not pr == None:
922                self.pr = pr
923                self.perform_inversion()
924        except:
925            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
926
927    def estimate_plot_inversion(self, alpha, nfunc, d_max,
928                                q_min=None, q_max=None,
929                                bck=False, height=0, width=0):
930        """
931            Estimate parameters from plotted data
932        """
933        self.alpha = alpha
934        self.nfunc = nfunc
935        self.max_length = d_max
936        self.q_min = q_min
937        self.q_max = q_max
938        self.has_bck = bck
939        self.slit_height = height
940        self.slit_width = width
941
942        try:
943            pr = self._create_plot_pr()
944            if not pr == None:
945                self.pr = pr
946                self.perform_estimate()
947        except:
948            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
949
950    def _create_plot_pr(self, estimate=False):
951        """
952        Create and prepare invertor instance from
953        a plottable data set.
954
955        :param path: path of the file to read in
956
957        """
958        # Sanity check
959        if self.current_plottable is None:
960            msg = "Please load a valid data set before proceeding."
961            wx.PostEvent(self.parent, StatusEvent(status=msg))
962            return None
963
964        # Get the data from the chosen data set and perform inversion
965        pr = Invertor()
966        pr.d_max = self.max_length
967        pr.alpha = self.alpha
968        pr.q_min = self.q_min
969        pr.q_max = self.q_max
970        pr.x = self.current_plottable.x
971        pr.y = self.current_plottable.y
972        pr.has_bck = self.has_bck
973        pr.slit_height = self.slit_height
974        pr.slit_width = self.slit_width
975
976        # Keep track of the plot window title to ensure that
977        # we can overlay the plots
978        pr.info["plot_group_id"] = self.current_plottable.group_id
979
980        # Fill in errors if none were provided
981        err = self.current_plottable.dy
982        all_zeros = True
983        if err == None:
984            err = numpy.zeros(len(pr.y))
985        else:
986            for i in range(len(err)):
987                if err[i] > 0:
988                    all_zeros = False
989
990        if all_zeros:
991            scale = None
992            min_err = 0.0
993            for i in range(len(pr.y)):
994                # Scale the error so that we can fit over several decades of Q
995                if scale == None:
996                    scale = 0.05 * math.sqrt(pr.y[i])
997                    min_err = 0.01 * pr.y[i]
998                err[i] = scale * math.sqrt(math.fabs(pr.y[i])) + min_err
999            message = "The loaded file had no error bars, "
1000            message += "statistical errors are assumed."
1001            wx.PostEvent(self.parent, StatusEvent(status=message))
1002
1003        pr.err = err
1004
1005        return pr
1006
1007    def setup_file_inversion(self, alpha, nfunc, d_max, data,
1008                             path=None, q_min=None, q_max=None,
1009                             bck=False, height=0, width=0):
1010        """
1011            Set up inversion
1012        """
1013        self.alpha = alpha
1014        self.nfunc = nfunc
1015        self.max_length = d_max
1016        self.q_min = q_min
1017        self.q_max = q_max
1018        self.has_bck = bck
1019        self.slit_height = height
1020        self.slit_width = width
1021
1022        try:
1023            pr = self._create_file_pr(data)
1024            if not pr == None:
1025                self.pr = pr
1026                self.perform_inversion()
1027        except:
1028            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1029
1030    def estimate_file_inversion(self, alpha, nfunc, d_max, data,
1031                                path=None, q_min=None, q_max=None,
1032                                bck=False, height=0, width=0):
1033        """
1034            Estimate parameters for inversion
1035        """
1036        self.alpha = alpha
1037        self.nfunc = nfunc
1038        self.max_length = d_max
1039        self.q_min = q_min
1040        self.q_max = q_max
1041        self.has_bck = bck
1042        self.slit_height = height
1043        self.slit_width = width
1044
1045        try:
1046            pr = self._create_file_pr(data)
1047            if not pr is None:
1048                self.pr = pr
1049                self.perform_estimate()
1050        except:
1051            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1052
1053    def _create_file_pr(self, data):
1054        """
1055        Create and prepare invertor instance from
1056        a file data set.
1057
1058        :param path: path of the file to read in
1059
1060        """
1061        # Reset the status bar so that we don't get mixed up
1062        # with old messages.
1063        #TODO: refactor this into a proper status handling
1064        wx.PostEvent(self.parent, StatusEvent(status=''))
1065        try:
1066            class FileData(object):
1067                x = None
1068                y = None
1069                err = None
1070                path = None
1071                def __init__(self, path):
1072                    self.path = path
1073
1074            self._current_file_data = FileData(data.path)
1075            self._current_file_data.x = data.x
1076            self._current_file_data.y = data.y
1077            self._current_file_data.err = data.dy
1078            x, y, err = data.x, data.y, data.dy
1079        except:
1080            load_error(sys.exc_value)
1081            return None
1082
1083        # If the file contains no data, just return
1084        if x is None or len(x) == 0:
1085            load_error("The loaded file contains no data")
1086            return None
1087
1088        # If we have not errors, add statistical errors
1089        if y is not None:
1090            if err == None or numpy.all(err) == 0:
1091                err = numpy.zeros(len(y))
1092                scale = None
1093                min_err = 0.0
1094                for i in range(len(y)):
1095                    # Scale the error so that we can fit over several decades of Q
1096                    if scale == None:
1097                        scale = 0.05 * math.sqrt(y[i])
1098                        min_err = 0.01 * y[i]
1099                    err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
1100                message = "The loaded file had no error bars, "
1101                message += "statistical errors are assumed."
1102                wx.PostEvent(self.parent, StatusEvent(status=message))
1103
1104        try:
1105            # Get the data from the chosen data set and perform inversion
1106            pr = Invertor()
1107            pr.d_max = self.max_length
1108            pr.alpha = self.alpha
1109            pr.q_min = self.q_min
1110            pr.q_max = self.q_max
1111            pr.x = x
1112            pr.y = y
1113            pr.err = err
1114            pr.has_bck = self.has_bck
1115            pr.slit_height = self.slit_height
1116            pr.slit_width = self.slit_width
1117            return pr
1118        except:
1119            load_error(sys.exc_value)
1120        return None
1121
1122    def perform_estimate(self):
1123        """
1124            Perform parameter estimation
1125        """
1126        from pr_thread import EstimatePr
1127
1128        # If a thread is already started, stop it
1129        if self.estimation_thread != None and \
1130            self.estimation_thread.isrunning():
1131            self.estimation_thread.stop()
1132            ## stop just raises the flag -- the thread is supposed to
1133            ## then kill itself. In August 2014 it was shown that this is
1134            ## incorrectly handled by fitting.py and a fix implemented.
1135            ## It is not clear whether it is properly used here, but the
1136            ## "fix" of waiting for the previous thread to end breaks the
1137            ## pr perspective completely as it causes an infinite loop.
1138            ## Thus it is likely the threading is bing properly handled.
1139            ## While the "fix" is no longer implemented the comment is
1140            ## left here till somebody ascertains that in fact the threads
1141            ## are being properly handled.
1142            ##
1143            ##    -PDB January 25, 2015
1144        pr = self.pr.clone()
1145        self.estimation_thread = EstimatePr(pr, self.nfunc,
1146                                            error_func=self._thread_error,
1147                                            completefn=self._estimate_completed,
1148                                            updatefn=None)
1149        self.estimation_thread.queue()
1150        self.estimation_thread.ready(2.5)
1151
1152    def perform_estimateNT(self):
1153        """
1154            Perform parameter estimation
1155        """
1156        from pr_thread import EstimateNT
1157
1158        # If a thread is already started, stop it
1159        if self.estimation_thread != None and self.estimation_thread.isrunning():
1160            self.estimation_thread.stop()
1161            ## stop just raises the flag -- the thread is supposed to
1162            ## then kill itself. In August 2014 it was shown that this is
1163            ## incorrectly handled by fitting.py and a fix implemented.
1164            ## It is not clear whether it is properly used here, but the
1165            ## "fix" of waiting for the previous thread to end breaks the
1166            ## pr perspective completely as it causes an infinite loop.
1167            ## Thus it is likely the threading is bing properly handled.
1168            ## While the "fix" is no longer implemented the comment is
1169            ## left here till somebody ascertains that in fact the threads
1170            ## are being properly handled.
1171            ##
1172            ##    -PDB January 25, 2015
1173        pr = self.pr.clone()
1174        # Skip the slit settings for the estimation
1175        # It slows down the application and it doesn't change the estimates
1176        pr.slit_height = 0.0
1177        pr.slit_width = 0.0
1178        self.estimation_thread = EstimateNT(pr, self.nfunc,
1179                                            error_func=self._thread_error,
1180                                            completefn=self._estimateNT_completed,
1181                                            updatefn=None)
1182        self.estimation_thread.queue()
1183        self.estimation_thread.ready(2.5)
1184
1185    def perform_inversion(self):
1186        """
1187            Perform inversion
1188        """
1189        self.start_thread()
1190
1191    def _on_context_inversion(self, event):
1192        """
1193            Call-back method for plot context menu
1194        """
1195        panel = event.GetEventObject()
1196        Plugin.on_perspective(self, event=event)
1197
1198        # If we have more than one displayed plot, make the user choose
1199        if len(panel.plots) >= 1 and \
1200            panel.graph.selected_plottable in panel.plots:
1201            dataset = panel.plots[panel.graph.selected_plottable].name
1202        else:
1203            logging.info("Prview Error: No data is available")
1204            return
1205
1206        # Store a reference to the current plottable
1207        # If we have a suggested value, use it.
1208        try:
1209            estimate = float(self.control_panel.alpha_estimate)
1210            self.control_panel.alpha = estimate
1211        except:
1212            self.control_panel.alpha = self.alpha
1213            logging.info("Prview :Alpha Not estimate yet")
1214        try:
1215            estimate = int(self.control_panel.nterms_estimate)
1216            self.control_panel.nfunc = estimate
1217        except:
1218            self.control_panel.nfunc = self.nfunc
1219            logging.info("Prview : ntemrs Not estimate yet")
1220
1221        self.current_plottable = panel.plots[panel.graph.selected_plottable]
1222        self.set_data([self.current_plottable])
1223        self.control_panel.plotname = dataset
1224        #self.control_panel.nfunc = self.nfunc
1225        self.control_panel.d_max = self.max_length
1226        #self.parent.set_perspective(self.perspective)
1227        self.control_panel._on_invert(None)
1228
1229    def get_panels(self, parent):
1230        """
1231            Create and return a list of panel objects
1232        """
1233        from inversion_panel import InversionControl
1234
1235        self.parent = parent
1236        self.frame = MDIFrame(self.parent, None, 'None', (100, 200))
1237        self.control_panel = InversionControl(self.frame, -1,
1238                                              style=wx.RAISED_BORDER)
1239        self.frame.set_panel(self.control_panel)
1240        self._frame_set_helper()
1241        self.control_panel.set_manager(self)
1242        self.control_panel.nfunc = self.nfunc
1243        self.control_panel.d_max = self.max_length
1244        self.control_panel.alpha = self.alpha
1245        self.perspective = []
1246        self.perspective.append(self.control_panel.window_name)
1247
1248        return [self.control_panel]
1249
1250    def set_data(self, data_list=None):
1251        """
1252        receive a list of data to compute pr
1253        """
1254        if data_list is None:
1255            data_list = []
1256        if len(data_list) >= 1:
1257            if len(data_list) == 1:
1258                data = data_list[0]
1259            else:
1260                data_1d_list = []
1261                data_2d_list = []
1262                error_msg = ""
1263                # separate data into data1d and data2d list
1264                for data in data_list:
1265                    if data is not None:
1266                        if issubclass(data.__class__, Data1D):
1267                            data_1d_list.append(data)
1268                        else:
1269                            error_msg += " %s type %s \n" % (str(data.name),
1270                                                             str(data.__class__.__name__))
1271                            data_2d_list.append(data)
1272                if len(data_2d_list) > 0:
1273                    msg = "PrView does not support the following data types:\n"
1274                    msg += error_msg
1275                if len(data_1d_list) == 0:
1276                    wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
1277                    return
1278                msg += "Prview does not allow multiple data!\n"
1279                msg += "Please select one.\n"
1280                if len(data_list) > 1:
1281                    from pr_widgets import DataDialog
1282                    dlg = DataDialog(data_list=data_1d_list, text=msg)
1283                    if dlg.ShowModal() == wx.ID_OK:
1284                        data = dlg.get_data()
1285                    else:
1286                        data = None
1287                    dlg.Destroy()
1288            if data is None:
1289                msg += "PrView receives no data. \n"
1290                wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
1291                return
1292            if issubclass(data.__class__, Data1D):
1293                try:
1294                    wx.PostEvent(self.parent,
1295                                 NewPlotEvent(action='remove',
1296                                              group_id=GROUP_ID_IQ_DATA,
1297                                              id=self.data_id))
1298                    self.data_id = data.id
1299                    self.control_panel._change_file(evt=None, data=data)
1300                except:
1301                    msg = "Prview Set_data: " + str(sys.exc_value)
1302                    wx.PostEvent(self.parent, StatusEvent(status=msg, info="error"))
1303            else:
1304                msg = "Pr cannot be computed for data of "
1305                msg += "type %s" % (data_list[0].__class__.__name__)
1306                wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
1307        else:
1308            msg = "Pr contain no data"
1309            wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
1310
1311    def post_init(self):
1312        """
1313            Post initialization call back to close the loose ends
1314            [Somehow openGL needs this call]
1315        """
1316        pass
Note: See TracBrowser for help on using the repository browser.