def trapez_moms(x_arg, y_arg, mom):
    """ Returns numerically integrated moments """
    if mom == 0:
        ret = trapz(y_arg, x_arg)
    else:
        ret = trapz(y_arg * x_arg**mom, x_arg)
    return ret

def dbltrapz(f, x, y):
    """ Perform a double trapezon integration of a vectorized function

        f(x,y) where x is an array of x values and y is an array of y 
        values to perform the integral over.
    
        Implementation taken from: 
            http://mail.scipy.org/pipermail/scipy-user/2011-February/028592.html

        I don't really understand it, but it seems to work...
        
        For example, if we wanted to integrate some really random function 

            f(x,y) = e^(-x^2*y)

        from 1<x<5 and 1<y<10.

        Using matheamtica, we find that the integral is ~0.09:
        
            N[Integrate[Integrate[Exp[-(x^2*y)], {x, 1, 5}], {y, 1, 10}], 20]

        Using our new function.

        >>> import numpy as np
        >>> f=lambda x,y: np.exp(-(x**2*y))
        >>> int=dbltrapz(f, np.linspace(1,5,1000), np.linspace(1,10,1000))
        >>> print np.allclose(int,0.089071862226039234609, rtol=1e-5, atol=1e-5)
        True
    """
    yy = y[:,np.newaxis]
    xx = x[np.newaxis,:]
    integrand=f(xx,yy)

    return integrate.trapz(integrate.trapz(integrand, yy, axis=0), x, axis=0)

 def cPofZ(self, arr, zx):
     ## hardcoding zmax for the time being, should fix it
     zmax = 1.5
     Ng = len(arr)
     dz = 0.001
     if not hasattr(self, "cnorm"):
         Nx = int(zmax / dz)
         xar = np.linspace(0, zmax, Nx)
         rect = np.zeros((Ng, Nx))
         for i, z in enumerate(xar):
             rect[:, i] = self.PofZ(arr, float(z), dz) / dz
         self.cnorm = trapz(rect, xar, axis=1)
     # for floats
     if type(zx) == type(0.1):
         Nx = int(zx / dz)
         xar = np.linspace(0, zx, Nx)
         rect = np.zeros((Ng, Nx))
         for i, z in enumerate(xar):
             rect[:, i] = self.PofZ(arr, float(z), dz) / dz
         unnormC = trapz(rect, xar, axis=1)
     else:
         # for arrays
         zxm = zx.max()
         Nx = int(zxm / dz)
         xar = np.linspace(0, zxm, Nx)
         rect = np.zeros((Ng, Nx))
         for i, z in enumerate(xar):
             rect[:, i] = self.PofZ(arr, float(z), dz) / dz
             rect[np.where(zx > z), i] = 0.0
         unnormC = trapz(rect, xar, axis=1)
     return unnormC / self.cnorm

    def get_SZ(self, psd, geometry):
        """
        Compute the scattering matrices for the given PSD and geometries.

        Returns:
            The new amplitude (S) and phase (Z) matrices.
        """
        if (self._S_table is None) or (self._Z_table is None):
            raise AttributeError(
                "Initialize or load the scattering table first.")

        if (not isinstance(psd, PSD)) or self._previous_psd != psd:
            self._S_dict = {}
            self._Z_dict = {}
            psd_w = psd(self._psd_D)

            for geom in self.geometries:
                self._S_dict[geom] = \
                    trapz(self._S_table[geom] * psd_w, self._psd_D)
                self._Z_dict[geom] = \
                    trapz(self._Z_table[geom] * psd_w, self._psd_D)

            self._previous_psd = psd

        return (self._S_dict[geometry], self._Z_dict[geometry])

