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

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since c972e3a was fa81e94, checked in by Piotr Rozyczko <rozyczko@…>, 7 years ago

Initial commit of the P(r) inversion perspective.
Code merged from Jeff Krzywon's ESS_GUI_Pr branch.
Also, minor 2to3 mods to sascalc/sasgui to enble error free setup.
  • Property mode set to 100755
File size: 46.1 KB
Line 
1"""
2    P(r) perspective for SasView
3"""
4################################################################################
5#This software was developed by the University of Tennessee as part of the
6#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
7#project funded by the US National Science Foundation.
8#
9#See the license text in license.txt
10#
11#copyright 2009, University of Tennessee
12################################################################################
13
14
15# Make sure the option of saving each curve is available
16# Use the I(q) curve as input and compare the output to P(r)
17
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_info()[1])
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_info()[1])
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_info()[1])
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_info()[1]
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_info()[1]))
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_info()[1]))
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_info()[1]))
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_info()[1]))
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_info()[1])
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_info()[1])
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_info()[1])
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.