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

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 6f73a08 was 14d05ba, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

minor update
  • Property mode set to 100644
File size: 39.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            raise RuntimeError, "This tool can only read 1D data"
369       
370        self._current_file_data.x = x
371        self._current_file_data.y = y
372        self._current_file_data.err = err
373        return x, y, err
374               
375    def load_columns(self, path = "sphere_60_q0_2.txt"):
376        """
377            Load 2- or 3- column ascii
378        """
379        import numpy, math, sys
380        # Read the data from the data file
381        data_x   = numpy.zeros(0)
382        data_y   = numpy.zeros(0)
383        data_err = numpy.zeros(0)
384        scale    = None
385        min_err  = 0.0
386        if not path == None:
387            input_f = open(path,'r')
388            buff    = input_f.read()
389            lines   = buff.split('\n')
390            for line in lines:
391                try:
392                    toks = line.split()
393                    x = float(toks[0])
394                    y = float(toks[1])
395                    if len(toks)>2:
396                        err = float(toks[2])
397                    else:
398                        if scale==None:
399                            scale = 0.05*math.sqrt(y)
400                            #scale = 0.05/math.sqrt(y)
401                            min_err = 0.01*y
402                        err = scale*math.sqrt(y)+min_err
403                        #err = 0
404                       
405                    data_x = numpy.append(data_x, x)
406                    data_y = numpy.append(data_y, y)
407                    data_err = numpy.append(data_err, err)
408                except:
409                    pass
410                   
411        if not scale==None:
412            message = "The loaded file had no error bars, statistical errors are assumed."
413            wx.PostEvent(self.parent, StatusEvent(status=message))
414        else:
415            wx.PostEvent(self.parent, StatusEvent(status=''))
416                       
417        return data_x, data_y, data_err     
418       
419    def load_abs(self, path):
420        """
421            Load an IGOR .ABS reduced file
422            @param path: file path
423            @return: x, y, err vectors
424        """
425        import numpy, math, sys
426        # Read the data from the data file
427        data_x   = numpy.zeros(0)
428        data_y   = numpy.zeros(0)
429        data_err = numpy.zeros(0)
430        scale    = None
431        min_err  = 0.0
432       
433        data_started = False
434        if not path == None:
435            input_f = open(path,'r')
436            buff    = input_f.read()
437            lines   = buff.split('\n')
438            for line in lines:
439                if data_started==True:
440                    try:
441                        toks = line.split()
442                        x = float(toks[0])
443                        y = float(toks[1])
444                        if len(toks)>2:
445                            err = float(toks[2])
446                        else:
447                            if scale==None:
448                                scale = 0.05*math.sqrt(y)
449                                #scale = 0.05/math.sqrt(y)
450                                min_err = 0.01*y
451                            err = scale*math.sqrt(y)+min_err
452                            #err = 0
453                           
454                        data_x = numpy.append(data_x, x)
455                        data_y = numpy.append(data_y, y)
456                        data_err = numpy.append(data_err, err)
457                    except:
458                        pass
459                elif line.find("The 6 columns")>=0:
460                    data_started = True     
461                   
462        if not scale==None:
463            message = "The loaded file had no error bars, statistical errors are assumed."
464            wx.PostEvent(self.parent, StatusEvent(status=message))
465        else:
466            wx.PostEvent(self.parent, StatusEvent(status=''))
467                       
468        return data_x, data_y, data_err     
469       
470       
471       
472    def pr_theory(self, r, R):
473        """
474           
475        """
476        if r<=2*R:
477            return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
478        else:
479            return 0.0
480
481    def get_context_menu(self, graph=None):
482        """
483            Get the context menu items available for P(r)
484            @param graph: the Graph object to which we attach the context menu
485            @return: a list of menu items with call-back function
486        """
487        # Look whether this Graph contains P(r) data
488        #if graph.selected_plottable==IQ_DATA_LABEL:
489        for item in graph.plottables:
490            if item.name==PR_FIT_LABEL:
491                m_list = [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
492                       ["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
493
494                if self._scale_output_unity==True or self._normalize_output==True:
495                    m_list.append(["Disable P(r) scaling", 
496                                   "Let the output P(r) keep the scale of the data", 
497                                   self._on_disable_scaling])
498               
499                if self._scale_output_unity==False:
500                    m_list.append(["Scale P_max(r) to unity", 
501                                   "Scale P(r) so that its maximum is 1", 
502                                   self._on_scale_unity])
503                   
504                if self._normalize_output==False:
505                    m_list.append(["Normalize P(r) to unity", 
506                                   "Normalize the integral of P(r) to 1", 
507                                   self._on_normalize])
508                   
509                return m_list
510                #return [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
511                #       ["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
512
513            elif item.name==graph.selected_plottable:
514                return [["Compute P(r)", "Compute P(r) from distribution", self._on_context_inversion]]     
515               
516        return []
517
518    def _on_disable_scaling(self, evt):
519        """
520            Disable P(r) scaling
521            @param evt: Menu event
522        """
523        self._normalize_output = False
524        self._scale_output_unity = False
525        self.show_pr(self._last_out, self._last_pr, self._last_cov)
526       
527        # Now replot the original added data
528        for plot in self._added_plots:
529            self._added_plots[plot].y = numpy.copy(self._default_Iq[plot])
530            wx.PostEvent(self.parent, NewPlotEvent(plot=self._added_plots[plot], 
531                                                   title=self._added_plots[plot].name,
532                                                   update=True))       
533       
534        # Need the update flag in the NewPlotEvent to protect against
535        # the plot no longer being there...
536       
537    def _on_normalize(self, evt):
538        """
539            Normalize the area under the P(r) curve to 1.
540            This operation is done for all displayed plots.
541           
542            @param evt: Menu event
543        """
544        self._normalize_output = True
545        self._scale_output_unity = False
546           
547        self.show_pr(self._last_out, self._last_pr, self._last_cov)
548       
549        # Now scale the added plots too
550        for plot in self._added_plots:
551            sum = numpy.sum(self._added_plots[plot].y)
552            npts = len(self._added_plots[plot].x)
553            sum *= self._added_plots[plot].x[npts-1]/npts
554            y = self._added_plots[plot].y/sum
555           
556            new_plot = Theory1D(self._added_plots[plot].x, y)
557            new_plot.name = self._added_plots[plot].name
558            new_plot.xaxis("\\rm{r}", 'A')
559            new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
560           
561            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, update=True,
562                                                   title=self._added_plots[plot].name))
563               
564       
565       
566    def _on_scale_unity(self, evt):
567        """
568            Scale the maximum P(r) value on each displayed plot to 1.
569           
570            @param evt: Menu event
571        """
572        self._scale_output_unity = True
573        self._normalize_output = False
574           
575        self.show_pr(self._last_out, self._last_pr, self._last_cov)
576       
577        # Now scale the added plots too
578        for plot in self._added_plots:
579            _max = 0
580            for y in self._added_plots[plot].y:
581                if y>_max: 
582                    _max = y
583            y = self._added_plots[plot].y/_max
584           
585            new_plot = Theory1D(self._added_plots[plot].x, y)
586            new_plot.name = self._added_plots[plot].name
587            new_plot.xaxis("\\rm{r}", 'A')
588            new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
589           
590            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, update=True,
591                                                   title=self._added_plots[plot].name))       
592       
593       
594    def _on_add_data(self, evt):
595        """
596            Add a data curve to the plot
597            WARNING: this will be removed once guiframe.plotting has its full functionality
598        """
599        path = self.choose_file()
600        if path==None:
601            return
602       
603        #x, y, err = self.parent.load_ascii_1D(path)
604        # Use data loader to load file
605        try:
606            dataread = Loader().load(path)
607            x = None
608            y = None
609            err = None
610            if dataread.__class__.__name__ == 'Data1D':
611                x = dataread.x
612                y = dataread.y
613                err = dataread.dy
614            else:
615                wx.PostEvent(self.parent, StatusEvent(status="This tool can only read 1D data"))
616                return
617        except:
618            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
619            return
620       
621        filename = os.path.basename(path)
622       
623        #new_plot = Data1D(x, y, dy=err)
624        new_plot = Theory1D(x, y)
625        new_plot.name = filename
626        new_plot.xaxis("\\rm{r}", 'A')
627        new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
628           
629        # Store a ref to the plottable for later use
630        self._added_plots[filename] = new_plot
631        self._default_Iq[filename]  = numpy.copy(y)
632       
633        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=filename))
634       
635       
636
637    def start_thread(self):
638        from pr_thread import CalcPr
639        from copy import deepcopy
640       
641        # If a thread is already started, stop it
642        if self.calc_thread != None and self.calc_thread.isrunning():
643            self.calc_thread.stop()
644               
645        pr = self.pr.clone()
646        self.calc_thread = CalcPr(pr, self.nfunc, error_func=self._thread_error, completefn=self._completed, updatefn=None)
647        self.calc_thread.queue()
648        self.calc_thread.ready(2.5)
649   
650    def _thread_error(self, error):
651        wx.PostEvent(self.parent, StatusEvent(status=error))
652   
653    def _estimate_completed(self, alpha, message, elapsed):
654        """
655            Parameter estimation completed,
656            display the results to the user
657            @param alpha: estimated best alpha
658            @param elapsed: computation time
659        """
660        # Save useful info
661        self.elapsed = elapsed
662        self.control_panel.alpha_estimate = alpha
663        if not message==None:
664            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
665           
666        self.perform_estimateNT()
667   
668
669   
670    def _estimateNT_completed(self, nterms, alpha, message, elapsed):
671        """
672            Parameter estimation completed,
673            display the results to the user
674            @param alpha: estimated best alpha
675            @param nterms: estimated number of terms
676            @param elapsed: computation time
677        """
678        # Save useful info
679        self.elapsed = elapsed
680        self.control_panel.nterms_estimate = nterms
681        self.control_panel.alpha_estimate = alpha
682        if not message==None:
683            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
684   
685    def _completed(self, out, cov, pr, elapsed):
686        """
687            Method called with the results when the inversion
688            is done
689           
690            @param out: output coefficient for the base functions
691            @param cov: covariance matrix
692            @param pr: Invertor instance
693            @param elapsed: time spent computing
694        """
695        from copy import deepcopy
696        # Save useful info
697        self.elapsed = elapsed
698        # Keep a copy of the last result
699        self._last_pr  = pr.clone()
700        self._last_out = out
701        self._last_cov = cov
702       
703        # Save Pr invertor
704        self.pr = pr
705       
706        #message = "Computation completed in %g seconds [chi2=%g]" % (elapsed, pr.chi2)
707        #wx.PostEvent(self.parent, StatusEvent(status=message))
708
709        cov = numpy.ascontiguousarray(cov)
710
711        # Show result on control panel
712        self.control_panel.chi2 = pr.chi2
713        self.control_panel.elapsed = elapsed
714        self.control_panel.oscillation = pr.oscillations(out)
715        #print "OSCILL", pr.oscillations(out)
716        #print "PEAKS:", pr.get_peaks(out)
717        self.control_panel.positive = pr.get_positive(out)
718        self.control_panel.pos_err  = pr.get_pos_err(out, cov)
719        self.control_panel.rg = pr.rg(out)
720        self.control_panel.iq0 = pr.iq0(out)
721        self.control_panel.bck = pr.background
722       
723        if False:
724            for i in range(len(out)):
725                try:
726                    print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
727                except: 
728                    print sys.exc_value
729                    print "%d: %g +- ?" % (i, out[i])       
730       
731            # Make a plot of I(q) data
732            new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
733            new_plot.name = IQ_DATA_LABEL
734            new_plot.xaxis("\\rm{Q}", 'A^{-1}')
735            new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
736            #new_plot.group_id = "test group"
737            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
738               
739        # Show I(q) fit
740        self.show_iq(out, self.pr)
741       
742        # Show P(r) fit
743        x_values, x_range = self.show_pr(out, self.pr, cov) 
744       
745        # Popup result panel
746        #result_panel = InversionResults(self.parent, -1, style=wx.RAISED_BORDER)
747       
748    def show_data(self, path=None, reset=False):
749        """
750            Show data read from a file
751            @param path: file path
752            @param reset: if True all other plottables will be cleared
753        """
754       
755       
756        if path is not None:
757            try:
758                pr = self._create_file_pr(path)
759                self.pr = pr
760            except: 
761                wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
762                return
763             
764        # Make a plot of I(q) data
765        if self.pr.err==None:
766            new_plot = Theory1D(self.pr.x, self.pr.y)
767        else:
768            new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
769        new_plot.name = IQ_DATA_LABEL
770        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
771        new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
772        new_plot.interactive = True
773        #new_plot.group_id = "test group"
774        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
775       
776        # Get Q range
777        self.control_panel.q_min = self.pr.x.min()
778        self.control_panel.q_max = self.pr.x.max()
779           
780
781       
782    def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None, 
783                             bck=False, height=0, width=0):
784        self.alpha = alpha
785        self.nfunc = nfunc
786        self.max_length = d_max
787        self.q_min = q_min
788        self.q_max = q_max
789        self.has_bck = bck
790        self.slit_height = height
791        self.slit_width  = width
792       
793        try:
794            pr = self._create_plot_pr()
795            if not pr==None:
796                self.pr = pr
797                self.perform_inversion()
798        except:
799            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
800
801    def estimate_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None, 
802                                bck=False, height=0, width=0):
803        self.alpha = alpha
804        self.nfunc = nfunc
805        self.max_length = d_max
806        self.q_min = q_min
807        self.q_max = q_max
808        self.has_bck = bck
809        self.slit_height = height
810        self.slit_width  = width
811       
812        try:
813            pr = self._create_plot_pr()
814            if not pr==None:
815                self.pr = pr
816                self.perform_estimate()
817        except:
818            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))           
819
820    def _create_plot_pr(self, estimate=False):
821        """
822            Create and prepare invertor instance from
823            a plottable data set.
824            @param path: path of the file to read in
825        """
826        # Get the data from the chosen data set and perform inversion
827        pr = Invertor()
828        pr.d_max = self.max_length
829        pr.alpha = self.alpha
830        pr.q_min = self.q_min
831        pr.q_max = self.q_max
832        pr.x = self.current_plottable.x
833        pr.y = self.current_plottable.y
834        pr.has_bck = self.has_bck
835        pr.slit_height = self.slit_height
836        pr.slit_width = self.slit_width
837       
838        # Fill in errors if none were provided
839        if self.current_plottable.dy == None:
840            print "no error", self.current_plottable.name
841            y = numpy.zeros(len(pr.y))
842            for i in range(len(pr.y)):
843                y[i] = math.sqrt(pr.y[i])
844            pr.err = y
845        else:
846            pr.err = self.current_plottable.dy
847       
848        #self.pr = pr
849        return pr
850
851         
852    def setup_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None, 
853                             bck=False, height=0, width=0):
854        self.alpha = alpha
855        self.nfunc = nfunc
856        self.max_length = d_max
857        self.q_min = q_min
858        self.q_max = q_max
859        self.has_bck = bck
860        self.slit_height = height
861        self.slit_width  = width
862       
863        try:
864            pr = self._create_file_pr(path)
865            if not pr==None:
866                self.pr = pr
867                self.perform_inversion()
868        except:
869            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
870         
871    def estimate_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None, 
872                                bck=False, height=0, width=0):
873        self.alpha = alpha
874        self.nfunc = nfunc
875        self.max_length = d_max
876        self.q_min = q_min
877        self.q_max = q_max
878        self.has_bck = bck
879        self.slit_height = height
880        self.slit_width  = width
881       
882        try:
883            pr = self._create_file_pr(path)
884            if not pr==None:
885                self.pr = pr
886                self.perform_estimate()
887        except:
888            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
889               
890         
891    def _create_file_pr(self, path):
892        """
893            Create and prepare invertor instance from
894            a file data set.
895            @param path: path of the file to read in
896        """
897        # Load data
898        if os.path.isfile(path):
899           
900            if self._current_file_data is not None \
901                and self._current_file_data.path==path:
902                # Protect against corrupted data from
903                # previous failed load attempt
904                if self._current_file_data.x is None:
905                    return None
906                x = self._current_file_data.x
907                y = self._current_file_data.y
908                err = self._current_file_data.err
909            else:
910                x, y, err = self.load(path)
911           
912            # If we have not errors, add statistical errors
913            if err==None and y is not None:
914                err = numpy.zeros(len(y))
915                scale = None
916                min_err = 0.0
917                for i in range(len(y)):
918                    # Scale the error so that we can fit over several decades of Q
919                    if scale==None:
920                        scale = 0.05*math.sqrt(y[i])
921                        min_err = 0.01*y[i]
922                    err[i] = scale*math.sqrt( math.fabs(y[i]) ) + min_err
923                message = "The loaded file had no error bars, statistical errors are assumed."
924                wx.PostEvent(self.parent, StatusEvent(status=message))
925            else:
926                wx.PostEvent(self.parent, StatusEvent(status=''))
927           
928            try:
929                # Get the data from the chosen data set and perform inversion
930                pr = Invertor()
931                pr.d_max = self.max_length
932                pr.alpha = self.alpha
933                pr.q_min = self.q_min
934                pr.q_max = self.q_max
935                pr.x = x
936                pr.y = y
937                pr.err = err
938                pr.has_bck = self.has_bck
939                pr.slit_height = self.slit_height
940                pr.slit_width = self.slit_width
941                return pr
942            except:
943                wx.PostEvent(self.parent, StatusEvent(status="Problem reading data: %s" % sys.exc_value))
944        return None
945       
946    def perform_estimate(self):
947        from pr_thread import EstimatePr
948        from copy import deepcopy
949       
950        # If a thread is already started, stop it
951        if self.estimation_thread != None and self.estimation_thread.isrunning():
952            self.estimation_thread.stop()
953               
954        pr = self.pr.clone()
955        self.estimation_thread = EstimatePr(pr, self.nfunc, error_func=self._thread_error, 
956                                            completefn = self._estimate_completed, 
957                                            updatefn   = None)
958        self.estimation_thread.queue()
959        self.estimation_thread.ready(2.5)
960   
961    def perform_estimateNT(self):
962        from pr_thread import EstimateNT
963        from copy import deepcopy
964       
965        # If a thread is already started, stop it
966        if self.estimation_thread != None and self.estimation_thread.isrunning():
967            self.estimation_thread.stop()
968               
969        pr = self.pr.clone()
970        # Skip the slit settings for the estimation
971        # It slows down the application and it doesn't change the estimates
972        pr.slit_height = 0.0
973        pr.slit_width  = 0.0
974        self.estimation_thread = EstimateNT(pr, self.nfunc, error_func=self._thread_error, 
975                                            completefn = self._estimateNT_completed, 
976                                            updatefn   = None)
977        self.estimation_thread.queue()
978        self.estimation_thread.ready(2.5)
979       
980    def perform_inversion(self):
981       
982        # Time estimate
983        #estimated = self.elapsed*self.nfunc**2
984        #message = "Computation time may take up to %g seconds" % self.elapsed
985        #wx.PostEvent(self.parent, StatusEvent(status=message))
986       
987        # Start inversion thread
988        self.start_thread()
989        return
990       
991        out, cov = self.pr.lstsq(self.nfunc)
992       
993        # Save useful info
994        self.elapsed = self.pr.elapsed
995       
996        for i in range(len(out)):
997            try:
998                print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
999            except: 
1000                print "%d: %g +- ?" % (i, out[i])       
1001       
1002       
1003       
1004        # Make a plot of I(q) data
1005        new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1006        new_plot.name = "I_{obs}(q)"
1007        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1008        new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1009        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1010               
1011        # Show I(q) fit
1012        self.show_iq(out, self.pr)
1013       
1014        # Show P(r) fit
1015        x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1016       
1017       
1018         
1019    def _on_context_inversion(self, event):
1020        panel = event.GetEventObject()
1021
1022        from inversion_panel import InversionDlg
1023       
1024        # If we have more than one displayed plot, make the user choose
1025        if len(panel.plots)>1 and panel.graph.selected_plottable in panel.plots:
1026            dataset = panel.graph.selected_plottable
1027            if False:
1028                dialog = InversionDlg(None, -1, "P(r) Inversion", panel.plots, pars=False)
1029                dialog.set_content(self.last_data, self.nfunc, self.alpha, self.max_length)
1030                if dialog.ShowModal() == wx.ID_OK:
1031                    dataset = dialog.get_content()
1032                    dialog.Destroy()
1033                else:
1034                    dialog.Destroy()
1035                    return
1036        elif len(panel.plots)==1:
1037            dataset = panel.plots.keys()[0]
1038        else:
1039            print "Error: No data is available"
1040            return
1041       
1042        # Store a reference to the current plottable
1043        # If we have a suggested value, use it.
1044        try:
1045            estimate = float(self.control_panel.alpha_estimate)
1046            self.control_panel.alpha = estimate
1047        except:
1048            self.control_panel.alpha = self.alpha
1049            print "No estimate yet"
1050            pass
1051        try:
1052            estimate = int(self.control_panel.nterms_estimate)
1053            self.control_panel.nfunc = estimate
1054        except:
1055            self.control_panel.nfunc = self.nfunc
1056            print "No estimate yet"
1057            pass
1058       
1059        self.current_plottable = panel.plots[dataset]
1060        self.control_panel.plotname = dataset
1061        #self.control_panel.nfunc = self.nfunc
1062        self.control_panel.d_max = self.max_length
1063        self.parent.set_perspective(self.perspective)
1064        self.control_panel._on_invert(None)
1065           
1066    def get_panels(self, parent):
1067        """
1068            Create and return a list of panel objects
1069        """
1070        from inversion_panel import InversionControl
1071       
1072        self.parent = parent
1073        self.control_panel = InversionControl(self.parent, -1, 
1074                                              style=wx.RAISED_BORDER,
1075                                              standalone=self.standalone)
1076        self.control_panel.set_manager(self)
1077        self.control_panel.nfunc = self.nfunc
1078        self.control_panel.d_max = self.max_length
1079        self.control_panel.alpha = self.alpha
1080       
1081        self.perspective = []
1082        self.perspective.append(self.control_panel.window_name)
1083       
1084        self.parent.Bind(EVT_PR_FILE, self._on_new_file)
1085       
1086        return [self.control_panel]
1087   
1088    def _on_new_file(self, evt):
1089        """
1090            Called when the application manager posted an
1091            EVT_PR_FILE event. Just prompt the control
1092            panel to load a new data file.
1093        """
1094        self.control_panel._change_file(None)
1095   
1096    def get_perspective(self):
1097        """
1098            Get the list of panel names for this perspective
1099        """
1100        return self.perspective
1101   
1102    def on_perspective(self, event):
1103        """
1104            Call back function for the perspective menu item.
1105            We notify the parent window that the perspective
1106            has changed.
1107        """
1108        self.parent.set_perspective(self.perspective)
1109   
1110    def post_init(self):
1111        """
1112            Post initialization call back to close the loose ends
1113            [Somehow openGL needs this call]
1114        """
1115        self.parent.set_perspective(self.perspective)
1116 
1117if __name__ == "__main__":
1118    i = Plugin()
1119    print i.perform_estimateNT()
1120   
1121   
1122   
1123   
Note: See TracBrowser for help on using the repository browser.