def SF_SD(m, wavelength, dp, ndp, minAngle=0, maxAngle=180, angularResolution=0.5, space='theta', angleMeasure='radians', normalization=None):
#  http://pymiescatt.readthedocs.io/en/latest/forward.html#SF_SD
  _steps = int(1+(maxAngle-minAngle)/angularResolution)
  ndp = coerceDType(ndp)
  dp = coerceDType(dp)
  SL = np.zeros(_steps)
  SR = np.zeros(_steps)
  SU = np.zeros(_steps)
  kwargs = {'minAngle':minAngle,
            'maxAngle':maxAngle,
            'angularResolution':angularResolution,
            'space':space,
            'normalization':None}
  for n,d in zip(ndp,dp):
    measure,l,r,u = ScatteringFunction(m,wavelength,d,**kwargs)
    SL += l*n
    SR += r*n
    SU += u*n
  if normalization in ['n','N','number','particles']:
    _n = trapz(ndp,dp)
    SL /= _n
    SR /= _n
    SU /= _n
  elif normalization in ['m','M','max','MAX']:
    SL /= np.max(SL)
    SR /= np.max(SR)
    SU /= np.max(SU)
  elif normalization in ['t','T','total','TOTAL']:
    SL /= trapz(SL,measure)
    SR /= trapz(SR,measure)
    SU /= trapz(SU,measure)
  return measure,SL,SR,SU

    def integrate(self, t0, tf, weights=None):
        """Integrates the data contained in the integrator from t0 to
        tf, inclusive. Applies weights to the data according to function
        weights(times - t0)
        """
        if t0 < self.times[0] or tf > self.times[-1]:
            return None

        interp = spi.interp1d(x=np.squeeze(self.times),
                              y=np.squeeze(self.vals), axis=0)

        # TODO Clean up using bisect
        istart = next(i for i, x in enumerate(self.times) if x > t0)
        ifinal = next((i for i, x in enumerate(self.times) if x > tf), -1)

        times = [t0] + self.times[istart:ifinal] + [tf]
        ref = [interp(t0)] + self.vals[istart:ifinal] + [interp(tf)]
        times = np.array(times)
        ref = np.asarray(ref)

        if weights is not None:
            w = weights(times - t0)
            ref = ref * w
            z = spt.trapz(y=w, x=times)
        else:
            z = 1.0

        return spt.trapz(y=ref, x=times) / z

