source: sasview/prview/perspectives/pr/pr.py @ 675e0ab

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 675e0ab was 410aad8, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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