source: sasview/src/sans/perspectives/pr/pr.py @ 5777106

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

Moving things around. Will definitely not build.

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