Index: nsinvariant/invariant.py
===================================================================
--- sansinvariant/invariant.py	(revision ad4053a35ebb020c5b26adb27098e74140045aca)
+++ 	(revision )
@@ -1,957 +1,0 @@
-#####################################################################
-#This software was developed by the University of Tennessee as part of the
-#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
-#project funded by the US National Science Foundation. 
-#See the license text in license.txt
-#copyright 2010, University of Tennessee
-######################################################################
-
-"""
-This module implements invariant and its related computations.
-
-:author: Gervaise B. Alina/UTK
-:author: Mathieu Doucet/UTK
-:author: Jae Cho/UTK
-
-"""
-import math 
-import numpy
-
-from DataLoader.data_info import Data1D as LoaderData1D
-
-# The minimum q-value to be used when extrapolating
-Q_MINIMUM  = 1e-5
-
-# The maximum q-value to be used when extrapolating
-Q_MAXIMUM  = 10
-
-# Number of steps in the extrapolation
-INTEGRATION_NSTEPS = 1000
-
-class Transform(object):
-    """
-    Define interface that need to compute a function or an inverse
-    function given some x, y 
-    """
-    
-    def linearize_data(self, data):
-        """
-        Linearize data so that a linear fit can be performed. 
-        Filter out the data that can't be transformed.
-        
-        :param data: LoadData1D instance
-        
-        """
-        # Check that the vector lengths are equal
-        assert(len(data.x)==len(data.y))
-        if data.dy is not None:
-            assert(len(data.x)==len(data.dy))
-            dy = data.dy
-        else:
-            dy = numpy.ones(len(data.y))
-            
-        # Transform the data
-        data_points = zip(data.x, data.y, dy)
-
-        output_points = [(self.linearize_q_value(p[0]),
-                          math.log(p[1]),
-                          p[2]/p[1]) for p in data_points if p[0]>0 and \
-                          p[1]>0 and p[2]>0]
-        
-        x_out, y_out, dy_out = zip(*output_points)
-       
-        # Create Data1D object
-        x_out = numpy.asarray(x_out)
-        y_out = numpy.asarray(y_out)
-        dy_out = numpy.asarray(dy_out)
-        linear_data = LoaderData1D(x=x_out, y=y_out, dy=dy_out)
-        
-        return linear_data
-    
-    def get_allowed_bins(self, data):
-        """
-        Goes through the data points and returns a list of boolean values
-        to indicate whether each points is allowed by the model or not.
-        
-        :param data: Data1D object
-        """
-        return [p[0]>0 and p[1]>0 and p[2]>0 for p in zip(data.x, data.y,
-                                                           data.dy)]
-        
-    def linearize_q_value(self, value):
-        """
-        Transform the input q-value for linearization
-        """
-        return NotImplemented
-
-    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
-        """
-        set private member
-        """
-        return NotImplemented
-     
-    def evaluate_model(self, x):
-        """
-        Returns an array f(x) values where f is the Transform function.
-        """
-        return NotImplemented
-    
-    def evaluate_model_errors(self, x):
-        """
-        Returns an array of I(q) errors
-        """
-        return NotImplemented
-    
-class Guinier(Transform):
-    """
-    class of type Transform that performs operations related to guinier 
-    function
-    """
-    def __init__(self, scale=1, radius=60):
-        Transform.__init__(self)
-        self.scale = scale
-        self.radius = radius
-        ## Uncertainty of scale parameter
-        self.dscale  = 0
-        ## Unvertainty of radius parameter
-        self.dradius = 0
-        
-    def linearize_q_value(self, value):
-        """
-        Transform the input q-value for linearization
-        
-        :param value: q-value
-        
-        :return: q*q
-        """
-        return value * value
-    
-    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
-        """
-	    assign new value to the scale and the radius
-    	"""
-        self.scale = math.exp(constant)
-        self.radius = math.sqrt(-3 * slope)
-        
-        # Errors
-        self.dscale = math.exp(constant)*dconstant
-        self.dradius = -3.0/2.0/math.sqrt(-3 * slope)*dslope
-        
-        return [self.radius, self.scale], [self.dradius, self.dscale]
-        
-    def evaluate_model(self, x):
-        """
-        return F(x)= scale* e-((radius*x)**2/3)
-        """
-        return self._guinier(x)
-             
-    def evaluate_model_errors(self, x):
-        """
-        Returns the error on I(q) for the given array of q-values
-        
-        :param x: array of q-values
-        """
-        p1 = numpy.array([self.dscale * math.exp(-((self.radius * q)**2/3)) \
-                          for q in x])
-        p2 = numpy.array([self.scale * math.exp(-((self.radius * q)**2/3))\
-                     * (-(q**2/3)) * 2 * self.radius * self.dradius for q in x])
-        diq2 = p1*p1 + p2*p2        
-        return numpy.array( [math.sqrt(err) for err in diq2] )
-             
-    def _guinier(self, x):
-        """
-        Retrive the guinier function after apply an inverse guinier function
-        to x
-        Compute a F(x) = scale* e-((radius*x)**2/3).
-        
-        :param x: a vector of q values 
-        :param scale: the scale value
-        :param radius: the guinier radius value
-        
-        :return: F(x)
-        """   
-        # transform the radius of coming from the inverse guinier function to a 
-        # a radius of a guinier function
-        if self.radius <= 0:
-            msg = "Rg expected positive value, but got %s"%self.radius
-            raise ValueError(msg)  
-        value = numpy.array([math.exp(-((self.radius * i)**2/3)) for i in x ]) 
-        return self.scale * value
-
-class PowerLaw(Transform):
-    """
-    class of type transform that perform operation related to power_law 
-    function
-    """
-    def __init__(self, scale=1, power=4):
-        Transform.__init__(self)
-        self.scale = scale
-        self.power = power
-   
-    def linearize_q_value(self, value):
-        """
-        Transform the input q-value for linearization
-        
-        :param value: q-value
-        
-        :return: log(q)
-        """
-        return math.log(value)
-    
-    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
-        """
-        Assign new value to the scale and the power 
-        """
-        self.power = -slope
-        self.scale = math.exp(constant)
-        
-        # Errors
-        self.dscale = math.exp(constant)*dconstant
-        self.dpower = -dslope
-        
-        return [self.power, self.scale], [self.dpower, self.dscale]
-        
-    def evaluate_model(self, x):
-        """
-        given a scale and a radius transform x, y using a power_law
-        function
-        """
-        return self._power_law(x)
-    
-    def evaluate_model_errors(self, x):
-        """
-        Returns the error on I(q) for the given array of q-values
-        :param x: array of q-values
-        """
-        p1 = numpy.array([self.dscale * math.pow(q, -self.power) for q in x])
-        p2 = numpy.array([self.scale * self.power * math.pow(q, -self.power-1)\
-                           * self.dpower for q in x])
-        diq2 = p1*p1 + p2*p2        
-        return numpy.array( [math.sqrt(err) for err in diq2] )
-       
-    def _power_law(self, x):
-        """
-        F(x) = scale* (x)^(-power)
-            when power= 4. the model is porod 
-            else power_law
-        The model has three parameters: ::
-            1. x: a vector of q values
-            2. power: power of the function
-            3. scale : scale factor value
-        
-        :param x: array
-        :return: F(x)
-        """
-        if self.power <= 0:
-            msg = "Power_law function expected positive power,"
-            msg += " but got %s"%self.power
-            raise ValueError(msg)
-        if self.scale <= 0:
-            msg = "scale expected positive value, but got %s"%self.scale
-            raise ValueError(msg) 
-       
-        value = numpy.array([ math.pow(i, -self.power) for i in x ])  
-        return self.scale * value
-
-class Extrapolator:
-    """
-    Extrapolate I(q) distribution using a given model
-    """
-    def __init__(self, data, model=None):
-        """
-        Determine a and b given a linear equation y = ax + b
-        
-        If a model is given, it will be used to linearize the data before 
-        the extrapolation is performed. If None,
-        a simple linear fit will be done.
-        
-        :param data: data containing x and y  such as  y = ax + b 
-        :param model: optional Transform object 
-        """
-        self.data  = data
-        self.model = model
-       
-        # Set qmin as the lowest non-zero value
-        self.qmin = Q_MINIMUM
-        for q_value in self.data.x:
-            if q_value > 0: 
-                self.qmin = q_value
-                break
-        self.qmax = max(self.data.x)
-             
-    def fit(self, power=None, qmin=None, qmax=None):
-        """
-        Fit data for y = ax + b  return a and b
-       
-        :param power: a fixed, otherwise None
-        :param qmin: Minimum Q-value
-        :param qmax: Maximum Q-value
-        """
-        if qmin is None:
-            qmin = self.qmin
-        if qmax is None:
-            qmax = self.qmax
-            
-        # Identify the bin range for the fit
-        idx = (self.data.x >= qmin) & (self.data.x <= qmax)
-        
-        fx = numpy.zeros(len(self.data.x))
-
-        # Uncertainty
-        if type(self.data.dy)==numpy.ndarray and \
-            len(self.data.dy)==len(self.data.x):
-            sigma = self.data.dy
-        else:
-            sigma = numpy.ones(len(self.data.x))
-            
-        # Compute theory data f(x)
-        fx[idx] = self.data.y[idx]
-        
-        # Linearize the data
-        if self.model is not None:
-            linearized_data = self.model.linearize_data(\
-                                            LoaderData1D(self.data.x[idx],
-                                                                fx[idx],
-                                                            dy = sigma[idx]))
-        else:
-            linearized_data = LoaderData1D(self.data.x[idx],
-                                           fx[idx],
-                                           dy = sigma[idx])
-        
-        ##power is given only for function = power_law    
-        if power != None:
-            sigma2 = linearized_data.dy * linearized_data.dy
-            a = -(power)
-            b = (numpy.sum(linearized_data.y/sigma2) \
-                 - a*numpy.sum(linearized_data.x/sigma2))/numpy.sum(1.0/sigma2)
-            
-            
-            deltas = linearized_data.x*a + \
-                    numpy.ones(len(linearized_data.x))*b-linearized_data.y
-            residuals = numpy.sum(deltas*deltas/sigma2)
-            
-            err = math.fabs(residuals) / numpy.sum(1.0/sigma2)
-            return [a, b], [0, math.sqrt(err)]
-        else:
-            A = numpy.vstack([ linearized_data.x/linearized_data.dy,
-                               1.0/linearized_data.dy]).T        
-            (p, residuals, rank, s) = numpy.linalg.lstsq(A,
-                                        linearized_data.y/linearized_data.dy)
-            
-            # Get the covariance matrix, defined as inv_cov = a_transposed * a
-            err = numpy.zeros(2)
-            try:
-                inv_cov = numpy.dot(A.transpose(), A)
-                cov = numpy.linalg.pinv(inv_cov)
-                err_matrix = math.fabs(residuals) * cov
-                err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])]
-            except:
-                err = [-1.0, -1.0]
-                
-            return p, err
-        
-
-class InvariantCalculator(object):
-    """
-    Compute invariant if data is given.
-    Can provide volume fraction and surface area if the user provides
-    Porod constant  and contrast values.
-    
-    :precondition:  the user must send a data of type DataLoader.Data1D
-                    the user provide background and scale values.
-                    
-    :note: Some computations depends on each others. 
-    """
-    def __init__(self, data, background=0, scale=1 ):
-        """
-        Initialize variables.
-        
-        :param data: data must be of type DataLoader.Data1D
-        :param background: Background value. The data will be corrected 
-            before processing
-        :param scale: Scaling factor for I(q). The data will be corrected
-            before processing
-        """
-        # Background and scale should be private data member if the only way to
-        # change them are by instantiating a new object.
-        self._background = background
-        self._scale = scale
-        # slit height for smeared data
-        self._smeared = None
-        # The data should be private
-        self._data = self._get_data(data)
-        # get the dxl if the data is smeared: This is done only once on init.
-        if self._data.dxl != None and self._data.dxl.all() >0:
-            # assumes constant dxl
-            self._smeared = self._data.dxl[0]
-      
-        # Since there are multiple variants of Q*, you should force the
-        # user to use the get method and keep Q* a private data member
-        self._qstar = None
-        
-        # You should keep the error on Q* so you can reuse it without
-        # recomputing the whole thing.
-        self._qstar_err = 0
-        
-        # Extrapolation parameters
-        self._low_extrapolation_npts = 4
-        self._low_extrapolation_function = Guinier()
-        self._low_extrapolation_power = None
-        self._low_extrapolation_power_fitted = None
-   
-        self._high_extrapolation_npts = 4
-        self._high_extrapolation_function = PowerLaw()
-        self._high_extrapolation_power = None
-        self._high_extrapolation_power_fitted = None
-        
-        # Extrapolation range
-        self._low_q_limit = Q_MINIMUM
-        
-    def _get_data(self, data):
-        """
-        :note: this function must be call before computing any type
-         of invariant
-         
-        :return: new data = self._scale *data - self._background
-        """
-        if not issubclass(data.__class__, LoaderData1D):
-            #Process only data that inherited from DataLoader.Data_info.Data1D
-            raise ValueError,"Data must be of type DataLoader.Data1D"
-        #from copy import deepcopy
-        new_data = (self._scale * data) - self._background   
-        
-        # Check that the vector lengths are equal
-        assert(len(new_data.x)==len(new_data.y))
-        
-        # Verify that the errors are set correctly
-        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
-            (min(new_data.dy)==0 and max(new_data.dy)==0):
-            new_data.dy = numpy.ones(len(new_data.x))  
-        return  new_data
-     
-    def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):
-        """
-        fit data with function using 
-        data = self._get_data()
-        fx = Functor(data , function)
-        y = data.y
-        slope, constant = linalg.lstsq(y,fx)
-        
-        :param qmin: data first q value to consider during the fit
-        :param qmax: data last q value to consider during the fit
-        :param power : power value to consider for power-law 
-        :param function: the function to use during the fit
-        
-        :return a: the scale of the function
-        :return b: the other parameter of the function for guinier will be radius
-                for power_law will be the power value
-        """
-        extrapolator = Extrapolator(data=self._data, model=model)
-        p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax) 
-        
-        return model.extract_model_parameters(constant=p[1], slope=p[0],
-                                              dconstant=dp[1], dslope=dp[0])
-    
-    def _get_qstar(self, data):
-        """
-        Compute invariant for pinhole data.
-        This invariant is given by: ::
-    
-            q_star = x0**2 *y0 *dx0 +x1**2 *y1 *dx1 
-                        + ..+ xn**2 *yn *dxn    for non smeared data
-                        
-            q_star = dxl0 *x0 *y0 *dx0 +dxl1 *x1 *y1 *dx1 
-                        + ..+ dlxn *xn *yn *dxn    for smeared data
-                        
-            where n >= len(data.x)-1
-            dxl = slit height dQl
-            dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
-            dx0 = (x1 - x0)/2
-            dxn = (xn - xn-1)/2
-            
-        :param data: the data to use to compute invariant.
-        
-        :return q_star: invariant value for pinhole data. q_star > 0
-        """
-        if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y):
-            msg =  "Length x and y must be equal"
-            msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
-                                                          len(data.y))
-            raise ValueError, msg
-        else:
-            # Take care of smeared data
-            if self._smeared is None:
-                gx = data.x * data.x
-            # assumes that len(x) == len(dxl).
-            else:               
-                gx = data.dxl * data.x 
-                
-            n = len(data.x)- 1
-            #compute the first delta q
-            dx0 = (data.x[1] - data.x[0])/2
-            #compute the last delta q
-            dxn = (data.x[n] - data.x[n-1])/2
-            sum = 0
-            sum += gx[0] * data.y[0] * dx0
-            sum += gx[n] * data.y[n] * dxn
-            
-            if len(data.x) == 2:
-                return sum
-            else:
-                #iterate between for element different
-                #from the first and the last
-                for i in xrange(1, n-1):
-                    dxi = (data.x[i+1] - data.x[i-1])/2
-                    sum += gx[i] * data.y[i] * dxi
-                return sum
-            
-    def _get_qstar_uncertainty(self, data):
-        """
-        Compute invariant uncertainty with with pinhole data.
-        This uncertainty is given as follow: ::
-        
-           dq_star = math.sqrt[(x0**2*(dy0)*dx0)**2 +
-                (x1**2 *(dy1)*dx1)**2 + ..+ (xn**2 *(dyn)*dxn)**2 ]
-        where n >= len(data.x)-1
-        dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
-        dx0 = (x1 - x0)/2
-        dxn = (xn - xn-1)/2
-        dyn: error on dy
-       
-        :param data:
-        :note: if data doesn't contain dy assume dy= math.sqrt(data.y)
-        """          
-        if len(data.x) <= 1 or len(data.y) <= 1 or \
-            len(data.x) != len(data.y) or \
-            (data.dy is not None and (len(data.dy) != len(data.y))):
-            msg = "Length of data.x and data.y must be equal"
-            msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
-                                                         len(data.y))
-            raise ValueError, msg
-        else:
-            #Create error for data without dy error
-            if data.dy is None:
-                dy = math.sqrt(y) 
-            else:
-                dy = data.dy
-            # Take care of smeared data
-            if self._smeared is None:
-                gx = data.x * data.x
-            # assumes that len(x) == len(dxl).
-            else:
-                gx = data.dxl * data.x
-  
-            n = len(data.x) - 1
-            #compute the first delta
-            dx0 = (data.x[1] - data.x[0])/2
-            #compute the last delta
-            dxn= (data.x[n] - data.x[n-1])/2
-            sum = 0
-            sum += (gx[0] * dy[0] * dx0)**2
-            sum += (gx[n] * dy[n] * dxn)**2
-            if len(data.x) == 2:
-                return math.sqrt(sum)
-            else:
-                #iterate between for element different
-                #from the first and the last
-                for i in xrange(1, n-1):
-                    dxi = (data.x[i+1] - data.x[i-1])/2
-                    sum += (gx[i] * dy[i] * dxi)**2
-                return math.sqrt(sum)
-        
-    def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS,
-                              q_start=Q_MINIMUM, q_end=Q_MAXIMUM):
-        """
-        :return: extrapolate data create from data
-        """
-        #create new Data1D to compute the invariant
-        q = numpy.linspace(start=q_start,
-                           stop=q_end,
-                           num=npts,
-                           endpoint=True)
-        iq = model.evaluate_model(q)
-        diq = model.evaluate_model_errors(q)
-         
-        result_data = LoaderData1D(x=q, y=iq, dy=diq)
-        if self._smeared != None:
-            result_data.dxl = self._smeared * numpy.ones(len(q))
-        return result_data
-    
-    def get_data(self):
-        """
-        :return: self._data
-        """
-        return self._data
-    
-    def get_extrapolation_power(self, range='high'):
-        """
-        :return: the fitted power for power law function for a given
-            extrapolation range
-        """
-        if range == 'low':
-            return self._low_extrapolation_power_fitted
-        return self._high_extrapolation_power_fitted
-    
-    def get_qstar_low(self):
-        """
-        Compute the invariant for extrapolated data at low q range.
-        
-        Implementation:
-            data = self._get_extra_data_low()
-            return self._get_qstar()
-            
-        :return q_star: the invariant for data extrapolated at low q.
-        """
-        # Data boundaries for fitting
-        qmin = self._data.x[0]
-        qmax = self._data.x[self._low_extrapolation_npts - 1]
-        
-        # Extrapolate the low-Q data
-        p, dp = self._fit(model=self._low_extrapolation_function,
-      		     	      qmin=qmin,
-                       	  qmax=qmax,
-                          power=self._low_extrapolation_power)
-        self._low_extrapolation_power_fitted = p[0]
-        
-        # Distribution starting point
-        self._low_q_limit = Q_MINIMUM
-        if Q_MINIMUM >= qmin:
-            self._low_q_limit = qmin/10
-        
-        data = self._get_extrapolated_data(\
-                                    model=self._low_extrapolation_function,
-                                            npts=INTEGRATION_NSTEPS,
-                                        q_start=self._low_q_limit, q_end=qmin)
-        
-        # Systematic error
-        # If we have smearing, the shape of the I(q) distribution at low Q will
-        # may not be a Guinier or simple power law. The following is 
-        # a conservative estimation for the systematic error.
-        err = qmin*qmin*math.fabs((qmin-self._low_q_limit)*\
-                                  (data.y[0] - data.y[INTEGRATION_NSTEPS-1]))
-        return self._get_qstar(data), self._get_qstar_uncertainty(data)+err
-        
-    def get_qstar_high(self):
-        """
-        Compute the invariant for extrapolated data at high q range.
-        
-        Implementation:
-            data = self._get_extra_data_high()
-            return self._get_qstar()
-            
-        :return q_star: the invariant for data extrapolated at high q.
-        """
-        # Data boundaries for fitting
-        x_len = len(self._data.x) - 1
-        qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)]
-        qmax = self._data.x[x_len]
-        
-        # fit the data with a model to get the appropriate parameters
-        p, dp = self._fit(model=self._high_extrapolation_function,
-                          qmin=qmin,
-                          qmax=qmax,
-                          power=self._high_extrapolation_power)
-        self._high_extrapolation_power_fitted = p[0]
-        
-        #create new Data1D to compute the invariant
-        data = self._get_extrapolated_data(\
-                                    model=self._high_extrapolation_function,
-                                           npts=INTEGRATION_NSTEPS,
-                                           q_start=qmax, q_end=Q_MAXIMUM)        
-        
-        return self._get_qstar(data), self._get_qstar_uncertainty(data)
-    
-    def get_extra_data_low(self, npts_in=None, q_start=None, npts=20):
-        """
-        Returns the extrapolated data used for the loew-Q invariant calculation.
-        By default, the distribution will cover the data points used for the 
-        extrapolation. The number of overlap points is a parameter (npts_in).
-        By default, the maximum q-value of the distribution will be  
-        the minimum q-value used when extrapolating for the purpose of the 
-        invariant calculation. 
-        
-        :param npts_in: number of data points for which
-            the extrapolated data overlap
-        :param q_start: is the minimum value to uses for extrapolated data
-        :param npts: the number of points in the extrapolated distribution
-           
-        """
-        # Get extrapolation range
-        if q_start is None:
-            q_start = self._low_q_limit
-            
-        if npts_in is None:
-            npts_in = self._low_extrapolation_npts
-        q_end = self._data.x[max(0, npts_in-1)]
-        
-        if q_start >= q_end:
-            return numpy.zeros(0), numpy.zeros(0)
-
-        return self._get_extrapolated_data(\
-                                    model=self._low_extrapolation_function,
-                                           npts=npts,
-                                           q_start=q_start, q_end=q_end)
-          
-    def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, npts=20):
-        """
-        Returns the extrapolated data used for the high-Q invariant calculation.
-        By default, the distribution will cover the data points used for the 
-        extrapolation. The number of overlap points is a parameter (npts_in).
-        By default, the maximum q-value of the distribution will be Q_MAXIMUM, 
-        the maximum q-value used when extrapolating for the purpose of the 
-        invariant calculation. 
-        
-        :param npts_in: number of data points for which the
-            extrapolated data overlap
-        :param q_end: is the maximum value to uses for extrapolated data
-        :param npts: the number of points in the extrapolated distribution
-        """
-        # Get extrapolation range
-        if npts_in is None:
-            npts_in = self._high_extrapolation_npts
-        _npts = len(self._data.x)
-        q_start = self._data.x[min(_npts, _npts-npts_in)]
-        
-        if q_start >= q_end:
-            return numpy.zeros(0), numpy.zeros(0)
-        
-        return self._get_extrapolated_data(\
-                                model=self._high_extrapolation_function,
-                                           npts=npts,
-                                           q_start=q_start, q_end=q_end)
-     
-    def set_extrapolation(self, range, npts=4, function=None, power=None):
-        """
-        Set the extrapolation parameters for the high or low Q-range.
-        Note that this does not turn extrapolation on or off.
-        
-        :param range: a keyword set the type of extrapolation . type string
-        :param npts: the numbers of q points of data to consider
-            for extrapolation
-        :param function: a keyword to select the function to use
-            for extrapolation.
-            of type string.
-        :param power: an power to apply power_low function
-                
-        """
-        range = range.lower()
-        if range not in ['high', 'low']:
-            raise ValueError, "Extrapolation range should be 'high' or 'low'"
-        function = function.lower()
-        if function not in ['power_law', 'guinier']:
-            msg = "Extrapolation function should be 'guinier' or 'power_law'"
-            raise ValueError, msg
-        
-        if range == 'high':
-            if function != 'power_law':
-                msg = "Extrapolation only allows a power law at high Q"
-                raise ValueError, msg
-            self._high_extrapolation_npts  = npts
-            self._high_extrapolation_power = power
-            self._high_extrapolation_power_fitted = power
-        else:
-            if function == 'power_law':
-                self._low_extrapolation_function = PowerLaw()
-            else:
-                self._low_extrapolation_function = Guinier()
-            self._low_extrapolation_npts  = npts
-            self._low_extrapolation_power = power
-            self._low_extrapolation_power_fitted = power
-        
-    def get_qstar(self, extrapolation=None):
-        """
-        Compute the invariant of the local copy of data.
-       
-        :param extrapolation: string to apply optional extrapolation 
-           
-        :return q_star: invariant of the data within data's q range
-        
-        :warning: When using setting data to Data1D ,
-            the user is responsible of
-            checking that the scale and the background are
-            properly apply to the data
-        
-        """
-        self._qstar = self._get_qstar(self._data)
-        self._qstar_err = self._get_qstar_uncertainty(self._data)
-        
-        if extrapolation is None:
-            return self._qstar
-        
-        # Compute invariant plus invariant of extrapolated data
-        extrapolation = extrapolation.lower()    
-        if extrapolation == "low":
-            qs_low, dqs_low = self.get_qstar_low()
-            qs_hi, dqs_hi   = 0, 0
-            
-        elif extrapolation == "high":
-            qs_low, dqs_low = 0, 0
-            qs_hi, dqs_hi   = self.get_qstar_high()
-            
-        elif extrapolation == "both":
-            qs_low, dqs_low = self.get_qstar_low()
-            qs_hi, dqs_hi   = self.get_qstar_high()
-            
-        self._qstar     += qs_low + qs_hi
-        self._qstar_err = math.sqrt(self._qstar_err*self._qstar_err \
-                                    + dqs_low*dqs_low + dqs_hi*dqs_hi)
-        
-        return self._qstar
-       
-    def get_surface(self, contrast, porod_const, extrapolation=None):
-        """
-        Compute the specific surface from the data.
-        
-        Implementation::
-        
-          V =  self.get_volume_fraction(contrast, extrapolation)
-    
-          Compute the surface given by:
-            surface = (2*pi *V(1- V)*porod_const)/ q_star
-           
-        :param contrast: contrast value to compute the volume
-        :param porod_const: Porod constant to compute the surface 
-        :param extrapolation: string to apply optional extrapolation
-        
-        :return: specific surface 
-        """
-        # Compute the volume
-        volume = self.get_volume_fraction(contrast, extrapolation)
-        return 2 * math.pi * volume *(1 - volume) * \
-            float(porod_const)/self._qstar
-        
-    def get_volume_fraction(self, contrast, extrapolation=None):
-        """
-        Compute volume fraction is deduced as follow: ::
-        
-            q_star = 2*(pi*contrast)**2* volume( 1- volume)
-            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2)
-            we get 2 values of volume:
-                 with   1 - 4 * k >= 0
-                 volume1 = (1- sqrt(1- 4*k))/2
-                 volume2 = (1+ sqrt(1- 4*k))/2
-           
-            q_star: the invariant value included extrapolation is applied
-                         unit  1/A^(3)*1/cm
-                    q_star = self.get_qstar()
-                    
-            the result returned will be 0 <= volume <= 1
-        
-        :param contrast: contrast value provides by the user of type float.
-                 contrast unit is 1/A^(2)= 10^(16)cm^(2)
-        :param extrapolation: string to apply optional extrapolation
-        
-        :return: volume fraction
-        
-        :note: volume fraction must have no unit
-        """
-        if contrast <= 0:
-            raise ValueError, "The contrast parameter must be greater than zero"  
-        
-        # Make sure Q star is up to date
-        self.get_qstar(extrapolation)
-        
-        if self._qstar <= 0:
-            msg = "Invalid invariant: Invariant Q* must be greater than zero"
-            raise RuntimeError, msg
-        
-        # Compute intermediate constant
-        k =  1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2)
-        # Check discriminant value
-        discrim = 1 - 4 * k
-        
-        # Compute volume fraction
-        if discrim < 0:
-            msg = "Could not compute the volume fraction: negative discriminant"
-            raise RuntimeError, msg
-        elif discrim == 0:
-            return 1/2
-        else:
-            volume1 = 0.5 * (1 - math.sqrt(discrim))
-            volume2 = 0.5 * (1 + math.sqrt(discrim))
-            
-            if 0 <= volume1 and volume1 <= 1:
-                return volume1
-            elif 0 <= volume2 and volume2 <= 1: 
-                return volume2 
-            msg = "Could not compute the volume fraction: inconsistent results"
-            raise RuntimeError, msg
-    
-    def get_qstar_with_error(self, extrapolation=None):
-        """
-        Compute the invariant uncertainty.
-        This uncertainty computation depends on whether or not the data is
-        smeared.
-        
-        :param extrapolation: string to apply optional extrapolation
-        
-        :return: invariant, the invariant uncertainty
-        """    
-        self.get_qstar(extrapolation)
-        return self._qstar, self._qstar_err
-    
-    def get_volume_fraction_with_error(self, contrast, extrapolation=None):
-        """
-        Compute uncertainty on volume value as well as the volume fraction
-        This uncertainty is given by the following equation: ::
-        
-            dV = 0.5 * (4*k* dq_star) /(2* math.sqrt(1-k* q_star))
-                                 
-            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2)
-            
-            q_star: the invariant value including extrapolated value if existing
-            dq_star: the invariant uncertainty
-            dV: the volume uncertainty
-        
-        The uncertainty will be set to -1 if it can't be computed.
-        
-        :param contrast: contrast value 
-        :param extrapolation: string to apply optional extrapolation
-        
-        :return: V, dV = volume fraction, error on volume fraction
-        """
-        volume = self.get_volume_fraction(contrast, extrapolation)
-        
-        # Compute error
-        k =  1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2)
-        # Check value inside the sqrt function
-        value = 1 - k * self._qstar
-        if (value) <= 0:
-            uncertainty = -1
-        # Compute uncertainty
-        uncertainty = math.fabs((0.5 * 4 * k * \
-                        self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar)))
-        
-        return volume, uncertainty
-    
-    def get_surface_with_error(self, contrast, porod_const, extrapolation=None):
-        """
-        Compute uncertainty of the surface value as well as the surface value.
-        The uncertainty is given as follow: ::
-        
-            dS = porod_const *2*pi[( dV -2*V*dV)/q_star
-                 + dq_star(v-v**2)
-                 
-            q_star: the invariant value
-            dq_star: the invariant uncertainty
-            V: the volume fraction value
-            dV: the volume uncertainty
-        
-        :param contrast: contrast value
-        :param porod_const: porod constant value 
-        :param extrapolation: string to apply optional extrapolation
-        
-        :return S, dS: the surface, with its uncertainty
-        """
-        # We get the volume fraction, with error
-        #   get_volume_fraction_with_error calls get_volume_fraction
-        #   get_volume_fraction calls get_qstar
-        #   which computes Qstar and dQstar
-        v, dv = self.get_volume_fraction_with_error(contrast, extrapolation)
-
-        s = self.get_surface(contrast=contrast, porod_const=porod_const, 
-                             extrapolation=extrapolation)
-        ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\
-                 + self._qstar_err * ( v - v**2))
-
-        return s, ds
Index: sansinvariant/src/sans/invariant/invariant.py
===================================================================
--- sansinvariant/src/sans/invariant/invariant.py	(revision c3af67c7fb93cbb849315ba9e9ba773db5e8b494)
+++ sansinvariant/src/sans/invariant/invariant.py	(revision c3af67c7fb93cbb849315ba9e9ba773db5e8b494)
@@ -0,0 +1,957 @@
+#####################################################################
+#This software was developed by the University of Tennessee as part of the
+#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+#project funded by the US National Science Foundation. 
+#See the license text in license.txt
+#copyright 2010, University of Tennessee
+######################################################################
+
+"""
+This module implements invariant and its related computations.
+
+:author: Gervaise B. Alina/UTK
+:author: Mathieu Doucet/UTK
+:author: Jae Cho/UTK
+
+"""
+import math 
+import numpy
+
+from DataLoader.data_info import Data1D as LoaderData1D
+
+# The minimum q-value to be used when extrapolating
+Q_MINIMUM  = 1e-5
+
+# The maximum q-value to be used when extrapolating
+Q_MAXIMUM  = 10
+
+# Number of steps in the extrapolation
+INTEGRATION_NSTEPS = 1000
+
+class Transform(object):
+    """
+    Define interface that need to compute a function or an inverse
+    function given some x, y 
+    """
+    
+    def linearize_data(self, data):
+        """
+        Linearize data so that a linear fit can be performed. 
+        Filter out the data that can't be transformed.
+        
+        :param data: LoadData1D instance
+        
+        """
+        # Check that the vector lengths are equal
+        assert(len(data.x)==len(data.y))
+        if data.dy is not None:
+            assert(len(data.x)==len(data.dy))
+            dy = data.dy
+        else:
+            dy = numpy.ones(len(data.y))
+            
+        # Transform the data
+        data_points = zip(data.x, data.y, dy)
+
+        output_points = [(self.linearize_q_value(p[0]),
+                          math.log(p[1]),
+                          p[2]/p[1]) for p in data_points if p[0]>0 and \
+                          p[1]>0 and p[2]>0]
+        
+        x_out, y_out, dy_out = zip(*output_points)
+       
+        # Create Data1D object
+        x_out = numpy.asarray(x_out)
+        y_out = numpy.asarray(y_out)
+        dy_out = numpy.asarray(dy_out)
+        linear_data = LoaderData1D(x=x_out, y=y_out, dy=dy_out)
+        
+        return linear_data
+    
+    def get_allowed_bins(self, data):
+        """
+        Goes through the data points and returns a list of boolean values
+        to indicate whether each points is allowed by the model or not.
+        
+        :param data: Data1D object
+        """
+        return [p[0]>0 and p[1]>0 and p[2]>0 for p in zip(data.x, data.y,
+                                                           data.dy)]
+        
+    def linearize_q_value(self, value):
+        """
+        Transform the input q-value for linearization
+        """
+        return NotImplemented
+
+    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
+        """
+        set private member
+        """
+        return NotImplemented
+     
+    def evaluate_model(self, x):
+        """
+        Returns an array f(x) values where f is the Transform function.
+        """
+        return NotImplemented
+    
+    def evaluate_model_errors(self, x):
+        """
+        Returns an array of I(q) errors
+        """
+        return NotImplemented
+    
+class Guinier(Transform):
+    """
+    class of type Transform that performs operations related to guinier 
+    function
+    """
+    def __init__(self, scale=1, radius=60):
+        Transform.__init__(self)
+        self.scale = scale
+        self.radius = radius
+        ## Uncertainty of scale parameter
+        self.dscale  = 0
+        ## Unvertainty of radius parameter
+        self.dradius = 0
+        
+    def linearize_q_value(self, value):
+        """
+        Transform the input q-value for linearization
+        
+        :param value: q-value
+        
+        :return: q*q
+        """
+        return value * value
+    
+    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
+        """
+	    assign new value to the scale and the radius
+    	"""
+        self.scale = math.exp(constant)
+        self.radius = math.sqrt(-3 * slope)
+        
+        # Errors
+        self.dscale = math.exp(constant)*dconstant
+        self.dradius = -3.0/2.0/math.sqrt(-3 * slope)*dslope
+        
+        return [self.radius, self.scale], [self.dradius, self.dscale]
+        
+    def evaluate_model(self, x):
+        """
+        return F(x)= scale* e-((radius*x)**2/3)
+        """
+        return self._guinier(x)
+             
+    def evaluate_model_errors(self, x):
+        """
+        Returns the error on I(q) for the given array of q-values
+        
+        :param x: array of q-values
+        """
+        p1 = numpy.array([self.dscale * math.exp(-((self.radius * q)**2/3)) \
+                          for q in x])
+        p2 = numpy.array([self.scale * math.exp(-((self.radius * q)**2/3))\
+                     * (-(q**2/3)) * 2 * self.radius * self.dradius for q in x])
+        diq2 = p1*p1 + p2*p2        
+        return numpy.array( [math.sqrt(err) for err in diq2] )
+             
+    def _guinier(self, x):
+        """
+        Retrive the guinier function after apply an inverse guinier function
+        to x
+        Compute a F(x) = scale* e-((radius*x)**2/3).
+        
+        :param x: a vector of q values 
+        :param scale: the scale value
+        :param radius: the guinier radius value
+        
+        :return: F(x)
+        """   
+        # transform the radius of coming from the inverse guinier function to a 
+        # a radius of a guinier function
+        if self.radius <= 0:
+            msg = "Rg expected positive value, but got %s"%self.radius
+            raise ValueError(msg)  
+        value = numpy.array([math.exp(-((self.radius * i)**2/3)) for i in x ]) 
+        return self.scale * value
+
+class PowerLaw(Transform):
+    """
+    class of type transform that perform operation related to power_law 
+    function
+    """
+    def __init__(self, scale=1, power=4):
+        Transform.__init__(self)
+        self.scale = scale
+        self.power = power
+   
+    def linearize_q_value(self, value):
+        """
+        Transform the input q-value for linearization
+        
+        :param value: q-value
+        
+        :return: log(q)
+        """
+        return math.log(value)
+    
+    def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0):
+        """
+        Assign new value to the scale and the power 
+        """
+        self.power = -slope
+        self.scale = math.exp(constant)
+        
+        # Errors
+        self.dscale = math.exp(constant)*dconstant
+        self.dpower = -dslope
+        
+        return [self.power, self.scale], [self.dpower, self.dscale]
+        
+    def evaluate_model(self, x):
+        """
+        given a scale and a radius transform x, y using a power_law
+        function
+        """
+        return self._power_law(x)
+    
+    def evaluate_model_errors(self, x):
+        """
+        Returns the error on I(q) for the given array of q-values
+        :param x: array of q-values
+        """
+        p1 = numpy.array([self.dscale * math.pow(q, -self.power) for q in x])
+        p2 = numpy.array([self.scale * self.power * math.pow(q, -self.power-1)\
+                           * self.dpower for q in x])
+        diq2 = p1*p1 + p2*p2        
+        return numpy.array( [math.sqrt(err) for err in diq2] )
+       
+    def _power_law(self, x):
+        """
+        F(x) = scale* (x)^(-power)
+            when power= 4. the model is porod 
+            else power_law
+        The model has three parameters: ::
+            1. x: a vector of q values
+            2. power: power of the function
+            3. scale : scale factor value
+        
+        :param x: array
+        :return: F(x)
+        """
+        if self.power <= 0:
+            msg = "Power_law function expected positive power,"
+            msg += " but got %s"%self.power
+            raise ValueError(msg)
+        if self.scale <= 0:
+            msg = "scale expected positive value, but got %s"%self.scale
+            raise ValueError(msg) 
+       
+        value = numpy.array([ math.pow(i, -self.power) for i in x ])  
+        return self.scale * value
+
+class Extrapolator:
+    """
+    Extrapolate I(q) distribution using a given model
+    """
+    def __init__(self, data, model=None):
+        """
+        Determine a and b given a linear equation y = ax + b
+        
+        If a model is given, it will be used to linearize the data before 
+        the extrapolation is performed. If None,
+        a simple linear fit will be done.
+        
+        :param data: data containing x and y  such as  y = ax + b 
+        :param model: optional Transform object 
+        """
+        self.data  = data
+        self.model = model
+       
+        # Set qmin as the lowest non-zero value
+        self.qmin = Q_MINIMUM
+        for q_value in self.data.x:
+            if q_value > 0: 
+                self.qmin = q_value
+                break
+        self.qmax = max(self.data.x)
+             
+    def fit(self, power=None, qmin=None, qmax=None):
+        """
+        Fit data for y = ax + b  return a and b
+       
+        :param power: a fixed, otherwise None
+        :param qmin: Minimum Q-value
+        :param qmax: Maximum Q-value
+        """
+        if qmin is None:
+            qmin = self.qmin
+        if qmax is None:
+            qmax = self.qmax
+            
+        # Identify the bin range for the fit
+        idx = (self.data.x >= qmin) & (self.data.x <= qmax)
+        
+        fx = numpy.zeros(len(self.data.x))
+
+        # Uncertainty
+        if type(self.data.dy)==numpy.ndarray and \
+            len(self.data.dy)==len(self.data.x):
+            sigma = self.data.dy
+        else:
+            sigma = numpy.ones(len(self.data.x))
+            
+        # Compute theory data f(x)
+        fx[idx] = self.data.y[idx]
+        
+        # Linearize the data
+        if self.model is not None:
+            linearized_data = self.model.linearize_data(\
+                                            LoaderData1D(self.data.x[idx],
+                                                                fx[idx],
+                                                            dy = sigma[idx]))
+        else:
+            linearized_data = LoaderData1D(self.data.x[idx],
+                                           fx[idx],
+                                           dy = sigma[idx])
+        
+        ##power is given only for function = power_law    
+        if power != None:
+            sigma2 = linearized_data.dy * linearized_data.dy
+            a = -(power)
+            b = (numpy.sum(linearized_data.y/sigma2) \
+                 - a*numpy.sum(linearized_data.x/sigma2))/numpy.sum(1.0/sigma2)
+            
+            
+            deltas = linearized_data.x*a + \
+                    numpy.ones(len(linearized_data.x))*b-linearized_data.y
+            residuals = numpy.sum(deltas*deltas/sigma2)
+            
+            err = math.fabs(residuals) / numpy.sum(1.0/sigma2)
+            return [a, b], [0, math.sqrt(err)]
+        else:
+            A = numpy.vstack([ linearized_data.x/linearized_data.dy,
+                               1.0/linearized_data.dy]).T        
+            (p, residuals, rank, s) = numpy.linalg.lstsq(A,
+                                        linearized_data.y/linearized_data.dy)
+            
+            # Get the covariance matrix, defined as inv_cov = a_transposed * a
+            err = numpy.zeros(2)
+            try:
+                inv_cov = numpy.dot(A.transpose(), A)
+                cov = numpy.linalg.pinv(inv_cov)
+                err_matrix = math.fabs(residuals) * cov
+                err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])]
+            except:
+                err = [-1.0, -1.0]
+                
+            return p, err
+        
+
+class InvariantCalculator(object):
+    """
+    Compute invariant if data is given.
+    Can provide volume fraction and surface area if the user provides
+    Porod constant  and contrast values.
+    
+    :precondition:  the user must send a data of type DataLoader.Data1D
+                    the user provide background and scale values.
+                    
+    :note: Some computations depends on each others. 
+    """
+    def __init__(self, data, background=0, scale=1 ):
+        """
+        Initialize variables.
+        
+        :param data: data must be of type DataLoader.Data1D
+        :param background: Background value. The data will be corrected 
+            before processing
+        :param scale: Scaling factor for I(q). The data will be corrected
+            before processing
+        """
+        # Background and scale should be private data member if the only way to
+        # change them are by instantiating a new object.
+        self._background = background
+        self._scale = scale
+        # slit height for smeared data
+        self._smeared = None
+        # The data should be private
+        self._data = self._get_data(data)
+        # get the dxl if the data is smeared: This is done only once on init.
+        if self._data.dxl != None and self._data.dxl.all() >0:
+            # assumes constant dxl
+            self._smeared = self._data.dxl[0]
+      
+        # Since there are multiple variants of Q*, you should force the
+        # user to use the get method and keep Q* a private data member
+        self._qstar = None
+        
+        # You should keep the error on Q* so you can reuse it without
+        # recomputing the whole thing.
+        self._qstar_err = 0
+        
+        # Extrapolation parameters
+        self._low_extrapolation_npts = 4
+        self._low_extrapolation_function = Guinier()
+        self._low_extrapolation_power = None
+        self._low_extrapolation_power_fitted = None
+   
+        self._high_extrapolation_npts = 4
+        self._high_extrapolation_function = PowerLaw()
+        self._high_extrapolation_power = None
+        self._high_extrapolation_power_fitted = None
+        
+        # Extrapolation range
+        self._low_q_limit = Q_MINIMUM
+        
+    def _get_data(self, data):
+        """
+        :note: this function must be call before computing any type
+         of invariant
+         
+        :return: new data = self._scale *data - self._background
+        """
+        if not issubclass(data.__class__, LoaderData1D):
+            #Process only data that inherited from DataLoader.Data_info.Data1D
+            raise ValueError,"Data must be of type DataLoader.Data1D"
+        #from copy import deepcopy
+        new_data = (self._scale * data) - self._background   
+        
+        # Check that the vector lengths are equal
+        assert(len(new_data.x)==len(new_data.y))
+        
+        # Verify that the errors are set correctly
+        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
+            (min(new_data.dy)==0 and max(new_data.dy)==0):
+            new_data.dy = numpy.ones(len(new_data.x))  
+        return  new_data
+     
+    def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None):
+        """
+        fit data with function using 
+        data = self._get_data()
+        fx = Functor(data , function)
+        y = data.y
+        slope, constant = linalg.lstsq(y,fx)
+        
+        :param qmin: data first q value to consider during the fit
+        :param qmax: data last q value to consider during the fit
+        :param power : power value to consider for power-law 
+        :param function: the function to use during the fit
+        
+        :return a: the scale of the function
+        :return b: the other parameter of the function for guinier will be radius
+                for power_law will be the power value
+        """
+        extrapolator = Extrapolator(data=self._data, model=model)
+        p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax) 
+        
+        return model.extract_model_parameters(constant=p[1], slope=p[0],
+                                              dconstant=dp[1], dslope=dp[0])
+    
+    def _get_qstar(self, data):
+        """
+        Compute invariant for pinhole data.
+        This invariant is given by: ::
+    
+            q_star = x0**2 *y0 *dx0 +x1**2 *y1 *dx1 
+                        + ..+ xn**2 *yn *dxn    for non smeared data
+                        
+            q_star = dxl0 *x0 *y0 *dx0 +dxl1 *x1 *y1 *dx1 
+                        + ..+ dlxn *xn *yn *dxn    for smeared data
+                        
+            where n >= len(data.x)-1
+            dxl = slit height dQl
+            dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
+            dx0 = (x1 - x0)/2
+            dxn = (xn - xn-1)/2
+            
+        :param data: the data to use to compute invariant.
+        
+        :return q_star: invariant value for pinhole data. q_star > 0
+        """
+        if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y):
+            msg =  "Length x and y must be equal"
+            msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
+                                                          len(data.y))
+            raise ValueError, msg
+        else:
+            # Take care of smeared data
+            if self._smeared is None:
+                gx = data.x * data.x
+            # assumes that len(x) == len(dxl).
+            else:               
+                gx = data.dxl * data.x 
+                
+            n = len(data.x)- 1
+            #compute the first delta q
+            dx0 = (data.x[1] - data.x[0])/2
+            #compute the last delta q
+            dxn = (data.x[n] - data.x[n-1])/2
+            sum = 0
+            sum += gx[0] * data.y[0] * dx0
+            sum += gx[n] * data.y[n] * dxn
+            
+            if len(data.x) == 2:
+                return sum
+            else:
+                #iterate between for element different
+                #from the first and the last
+                for i in xrange(1, n-1):
+                    dxi = (data.x[i+1] - data.x[i-1])/2
+                    sum += gx[i] * data.y[i] * dxi
+                return sum
+            
+    def _get_qstar_uncertainty(self, data):
+        """
+        Compute invariant uncertainty with with pinhole data.
+        This uncertainty is given as follow: ::
+        
+           dq_star = math.sqrt[(x0**2*(dy0)*dx0)**2 +
+                (x1**2 *(dy1)*dx1)**2 + ..+ (xn**2 *(dyn)*dxn)**2 ]
+        where n >= len(data.x)-1
+        dxi = 1/2*(xi+1 - xi) + (xi - xi-1)
+        dx0 = (x1 - x0)/2
+        dxn = (xn - xn-1)/2
+        dyn: error on dy
+       
+        :param data:
+        :note: if data doesn't contain dy assume dy= math.sqrt(data.y)
+        """          
+        if len(data.x) <= 1 or len(data.y) <= 1 or \
+            len(data.x) != len(data.y) or \
+            (data.dy is not None and (len(data.dy) != len(data.y))):
+            msg = "Length of data.x and data.y must be equal"
+            msg += " and greater than 1; got x=%s, y=%s"%(len(data.x),
+                                                         len(data.y))
+            raise ValueError, msg
+        else:
+            #Create error for data without dy error
+            if data.dy is None:
+                dy = math.sqrt(y) 
+            else:
+                dy = data.dy
+            # Take care of smeared data
+            if self._smeared is None:
+                gx = data.x * data.x
+            # assumes that len(x) == len(dxl).
+            else:
+                gx = data.dxl * data.x
+  
+            n = len(data.x) - 1
+            #compute the first delta
+            dx0 = (data.x[1] - data.x[0])/2
+            #compute the last delta
+            dxn= (data.x[n] - data.x[n-1])/2
+            sum = 0
+            sum += (gx[0] * dy[0] * dx0)**2
+            sum += (gx[n] * dy[n] * dxn)**2
+            if len(data.x) == 2:
+                return math.sqrt(sum)
+            else:
+                #iterate between for element different
+                #from the first and the last
+                for i in xrange(1, n-1):
+                    dxi = (data.x[i+1] - data.x[i-1])/2
+                    sum += (gx[i] * dy[i] * dxi)**2
+                return math.sqrt(sum)
+        
+    def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS,
+                              q_start=Q_MINIMUM, q_end=Q_MAXIMUM):
+        """
+        :return: extrapolate data create from data
+        """
+        #create new Data1D to compute the invariant
+        q = numpy.linspace(start=q_start,
+                           stop=q_end,
+                           num=npts,
+                           endpoint=True)
+        iq = model.evaluate_model(q)
+        diq = model.evaluate_model_errors(q)
+         
+        result_data = LoaderData1D(x=q, y=iq, dy=diq)
+        if self._smeared != None:
+            result_data.dxl = self._smeared * numpy.ones(len(q))
+        return result_data
+    
+    def get_data(self):
+        """
+        :return: self._data
+        """
+        return self._data
+    
+    def get_extrapolation_power(self, range='high'):
+        """
+        :return: the fitted power for power law function for a given
+            extrapolation range
+        """
+        if range == 'low':
+            return self._low_extrapolation_power_fitted
+        return self._high_extrapolation_power_fitted
+    
+    def get_qstar_low(self):
+        """
+        Compute the invariant for extrapolated data at low q range.
+        
+        Implementation:
+            data = self._get_extra_data_low()
+            return self._get_qstar()
+            
+        :return q_star: the invariant for data extrapolated at low q.
+        """
+        # Data boundaries for fitting
+        qmin = self._data.x[0]
+        qmax = self._data.x[self._low_extrapolation_npts - 1]
+        
+        # Extrapolate the low-Q data
+        p, dp = self._fit(model=self._low_extrapolation_function,
+      		     	      qmin=qmin,
+                       	  qmax=qmax,
+                          power=self._low_extrapolation_power)
+        self._low_extrapolation_power_fitted = p[0]
+        
+        # Distribution starting point
+        self._low_q_limit = Q_MINIMUM
+        if Q_MINIMUM >= qmin:
+            self._low_q_limit = qmin/10
+        
+        data = self._get_extrapolated_data(\
+                                    model=self._low_extrapolation_function,
+                                            npts=INTEGRATION_NSTEPS,
+                                        q_start=self._low_q_limit, q_end=qmin)
+        
+        # Systematic error
+        # If we have smearing, the shape of the I(q) distribution at low Q will
+        # may not be a Guinier or simple power law. The following is 
+        # a conservative estimation for the systematic error.
+        err = qmin*qmin*math.fabs((qmin-self._low_q_limit)*\
+                                  (data.y[0] - data.y[INTEGRATION_NSTEPS-1]))
+        return self._get_qstar(data), self._get_qstar_uncertainty(data)+err
+        
+    def get_qstar_high(self):
+        """
+        Compute the invariant for extrapolated data at high q range.
+        
+        Implementation:
+            data = self._get_extra_data_high()
+            return self._get_qstar()
+            
+        :return q_star: the invariant for data extrapolated at high q.
+        """
+        # Data boundaries for fitting
+        x_len = len(self._data.x) - 1
+        qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)]
+        qmax = self._data.x[x_len]
+        
+        # fit the data with a model to get the appropriate parameters
+        p, dp = self._fit(model=self._high_extrapolation_function,
+                          qmin=qmin,
+                          qmax=qmax,
+                          power=self._high_extrapolation_power)
+        self._high_extrapolation_power_fitted = p[0]
+        
+        #create new Data1D to compute the invariant
+        data = self._get_extrapolated_data(\
+                                    model=self._high_extrapolation_function,
+                                           npts=INTEGRATION_NSTEPS,
+                                           q_start=qmax, q_end=Q_MAXIMUM)        
+        
+        return self._get_qstar(data), self._get_qstar_uncertainty(data)
+    
+    def get_extra_data_low(self, npts_in=None, q_start=None, npts=20):
+        """
+        Returns the extrapolated data used for the loew-Q invariant calculation.
+        By default, the distribution will cover the data points used for the 
+        extrapolation. The number of overlap points is a parameter (npts_in).
+        By default, the maximum q-value of the distribution will be  
+        the minimum q-value used when extrapolating for the purpose of the 
+        invariant calculation. 
+        
+        :param npts_in: number of data points for which
+            the extrapolated data overlap
+        :param q_start: is the minimum value to uses for extrapolated data
+        :param npts: the number of points in the extrapolated distribution
+           
+        """
+        # Get extrapolation range
+        if q_start is None:
+            q_start = self._low_q_limit
+            
+        if npts_in is None:
+            npts_in = self._low_extrapolation_npts
+        q_end = self._data.x[max(0, npts_in-1)]
+        
+        if q_start >= q_end:
+            return numpy.zeros(0), numpy.zeros(0)
+
+        return self._get_extrapolated_data(\
+                                    model=self._low_extrapolation_function,
+                                           npts=npts,
+                                           q_start=q_start, q_end=q_end)
+          
+    def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, npts=20):
+        """
+        Returns the extrapolated data used for the high-Q invariant calculation.
+        By default, the distribution will cover the data points used for the 
+        extrapolation. The number of overlap points is a parameter (npts_in).
+        By default, the maximum q-value of the distribution will be Q_MAXIMUM, 
+        the maximum q-value used when extrapolating for the purpose of the 
+        invariant calculation. 
+        
+        :param npts_in: number of data points for which the
+            extrapolated data overlap
+        :param q_end: is the maximum value to uses for extrapolated data
+        :param npts: the number of points in the extrapolated distribution
+        """
+        # Get extrapolation range
+        if npts_in is None:
+            npts_in = self._high_extrapolation_npts
+        _npts = len(self._data.x)
+        q_start = self._data.x[min(_npts, _npts-npts_in)]
+        
+        if q_start >= q_end:
+            return numpy.zeros(0), numpy.zeros(0)
+        
+        return self._get_extrapolated_data(\
+                                model=self._high_extrapolation_function,
+                                           npts=npts,
+                                           q_start=q_start, q_end=q_end)
+     
+    def set_extrapolation(self, range, npts=4, function=None, power=None):
+        """
+        Set the extrapolation parameters for the high or low Q-range.
+        Note that this does not turn extrapolation on or off.
+        
+        :param range: a keyword set the type of extrapolation . type string
+        :param npts: the numbers of q points of data to consider
+            for extrapolation
+        :param function: a keyword to select the function to use
+            for extrapolation.
+            of type string.
+        :param power: an power to apply power_low function
+                
+        """
+        range = range.lower()
+        if range not in ['high', 'low']:
+            raise ValueError, "Extrapolation range should be 'high' or 'low'"
+        function = function.lower()
+        if function not in ['power_law', 'guinier']:
+            msg = "Extrapolation function should be 'guinier' or 'power_law'"
+            raise ValueError, msg
+        
+        if range == 'high':
+            if function != 'power_law':
+                msg = "Extrapolation only allows a power law at high Q"
+                raise ValueError, msg
+            self._high_extrapolation_npts  = npts
+            self._high_extrapolation_power = power
+            self._high_extrapolation_power_fitted = power
+        else:
+            if function == 'power_law':
+                self._low_extrapolation_function = PowerLaw()
+            else:
+                self._low_extrapolation_function = Guinier()
+            self._low_extrapolation_npts  = npts
+            self._low_extrapolation_power = power
+            self._low_extrapolation_power_fitted = power
+        
+    def get_qstar(self, extrapolation=None):
+        """
+        Compute the invariant of the local copy of data.
+       
+        :param extrapolation: string to apply optional extrapolation 
+           
+        :return q_star: invariant of the data within data's q range
+        
+        :warning: When using setting data to Data1D ,
+            the user is responsible of
+            checking that the scale and the background are
+            properly apply to the data
+        
+        """
+        self._qstar = self._get_qstar(self._data)
+        self._qstar_err = self._get_qstar_uncertainty(self._data)
+        
+        if extrapolation is None:
+            return self._qstar
+        
+        # Compute invariant plus invariant of extrapolated data
+        extrapolation = extrapolation.lower()    
+        if extrapolation == "low":
+            qs_low, dqs_low = self.get_qstar_low()
+            qs_hi, dqs_hi   = 0, 0
+            
+        elif extrapolation == "high":
+            qs_low, dqs_low = 0, 0
+            qs_hi, dqs_hi   = self.get_qstar_high()
+            
+        elif extrapolation == "both":
+            qs_low, dqs_low = self.get_qstar_low()
+            qs_hi, dqs_hi   = self.get_qstar_high()
+            
+        self._qstar     += qs_low + qs_hi
+        self._qstar_err = math.sqrt(self._qstar_err*self._qstar_err \
+                                    + dqs_low*dqs_low + dqs_hi*dqs_hi)
+        
+        return self._qstar
+       
+    def get_surface(self, contrast, porod_const, extrapolation=None):
+        """
+        Compute the specific surface from the data.
+        
+        Implementation::
+        
+          V =  self.get_volume_fraction(contrast, extrapolation)
+    
+          Compute the surface given by:
+            surface = (2*pi *V(1- V)*porod_const)/ q_star
+           
+        :param contrast: contrast value to compute the volume
+        :param porod_const: Porod constant to compute the surface 
+        :param extrapolation: string to apply optional extrapolation
+        
+        :return: specific surface 
+        """
+        # Compute the volume
+        volume = self.get_volume_fraction(contrast, extrapolation)
+        return 2 * math.pi * volume *(1 - volume) * \
+            float(porod_const)/self._qstar
+        
+    def get_volume_fraction(self, contrast, extrapolation=None):
+        """
+        Compute volume fraction is deduced as follow: ::
+        
+            q_star = 2*(pi*contrast)**2* volume( 1- volume)
+            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2)
+            we get 2 values of volume:
+                 with   1 - 4 * k >= 0
+                 volume1 = (1- sqrt(1- 4*k))/2
+                 volume2 = (1+ sqrt(1- 4*k))/2
+           
+            q_star: the invariant value included extrapolation is applied
+                         unit  1/A^(3)*1/cm
+                    q_star = self.get_qstar()
+                    
+            the result returned will be 0 <= volume <= 1
+        
+        :param contrast: contrast value provides by the user of type float.
+                 contrast unit is 1/A^(2)= 10^(16)cm^(2)
+        :param extrapolation: string to apply optional extrapolation
+        
+        :return: volume fraction
+        
+        :note: volume fraction must have no unit
+        """
+        if contrast <= 0:
+            raise ValueError, "The contrast parameter must be greater than zero"  
+        
+        # Make sure Q star is up to date
+        self.get_qstar(extrapolation)
+        
+        if self._qstar <= 0:
+            msg = "Invalid invariant: Invariant Q* must be greater than zero"
+            raise RuntimeError, msg
+        
+        # Compute intermediate constant
+        k =  1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2)
+        # Check discriminant value
+        discrim = 1 - 4 * k
+        
+        # Compute volume fraction
+        if discrim < 0:
+            msg = "Could not compute the volume fraction: negative discriminant"
+            raise RuntimeError, msg
+        elif discrim == 0:
+            return 1/2
+        else:
+            volume1 = 0.5 * (1 - math.sqrt(discrim))
+            volume2 = 0.5 * (1 + math.sqrt(discrim))
+            
+            if 0 <= volume1 and volume1 <= 1:
+                return volume1
+            elif 0 <= volume2 and volume2 <= 1: 
+                return volume2 
+            msg = "Could not compute the volume fraction: inconsistent results"
+            raise RuntimeError, msg
+    
+    def get_qstar_with_error(self, extrapolation=None):
+        """
+        Compute the invariant uncertainty.
+        This uncertainty computation depends on whether or not the data is
+        smeared.
+        
+        :param extrapolation: string to apply optional extrapolation
+        
+        :return: invariant, the invariant uncertainty
+        """    
+        self.get_qstar(extrapolation)
+        return self._qstar, self._qstar_err
+    
+    def get_volume_fraction_with_error(self, contrast, extrapolation=None):
+        """
+        Compute uncertainty on volume value as well as the volume fraction
+        This uncertainty is given by the following equation: ::
+        
+            dV = 0.5 * (4*k* dq_star) /(2* math.sqrt(1-k* q_star))
+                                 
+            for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2)
+            
+            q_star: the invariant value including extrapolated value if existing
+            dq_star: the invariant uncertainty
+            dV: the volume uncertainty
+        
+        The uncertainty will be set to -1 if it can't be computed.
+        
+        :param contrast: contrast value 
+        :param extrapolation: string to apply optional extrapolation
+        
+        :return: V, dV = volume fraction, error on volume fraction
+        """
+        volume = self.get_volume_fraction(contrast, extrapolation)
+        
+        # Compute error
+        k =  1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2)
+        # Check value inside the sqrt function
+        value = 1 - k * self._qstar
+        if (value) <= 0:
+            uncertainty = -1
+        # Compute uncertainty
+        uncertainty = math.fabs((0.5 * 4 * k * \
+                        self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar)))
+        
+        return volume, uncertainty
+    
+    def get_surface_with_error(self, contrast, porod_const, extrapolation=None):
+        """
+        Compute uncertainty of the surface value as well as the surface value.
+        The uncertainty is given as follow: ::
+        
+            dS = porod_const *2*pi[( dV -2*V*dV)/q_star
+                 + dq_star(v-v**2)
+                 
+            q_star: the invariant value
+            dq_star: the invariant uncertainty
+            V: the volume fraction value
+            dV: the volume uncertainty
+        
+        :param contrast: contrast value
+        :param porod_const: porod constant value 
+        :param extrapolation: string to apply optional extrapolation
+        
+        :return S, dS: the surface, with its uncertainty
+        """
+        # We get the volume fraction, with error
+        #   get_volume_fraction_with_error calls get_volume_fraction
+        #   get_volume_fraction calls get_qstar
+        #   which computes Qstar and dQstar
+        v, dv = self.get_volume_fraction_with_error(contrast, extrapolation)
+
+        s = self.get_surface(contrast=contrast, porod_const=porod_const, 
+                             extrapolation=extrapolation)
+        ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\
+                 + self._qstar_err * ( v - v**2))
+
+        return s, ds
