source: sasview/src/sans/perspectives/pr/inversion_state.py @ 3600c12

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

Moving things around. Will definitely not build.

  • Property mode set to 100644
File size: 19.7 KB
Line 
1"""
2    Handling of P(r) inversion states
3"""
4################################################################################
5#This software was developed by the University of Tennessee as part of the
6#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
7#project funded by the US National Science Foundation.
8#
9#See the license text in license.txt
10#
11#copyright 2009, University of Tennessee
12################################################################################
13
14
15import time
16import os
17import sys
18import logging
19from lxml import etree
20from sans.guiframe.dataFitting import Data1D
21from sans.dataloader.readers.cansas_reader import Reader as CansasReader
22from sans.dataloader.readers.cansas_reader import get_content
23
24PRNODE_NAME = 'pr_inversion'
25CANSAS_NS = "cansas1d/1.0"
26
27# Translation of names between stored and object data
28## List of P(r) inversion inputs
29in_list =  [["nterms",       "nfunc"],
30           ["d_max",        "d_max"],
31           ["alpha",        "alpha"],
32           ["slit_width",   "width"],
33           ["slit_height",  "height"],
34           ["qmin",         "qmin"],
35           ["qmax",         "qmax"],
36           ["estimate_bck", "estimate_bck"]]                     
37
38## List of P(r) inversion outputs
39out_list = [["elapsed", "elapsed"],
40           ["rg",      "rg"],
41           ["iq0",     "iq0"],
42           ["bck",     "bck"],
43           ["chi2",    "chi2"],
44           ["osc",     "osc"],
45           ["pos",     "pos"],
46           ["pos_err", "pos_err"],
47           ["alpha_estimate", "alpha_estimate"],
48           ["nterms_estimate", "nterms_estimate"]]
49
50class InversionState(object):
51    """
52    Class to hold the state information of the InversionControl panel.
53    """
54    def __init__(self):
55        """
56        Default values
57        """
58        # Input
59        self.file  = None
60        self.estimate_bck = False
61        self.timestamp = time.time()
62       
63        # Inversion parameters
64        self.nfunc = None
65        self.d_max = None
66        self.alpha = None
67       
68        # Slit parameters
69        self.height = None
70        self.width  = None
71       
72        # Q range
73        self.qmin  = None
74        self.qmax  = None
75       
76        # Outputs
77        self.elapsed = None
78        self.rg    = None
79        self.iq0   = None
80        self.bck   = None
81        self.chi2  = None
82        self.osc   = None
83        self.pos   = None
84        self.pos_err = None
85       
86        # Estimates
87        self.alpha_estimate = None
88        self.nterms_estimate = None
89       
90        # Data
91        self.q       = None
92        self.iq_obs  = None
93        self.iq_calc = None
94       
95        # Coefficients
96        self.coefficients = None
97        self.covariance = None
98   
99    def __str__(self):
100        """
101        Pretty print
102       
103        :return: string representing the state
104       
105        """
106        state  = "File:         %s\n" % self.file
107        state += "Timestamp:    %s\n" % self.timestamp
108        state += "Estimate bck: %s\n" % str(self.estimate_bck)
109        state += "No. terms:    %s\n" % str(self.nfunc)
110        state += "D_max:        %s\n" % str(self.d_max)
111        state += "Alpha:        %s\n" % str(self.alpha)
112       
113        state += "Slit height:  %s\n" % str(self.height)
114        state += "Slit width:   %s\n" % str(self.width)
115       
116        state += "Qmin:         %s\n" % str(self.qmin)
117        state += "Qmax:         %s\n" % str(self.qmax)
118       
119        state += "\nEstimates:\n"
120        state += "  Alpha:      %s\n" % str(self.alpha_estimate)
121        state += "  Nterms:     %s\n" % str(self.nterms_estimate)
122       
123        state += "\nOutputs:\n"
124        state += "  Elapsed:    %s\n" % str(self.elapsed)
125        state += "  Rg:         %s\n" % str(self.rg)
126        state += "  I(q=0):     %s\n" % str(self.iq0)
127        state += "  Bck:        %s\n" % str(self.bck)
128        state += "  Chi^2:      %s\n" % str(self.chi2)
129        state += "  Oscillation:%s\n" % str(self.osc)
130        state += "  Positive:   %s\n" % str(self.pos)
131        state += "  1-sigma pos:%s\n" % str(self.pos_err)
132       
133        return state
134       
135    def toXML(self, file="pr_state.prv", doc=None, entry_node=None):
136        """
137        Writes the state of the InversionControl panel to file, as XML.
138       
139        Compatible with standalone writing, or appending to an
140        already existing XML document. In that case, the XML document
141        is required. An optional entry node in the XML document
142        may also be given.
143       
144        :param file: file to write to
145        :param doc: XML document object [optional]
146        :param entry_node: XML node within the XML document at which
147            we will append the data [optional]
148       
149        """
150        from xml.dom.minidom import getDOMImplementation
151
152        # Check whether we have to write a standalone XML file
153        if doc is None:
154            impl = getDOMImplementation()
155       
156            doc_type = impl.createDocumentType(PRNODE_NAME, "1.0", "1.0")     
157       
158            newdoc = impl.createDocument(None, PRNODE_NAME, doc_type)
159            top_element = newdoc.documentElement
160        else:
161            # We are appending to an existing document
162            newdoc = doc
163            top_element = newdoc.createElement(PRNODE_NAME)
164            if entry_node is None:
165                newdoc.documentElement.appendChild(top_element)
166            else:
167                entry_node.appendChild(top_element)
168           
169        attr = newdoc.createAttribute("version")
170        attr.nodeValue = '1.0'
171        top_element.setAttributeNode(attr)
172       
173        # File name
174        element = newdoc.createElement("filename")
175        if self.file is not None:
176            element.appendChild(newdoc.createTextNode(str(self.file)))
177        else:
178            element.appendChild(newdoc.createTextNode(str(file)))
179        top_element.appendChild(element)
180       
181        element = newdoc.createElement("timestamp")
182        element.appendChild(newdoc.createTextNode(time.ctime(self.timestamp)))
183        attr = newdoc.createAttribute("epoch")
184        attr.nodeValue = str(self.timestamp)
185        element.setAttributeNode(attr)
186        top_element.appendChild(element)
187       
188        # Inputs
189        inputs = newdoc.createElement("inputs")
190        top_element.appendChild(inputs)
191       
192        for item in in_list:
193            element = newdoc.createElement(item[0])
194            cmd = "element.appendChild(newdoc.createTextNode(str(self.%s)))"
195            exec  cmd % item[1]
196            inputs.appendChild(element)
197             
198        # Outputs
199        outputs = newdoc.createElement("outputs")
200        top_element.appendChild(outputs)
201       
202        for item in out_list:
203            element = newdoc.createElement(item[0])
204            cmd = "element.appendChild(newdoc.createTextNode(str(self.%s)))"
205            exec  cmd % item[1]
206            outputs.appendChild(element)
207                   
208        # Save output coefficients and its covariance matrix
209        element = newdoc.createElement("coefficients")
210        element.appendChild(newdoc.createTextNode(str(self.coefficients)))
211        outputs.appendChild(element)
212        element = newdoc.createElement("covariance")
213        element.appendChild(newdoc.createTextNode(str(self.covariance)))
214        outputs.appendChild(element)
215                   
216        # Save the file
217        if doc is None:
218            fd = open(file, 'w')
219            fd.write(newdoc.toprettyxml())
220            fd.close()
221            return None
222        else:
223            return newdoc.toprettyxml()
224
225    def fromXML(self, file=None, node=None):
226        """
227        Load a P(r) inversion state from a file
228       
229        :param file: .prv file
230        :param node: node of a XML document to read from
231       
232        """
233        if file is not None:
234            msg = "InversionState no longer supports non-CanSAS"
235            msg += " format for P(r) files"
236            raise RuntimeError, msg
237           
238        if node.get('version') and node.get('version') == '1.0':
239           
240            # Get file name
241            entry = get_content('ns:filename', node)
242            if entry is not None:
243                self.file = entry.text.strip()
244           
245            # Get time stamp
246            entry = get_content('ns:timestamp', node)
247            if entry is not None and entry.get('epoch'):
248                try:
249                    self.timestamp = float(entry.get('epoch'))
250                except:
251                    msg = "InversionState.fromXML: Could not read "
252                    msg += "timestamp\n %s" % sys.exc_value
253                    logging.error(msg)
254           
255            # Parse inversion inputs
256            entry = get_content('ns:inputs', node)
257            if entry is not None:
258                for item in in_list:
259                    input_field = get_content('ns:%s' % item[0], entry)
260                    if input_field is not None:
261                        try:
262                            cmd = 'self.%s = float(input_field.text.strip())' 
263                            exec  cmd % item[1]
264                        except:
265                            exec 'self.%s = None' % item[1]
266                input_field = get_content('ns:estimate_bck', entry)
267                if input_field is not None:
268                    try:
269                        self.estimate_bck = input_field.text.strip()=='True'
270                    except:
271                        self.estimate_bck = False
272                   
273            # Parse inversion outputs
274            entry = get_content('ns:outputs', node)
275            if entry is not None:
276                # Output parameters (scalars)
277                for item in out_list:
278                    input_field = get_content('ns:%s' % item[0], entry)
279                    if input_field is not None:
280                        try:
281                            cmd = 'self.%s = float(input_field.text.strip())'
282                            exec  cmd % item[1]
283                        except:
284                            exec 'self.%s = None' % item[1]
285           
286                # Look for coefficients
287                # Format is [value, value, value, value]
288                coeff = get_content('ns:coefficients', entry)
289                if coeff is not None:
290                    # Remove brackets
291                    c_values = coeff.text.strip().replace('[','')
292                    c_values = c_values.replace(']','')
293                    toks = c_values.split()
294                    self.coefficients = []
295                    for c in toks:
296                        try:
297                            self.coefficients.append(float(c))
298                        except:
299                            # Bad data, skip. We will count the number of
300                            # coefficients at the very end and deal with
301                            # inconsistencies then.
302                            pass
303                    # Sanity check
304                    if not len(self.coefficients) == self.nfunc:
305                        # Inconsistent number of coefficients.
306                        # Don't keep the data.
307                        err_msg = "InversionState.fromXML: inconsistant "
308                        err_msg += "number of coefficients: "
309                        err_msg += "%d %d" % (len(self.coefficients),
310                                              self.nfunc)
311                        logging.error(err_msg)
312                        self.coefficients = None
313               
314                # Look for covariance matrix
315                # Format is [ [value, value], [value, value] ]
316                coeff = get_content('ns:covariance', entry)
317                if coeff is not None:
318                    # Parse rows
319                    rows = coeff.text.strip().split('[')
320                    self.covariance = []
321                    for row in rows:
322                        row = row.strip()
323                        if len(row) == 0: continue
324                        # Remove end bracket
325                        row = row.replace(']','')
326                        c_values = row.split()
327                        cov_row = []
328                        for c in c_values:
329                            try:
330                                cov_row.append(float(c))
331                            except:
332                                # Bad data, skip. We will count the number of
333                                # coefficients at the very end and deal with
334                                # inconsistencies then.
335                                pass
336                        # Sanity check: check the number of entries in the row
337                        if len(cov_row) == self.nfunc:
338                            self.covariance.append(cov_row)
339                    # Sanity check: check the number of rows in the covariance
340                    # matrix
341                    if not len(self.covariance) == self.nfunc:
342                        # Inconsistent dimensions of the covariance matrix.
343                        # Don't keep the data.
344                        err_msg = "InversionState.fromXML: "
345                        err_msg += "inconsistant dimensions of the "
346                        err_msg += " covariance matrix: "
347                        err_msg += "%d %d" % (len(self.covariance), self.nfunc)
348                        logging.error(err_msg)
349                        self.covariance = None
350   
351class Reader(CansasReader):
352    """
353    Class to load a .prv P(r) inversion file
354    """
355    ## File type
356    type_name = "P(r)"
357   
358    ## Wildcards
359    type = ["P(r) files (*.prv)|*.prv",
360            "SANSView files (*.svs)|*.svs"]
361    ## List of allowed extensions
362    ext = ['.prv', '.PRV', '.svs', '.SVS'] 
363   
364    def __init__(self, call_back, cansas=True):
365        """
366        Initialize the call-back method to be called
367        after we load a file
368       
369        :param call_back: call-back method
370        :param cansas:  True = files will be written/read in CanSAS format
371                        False = write CanSAS format
372           
373        """
374        ## Call back method to be executed after a file is read
375        self.call_back = call_back
376        ## CanSAS format flag
377        self.cansas = cansas
378        self.state = None
379       
380    def read(self, path):
381        """
382        Load a new P(r) inversion state from file
383       
384        :param path: file path
385       
386        :return: None
387       
388        """
389        if self.cansas==True:
390            return self._read_cansas(path)
391        else:
392            return self._read_standalone(path)
393       
394    def _read_standalone(self, path):
395        """
396        Load a new P(r) inversion state from file.
397        The P(r) node is assumed to be the top element.
398       
399        :param path: file path
400       
401        :return: None
402       
403        """
404        # Read the new state from file
405        state = InversionState()
406        state.fromXML(file=path)
407       
408        # Call back to post the new state
409        self.state = state
410        #self.call_back(state)
411        return None
412   
413    def _parse_prstate(self, entry):
414        """
415        Read a p(r) inversion result from an XML node
416       
417        :param entry: XML node to read from
418       
419        :return: InversionState object
420       
421        """
422        state = None
423       
424        # Locate the P(r) node
425        try:
426            nodes = entry.xpath('ns:%s' % PRNODE_NAME,
427                                namespaces={'ns': CANSAS_NS})
428            if nodes !=[]:
429                # Create an empty state
430                state =  InversionState()
431                state.fromXML(node=nodes[0])
432        except:
433            msg = "XML document does not contain P(r) "
434            msg += "information.\n %s" % sys.exc_value
435            logging.info(msg)
436           
437        return state
438   
439    def _read_cansas(self, path):
440        """
441        Load data and P(r) information from a CanSAS XML file.
442       
443        :param path: file path
444       
445        :return: Data1D object if a single SASentry was found,
446                    or a list of Data1D objects if multiple entries were found,
447                    or None of nothing was found
448                   
449        :raise RuntimeError: when the file can't be opened
450        :raise ValueError: when the length of the data vectors are inconsistent
451       
452        """
453        output = []
454       
455        if os.path.isfile(path):
456            basename  = os.path.basename(path)
457            root, extension = os.path.splitext(basename)
458            #TODO: eventually remove the check for .xml once
459            # the P(r) writer/reader is truly complete.
460            if  extension.lower() in self.ext or extension.lower() == '.xml':
461               
462                tree = etree.parse(path, parser=etree.ETCompatXMLParser())
463                # Check the format version number
464                # Specifying the namespace will take care of the file
465                #format version
466                root = tree.getroot()
467               
468                entry_list = root.xpath('/ns:SASroot/ns:SASentry',
469                                        namespaces={'ns': CANSAS_NS})
470
471                for entry in entry_list:
472                    sas_entry = self._parse_entry(entry)
473                    prstate = self._parse_prstate(entry)
474                    #prstate could be None when .svs file is loaded
475                    #in this case, skip appending to output
476                    if prstate != None:
477                        sas_entry.meta_data['prstate'] = prstate
478                        sas_entry.filename = prstate.file
479                        output.append(sas_entry)
480        else:
481            raise RuntimeError, "%s is not a file" % path
482       
483        # Return output consistent with the loader's api
484        if len(output) == 0:
485            return None
486        elif len(output) == 1:
487            # Call back to post the new state
488            self.call_back(output[0].meta_data['prstate'], datainfo = output[0])
489            #self.state = output[0].meta_data['prstate']
490            return output[0]
491        else:
492            return output               
493   
494   
495    def write(self, filename, datainfo=None, prstate=None):
496        """
497        Write the content of a Data1D as a CanSAS XML file
498       
499        :param filename: name of the file to write
500        :param datainfo: Data1D object
501        :param prstate: InversionState object
502       
503        """
504        # Sanity check
505        if self.cansas == True:
506            doc =self.write_toXML(datainfo, prstate)       
507            # Write the XML document
508            fd = open(filename, 'w')
509            fd.write(doc.toprettyxml())
510            fd.close()
511        else:
512            prstate.toXML(file=filename)
513       
514    def write_toXML(self, datainfo=None, state=None):
515        """
516        Write toXML, a helper for write()
517       
518        : return: xml doc
519        """
520        if datainfo is None:
521            datainfo = Data1D(x=[], y=[])   
522        elif not issubclass(datainfo.__class__, Data1D):
523            msg = "The cansas writer expects a Data1D "
524            msg += "instance: %s" % str(datainfo.__class__.__name__)
525            raise RuntimeError, msg
526   
527        # Create basic XML document
528        doc, sasentry = self._to_xml_doc(datainfo)
529   
530        # Add the invariant information to the XML document
531        if state is not None:
532            state.toXML(doc=doc, entry_node=sasentry)
533           
534        return doc
535   
Note: See TracBrowser for help on using the repository browser.