source: sasview/prview/perspectives/pr/pr.py @ 6fea919

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 6fea919 was ceaf16e, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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