source: sasview/sansview/perspectives/pr/pr.py @ 111acd2

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 111acd2 was 4f20b27, checked in by Jae Cho <jhjcho@…>, 15 years ago

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