source: sasview/prview/perspectives/pr/pr.py @ 7f84e22

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 7f84e22 was 3e41f43, checked in by Gervaise Alina <gervyh@…>, 14 years ago

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