source: sasview/prview/perspectives/pr/pr.py @ 0b12abb5

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 0b12abb5 was 7116b6e0, checked in by Gervaise Alina <gervyh@…>, 15 years ago

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