source: sasview/src/sas/perspectives/pr/inversion_state.py @ b9a5f0e

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 b9a5f0e was b9a5f0e, checked in by krzywon, 9 years ago

90% complete with the conversion.

  • 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 sas.guiframe.dataFitting import Data1D
21from sas.dataloader.readers.cansas_reader import Reader as CansasReader
22from sas.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        #TODO: Get this to work
151        from xml.dom.minidom import getDOMImplementation
152
153        # Check whether we have to write a standalone XML file
154        if doc is None:
155            impl = getDOMImplementation()
156       
157            doc_type = impl.createDocumentType(PRNODE_NAME, "1.0", "1.0")     
158       
159            newdoc = impl.createDocument(None, PRNODE_NAME, doc_type)
160            top_element = newdoc.documentElement
161        else:
162            # We are appending to an existing document
163            newdoc = doc
164            top_element = newdoc.createElement(PRNODE_NAME)
165            if entry_node is None:
166                newdoc.documentElement.appendChild(top_element)
167            else:
168                entry_node.appendChild(top_element)
169           
170        attr = newdoc.createAttribute("version")
171        attr.nodeValue = '1.0'
172        top_element.setAttributeNode(attr)
173       
174        # File name
175        element = newdoc.createElement("filename")
176        if self.file is not None:
177            element.appendChild(newdoc.createTextNode(str(self.file)))
178        else:
179            element.appendChild(newdoc.createTextNode(str(file)))
180        top_element.appendChild(element)
181       
182        element = newdoc.createElement("timestamp")
183        element.appendChild(newdoc.createTextNode(time.ctime(self.timestamp)))
184        attr = newdoc.createAttribute("epoch")
185        attr.nodeValue = str(self.timestamp)
186        element.setAttributeNode(attr)
187        top_element.appendChild(element)
188       
189        # Inputs
190        inputs = newdoc.createElement("inputs")
191        top_element.appendChild(inputs)
192       
193        for item in in_list:
194            element = newdoc.createElement(item[0])
195            cmd = "element.appendChild(newdoc.createTextNode(str(self.%s)))"
196            exec  cmd % item[1]
197            inputs.appendChild(element)
198             
199        # Outputs
200        outputs = newdoc.createElement("outputs")
201        top_element.appendChild(outputs)
202       
203        for item in out_list:
204            element = newdoc.createElement(item[0])
205            cmd = "element.appendChild(newdoc.createTextNode(str(self.%s)))"
206            exec  cmd % item[1]
207            outputs.appendChild(element)
208                   
209        # Save output coefficients and its covariance matrix
210        element = newdoc.createElement("coefficients")
211        element.appendChild(newdoc.createTextNode(str(self.coefficients)))
212        outputs.appendChild(element)
213        element = newdoc.createElement("covariance")
214        element.appendChild(newdoc.createTextNode(str(self.covariance)))
215        outputs.appendChild(element)
216                   
217        # Save the file
218        if doc is None:
219            fd = open(file, 'w')
220            fd.write(newdoc.toprettyxml())
221            fd.close()
222            return None
223        else:
224            return newdoc
225
226    def fromXML(self, file=None, node=None):
227        """
228        Load a P(r) inversion state from a file
229       
230        :param file: .prv file
231        :param node: node of a XML document to read from
232       
233        """
234        if file is not None:
235            msg = "InversionState no longer supports non-CanSAS"
236            msg += " format for P(r) files"
237            raise RuntimeError, msg
238           
239        if node.get('version') and node.get('version') == '1.0':
240           
241            # Get file name
242            entry = get_content('ns:filename', node)
243            if entry is not None:
244                self.file = entry.text.strip()
245           
246            # Get time stamp
247            entry = get_content('ns:timestamp', node)
248            if entry is not None and entry.get('epoch'):
249                try:
250                    self.timestamp = float(entry.get('epoch'))
251                except:
252                    msg = "InversionState.fromXML: Could not read "
253                    msg += "timestamp\n %s" % sys.exc_value
254                    logging.error(msg)
255           
256            # Parse inversion inputs
257            entry = get_content('ns:inputs', node)
258            if entry is not None:
259                for item in in_list:
260                    input_field = get_content('ns:%s' % item[0], entry)
261                    if input_field is not None:
262                        try:
263                            cmd = 'self.%s = float(input_field.text.strip())' 
264                            exec  cmd % item[1]
265                        except:
266                            exec 'self.%s = None' % item[1]
267                input_field = get_content('ns:estimate_bck', entry)
268                if input_field is not None:
269                    try:
270                        self.estimate_bck = input_field.text.strip()=='True'
271                    except:
272                        self.estimate_bck = False
273                   
274            # Parse inversion outputs
275            entry = get_content('ns:outputs', node)
276            if entry is not None:
277                # Output parameters (scalars)
278                for item in out_list:
279                    input_field = get_content('ns:%s' % item[0], entry)
280                    if input_field is not None:
281                        try:
282                            cmd = 'self.%s = float(input_field.text.strip())'
283                            exec  cmd % item[1]
284                        except:
285                            exec 'self.%s = None' % item[1]
286           
287                # Look for coefficients
288                # Format is [value, value, value, value]
289                coeff = get_content('ns:coefficients', entry)
290                if coeff is not None:
291                    # Remove brackets
292                    c_values = coeff.text.strip().replace('[','')
293                    c_values = c_values.replace(']','')
294                    toks = c_values.split()
295                    self.coefficients = []
296                    for c in toks:
297                        try:
298                            self.coefficients.append(float(c))
299                        except:
300                            # Bad data, skip. We will count the number of
301                            # coefficients at the very end and deal with
302                            # inconsistencies then.
303                            pass
304                    # Sanity check
305                    if not len(self.coefficients) == self.nfunc:
306                        # Inconsistent number of coefficients.
307                        # Don't keep the data.
308                        err_msg = "InversionState.fromXML: inconsistant "
309                        err_msg += "number of coefficients: "
310                        err_msg += "%d %d" % (len(self.coefficients),
311                                              self.nfunc)
312                        logging.error(err_msg)
313                        self.coefficients = None
314               
315                # Look for covariance matrix
316                # Format is [ [value, value], [value, value] ]
317                coeff = get_content('ns:covariance', entry)
318                if coeff is not None:
319                    # Parse rows
320                    rows = coeff.text.strip().split('[')
321                    self.covariance = []
322                    for row in rows:
323                        row = row.strip()
324                        if len(row) == 0: continue
325                        # Remove end bracket
326                        row = row.replace(']','')
327                        c_values = row.split()
328                        cov_row = []
329                        for c in c_values:
330                            try:
331                                cov_row.append(float(c))
332                            except:
333                                # Bad data, skip. We will count the number of
334                                # coefficients at the very end and deal with
335                                # inconsistencies then.
336                                pass
337                        # Sanity check: check the number of entries in the row
338                        if len(cov_row) == self.nfunc:
339                            self.covariance.append(cov_row)
340                    # Sanity check: check the number of rows in the covariance
341                    # matrix
342                    if not len(self.covariance) == self.nfunc:
343                        # Inconsistent dimensions of the covariance matrix.
344                        # Don't keep the data.
345                        err_msg = "InversionState.fromXML: "
346                        err_msg += "inconsistant dimensions of the "
347                        err_msg += " covariance matrix: "
348                        err_msg += "%d %d" % (len(self.covariance), self.nfunc)
349                        logging.error(err_msg)
350                        self.covariance = None
351   
352class Reader(CansasReader):
353    """
354    Class to load a .prv P(r) inversion file
355    """
356    ## File type
357    type_name = "P(r)"
358   
359    ## Wildcards
360    type = ["P(r) files (*.prv)|*.prv",
361            "SASView files (*.svs)|*.svs"]
362    ## List of allowed extensions
363    ext = ['.prv', '.PRV', '.svs', '.SVS'] 
364   
365    def __init__(self, call_back, cansas=True):
366        """
367        Initialize the call-back method to be called
368        after we load a file
369       
370        :param call_back: call-back method
371        :param cansas:  True = files will be written/read in CanSAS format
372                        False = write CanSAS format
373           
374        """
375        ## Call back method to be executed after a file is read
376        self.call_back = call_back
377        ## CanSAS format flag
378        self.cansas = cansas
379        self.state = None
380       
381    def read(self, path):
382        """
383        Load a new P(r) inversion state from file
384       
385        :param path: file path
386       
387        :return: None
388       
389        """
390        if self.cansas==True:
391            return self._read_cansas(path)
392        else:
393            return self._read_standalone(path)
394       
395    def _read_standalone(self, path):
396        """
397        Load a new P(r) inversion state from file.
398        The P(r) node is assumed to be the top element.
399       
400        :param path: file path
401       
402        :return: None
403       
404        """
405        # Read the new state from file
406        state = InversionState()
407        state.fromXML(file=path)
408       
409        # Call back to post the new state
410        self.state = state
411        #self.call_back(state)
412        return None
413   
414    def _parse_prstate(self, entry):
415        """
416        Read a p(r) inversion result from an XML node
417       
418        :param entry: XML node to read from
419       
420        :return: InversionState object
421       
422        """
423        state = None
424       
425        # Locate the P(r) node
426        try:
427            nodes = entry.xpath('ns:%s' % PRNODE_NAME,
428                                namespaces={'ns': CANSAS_NS})
429            if nodes !=[]:
430                # Create an empty state
431                state =  InversionState()
432                state.fromXML(node=nodes[0])
433        except:
434            msg = "XML document does not contain P(r) "
435            msg += "information.\n %s" % sys.exc_value
436            logging.info(msg)
437           
438        return state
439   
440    def _read_cansas(self, path):
441        """
442        Load data and P(r) information from a CanSAS XML file.
443       
444        :param path: file path
445       
446        :return: Data1D object if a single SASentry was found,
447                    or a list of Data1D objects if multiple entries were found,
448                    or None of nothing was found
449                   
450        :raise RuntimeError: when the file can't be opened
451        :raise ValueError: when the length of the data vectors are inconsistent
452       
453        """
454        output = []
455       
456        if os.path.isfile(path):
457            basename  = os.path.basename(path)
458            root, extension = os.path.splitext(basename)
459            #TODO: eventually remove the check for .xml once
460            # the P(r) writer/reader is truly complete.
461            if  extension.lower() in self.ext or extension.lower() == '.xml':
462               
463                tree = etree.parse(path, parser=etree.ETCompatXMLParser())
464                # Check the format version number
465                # Specifying the namespace will take care of the file
466                #format version
467                root = tree.getroot()
468               
469                entry_list = root.xpath('/ns:SASroot/ns:SASentry',
470                                        namespaces={'ns': CANSAS_NS})
471
472                for entry in entry_list:
473                    sas_entry, _ = self._parse_entry(entry)
474                    prstate = self._parse_prstate(entry)
475                    #prstate could be None when .svs file is loaded
476                    #in this case, skip appending to output
477                    if prstate != None:
478                        sas_entry.meta_data['prstate'] = prstate
479                        sas_entry.filename = prstate.file
480                        output.append(sas_entry)
481        else:
482            raise RuntimeError, "%s is not a file" % path
483       
484        # Return output consistent with the loader's api
485        if len(output) == 0:
486            return None
487        elif len(output) == 1:
488            # Call back to post the new state
489            self.call_back(output[0].meta_data['prstate'], datainfo = output[0])
490            #self.state = output[0].meta_data['prstate']
491            return output[0]
492        else:
493            return output               
494   
495   
496    def write(self, filename, datainfo=None, prstate=None):
497        """
498        Write the content of a Data1D as a CanSAS XML file
499       
500        :param filename: name of the file to write
501        :param datainfo: Data1D object
502        :param prstate: InversionState object
503       
504        """
505        # Sanity check
506        if self.cansas == True:
507            doc = self.write_toXML(datainfo, prstate)       
508            # Write the XML document
509            fd = open(filename, 'w')
510            fd.write(doc.toprettyxml())
511            fd.close()
512        else:
513            prstate.toXML(file=filename)
514       
515    def write_toXML(self, datainfo=None, state=None):
516        """
517        Write toXML, a helper for write()
518       
519        : return: xml doc
520        """
521        if datainfo is None:
522            datainfo = Data1D(x=[], y=[])   
523        elif not issubclass(datainfo.__class__, Data1D):
524            msg = "The cansas writer expects a Data1D "
525            msg += "instance: %s" % str(datainfo.__class__.__name__)
526            raise RuntimeError, msg
527   
528        # Create basic XML document
529        doc, sasentry = self._to_xml_doc(datainfo)
530       
531        # Add the invariant information to the XML document
532        if state is not None:
533            doc = state.toXML(doc=doc, entry_node=sasentry)
534           
535        return doc
536   
Note: See TracBrowser for help on using the repository browser.