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

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 cacbd7d was 75a7ece, checked in by Gervaise Alina <gervyh@…>, 14 years ago

allow selection of multiple file for prview at laoding time

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