source: sasview/prview/perspectives/pr/pr.py @ d4ccbb1

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 d4ccbb1 was 302510b, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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