source: sasview/prview/perspectives/pr/pr.py @ 4e74e13

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 4e74e13 was f0ff8d65, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

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