source: sasview/prview/perspectives/pr/inversion_state.py @ ba69349

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 ba69349 was 80d2872, checked in by Mathieu Doucet <doucetm@…>, 15 years ago

prview: ensure that data selected for inversion are of the right class.

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