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

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.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since c1d5aea was c1d5aea, checked in by andyfaff, 7 years ago

MAINT: use 'x is not None', instead of 'not x is None'.

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