def compute_expected_coordinates1D(grid, x0):
    """
    Computes expected value and variance of a 1D grid along its (only) axis.

    :param grid: 1D grid to integrate
    :param x0: coordinates of all points along the integration axis

    :type grid: numpy array
    :type x0: numpy array

    :rtype: float
    :returns:
        * exp_x0 : expected value
        * var_x0 : variance

    """

    # normalize 1D grid
    grid_integral = compute_integral1D(grid, x0)
    prob_x0 = grid / grid_integral

    # get 1D marginals, expected values and variances
    exp_x0 = si.trapz(x0*prob_x0, x=x0)
    var_x0 = si.trapz((x0-exp_x0)*(x0-exp_x0)*prob_x0, x=x0)

    return exp_x0, var_x0

    def get_angular_integrated(self, psd, geometry, property_name):
        if self._angular_table is None:
            raise AttributeError(
                "Initialize or load the table of angular-integrated " + 
                "quantities first."
            )

        psd_w = psd(self._psd_D)

        def sca_xsect(geom):
            return trapz(
                self._angular_table["sca_xsect"][geom] * psd_w, 
                self._psd_D
            )
    
        if property_name == "sca_xsect":
            sca_prop = sca_xsect(geometry)
        elif property_name == "ext_xsect":
            sca_prop = trapz(
                self._angular_table["ext_xsect"][geometry] * psd_w, 
                self._psd_D
            )
        elif property_name == "asym":
            sca_xsect_int = sca_xsect(geometry)
            if sca_xsect_int > 0:
                sca_prop = trapz(
                    self._angular_table["asym"][geometry] * \
                    self._angular_table["sca_xsect"][geometry] * psd_w,  
                    self._psd_D
                )
                sca_prop /= sca_xsect_int
            else:
                sca_prop = 0.0

        return sca_prop

    def _ZZlk(self):
        hist = np.zeros(self.result_bins)
        histLB = np.zeros(self.result_bins)
        self.ZZ_phist = [ np.zeros(self.prob_pts) for i in range(4) ]

        for i in range(self.blocks):
            pp00,pp01,pp10 = \
                    np.mgrid[i*self.prob_pts/self.blocks:(i+1)*self.prob_pts/self.blocks,\
                    0:self.prob_pts,0:self.prob_pts]/float(self.prob_pts)
            plk = self.diag_likelihood(self.N_ZZ, pp00,pp01,pp10)
            hist += np.histogram(pp01+pp10, bins=self.result_bins,
                    range=(0,1), density=False, weights=plk)[0]
            histLB += np.histogram(pp01+pp10-2*np.sqrt(pp00*(1.-pp00-pp01-pp10)), 
                    bins=self.result_bins, range=(0,1), density=False, weights=plk)[0]
            
            for i,pp in enumerate([pp00,pp01,pp10,1.-pp00-pp01-pp10]):
                self.ZZ_phist[i] += np.histogram(pp, bins=self.prob_pts,
                        range=(0,1), density=False, weights=plk)[0]

        I = integrate.trapz(hist, dx=self.ZZ_dx)
        hist /= I
        ILB = integrate.trapz(histLB, dx=self.ZZ_dx)
        histLB /= ILB
        
        for h in self.ZZ_phist:
            Ip = integrate.trapz(h, dx=1./self.prob_pts)
            h /= Ip        
        
        return hist, histLB

    def draw_reabsorption(self, columns, reabsorption=1):
        try:
            f = plt.figure(self.graphnumber_reab)
            plt.xlabel('Wavelengths, nm')
            plt.ylabel('Intensity, a.u.')

            Xe = self.data['Xe']
            Ec = self.data['Ec']
            Ec = Ec / integrate.trapz(Ec, dx=self.step_E)
            Ed = self.data['Ed']
            Ed = Ed / integrate.trapz(Ed, dx=self.step_E) * reabsorption

            if len(Xe) > len(Ec):
                Ec = np.append(Ec, np.zeros(len(Xe)-len(Ec)))
            if len(Xe) > len(Ed):
                Ed = np.append(Ed, np.zeros(len(Xe)-len(Ed)))

            plt.plot(Xe, Ec)
            plt.plot(Xe, Ed)

            f.canvas.set_window_title('Reabsorption (' + str(reabsorption) + '; Xe: Ec, Ed)')
            plt.title('Xe: Ec, Ed')

            f.show()
            self.graphnumber_reab += 1
            self.logger.info('Reabsorption calculation finished')
        except:
            self.logger.info('No reabsortion graph could be drawn')

    def integration(self,is_quadrant=0,use_flux_per_mrad2_or_mm2=1):
        if (self.X is None or self.Y is None):
            raise Exception(" X and Y must be array for integration")
        if self.X.shape != self.Y.shape:
            raise Exception(" X and Y must have the same shape")

        if len(self.X.shape)==2 :
            X = np.linspace(self.X.min(), self.X.max(), self.X.shape[0])
            Y = np.linspace(self.Y.min(), self.Y.max(), self.Y.shape[1])

            res1=integrate.trapz(self.intensity, X)
            res = integrate.trapz(res1, Y)

            # res = integrate.simps(integrate.simps(self.intensity, X), Y)

            # res = self.intensity.sum() * (self.X[1,0] - self.X[0,0]) * (self.Y[0,1] - self.X[0,0]) # regular grid only
        else : # X and Y are 1d array
            if len(self.X) == 1:
                res = self.intensity[0]
            else: # choix arbitraire
                XY = np.zeros_like(self.X)
                for i in range(1, len(self.X)):
                    XY[i] = XY[i-1]+np.sqrt((self.X[i] - self.X[i-1]) * 2 + (self.Y[i] - self.Y[i-1]) ** 2)
                res = np.trapz(self.intensity, XY)

        # Note that the value of flux is in phot/s/0.1%bw/mrad2 (or .../mm2) and
        # our grid is in rad (or m), therefore we must account for this:
        if use_flux_per_mrad2_or_mm2:
            res *= 1e6

        # in case the calculation is for a quadrant, the integral is four times the calculated value
        if is_quadrant:
            res *= 4
        return res

