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

Last change on this file since d620d03c was 5251ec6, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

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