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

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 a1b8fee was a1b8fee, checked in by andyfaff, 7 years ago

MAINT: from future import print_function

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