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

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 cb62bd5 was cb62bd5, checked in by lewis, 5 years ago

Use manually inputted background level in Pr calculation

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