source: sasview/src/sans/perspectives/pr/pr.py @ f4167637

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.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since f4167637 was f4167637, checked in by butler, 10 years ago

update Pr and simulations to fix thread problem. note: I don't think simulations is used at this point but figure it should be fixed if in repo.

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