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

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

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