source: sasview/DataLoader/data_info.py @ 3f36675

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 3f36675 was 99abf5d, checked in by Gervaise Alina <gervyh@…>, 14 years ago

add way to edit process on data
  • Property mode set to 100644
File size: 28.5 KB
Line 
1"""
2    Module that contains classes to hold information read from
3    reduced data files.
4   
5    A good description of the data members can be found in
6    the CanSAS 1D XML data format:
7   
8    http://www.smallangles.net/wgwiki/index.php/cansas1d_documentation
9"""
10
11
12#This software was developed by the University of Tennessee as part of the
13#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
14#project funded by the US National Science Foundation.
15#If you use DANSE applications to do scientific research that leads to
16#publication, we ask that you acknowledge the use of the software with the
17#following sentence:
18#This work benefited from DANSE software developed under NSF award DMR-0520547.
19#copyright 2008, University of Tennessee
20
21
22#TODO: Keep track of data manipulation in the 'process' data structure.
23#TODO: This module should be independent of plottables. We should write
24#        an adapter class for plottables when needed.
25
26#from sans.guitools.plottables import Data1D as plottable_1D
27from data_util.uncertainty import Uncertainty
28import numpy
29import math
30
31class plottable_1D:
32    """
33    Data1D is a place holder for 1D plottables.
34    """
35    # The presence of these should be mutually
36    # exclusive with the presence of Qdev (dx)
37    x = None
38    y = None
39    dx = None
40    dy = None
41    ## Slit smearing length
42    dxl = None
43    ## Slit smearing width
44    dxw = None
45   
46    # Units
47    _xaxis = ''
48    _xunit = ''
49    _yaxis = ''
50    _yunit = ''
51   
52    def __init__(self, x, y, dx=None, dy=None, dxl=None, dxw=None):
53        self.x = numpy.asarray(x)
54        self.y = numpy.asarray(y)
55        if dx is not None:
56             self.dx = numpy.asarray(dx)
57        if dy is not None: 
58            self.dy = numpy.asarray(dy)
59        if dxl is not None: 
60            self.dxl = numpy.asarray(dxl)
61        if dxw is not None: 
62            self.dxw = numpy.asarray(dxw)
63
64    def xaxis(self, label, unit):
65        """
66        set the x axis label and unit
67        """
68        self._xaxis = label
69        self._xunit = unit
70       
71    def yaxis(self, label, unit):
72        """
73        set the y axis label and unit
74        """
75        self._yaxis = label
76        self._yunit = unit
77
78class plottable_2D:
79    """
80    Data2D is a place holder for 2D plottables.
81    """
82    xmin = None
83    xmax = None
84    ymin = None
85    ymax = None
86    data = None
87    qx_data     = None
88    qy_data     = None
89    q_data      = None
90    err_data    = None
91    dqx_data    = None
92    dqy_data    = None
93    mask        = None
94   
95    # Units
96    _xaxis = ''
97    _xunit = ''
98    _yaxis = ''
99    _yunit = ''
100    _zaxis = ''
101    _zunit = ''
102   
103    def __init__(self, data=None, err_data=None, qx_data=None,
104                  qy_data=None, q_data=None, mask=None,
105                   dqx_data=None, dqy_data=None):
106        self.data = numpy.asarray(data)
107        self.qx_data = numpy.asarray(qx_data)
108        self.qy_data = numpy.asarray(qy_data)
109        self.q_data = numpy.asarray(q_data)
110        self.mask = numpy.asarray(mask)
111        self.err_data = numpy.asarray(err_data)
112        if dqx_data is not None: self.dqx_data = numpy.asarray(dqx_data) 
113        if dqy_data is not None: self.dqy_data = numpy.asarray(dqy_data) 
114               
115    def xaxis(self, label, unit):
116        """
117        set the x axis label and unit
118        """
119        self._xaxis = label
120        self._xunit = unit
121       
122    def yaxis(self, label, unit):
123        """
124        set the y axis label and unit
125        """
126        self._yaxis = label
127        self._yunit = unit
128           
129    def zaxis(self, label, unit):
130        """
131        set the z axis label and unit
132        """
133        self._zaxis = label
134        self._zunit = unit
135
136           
137class Vector:
138    """
139    Vector class to hold multi-dimensional objects
140    """
141    ## x component
142    x = None
143    ## y component
144    y = None
145    ## z component
146    z = None
147   
148    def __init__(self, x=None, y=None, z=None):
149        """
150        Initialization. Components that are not
151        set a set to None by default.
152       
153        :param x: x component
154        :param y: y component
155        :param z: z component
156       
157        """
158        self.x = x
159        self.y = y
160        self.z = z
161       
162    def __str__(self):
163        msg = "x = %s\ty = %s\tz = %s" % (str(self.x), str(self.y), str(self.z))
164        return msg
165       
166
167class Detector:
168    """
169    Class to hold detector information
170    """
171    ## Name of the instrument [string]
172    name = None
173    ## Sample to detector distance [float] [mm]
174    distance = None
175    distance_unit = 'mm'
176    ## Offset of this detector position in X, Y,
177    #(and Z if necessary) [Vector] [mm]
178    offset = None
179    offset_unit = 'm'
180    ## Orientation (rotation) of this detector in roll,
181    # pitch, and yaw [Vector] [degrees]
182    orientation = None
183    orientation_unit = 'degree'
184    ## Center of the beam on the detector in X and Y
185    #(and Z if necessary) [Vector] [mm]
186    beam_center = None
187    beam_center_unit = 'mm'
188    ## Pixel size in X, Y, (and Z if necessary) [Vector] [mm]
189    pixel_size = None
190    pixel_size_unit = 'mm'
191    ## Slit length of the instrument for this detector.[float] [mm]
192    slit_length = None
193    slit_length_unit = 'mm'
194   
195    def __init__(self):
196        """
197       
198        Initialize class attribute that are objects...
199       
200        """
201        self.offset      = Vector()
202        self.orientation = Vector()
203        self.beam_center = Vector()
204        self.pixel_size  = Vector()
205       
206   
207    def __str__(self):
208        _str  = "Detector:\n"
209        _str += "   Name:         %s\n" % self.name
210        _str += "   Distance:     %s [%s]\n" % \
211            (str(self.distance), str(self.distance_unit))
212        _str += "   Offset:       %s [%s]\n" % \
213            (str(self.offset), str(self.offset_unit))
214        _str += "   Orientation:  %s [%s]\n" % \
215            (str(self.orientation), str(self.orientation_unit))
216        _str += "   Beam center:  %s [%s]\n" % \
217            (str(self.beam_center), str(self.beam_center_unit))
218        _str += "   Pixel size:   %s [%s]\n" % \
219            (str(self.pixel_size), str(self.pixel_size_unit))
220        _str += "   Slit length:  %s [%s]\n" % \
221            (str(self.slit_length), str(self.slit_length_unit))
222        return _str
223
224class Aperture:
225    ## Name
226    name = None
227    ## Type
228    type = None
229    ## Size name
230    size_name = None
231    ## Aperture size [Vector]
232    size = None
233    size_unit = 'mm'
234    ## Aperture distance [float]
235    distance = None
236    distance_unit = 'mm'
237   
238    def __init__(self):
239        self.size = Vector()
240   
241class Collimation:
242    """
243    Class to hold collimation information
244    """
245    ## Name
246    name = None
247    ## Length [float] [mm]
248    length = None
249    length_unit = 'mm'
250    ## Aperture
251    aperture = None
252   
253    def __init__(self):
254        self.aperture = []
255   
256    def __str__(self):
257        _str = "Collimation:\n"
258        _str += "   Length:       %s [%s]\n" % \
259            (str(self.length), str(self.length_unit))
260        for item in self.aperture:
261            _str += "   Aperture size:%s [%s]\n" % \
262                (str(item.size), str(item.size_unit))
263            _str += "   Aperture_dist:%s [%s]\n" % \
264                (str(item.distance), str(item.distance_unit))
265        return _str
266
267class Source:
268    """
269    Class to hold source information
270    """ 
271    ## Name
272    name = None
273    ## Radiation type [string]
274    radiation = None
275    ## Beam size name
276    beam_size_name = None
277    ## Beam size [Vector] [mm]
278    beam_size = None
279    beam_size_unit = 'mm'
280    ## Beam shape [string]
281    beam_shape = None
282    ## Wavelength [float] [Angstrom]
283    wavelength = None
284    wavelength_unit = 'A'
285    ## Minimum wavelength [float] [Angstrom]
286    wavelength_min = None
287    wavelength_min_unit = 'nm'
288    ## Maximum wavelength [float] [Angstrom]
289    wavelength_max = None
290    wavelength_max_unit = 'nm'
291    ## Wavelength spread [float] [Angstrom]
292    wavelength_spread = None
293    wavelength_spread_unit = 'percent'
294   
295    def __init__(self):
296        self.beam_size = Vector()
297       
298   
299    def __str__(self):
300        _str  = "Source:\n"
301        _str += "   Radiation:    %s\n" % str(self.radiation)
302        _str += "   Shape:        %s\n" % str(self.beam_shape)
303        _str += "   Wavelength:   %s [%s]\n" % \
304            (str(self.wavelength), str(self.wavelength_unit))
305        _str += "   Waveln_min:   %s [%s]\n" % \
306            (str(self.wavelength_min), str(self.wavelength_min_unit))
307        _str += "   Waveln_max:   %s [%s]\n" % \
308            (str(self.wavelength_max), str(self.wavelength_max_unit))
309        _str += "   Waveln_spread:%s [%s]\n" % \
310            (str(self.wavelength_spread), str(self.wavelength_spread_unit))
311        _str += "   Beam_size:    %s [%s]\n" % \
312            (str(self.beam_size), str(self.beam_size_unit))
313        return _str
314   
315   
316"""
317Definitions of radiation types
318"""
319NEUTRON  = 'neutron'
320XRAY     = 'x-ray'
321MUON     = 'muon'
322ELECTRON = 'electron'
323   
324class Sample:
325    """
326    Class to hold the sample description
327    """
328    ## Short name for sample
329    name = ''
330    ## ID
331    ID = ''
332    ## Thickness [float] [mm]
333    thickness = None
334    thickness_unit = 'mm'
335    ## Transmission [float] [fraction]
336    transmission = None
337    ## Temperature [float] [C]
338    temperature = None
339    temperature_unit = 'C'
340    ## Position [Vector] [mm]
341    position = None
342    position_unit = 'mm'
343    ## Orientation [Vector] [degrees]
344    orientation = None
345    orientation_unit = 'degree'
346    ## Details
347    details = None
348   
349    def __init__(self):
350        self.position    = Vector()
351        self.orientation = Vector()
352        self.details     = []
353   
354    def __str__(self):
355        _str  = "Sample:\n"
356        _str += "   ID:           %s\n" % str(self.ID)
357        _str += "   Transmission: %s\n" % str(self.transmission)
358        _str += "   Thickness:    %s [%s]\n" % \
359            (str(self.thickness), str(self.thickness_unit))
360        _str += "   Temperature:  %s [%s]\n" % \
361            (str(self.temperature), str(self.temperature_unit))
362        _str += "   Position:     %s [%s]\n" % \
363            (str(self.position), str(self.position_unit))
364        _str += "   Orientation:  %s [%s]\n" % \
365            (str(self.orientation), str(self.orientation_unit))
366       
367        _str += "   Details:\n"
368        for item in self.details:
369            _str += "      %s\n" % item
370           
371        return _str
372 
373class Process:
374    """
375    Class that holds information about the processes
376    performed on the data.
377    """
378    name = ''
379    date = ''
380    description = ''
381    term = None
382    notes = None
383   
384    def __init__(self):
385        self.term = []
386        self.notes = []
387   
388    def __str__(self):
389        _str  = "Process:\n"
390        _str += "   Name:         %s\n" % self.name
391        _str += "   Date:         %s\n" % self.date
392        _str += "   Description:  %s\n" % self.description
393        for item in self.term:
394            _str += "   Term:         %s\n" % item
395        for item in self.notes:
396            _str += "   Note:         %s\n" % item
397        return _str
398   
399 
400class DataInfo:
401    """
402    Class to hold the data read from a file.
403    It includes four blocks of data for the
404    instrument description, the sample description,
405    the data itself and any other meta data.
406    """
407    ## Title
408    title      = ''
409    ## Run number
410    run        = None
411    ## Run name
412    run_name   = None
413    ## File name
414    filename   = ''
415    ## Notes
416    notes      = None
417    ## Processes (Action on the data)
418    process    = None
419    ## Instrument name
420    instrument = ''
421    ## Detector information
422    detector   = None
423    ## Sample information
424    sample     = None
425    ## Source information
426    source     = None
427    ## Collimation information
428    collimation = None
429    ## Additional meta-data
430    meta_data  = None
431    ## Loading errors
432    errors = None
433           
434    def __init__(self):
435        """
436        Initialization
437        """
438        ## Title
439        self.title      = ''
440        ## Run number
441        self.run        = []
442        self.run_name   = {}
443        ## File name
444        self.filename   = ''
445        ## Notes
446        self.notes      = []
447        ## Processes (Action on the data)
448        self.process    = []
449        ## Instrument name
450        self.instrument = ''
451        ## Detector information
452        self.detector   = []
453        ## Sample information
454        self.sample     = Sample()
455        ## Source information
456        self.source     = Source()
457        ## Collimation information
458        self.collimation = []
459        ## Additional meta-data
460        self.meta_data  = {}
461        ## Loading errors
462        self.errors = []       
463       
464    def append_empty_process(self):
465        """
466        """
467        self.process.append(Process())
468       
469    def add_notes(self, message=""):
470        """
471        Add notes to datainfo
472        """
473        self.notes.append(message)
474       
475    def __str__(self):
476        """
477        Nice printout
478        """
479        _str =  "File:            %s\n" % self.filename
480        _str += "Title:           %s\n" % self.title
481        _str += "Run:             %s\n" % str(self.run)
482        _str += "Instrument:      %s\n" % str(self.instrument)
483        _str += "%s\n" % str(self.sample)
484        _str += "%s\n" % str(self.source)
485        for item in self.detector:
486            _str += "%s\n" % str(item)
487        for item in self.collimation:
488            _str += "%s\n" % str(item)
489        for item in self.process:
490            _str += "%s\n" % str(item)
491        for item in self.notes:
492            _str += "%s\n" % str(item)
493
494        return _str
495           
496    # Private method to perform operation. Not implemented for DataInfo,
497    # but should be implemented for each data class inherited from DataInfo
498    # that holds actual data (ex.: Data1D)
499    def _perform_operation(self, other, operation): 
500        """
501        Private method to perform operation. Not implemented for DataInfo,
502        but should be implemented for each data class inherited from DataInfo
503        that holds actual data (ex.: Data1D)
504        """
505        return NotImplemented
506
507    def __add__(self, other):
508        """
509        Add two data sets
510       
511        :param other: data set to add to the current one
512        :return: new data set
513        :raise ValueError: raised when two data sets are incompatible
514        """
515        def operation(a, b):
516            return a + b
517        return self._perform_operation(other, operation)
518       
519    def __radd__(self, other):
520        """
521        Add two data sets
522       
523        :param other: data set to add to the current one
524       
525        :return: new data set
526       
527        :raise ValueError: raised when two data sets are incompatible
528       
529        """
530        def operation(a, b):
531            return b + a
532        return self._perform_operation(other, operation)
533       
534    def __sub__(self, other):
535        """
536        Subtract two data sets
537       
538        :param other: data set to subtract from the current one
539       
540        :return: new data set
541       
542        :raise ValueError: raised when two data sets are incompatible
543       
544        """
545        def operation(a, b):
546            return a - b
547        return self._perform_operation(other, operation)
548       
549    def __rsub__(self, other):
550        """
551        Subtract two data sets
552       
553        :param other: data set to subtract from the current one
554       
555        :return: new data set
556       
557        :raise ValueError: raised when two data sets are incompatible
558       
559        """
560        def operation(a, b):
561            return b - a
562        return self._perform_operation(other, operation)
563       
564    def __mul__(self, other):
565        """
566        Multiply two data sets
567       
568        :param other: data set to subtract from the current one
569       
570        :return: new data set
571       
572        :raise ValueError: raised when two data sets are incompatible
573       
574        """
575        def operation(a, b):
576            return a * b
577        return self._perform_operation(other, operation)
578       
579    def __rmul__(self, other):
580        """
581        Multiply two data sets
582       
583        :param other: data set to subtract from the current one
584       
585        :return: new data set
586       
587        :raise ValueError: raised when two data sets are incompatible
588        """
589        def operation(a, b):
590            return b * a
591        return self._perform_operation(other, operation)
592       
593    def __div__(self, other):
594        """
595        Divided a data set by another
596       
597        :param other: data set that the current one is divided by
598       
599        :return: new data set
600       
601        :raise ValueError: raised when two data sets are incompatible
602       
603        """
604        def operation(a, b):
605            return a/b
606        return self._perform_operation(other, operation)
607       
608    def __rdiv__(self, other):
609        """
610        Divided a data set by another
611       
612        :param other: data set that the current one is divided by
613       
614        :return: new data set
615       
616        :raise ValueError: raised when two data sets are incompatible
617       
618        """
619        def operation(a, b):
620            return b/a
621        return self._perform_operation(other, operation)           
622           
623class Data1D(plottable_1D, DataInfo):
624    """
625    1D data class
626    """
627    x_unit = '1/A'
628    y_unit = '1/cm'
629   
630    def __init__(self, x, y, dx=None, dy=None):
631        DataInfo.__init__(self)
632        plottable_1D.__init__(self, x, y, dx, dy)
633       
634       
635    def __str__(self):
636        """
637        Nice printout
638        """
639        _str =  "%s\n" % DataInfo.__str__(self)
640   
641        _str += "Data:\n"
642        _str += "   Type:         %s\n" % self.__class__.__name__
643        _str += "   X-axis:       %s\t[%s]\n" % (self._xaxis, self._xunit)
644        _str += "   Y-axis:       %s\t[%s]\n" % (self._yaxis, self._yunit)
645        _str += "   Length:       %g\n" % len(self.x)
646
647        return _str
648
649    def is_slit_smeared(self):
650        """
651        Check whether the data has slit smearing information
652       
653        :return: True is slit smearing info is present, False otherwise
654       
655        """
656        def _check(v):           
657            if (v.__class__ == list or v.__class__ == numpy.ndarray) \
658                and len(v) > 0 and min(v) > 0:
659                return True
660           
661            return False
662       
663        return _check(self.dxl) or _check(self.dxw)
664       
665    def clone_without_data(self, length=0, clone=None):
666        """
667        Clone the current object, without copying the data (which
668        will be filled out by a subsequent operation).
669        The data arrays will be initialized to zero.
670       
671        :param length: length of the data array to be initialized
672        :param clone: if provided, the data will be copied to clone
673        """
674        from copy import deepcopy
675       
676        if clone is None or not issubclass(clone.__class__, Data1D):
677            x  = numpy.zeros(length) 
678            dx = numpy.zeros(length) 
679            y  = numpy.zeros(length) 
680            dy = numpy.zeros(length) 
681            clone = Data1D(x, y, dx=dx, dy=dy)
682       
683        clone.title       = self.title
684        clone.run         = self.run
685        clone.filename    = self.filename
686        clone.notes       = deepcopy(self.notes) 
687        clone.process     = deepcopy(self.process) 
688        clone.detector    = deepcopy(self.detector) 
689        clone.sample      = deepcopy(self.sample) 
690        clone.source      = deepcopy(self.source) 
691        clone.collimation = deepcopy(self.collimation) 
692        clone.meta_data   = deepcopy(self.meta_data) 
693        clone.errors      = deepcopy(self.errors) 
694       
695        return clone
696
697    def _validity_check(self, other):
698        """
699        Checks that the data lengths are compatible.
700        Checks that the x vectors are compatible.
701        Returns errors vectors equal to original
702        errors vectors if they were present or vectors
703        of zeros when none was found.
704       
705        :param other: other data set for operation
706       
707        :return: dy for self, dy for other [numpy arrays]
708       
709        :raise ValueError: when lengths are not compatible
710       
711        """
712        dy_other = None
713        if isinstance(other, Data1D):
714            # Check that data lengths are the same
715            if len(self.x) != len(other.x) or \
716                len(self.y) != len(other.y):
717                msg =  "Unable to perform operation: data length are not equal"
718                raise ValueError, msg
719           
720            # Here we could also extrapolate between data points
721            for i in range(len(self.x)):
722                if self.x[i] != other.x[i]:
723                    msg = "Incompatible data sets: x-values do not match"
724                    raise ValueError, msg
725           
726            # Check that the other data set has errors, otherwise
727            # create zero vector
728            dy_other = other.dy
729            if other.dy == None or (len(other.dy) != len(other.y)):
730                dy_other = numpy.zeros(len(other.y))
731           
732        # Check that we have errors, otherwise create zero vector
733        dy = self.dy
734        if self.dy == None or (len(self.dy) != len(self.y)):
735            dy = numpy.zeros(len(self.y))           
736           
737        return dy, dy_other
738
739    def _perform_operation(self, other, operation):
740        """
741        """
742        # First, check the data compatibility
743        dy, dy_other = self._validity_check(other)
744        result = self.clone_without_data(len(self.x))
745       
746        for i in range(len(self.x)):
747            result.x[i] = self.x[i]
748            if self.dx is not None and len(self.x) == len(self.dx):
749                result.dx[i] = self.dx[i]
750           
751            a = Uncertainty(self.y[i], dy[i]**2)
752            if isinstance(other, Data1D):
753                b = Uncertainty(other.y[i], dy_other[i]**2)
754            else:
755                b = other
756           
757            output = operation(a, b)
758            result.y[i] = output.x
759            result.dy[i] = math.sqrt(math.fabs(output.variance))
760        return result
761       
762class Data2D(plottable_2D, DataInfo):
763    """
764    2D data class
765    """
766    ## Units for Q-values
767    Q_unit = '1/A'
768   
769    ## Units for I(Q) values
770    I_unit = '1/cm'
771   
772    ## Vector of Q-values at the center of each bin in x
773    x_bins = None
774   
775    ## Vector of Q-values at the center of each bin in y
776    y_bins = None
777   
778   
779    def __init__(self, data=None, err_data=None, qx_data=None,
780                  qy_data=None, q_data=None, mask=None, 
781                  dqx_data=None, dqy_data=None):
782        self.y_bins = []
783        self.x_bins = []
784        DataInfo.__init__(self)
785        plottable_2D.__init__(self, data, err_data, qx_data,
786                              qy_data, q_data,mask, dqx_data, dqy_data)
787        if len(self.detector) > 0:
788            raise RuntimeError, "Data2D: Detector bank already filled at init"
789
790    def __str__(self):
791        _str =  "%s\n" % DataInfo.__str__(self)
792       
793        _str += "Data:\n"
794        _str += "   Type:         %s\n" % self.__class__.__name__
795        _str += "   X- & Y-axis:  %s\t[%s]\n" % (self._yaxis, self._yunit)
796        _str += "   Z-axis:       %s\t[%s]\n" % (self._zaxis, self._zunit)
797        #leny = 0
798        #if len(self.data) > 0:
799        #    leny = len(self.data)
800        _str += "   Length:       %g \n" % (len(self.data))
801       
802        return _str
803 
804    def clone_without_data(self, length=0, clone=None):
805        """
806        Clone the current object, without copying the data (which
807        will be filled out by a subsequent operation).
808        The data arrays will be initialized to zero.
809       
810        :param length: length of the data array to be initialized
811        :param clone: if provided, the data will be copied to clone
812        """
813        from copy import deepcopy
814       
815        if clone is None or not issubclass(clone.__class__, Data2D): 
816            data     = numpy.zeros(length) 
817            err_data = numpy.zeros(length) 
818            qx_data = numpy.zeros(length)
819            qy_data = numpy.zeros(length)
820            q_data = numpy.zeros(length)
821            mask = numpy.zeros(length)
822            dqx_data = None
823            dqy_data = None
824            clone = Data2D(data, err_data, qx_data, qy_data,
825                            q_data,mask, dqx_data=dqx_data, dqy_data=dqy_data)
826
827        clone.title       = self.title
828        clone.run         = self.run
829        clone.filename    = self.filename
830        clone.notes       = deepcopy(self.notes) 
831        clone.process     = deepcopy(self.process) 
832        clone.detector    = deepcopy(self.detector) 
833        clone.sample      = deepcopy(self.sample) 
834        clone.source      = deepcopy(self.source) 
835        clone.collimation = deepcopy(self.collimation) 
836        clone.meta_data   = deepcopy(self.meta_data) 
837        clone.errors      = deepcopy(self.errors) 
838       
839        return clone
840 
841 
842    def _validity_check(self, other):
843        """
844        Checks that the data lengths are compatible.
845        Checks that the x vectors are compatible.
846        Returns errors vectors equal to original
847        errors vectors if they were present or vectors
848        of zeros when none was found.
849       
850        :param other: other data set for operation
851       
852        :return: dy for self, dy for other [numpy arrays]
853       
854        :raise ValueError: when lengths are not compatible
855       
856        """
857        err_other = None
858        if isinstance(other, Data2D):
859            # Check that data lengths are the same
860            if numpy.size(self.data) != numpy.size(other.data):
861                msg = "Unable to perform operation: data length are not equal"
862                raise ValueError, msg
863               
864            # Check that the scales match
865            #TODO: matching scales?     
866           
867            # Check that the other data set has errors, otherwise
868            # create zero vector
869            #TODO: test this
870            err_other = other.err_data
871            if other.err_data == None or \
872                (numpy.size(other.err_data) != numpy.size(other.data)):
873                err_other = numpy.zeros([numpy.size(other.data, 0),
874                                          numpy.size(other.data, 1)])
875           
876        # Check that we have errors, otherwise create zero vector
877        err = self.err_data
878        if self.err_data == None or \
879            (numpy.size(self.err_data) != numpy.size(self.data)):
880            err = numpy.zeros([numpy.size(self.data, 0),
881                                numpy.size(self.data, 1)])
882           
883        return err, err_other
884 
885 
886    def _perform_operation(self, other, operation):
887        """
888        Perform 2D operations between data sets
889       
890        :param other: other data set
891        :param operation: function defining the operation
892       
893        """
894        # First, check the data compatibility
895        dy, dy_other = self._validity_check(other)
896   
897        result = self.clone_without_data([numpy.size(self.data, 0), 
898                                          numpy.size(self.data, 1)])
899       
900        for i in range(numpy.size(self.data, 0)):
901            for j in range(numpy.size(self.data, 1)):
902                result.data[i][j] = self.data[i][j]
903                if self.err_data is not None and \
904                    numpy.size(self.data) == numpy.size(self.err_data):
905                    result.err_data[i][j] = self.err_data[i][j]
906               
907                a = Uncertainty(self.data[i][j], dy[i][j]**2)
908                if isinstance(other, Data2D):
909                    b = Uncertainty(other.data[i][j], dy_other[i][j]**2)
910                else:
911                    b = other
912               
913                output = operation(a, b)
914                result.data[i][j] = output.x
915                result.err_data[i][j] = math.sqrt(math.fabs(output.variance))
916        return result
917   
918 
919 
Note: See TracBrowser for help on using the repository browser.