def calculate_variance(beta):
    """
    This function calculates variance of curve beta

    :param beta: numpy ndarray of shape (2,M) of M samples

    :rtype: numpy ndarray
    :return variance: variance

    """
    n, T = beta.shape
    betadot = gradient(beta, 1. / (T - 1))
    betadot = betadot[1]
    normbetadot = zeros(T)
    centroid = calculatecentroid(beta)
    integrand = zeros((n, n, T))
    t = linspace(0, 1, T)
    for i in range(0, T):
        normbetadot[i] = norm(betadot[:, i])
        a1 = (beta[:, i] - centroid)
        a1 = a1.reshape((n, 1))
        integrand[:, :, i] = a1.dot(a1.T) * normbetadot[i]

    l = trapz(normbetadot, t)
    variance = trapz(integrand, t, axis=2)
    variance /= l

    return (variance)

def calculate_energy(alphadot, T=100, k=5):
    """
    calculates energy along path

    :param alphadot: numpy ndarray of shape (2,M) of M samples
    :param T: Number of samples of curve (Default = 100)
    :param k: number of samples along path (Default = 5)

    :rtype: numpy scalar
    :return E: energy

    """
    integrand1 = zeros((k, T))
    integrand2 = zeros(k)

    for i in range(0, k):
        for j in range(1, T):
            tmp = alphadot[:, j, i].T
            integrand1[i, j] = tmp.dot(alphadot[:, j, i])

        integrand2[i] = trapz(integrand1[i, :], linspace(0, 1, T))

    E = 0.5 * trapz(integrand2, linspace(0, 1, k))

    return E

    def read_intensity(self, freq=[0.0], t_start=0.0, t_end=100.0):
        '''
        Calculates the intensity of a certain frequency in the FID. It calculates 
        the Fourier coefficents for that specific frequency.
        
        Parameters
        ----------
        freq (list, float): List of frequencies for which the intensity should be 
        calculated (Frequencies in MHz)

        Returns
        -------
        intensity (1D-array, float): Array of the intensities
        
        Notes
        -----
        none
        '''
        if len(self.fid) > 0:
            fid = slice_fid(self.fid, t_start, t_end)
            intensity = np.array([])
            for x in freq:
                cos2 = np.cos(fid[:,0]*x*1.E6*2.*np.pi)
                sin2 = np.sin(fid[:,0]*x*1.E6*2.*np.pi)
        
                a = inte.trapz(cos2*fid[:,1], fid[:,0]) / abs(fid[-1,0]) * 2.
                b = inte.trapz(sin2*fid[:,1], fid[:,0]) / abs(fid[-1,0]) * 2.
       
                intensity = np.hstack((intensity, np.sqrt(a**2 + b**2)))
            
            return intensity
                
        else:
            return np.array([])

