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

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 ee6f84c was 6f1f129, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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