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

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 f8b79d6 was 0ab485b, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

prview: check that the data being sent to the prview writer is of the right type.
  • Property mode set to 100644
File size: 45.5 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            #new_plot.group_id = "test group"
329            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)"))
330       
331       
332    def _on_pr_npts(self, evt):
333        """
334            Redisplay P(r) with a different number of points
335        """   
336        from inversion_panel import PrDistDialog
337        dialog = PrDistDialog(None, -1)
338        dialog.set_content(self._pr_npts)
339        if dialog.ShowModal() == wx.ID_OK:
340            self._pr_npts= dialog.get_content()
341            dialog.Destroy()
342            self.show_pr(self._last_out, self._last_pr, self._last_cov)
343        else:
344            dialog.Destroy()
345       
346       
347    def show_pr(self, out, pr, cov=None):
348        import numpy
349        import pylab
350        import math
351        from sans.guicomm.events import NewPlotEvent           
352        from danse.common.plottools import Data1D, Theory1D
353       
354        # Show P(r)
355        x = pylab.arange(0.0, pr.d_max, pr.d_max/self._pr_npts)
356   
357        y = numpy.zeros(len(x))
358        dy = numpy.zeros(len(x))
359        y_true = numpy.zeros(len(x))
360
361        sum = 0.0
362        pmax = 0.0
363        cov2 = numpy.ascontiguousarray(cov)
364       
365        for i in range(len(x)):
366            if cov2==None:
367                value = pr.pr(out, x[i])
368            else:
369                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
370            sum += value*pr.d_max/len(x)
371           
372            # keep track of the maximum P(r) value
373            if value>pmax:
374                pmax = value
375               
376            y[i] = value
377               
378        if self._normalize_output==True:
379            y = y/sum
380            dy = dy/sum
381        elif self._scale_output_unity==True:
382            y = y/pmax
383            dy = dy/pmax
384       
385        if cov2==None:
386            new_plot = Theory1D(x, y)
387        else:
388            new_plot = Data1D(x, y, dy=dy)
389        new_plot.name = PR_FIT_LABEL
390        new_plot.xaxis("\\rm{r}", 'A')
391        new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
392        # Make sure that the plot is linear
393        new_plot.xtransform="x"
394        new_plot.ytransform="y"                 
395        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
396       
397        return x, pr.d_max
398       
399       
400    def choose_file(self, path=None):
401        """
402       
403        """
404        #TODO: this should be in a common module
405        return self.parent.choose_file(path=path)
406               
407               
408    def load(self, path):
409        """
410            Load data. This will eventually be replaced
411            by our standard DataLoader class.
412        """
413        class FileData:
414            x = None
415            y = None
416            err = None
417            path = None
418           
419            def __init__(self, path):
420                self.path = path
421               
422        self._current_file_data = FileData(path)
423       
424        # Use data loader to load file
425        dataread = Loader().load(path)
426       
427        # Notify the user if we could not read the file
428        if dataread is None:
429            raise RuntimeError, "Invalid data"
430           
431        x = None
432        y = None
433        err = None
434        if dataread.__class__.__name__ == 'Data1D':
435            x = dataread.x
436            y = dataread.y
437            err = dataread.dy
438        else:
439            if isinstance(dataread, list) and len(dataread)>0:
440                x = dataread[0].x
441                y = dataread[0].y
442                err = dataread[0].dy
443                msg = "PrView only allows a single data set at a time. "
444                msg += "Only the first data set was loaded." 
445                wx.PostEvent(self.parent, StatusEvent(status=msg))
446            else:
447                if dataread is None:
448                    return x, y, err
449                raise RuntimeError, "This tool can only read 1D data"
450       
451        self._current_file_data.x = x
452        self._current_file_data.y = y
453        self._current_file_data.err = err
454        return x, y, err
455               
456    def load_columns(self, path = "sphere_60_q0_2.txt"):
457        """
458            Load 2- or 3- column ascii
459        """
460        import numpy, math, sys
461        # Read the data from the data file
462        data_x   = numpy.zeros(0)
463        data_y   = numpy.zeros(0)
464        data_err = numpy.zeros(0)
465        scale    = None
466        min_err  = 0.0
467        if not path == None:
468            input_f = open(path,'r')
469            buff    = input_f.read()
470            lines   = buff.split('\n')
471            for line in lines:
472                try:
473                    toks = line.split()
474                    x = float(toks[0])
475                    y = float(toks[1])
476                    if len(toks)>2:
477                        err = float(toks[2])
478                    else:
479                        if scale==None:
480                            scale = 0.05*math.sqrt(y)
481                            #scale = 0.05/math.sqrt(y)
482                            min_err = 0.01*y
483                        err = scale*math.sqrt(y)+min_err
484                        #err = 0
485                       
486                    data_x = numpy.append(data_x, x)
487                    data_y = numpy.append(data_y, y)
488                    data_err = numpy.append(data_err, err)
489                except:
490                    pass
491                   
492        if not scale==None:
493            message = "The loaded file had no error bars, statistical errors are assumed."
494            wx.PostEvent(self.parent, StatusEvent(status=message))
495        else:
496            wx.PostEvent(self.parent, StatusEvent(status=''))
497                       
498        return data_x, data_y, data_err     
499       
500    def load_abs(self, path):
501        """
502            Load an IGOR .ABS reduced file
503            @param path: file path
504            @return: x, y, err vectors
505        """
506        import numpy, math, sys
507        # Read the data from the data file
508        data_x   = numpy.zeros(0)
509        data_y   = numpy.zeros(0)
510        data_err = numpy.zeros(0)
511        scale    = None
512        min_err  = 0.0
513       
514        data_started = False
515        if not path == None:
516            input_f = open(path,'r')
517            buff    = input_f.read()
518            lines   = buff.split('\n')
519            for line in lines:
520                if data_started==True:
521                    try:
522                        toks = line.split()
523                        x = float(toks[0])
524                        y = float(toks[1])
525                        if len(toks)>2:
526                            err = float(toks[2])
527                        else:
528                            if scale==None:
529                                scale = 0.05*math.sqrt(y)
530                                #scale = 0.05/math.sqrt(y)
531                                min_err = 0.01*y
532                            err = scale*math.sqrt(y)+min_err
533                            #err = 0
534                           
535                        data_x = numpy.append(data_x, x)
536                        data_y = numpy.append(data_y, y)
537                        data_err = numpy.append(data_err, err)
538                    except:
539                        pass
540                elif line.find("The 6 columns")>=0:
541                    data_started = True     
542                   
543        if not scale==None:
544            message = "The loaded file had no error bars, statistical errors are assumed."
545            wx.PostEvent(self.parent, StatusEvent(status=message))
546        else:
547            wx.PostEvent(self.parent, StatusEvent(status=''))
548                       
549        return data_x, data_y, data_err     
550       
551       
552       
553    def pr_theory(self, r, R):
554        """
555           
556        """
557        if r<=2*R:
558            return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
559        else:
560            return 0.0
561
562    def get_context_menu(self, graph=None):
563        """
564            Get the context menu items available for P(r)
565            @param graph: the Graph object to which we attach the context menu
566            @return: a list of menu items with call-back function
567        """
568        # Look whether this Graph contains P(r) data
569        #if graph.selected_plottable==IQ_DATA_LABEL:
570        for item in graph.plottables:
571            if item.name==PR_FIT_LABEL:
572                m_list = [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
573                       ["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
574
575                if self._scale_output_unity==True or self._normalize_output==True:
576                    m_list.append(["Disable P(r) scaling", 
577                                   "Let the output P(r) keep the scale of the data", 
578                                   self._on_disable_scaling])
579               
580                if self._scale_output_unity==False:
581                    m_list.append(["Scale P_max(r) to unity", 
582                                   "Scale P(r) so that its maximum is 1", 
583                                   self._on_scale_unity])
584                   
585                if self._normalize_output==False:
586                    m_list.append(["Normalize P(r) to unity", 
587                                   "Normalize the integral of P(r) to 1", 
588                                   self._on_normalize])
589                   
590                return m_list
591                #return [["Add P(r) data", "Load a data file and display it on this plot", self._on_add_data],
592                #       ["Change number of P(r) points", "Change the number of points on the P(r) output", self._on_pr_npts]]
593
594            elif item.name==graph.selected_plottable:
595                    #TODO: we might want to check that the units are consistent with I(q)
596                    #      before allowing this menu item
597                if issubclass(item.__class__, DataLoader.data_info.Data1D):
598                    return [["Compute P(r)", "Compute P(r) from distribution", self._on_context_inversion]]     
599               
600        return []
601
602    def _on_disable_scaling(self, evt):
603        """
604            Disable P(r) scaling
605            @param evt: Menu event
606        """
607        self._normalize_output = False
608        self._scale_output_unity = False
609        self.show_pr(self._last_out, self._last_pr, self._last_cov)
610       
611        # Now replot the original added data
612        for plot in self._added_plots:
613            self._added_plots[plot].y = numpy.copy(self._default_Iq[plot])
614            wx.PostEvent(self.parent, NewPlotEvent(plot=self._added_plots[plot], 
615                                                   title=self._added_plots[plot].name,
616                                                   update=True))       
617       
618        # Need the update flag in the NewPlotEvent to protect against
619        # the plot no longer being there...
620       
621    def _on_normalize(self, evt):
622        """
623            Normalize the area under the P(r) curve to 1.
624            This operation is done for all displayed plots.
625           
626            @param evt: Menu event
627        """
628        self._normalize_output = True
629        self._scale_output_unity = False
630           
631        self.show_pr(self._last_out, self._last_pr, self._last_cov)
632       
633        # Now scale the added plots too
634        for plot in self._added_plots:
635            sum = numpy.sum(self._added_plots[plot].y)
636            npts = len(self._added_plots[plot].x)
637            sum *= self._added_plots[plot].x[npts-1]/npts
638            y = self._added_plots[plot].y/sum
639           
640            new_plot = Theory1D(self._added_plots[plot].x, y)
641            new_plot.name = self._added_plots[plot].name
642            new_plot.xaxis("\\rm{r}", 'A')
643            new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
644           
645            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, update=True,
646                                                   title=self._added_plots[plot].name))
647               
648       
649       
650    def _on_scale_unity(self, evt):
651        """
652            Scale the maximum P(r) value on each displayed plot to 1.
653           
654            @param evt: Menu event
655        """
656        self._scale_output_unity = True
657        self._normalize_output = False
658           
659        self.show_pr(self._last_out, self._last_pr, self._last_cov)
660       
661        # Now scale the added plots too
662        for plot in self._added_plots:
663            _max = 0
664            for y in self._added_plots[plot].y:
665                if y>_max: 
666                    _max = y
667            y = self._added_plots[plot].y/_max
668           
669            new_plot = Theory1D(self._added_plots[plot].x, y)
670            new_plot.name = self._added_plots[plot].name
671            new_plot.xaxis("\\rm{r}", 'A')
672            new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
673           
674            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, update=True,
675                                                   title=self._added_plots[plot].name))       
676       
677       
678    def _on_add_data(self, evt):
679        """
680            Add a data curve to the plot
681            WARNING: this will be removed once guiframe.plotting has its full functionality
682        """
683        path = self.choose_file()
684        if path==None:
685            return
686       
687        #x, y, err = self.parent.load_ascii_1D(path)
688        # Use data loader to load file
689        try:
690            dataread = Loader().load(path)
691            x = None
692            y = None
693            err = None
694            if dataread.__class__.__name__ == 'Data1D':
695                x = dataread.x
696                y = dataread.y
697                err = dataread.dy
698            else:
699                if isinstance(dataread, list) and len(dataread)>0:
700                    x = dataread[0].x
701                    y = dataread[0].y
702                    err = dataread[0].dy
703                    msg = "PrView only allows a single data set at a time. "
704                    msg += "Only the first data set was loaded." 
705                    wx.PostEvent(self.parent, StatusEvent(status=msg))
706                else:
707                    wx.PostEvent(self.parent, StatusEvent(status="This tool can only read 1D data"))
708                    return
709           
710        except:
711            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
712            return
713       
714        filename = os.path.basename(path)
715       
716        #new_plot = Data1D(x, y, dy=err)
717        new_plot = Theory1D(x, y)
718        new_plot.name = filename
719        new_plot.xaxis("\\rm{r}", 'A')
720        new_plot.yaxis("\\rm{P(r)} ","cm^{-3}")
721           
722        # Store a ref to the plottable for later use
723        self._added_plots[filename] = new_plot
724        self._default_Iq[filename]  = numpy.copy(y)
725       
726        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=filename))
727       
728       
729
730    def start_thread(self):
731        from pr_thread import CalcPr
732        from copy import deepcopy
733       
734        # If a thread is already started, stop it
735        if self.calc_thread != None and self.calc_thread.isrunning():
736            self.calc_thread.stop()
737               
738        pr = self.pr.clone()
739        self.calc_thread = CalcPr(pr, self.nfunc, error_func=self._thread_error, completefn=self._completed, updatefn=None)
740        self.calc_thread.queue()
741        self.calc_thread.ready(2.5)
742   
743    def _thread_error(self, error):
744        wx.PostEvent(self.parent, StatusEvent(status=error))
745   
746    def _estimate_completed(self, alpha, message, elapsed):
747        """
748            Parameter estimation completed,
749            display the results to the user
750            @param alpha: estimated best alpha
751            @param elapsed: computation time
752        """
753        # Save useful info
754        self.elapsed = elapsed
755        self.control_panel.alpha_estimate = alpha
756        if not message==None:
757            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
758           
759        self.perform_estimateNT()
760   
761
762   
763    def _estimateNT_completed(self, nterms, alpha, message, elapsed):
764        """
765            Parameter estimation completed,
766            display the results to the user
767            @param alpha: estimated best alpha
768            @param nterms: estimated number of terms
769            @param elapsed: computation time
770        """
771        # Save useful info
772        self.elapsed = elapsed
773        self.control_panel.nterms_estimate = nterms
774        self.control_panel.alpha_estimate = alpha
775        if not message==None:
776            wx.PostEvent(self.parent, StatusEvent(status=str(message)))
777   
778    def _completed(self, out, cov, pr, elapsed):
779        """
780            Method called with the results when the inversion
781            is done
782           
783            @param out: output coefficient for the base functions
784            @param cov: covariance matrix
785            @param pr: Invertor instance
786            @param elapsed: time spent computing
787        """
788        from copy import deepcopy
789        # Save useful info
790        self.elapsed = elapsed
791        # Keep a copy of the last result
792        self._last_pr  = pr.clone()
793        self._last_out = out
794        self._last_cov = cov
795       
796        # Save Pr invertor
797        self.pr = pr
798       
799        #message = "Computation completed in %g seconds [chi2=%g]" % (elapsed, pr.chi2)
800        #wx.PostEvent(self.parent, StatusEvent(status=message))
801
802        cov = numpy.ascontiguousarray(cov)
803
804        # Show result on control panel
805        self.control_panel.chi2 = pr.chi2
806        self.control_panel.elapsed = elapsed
807        self.control_panel.oscillation = pr.oscillations(out)
808        #print "OSCILL", pr.oscillations(out)
809        #print "PEAKS:", pr.get_peaks(out)
810        self.control_panel.positive = pr.get_positive(out)
811        self.control_panel.pos_err  = pr.get_pos_err(out, cov)
812        self.control_panel.rg = pr.rg(out)
813        self.control_panel.iq0 = pr.iq0(out)
814        self.control_panel.bck = pr.background
815       
816        if False:
817            for i in range(len(out)):
818                try:
819                    print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
820                except: 
821                    print sys.exc_value
822                    print "%d: %g +- ?" % (i, out[i])       
823       
824            # Make a plot of I(q) data
825            new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
826            new_plot.name = IQ_DATA_LABEL
827            new_plot.xaxis("\\rm{Q}", 'A^{-1}')
828            new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
829            #new_plot.group_id = "test group"
830            wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
831               
832        # Show I(q) fit
833        self.show_iq(out, self.pr)
834       
835        # Show P(r) fit
836        x_values, x_range = self.show_pr(out, self.pr, cov) 
837       
838        # Popup result panel
839        #result_panel = InversionResults(self.parent, -1, style=wx.RAISED_BORDER)
840       
841    def show_data(self, path=None, reset=False):
842        """
843            Show data read from a file
844            @param path: file path
845            @param reset: if True all other plottables will be cleared
846        """
847        if path is not None:
848            try:
849                pr = self._create_file_pr(path)
850            except:
851                status="Problem reading data: %s" % sys.exc_value
852                wx.PostEvent(self.parent, StatusEvent(status=status))
853                raise RuntimeError, status
854               
855            # If the file contains nothing, just return
856            if pr is None:
857                raise RuntimeError, "Loaded data is invalid"
858           
859            self.pr = pr
860             
861        # Make a plot of I(q) data
862        if self.pr.err==None:
863            new_plot = Theory1D(self.pr.x, self.pr.y)
864        else:
865            new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
866        new_plot.name = IQ_DATA_LABEL
867        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
868        new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
869        new_plot.interactive = True
870        #new_plot.group_id = "test group"
871        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
872       
873        self.current_plottable = new_plot
874        self.current_plottable.group_id = IQ_DATA_LABEL
875       
876       
877        # Get Q range
878        self.control_panel.q_min = self.pr.x.min()
879        self.control_panel.q_max = self.pr.x.max()
880           
881    def save_data(self, filepath, prstate=None):
882        """
883            Save data in provided state object.
884            TODO: move the state code away from inversion_panel and move it here.
885                    Then remove the "prstate" input and make this method private.
886                   
887            @param filepath: path of file to write to
888            @param prstate: P(r) inversion state
889        """
890        #TODO: do we need this or can we use DataLoader.loader.save directly?
891       
892        # Add output data and coefficients to state
893        prstate.coefficients = self._last_out
894        prstate.covariance = self._last_cov
895       
896        # Write the output to file
897        # First, check that the data is of the right type
898        if issubclass(self.current_plottable.__class__, DataLoader.data_info.Data1D):
899            self.state_reader.write(filepath, self.current_plottable, prstate)
900        else:
901            raise RuntimeError, "pr.save_data: the data being saved is not a DataLoader.data_info.Data1D object" 
902       
903       
904    def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None, 
905                             bck=False, height=0, width=0):
906        self.alpha = alpha
907        self.nfunc = nfunc
908        self.max_length = d_max
909        self.q_min = q_min
910        self.q_max = q_max
911        self.has_bck = bck
912        self.slit_height = height
913        self.slit_width  = width
914       
915        try:
916            pr = self._create_plot_pr()
917            if not pr==None:
918                self.pr = pr
919                self.perform_inversion()
920        except:
921            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
922
923    def estimate_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None, 
924                                bck=False, height=0, width=0):
925        self.alpha = alpha
926        self.nfunc = nfunc
927        self.max_length = d_max
928        self.q_min = q_min
929        self.q_max = q_max
930        self.has_bck = bck
931        self.slit_height = height
932        self.slit_width  = width
933       
934        try:
935            pr = self._create_plot_pr()
936            if not pr==None:
937                self.pr = pr
938                self.perform_estimate()
939        except:
940            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))           
941
942    def _create_plot_pr(self, estimate=False):
943        """
944            Create and prepare invertor instance from
945            a plottable data set.
946            @param path: path of the file to read in
947        """
948        # Sanity check
949        if self.current_plottable is None:
950            msg = "Please load a valid data set before proceeding."
951            wx.PostEvent(self.parent, StatusEvent(status=msg)) 
952            return None   
953       
954        # Get the data from the chosen data set and perform inversion
955        pr = Invertor()
956        pr.d_max = self.max_length
957        pr.alpha = self.alpha
958        pr.q_min = self.q_min
959        pr.q_max = self.q_max
960        pr.x = self.current_plottable.x
961        pr.y = self.current_plottable.y
962        pr.has_bck = self.has_bck
963        pr.slit_height = self.slit_height
964        pr.slit_width = self.slit_width
965       
966        # Keep track of the plot window title to ensure that
967        # we can overlay the plots
968        if hasattr(self.current_plottable, "group_id"):
969            pr.info["plot_group_id"] = self.current_plottable.group_id
970       
971        # Fill in errors if none were provided
972        err = self.current_plottable.dy
973        all_zeros = True
974        if err == None:
975            err = numpy.zeros(len(pr.y)) 
976        else:   
977            for i in range(len(err)):
978                if err[i]>0:
979                    all_zeros = False
980       
981        if all_zeros:       
982            scale = None
983            min_err = 0.0
984            for i in range(len(pr.y)):
985                # Scale the error so that we can fit over several decades of Q
986                if scale==None:
987                    scale = 0.05*math.sqrt(pr.y[i])
988                    min_err = 0.01*pr.y[i]
989                err[i] = scale*math.sqrt( math.fabs(pr.y[i]) ) + min_err
990            message = "The loaded file had no error bars, statistical errors are assumed."
991            wx.PostEvent(self.parent, StatusEvent(status=message))
992
993        pr.err = err
994       
995        return pr
996
997         
998    def setup_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None, 
999                             bck=False, height=0, width=0):
1000        self.alpha = alpha
1001        self.nfunc = nfunc
1002        self.max_length = d_max
1003        self.q_min = q_min
1004        self.q_max = q_max
1005        self.has_bck = bck
1006        self.slit_height = height
1007        self.slit_width  = width
1008       
1009        try:
1010            pr = self._create_file_pr(path)
1011            if not pr==None:
1012                self.pr = pr
1013                self.perform_inversion()
1014        except:
1015            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1016         
1017    def estimate_file_inversion(self, alpha, nfunc, d_max, path, q_min=None, q_max=None, 
1018                                bck=False, height=0, width=0):
1019        self.alpha = alpha
1020        self.nfunc = nfunc
1021        self.max_length = d_max
1022        self.q_min = q_min
1023        self.q_max = q_max
1024        self.has_bck = bck
1025        self.slit_height = height
1026        self.slit_width  = width
1027       
1028        try:
1029            pr = self._create_file_pr(path)
1030            if not pr==None:
1031                self.pr = pr
1032                self.perform_estimate()
1033        except:
1034            wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
1035               
1036         
1037    def _create_file_pr(self, path):
1038        """
1039            Create and prepare invertor instance from
1040            a file data set.
1041            @param path: path of the file to read in
1042        """
1043        # Load data
1044        if os.path.isfile(path):
1045           
1046            if self._current_file_data is not None \
1047                and self._current_file_data.path==path:
1048                # Protect against corrupted data from
1049                # previous failed load attempt
1050                if self._current_file_data.x is None:
1051                    return None
1052                x = self._current_file_data.x
1053                y = self._current_file_data.y
1054                err = self._current_file_data.err
1055               
1056                message = "The data from this file has already been loaded."
1057                wx.PostEvent(self.parent, StatusEvent(status=message))
1058            else:
1059                # Reset the status bar so that we don't get mixed up
1060                # with old messages.
1061                #TODO: refactor this into a proper status handling
1062                wx.PostEvent(self.parent, StatusEvent(status=''))
1063                try:
1064                    x, y, err = self.load(path)
1065                except:
1066                    load_error(sys.exc_value)
1067                    return None
1068               
1069                # If the file contains no data, just return
1070                if x is None or len(x)==0:
1071                    load_error("The loaded file contains no data")
1072                    return None
1073           
1074            # If we have not errors, add statistical errors
1075            if err==None and y is not None:
1076                err = numpy.zeros(len(y))
1077                scale = None
1078                min_err = 0.0
1079                for i in range(len(y)):
1080                    # Scale the error so that we can fit over several decades of Q
1081                    if scale==None:
1082                        scale = 0.05*math.sqrt(y[i])
1083                        min_err = 0.01*y[i]
1084                    err[i] = scale*math.sqrt( math.fabs(y[i]) ) + min_err
1085                message = "The loaded file had no error bars, statistical errors are assumed."
1086                wx.PostEvent(self.parent, StatusEvent(status=message))
1087           
1088            try:
1089                # Get the data from the chosen data set and perform inversion
1090                pr = Invertor()
1091                pr.d_max = self.max_length
1092                pr.alpha = self.alpha
1093                pr.q_min = self.q_min
1094                pr.q_max = self.q_max
1095                pr.x = x
1096                pr.y = y
1097                pr.err = err
1098                pr.has_bck = self.has_bck
1099                pr.slit_height = self.slit_height
1100                pr.slit_width = self.slit_width
1101                return pr
1102            except:
1103                load_error(sys.exc_value)
1104        return None
1105       
1106    def perform_estimate(self):
1107        from pr_thread import EstimatePr
1108        from copy import deepcopy
1109       
1110        # If a thread is already started, stop it
1111        if self.estimation_thread != None and self.estimation_thread.isrunning():
1112            self.estimation_thread.stop()
1113               
1114        pr = self.pr.clone()
1115        self.estimation_thread = EstimatePr(pr, self.nfunc, error_func=self._thread_error, 
1116                                            completefn = self._estimate_completed, 
1117                                            updatefn   = None)
1118        self.estimation_thread.queue()
1119        self.estimation_thread.ready(2.5)
1120   
1121    def perform_estimateNT(self):
1122        from pr_thread import EstimateNT
1123        from copy import deepcopy
1124       
1125        # If a thread is already started, stop it
1126        if self.estimation_thread != None and self.estimation_thread.isrunning():
1127            self.estimation_thread.stop()
1128               
1129        pr = self.pr.clone()
1130        # Skip the slit settings for the estimation
1131        # It slows down the application and it doesn't change the estimates
1132        pr.slit_height = 0.0
1133        pr.slit_width  = 0.0
1134        self.estimation_thread = EstimateNT(pr, self.nfunc, error_func=self._thread_error, 
1135                                            completefn = self._estimateNT_completed, 
1136                                            updatefn   = None)
1137        self.estimation_thread.queue()
1138        self.estimation_thread.ready(2.5)
1139       
1140    def perform_inversion(self):
1141       
1142        # Time estimate
1143        #estimated = self.elapsed*self.nfunc**2
1144        #message = "Computation time may take up to %g seconds" % self.elapsed
1145        #wx.PostEvent(self.parent, StatusEvent(status=message))
1146       
1147        # Start inversion thread
1148        self.start_thread()
1149        return
1150       
1151        out, cov = self.pr.lstsq(self.nfunc)
1152       
1153        # Save useful info
1154        self.elapsed = self.pr.elapsed
1155       
1156        for i in range(len(out)):
1157            try:
1158                print "%d: %g +- %g" % (i, out[i], math.sqrt(math.fabs(cov[i][i])))
1159            except: 
1160                print "%d: %g +- ?" % (i, out[i])       
1161       
1162       
1163       
1164        # Make a plot of I(q) data
1165        new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
1166        new_plot.name = "I_{obs}(q)"
1167        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
1168        new_plot.yaxis("\\rm{Intensity} ","cm^{-1}")
1169        wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="Iq"))
1170               
1171        # Show I(q) fit
1172        self.show_iq(out, self.pr)
1173       
1174        # Show P(r) fit
1175        x_values, x_range = self.show_pr(out, self.pr, cov=cov)
1176       
1177       
1178         
1179    def _on_context_inversion(self, event):
1180        panel = event.GetEventObject()
1181
1182        # If we have more than one displayed plot, make the user choose
1183        if len(panel.plots)>1 and panel.graph.selected_plottable in panel.plots:
1184            dataset = panel.graph.selected_plottable
1185        elif len(panel.plots)==1:
1186            dataset = panel.plots.keys()[0]
1187        else:
1188            print "Error: No data is available"
1189            return
1190       
1191        # Store a reference to the current plottable
1192        # If we have a suggested value, use it.
1193        try:
1194            estimate = float(self.control_panel.alpha_estimate)
1195            self.control_panel.alpha = estimate
1196        except:
1197            self.control_panel.alpha = self.alpha
1198            print "No estimate yet"
1199            pass
1200        try:
1201            estimate = int(self.control_panel.nterms_estimate)
1202            self.control_panel.nfunc = estimate
1203        except:
1204            self.control_panel.nfunc = self.nfunc
1205            print "No estimate yet"
1206            pass
1207       
1208        self.current_plottable = panel.plots[dataset]
1209        self.control_panel.plotname = dataset
1210        #self.control_panel.nfunc = self.nfunc
1211        self.control_panel.d_max = self.max_length
1212        self.parent.set_perspective(self.perspective)
1213        self.control_panel._on_invert(None)
1214           
1215    def get_panels(self, parent):
1216        """
1217            Create and return a list of panel objects
1218        """
1219        from inversion_panel import InversionControl
1220       
1221        self.parent = parent
1222        self.control_panel = InversionControl(self.parent, -1, 
1223                                              style=wx.RAISED_BORDER,
1224                                              standalone=self.standalone)
1225        self.control_panel.set_manager(self)
1226        self.control_panel.nfunc = self.nfunc
1227        self.control_panel.d_max = self.max_length
1228        self.control_panel.alpha = self.alpha
1229       
1230        self.perspective = []
1231        self.perspective.append(self.control_panel.window_name)
1232       
1233        self.parent.Bind(EVT_PR_FILE, self._on_new_file)
1234       
1235        return [self.control_panel]
1236   
1237    def _on_new_file(self, evt):
1238        """
1239            Called when the application manager posted an
1240            EVT_PR_FILE event. Just prompt the control
1241            panel to load a new data file.
1242        """
1243        self.control_panel._change_file(None)
1244   
1245    def get_perspective(self):
1246        """
1247            Get the list of panel names for this perspective
1248        """
1249        return self.perspective
1250   
1251    def on_perspective(self, event):
1252        """
1253            Call back function for the perspective menu item.
1254            We notify the parent window that the perspective
1255            has changed.
1256        """
1257        self.parent.set_perspective(self.perspective)
1258   
1259    def post_init(self):
1260        """
1261            Post initialization call back to close the loose ends
1262            [Somehow openGL needs this call]
1263        """
1264        if self.standalone==True:
1265            self.parent.set_perspective(self.perspective)
1266 
1267if __name__ == "__main__":
1268    i = Plugin()
1269    print i.perform_estimateNT()
1270   
1271   
1272   
1273   
Note: See TracBrowser for help on using the repository browser.