source: sasview/prview/perspectives/pr/pr.py @ 3c44c66

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 3c44c66 was ed90edb, checked in by Mathieu Doucet <doucetm@…>, 14 years ago

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