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

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 c80c06f was e43c012, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

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