source: sasview/Invariant/invariant.py @ 25a2ee3

"""
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
    
    
    TODO: 
        - intro / documentation

"""
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:
            raise ValueError("Rg expected positive value, but got %s"%self.radius)  
        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: 
            @param x: a vector of q values
            @param power: power of the function
            @param scale : scale factor value
            @param F(x)
        """
        if self.power <= 0:
            raise ValueError("Power_law function expected positive power, but got %s"%self.power)
        if self.scale <= 0:
            raise ValueError("scale expected positive value, but got %s"%self.scale) 
       
        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
        
        # The data should be private
        self._data = self._get_data(data)
      
        # 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 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 
                            
            where n >= len(data.x)-1
            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:
            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 += data.x[0] * data.x[0] * data.y[0] * dx0
            sum += data.x[n] * data.x[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 += data.x[i] * data.x[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
                
            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 += (data.x[0] * data.x[0] * dy[0] * dx0)**2
            sum += (data.x[n] * data.x[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 += (data.x[i] * data.x[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)
        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']:
            raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'"
        
        if range == 'high':
            if function != 'power_law':
                raise ValueError, "Extrapolation only allows a power law at high Q"
            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:
            raise RuntimeError, "Invalid invariant: Invariant Q* must be greater than zero"
        
        # 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:
            raise RuntimeError, "Could not compute the volume fraction: negative discriminant"
        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 
            raise RuntimeError, "Could not compute the volume fraction: inconsistent results"
    
    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