def calculate_MOC_initial_condition(number_of_classes, U=0.177, D=0.03):
    data = genfromtxt("paper_data/re6400/upstream.csv")

    # Data provided are diameters in mm. We need to covert to volume
    mm2m = 0.001
    v_data = pi / 6 * (data[:, 0] * mm2m)**3
    p_data = data[:, 1]

    p_interpolated = interp1d(
        v_data, p_data, kind='nearest', fill_value=0, bounds_error=False)

    # From Galinat we know that phase volumetric fraction is 1.7%-3%
    alpha = (0.03 + 0.017) / 2
    A = pi / 4 * D**2
    Vc = 1.0 * A  # Assume unit height column

    # Estimate the total number from total concentration
    v = linspace(0, 1.2 * max(v_data), number_of_classes + 1)
    dv = v[1]
    p_empirical = p_interpolated(v) / trapz(p_interpolated(v), x=v)
    cdf_empirical = cumsum(p_empirical) * dv  # Cumulative distribution

    meanVd = trapz(v * p_empirical, x=v)  # Dispersed phase mean volume
    N0 = alpha * Vc / meanVd  # Number of drops in 1m column

    # print("Total number of drops in unit column {0:0.0f}".format(N0))
    # print("Particle influx will be: {0:0.2f} drops/s".format(U * N0))
    # print(U * N0 * meanVd / (U * A))

    return dv, clip(
        N0 * (cdf_empirical[1:] - cdf_empirical[:-1]),
        a_min=1e-4, a_max=None
    )

    def __init__(self, waveLengths, transmission):
        """Initialise the filter transmission function

        @param waveLengths  
        @param transmission

        """        
 
        # keep track of original range and resolution of transmission function 
        self.lamMin = waveLengths[0]
        self.lamMax = waveLengths[-1]
        self.nLam = len(waveLengths)   
		#this could come handy if we need to defer the construction of the interpolator
        self.wavelengths = waveLengths
        self.transmission = transmission
        #cache these values for later reuse. Note that this will fail if the filter
        #support is disconnected... need to build a safer interface?
        self.integral1 = integ.trapz(transmission*waveLengths,waveLengths)
        self.integral2 = integ.trapz(transmission/waveLengths,waveLengths)
        self.effLambda = np.sqrt(self.integral1/self.integral2)
        
        # check wavelength and transmission domains are valid
        if (self.lamMin<0.):
            raise ValueError("ERROR! wavelengths cannot be less than zero!")
            
        if np.any(transmission<0.):
            raise ValueError("ERROR! cannot have negative transmission")
           
        if np.any(transmission>1.000001): #why not just 1.?
            raise ValueError("ERROR! cannot have transmission > 1") 

        # filter is now represeted as an interpolation object (linear)
        self.filt = interp.InterpolatedUnivariateSpline(waveLengths, transmission, k=1)

 def computeColor(self, filtX, filtY, z=0):
     """Compute color (flux in filter X - filter Y) of SED at redshift z, return color in magnitudes
     
        @param filtX   lower wavelength filter
        @param filtY   higher wavelength filter
        @param z       redshift of SED
     """
     
     if filtX not in self.filterDict:
         emsg = "Filter " + filtX + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtY not in self.filterDict:
         emsg = "Filter " + filtY + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtX == filtY:
         raise ValueError("ERROR! cannot have color as difference between same filter")
     
     # define integral limits by filter extents
     aX, bX = self.filterDict[filtX].returnFilterRange()
     aY, bY = self.filterDict[filtY].returnFilterRange()
     
     
     # int S(lam_obs)*X(lam)*lam dlam
     lam = self.filterDict[filtX].wavelengths
     res = self.sed.getFlux(lam, z)*self.filterDict[filtX].getTrans(lam)*lam
     integrand = interp.InterpolatedUnivariateSpline(lam, res, k=1)
     if self.FAST_INTEG:
         x = np.linspace(aX, bX, self.INTEG_PREC)
         y = integrand(x)
         int1 = integ.trapz(y,x)
     else:
         int1 = integ.quad(integrand, aX, bX)[0]
     
     
     
     # int S(lam_obs)*Y(lam)*lam dlam
     lam = self.filterDict[filtY].wavelengths
     res = self.sed.getFlux(lam, z)*self.filterDict[filtY].getTrans(lam)*lam
     integrand = interp.InterpolatedUnivariateSpline(lam, res, k=1)
     if self.FAST_INTEG:
         x = np.linspace(aY, bY, self.INTEG_PREC)
         y = integrand(x)
         int2 = integ.trapz(y,x)
     else:
         int2 = integ.quad(integrand, aY, bY)[0]
     
     
     # Not sure about this zero-point term? But colors are totally off without it
     # .integral2 = int filter(lam)/lam dlam
     int3 = self.filterDict[filtX].integral2
     int4 = self.filterDict[filtY].integral2
     
     zp = -2.5*math.log10(int4/int3)
     
     
     # zero observed flux in either filter so color should be infinite
     if (int1==0. or int2==0.):
         return float("inf")
         
     return -2.5*math.log10(int1/int2) + zp

 def convolve(self,filt,returnMag=True):
     if not filt.type=='filter':
         raise TypeError('The "filter" you specified is invalid.')
     # Find the part of the spectrum where the FTC is defined
     startN,stopN = np.searchsorted(self.wavelengths.value,[filt.wavelengths[0].value,filt.wavelengths[-1].value])            
     if startN!=0:
         startN-=1
     if stopN!=(len(self.wavelengths)-1):
         stopN+=1
     ls = self.wavelengths[startN:stopN].value
     # Resample filter
     newFtc = np.interp(ls,filt.wavelengths.value,filt.ftc)
     num = interp1d(ls,newFtc * self.spec[startN:stopN])
     denom = interp1d(ls,newFtc)
     # Integrate
     numerator = trapz(num(filt.wavelengths[0:-1].value),x=filt.wavelengths[0:-1].value)
     denominator = trapz(denom(filt.wavelengths[0:-1].value),x=filt.wavelengths[0:-1].value)
    
     if returnMag:
         try:
             ABmag = -2.5 * np.log10(numerator/denominator) - 48.60
         except:
             global badFilt
             badFilt = filt
             
             global badspec
             badspec = self
             
             ABmag = -2.5 * np.log10(numerator/denominator) - 48.60
         if np.isnan(ABmag):
             print(self.params)
         return(ABmag*u.Unit('mag'))
     else:
         return(u.Unit('erg') * u.Unit('s')**-1 * u.Unit('cm')**-2 * u.Unit('Hz')**-1\
             *numerator/denominator) # erg s−1 cm−2 Hz−1

