source: sasview/prview/perspectives/pr/pr.py @ 1c8015f

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 1c8015f was d2ee6f6, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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