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

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 0d72cac was fa81e94, checked in by Piotr Rozyczko <rozyczko@…>, 7 years ago

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