def compute_integral4D(grid, x0, x1, x2, x3):
    """
    Computes norm of a 4D grid

    :param grid: 4D grid to integrate
    :param x0: coordinates of all points along axis=0
    :param x1: coordinates of all points along axis=1
    :param x2: coordinates of all points along axis=2
    :param x3: coordinates of all points along axis=3

    :type grid: numpy array
    :type x0: numpy array
    :type x1: numpy array
    :type x2: numpy array
    :type x3: numpy array

    :rtype: float
    :returns: grid norm

    """

    grid_norm = si.trapz(si.trapz(si.trapz(si.trapz(grid, x=x0, axis=0),
                                           x=x1, axis=0), x=x2, axis=0),
                         x=x3, axis=0)
    return grid_norm

def find_basis_normal(q):
    """
    Finds the basis normal to the srvf

    :param q1: numpy ndarray of shape (2,M) of M samples

    :rtype: list of numpy ndarray
    :return basis: list containing basis vectors

    """
    n, T = q.shape

    f1 = zeros((n, T))
    f2 = zeros((n, T))
    for i in range(0, T):
        f1[:, i] = q[0, i] * q[:, i] / norm(q[:, i]) + array([norm(q[:, i]),
                                                              0])
        f2[:, i] = q[1, i] * q[:, i] / norm(q[:, i]) + array([0,
                                                              norm(q[:, i])])

    h3 = f1
    h4 = f2
    integrandb3 = zeros(T)
    integrandb4 = zeros(T)
    for i in range(0, T):
        a = q[:, i].T
        integrandb3[i] = a.dot(h3[:, i])
        integrandb4[i] = a.dot(h4[:, i])

    b3 = h3 - q * trapz(integrandb3, linspace(0, 1, T))
    b4 = h4 - q * trapz(integrandb4, linspace(0, 1, T))

    basis = [b3, b4]

    return (basis)