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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since ded2ce3 was d85c194, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 9 years ago

Remaining modules refactored

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