Python scipy 模块，integrate() 实例源码

def schechter_cdf(m,A=1,beta=2,m0=100,mmin=10,mmax=None,npts=1e4):
        """
        Return the CDF value of a given mass for a set mmin,mmax
        mmax will default to 10 m0 if not specified

        Analytic integral of the Schechter function:
        http://www.wolframalpha.com/input/?i=integral%28x^-a+exp%28-x%2Fm%29+dx%29
        """
        if mmax is None:
            mmax = 10*m0

        # integrate the CDF from the minimum to maximum 
        # undefined posint = -m0 * mmax**-beta * (mmax/m0)**beta * scipy.special.gammainc(1-beta, mmax/m0)
        # undefined negint = -m0 * mmin**-beta * (mmin/m0)**beta * scipy.special.gammainc(1-beta, mmin/m0)
        posint = -mmax**(1-beta) * scipy.special.expn(beta, mmax/m0)
        negint = -mmin**(1-beta) * scipy.special.expn(beta, mmin/m0)
        tot = posint-negint

        # normalize by the integral
        # undefined ret = (-m0 * m**-beta * (m/m0)**beta * scipy.special.gammainc(1-beta, m/m0)) / tot
        ret = (-m**(1-beta) * scipy.special.expn(beta, m/m0) - negint)/ tot

        return ret

def pan_tompkins(self,data_buf):
        '''
        1) 3rd differential of the data_buffer (512 samples)
        2) square of that
        3) integrate -> integrate(data_buffer)
        4) take the pulse from that 
        '''
        order = 3
        diff = np.diff(data_buf,3)
        for i in range(order):
            diff = np.insert(diff,0,0)
        square = np.square(diff)
        window = 64
        integral = np.zeros((len(square)))
        for i in range(len(square)):
            integral[i] = np.trapz(square[i:i+window])
        # print(integral)
        return integral

def integrate(self, rmin=0, rmax=numpy.inf):
        """
        Calculate the 2D integral of the 1D surface brightness profile 
        (i.e, the flux) between rmin and rmax (elliptical radii). 

        Parameters:
        -----------
        rmin : minimum integration radius (deg)
        rmax : maximum integration radius (deg)

        Returns:
        --------
        integral : Solid angle integral (deg^2)
        """
        if rmin < 0: raise Exception('rmin must be >= 0')
        integrand = lambda r: self._pdf(r) * 2*numpy.pi * r
        return scipy.integrate.quad(integrand,rmin,rmax,full_output=True,epsabs=0)[0]

def int_imf_dm(m1,m2,m,imf,bywhat='bymass',integral='normal'):
    ''' 
    Integrate IMF between m1 and m2.

    Parameters
    ----------
    m1 : float
        Min mass
    m2 : float
        Max mass
    m : float
        Mass array
    imf : float
        IMF array
    bywhat : string, optional
        'bymass' integrates the mass that goes into stars of
        that mass interval; or 'bynumber' which integrates the number
        of stars in that mass interval.  The default is 'bymass'.
    integrate : string, optional
        'normal' uses sc.integrate.trapz; 'cum' returns cumulative
        trapezoidal integral.  The default is 'normal'.

    '''


    ind_m = (m >= min(m1,m2)) & (m <= max(m1,m2))
    if integral is 'normal':
        int_func = sc.integrate.trapz
    elif integral is 'cum':
        int_func = sc.integrate.cumtrapz
    else:
        print("Error in int_imf_dm: don't know how to integrate")
        return 0

    if bywhat is 'bymass':
        return int_func(m[ind_m]*imf[ind_m],m[ind_m])
    elif bywhat is 'bynumber':
        return int_func(imf[ind_m],m[ind_m])
    else:
        print("Error in int_imf_dm: don't know by what to integrate")
        return 0

def __init__(self, sensor, layer):

        # Set size of phase matrix: active needs an extended phase matrix
        if sensor.mode == 'P':
            self.npol = 2
            self.m_max = 0
        else:
            self.npol = 3
            self.m_max = 3

        # Bring layer and sensor properties into emmodel
        self.frac_volume = layer.frac_volume
        self.microstructure = layer.microstructure  # Do this here, so can pass FT of correlation fn to phase function
        self.e0 = layer.permittivity(0, sensor.frequency)  # background permittivity
        self.eps = layer.permittivity(1, sensor.frequency)  # scatterer permittivity
        self.k0 = 2 * np.pi * sensor.frequency / C_SPEED  # Wavenumber in free space

        # Calculate depolarization factors and iba_coefficient
        self.depol_xyz = depolarization_factors()
        self._effective_permittivity = self.effective_permittivity()
        self.iba_coeff = self.compute_iba_coeff()

        # Absorption coefficient for general lossy medium under assumption of low-loss medium.
        self.ka = self.compute_ka()

        # Calculate scattering coefficient: integrate p11+p12 over mu
        k = 6  # number of samples. This should be adaptative depending on the size/wavelength
        mu = np.linspace(1, -1, 2**k + 1)
        y = self.ks_integrand(mu)
        ks_int = scipy.integrate.romb(y, mu[0] - mu[1])  # integrate between 0 and pi (i.e. mu between -1 and 1)
        self.ks = ks_int / 4.  # Ding et al. (2010), normalised by (1/4pi)

        assert(self.ks >= 0)

def __init__(self, sensor, layer):
        IBA.__init__(self, sensor, layer)  # Gives all IBA parameters. Some need to be recalculated (effective permittivity, scattering and absorption coefficients)

        self._effective_permittivity = polder_van_santen(self.frac_volume)

        # Imaginary component for effective permittivity from Wiesmann and Matzler (1999)
        y2 = self.mean_sq_field_ratio(self.e0, self.eps)
        effective_permittivity_imag = self.frac_volume * self.eps.imag * y2 * np.sqrt(self._effective_permittivity)
        self._effective_permittivity = self._effective_permittivity + 1j * effective_permittivity_imag

        self.iba_coeff = self.compute_iba_coeff()
        ks_int, ks_err = scipy.integrate.quad(self._mm_integrand, 0, np.pi)
        self.ks = ks_int / 2.  # Matzler and Wiesmann, RSE, 1999, eqn (8)
        # General lossy medium under assumption of low-loss medium.
        self.ka = self.compute_ka()

def integrate(self, params, sel_dist, theta):
        """
        integration without re-normalizing the DFE. This assumes the
        portion of the DFE that is not integrated is not seen in your
        sample. 
        """
        #need to include tuple() here to make this function play nice
        #with numpy arrays
        sel_args = (self.gammas,) + tuple(params)
        #compute weights for each fs
        weights = sel_dist(*sel_args)

        #compute weight for the effectively neutral portion. not using
        #CDF function because I want this to be able to compute weight
        #for arbitrary mass functions
        weight_neu, err_neu = scipy.integrate.quad(sel_dist, self.gammas[-1], 
                                                   0, args=tuple(params))

        #function's adaptable for demographic models from 1-3 populations
        pops = len(self.neu_spec.shape)
        if pops == 1:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis]*self.spectra, self.gammas, axis=0)
        elif pops == 2:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis]*self.spectra, 
                self.gammas, axis=0)
        elif pops == 3:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis,numpy.newaxis]*self.spectra,
                self.gammas, axis=0)
        else:
            raise IndexError("Must have one to three populations")

        integrated_fs = Spectrum(integrated, extrap_x=self.extrap_x)

        #no normalization, allow lethal mutations to fall out
        return integrated_fs * theta

def integrate_norm(self, params, sel_dist, theta):
        """
        """
        #need to include tuple() here to make this function play nice
        #with numpy arrays
        #compute weights for each fs
        sel_args = (self.gammas,) + tuple(params)
        weights = sel_dist(*sel_args)

        #compute weight for the effectively neutral portion. not using
        #CDF function because I want this to be able to compute weight
        #for arbitrary mass functions
        weight_neu, err_neu = scipy.integrate.quad(sel_dist, self.gammas[-1], 
                                                   0, args=tuple(params))

        #function's adaptable for demographic models from 1-3
        #populations but this assumes the selection coefficient is the
        #same in both populations
        pops = len(self.neu_spec.shape)
        if pops == 1:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis]*self.spectra, self.gammas, axis=0)
        elif pops == 2:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis]*self.spectra, 
                self.gammas, axis=0)
        elif pops == 3:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis,numpy.newaxis]*self.spectra,
                self.gammas, axis=0)
        else:
            raise IndexError("Must have one to three populations")

        integrated_fs = Spectrum(integrated, extrap_x=self.extrap_x)

        #normalization
        dist_int = Numerics.trapz(weights, self.gammas) + weight_neu
        return integrated_fs/dist_int * theta


#define a bunch of default distributions just to make everything easier

def integrate(self, params, sel_dist, theta):
        """
        integration without re-normalizing the DFE. This assumes the
        portion of the DFE that is not integrated is not seen in your
        sample. 
        """
        #need to include tuple() here to make this function play nice
        #with numpy arrays
        sel_args = (self.gammas,) + tuple(params)
        #compute weights for each fs
        weights = sel_dist(*sel_args)

        #compute weight for the effectively neutral portion. not using
        #CDF function because I want this to be able to compute weight
        #for arbitrary mass functions
        weight_neu, err_neu = scipy.integrate.quad(sel_dist, self.gammas[-1], 
                                                   0, args=tuple(params))

        #function's adaptable for demographic models from 1-3 populations
        pops = len(self.neu_spec.shape)
        if pops == 1:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis]*self.spectra, self.gammas, axis=0)
        elif pops == 2:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis]*self.spectra, 
                self.gammas, axis=0)
        elif pops == 3:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis,numpy.newaxis]*self.spectra,
                self.gammas, axis=0)
        else:
            raise IndexError("Must have one to three populations")

        integrated_fs = Spectrum(integrated, extrap_x=self.extrap_x)

        #no normalization, allow lethal mutations to fall out
        return integrated_fs * theta

def integrate_norm(self, params, sel_dist, theta):
        """
        """
        #need to include tuple() here to make this function play nice
        #with numpy arrays
        #compute weights for each fs
        sel_args = (self.gammas,) + tuple(params)
        weights = sel_dist(*sel_args)

        #compute weight for the effectively neutral portion. not using
        #CDF function because I want this to be able to compute weight
        #for arbitrary mass functions
        weight_neu, err_neu = scipy.integrate.quad(sel_dist, self.gammas[-1], 
                                                   0, args=tuple(params))

        #function's adaptable for demographic models from 1-3
        #populations but this assumes the selection coefficient is the
        #same in both populations
        pops = len(self.neu_spec.shape)
        if pops == 1:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis]*self.spectra, self.gammas, axis=0)
        elif pops == 2:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis]*self.spectra, 
                self.gammas, axis=0)
        elif pops == 3:
            integrated = self.neu_spec*weight_neu + Numerics.trapz(
                weights[:,numpy.newaxis,numpy.newaxis,numpy.newaxis]*self.spectra,
                self.gammas, axis=0)
        else:
            raise IndexError("Must have one to three populations")

        integrated_fs = Spectrum(integrated, extrap_x=self.extrap_x)

        #normalization
        dist_int = Numerics.trapz(weights, self.gammas) + weight_neu
        return integrated_fs/dist_int * theta


#define a bunch of default distributions just to make everything easier

def integrate(self):
        if self._integrate is None:
            try:
                import scipy.integrate as integrate
            except ImportError:
                integrate = NotAModule(self._name)
            self._integrate = integrate
        return self._integrate

def _scipy_skip(self):
        if SCIPY_INT is None:  # pragma: NO COVER
            self.skipTest('SciPy not installed')

def integrate(self, mlow, mhigh, **kwargs):
        """
        Integrate the mass function over some range
        """
        import scipy.integrate
        return scipy.integrate.quad(self, mlow, mhigh, **kwargs)

def __init__(self, breaks, mmin, mmax):
        self.breaks = breaks
        self.normalization = self.integrate(mmin, mmax)[0]

def integrate(fn=kroupa, bins=np.logspace(-2,2,500)):
    xax = (bins[:-1]+bins[1:])/2.
    integral = (bins[1:]-bins[:-1]) * (fn(bins[:-1])+fn(bins[1:])) / 2.

    return xax,integral

def cumint(fn=kroupa, bins=np.logspace(-2,2,500)):
    xax,integral = integrate(fn,bins)
    return integral.cumsum() / integral.sum()

def integratedTimeOverlap(self, otherTimeNuclearWF):
        return scipy.integrate.simps(self.timeOverlap(otherTimeNuclearWF), dx = self.mySpace.dt)

def integrate(self, dt):
        "Integrates Amplitude to desired dt"
        DT = dt 
        tValues = np.arange(self.minTime(), self.maxTime() + DT, DT)
        Values = self.valueAtTime(tValues)
        return scipy.integrate.simps(Values, dx = DT)

def integrateInSimulationSpace(self):
        "Integrates Amplitude to mySpace.dt"
        return self.integrate(self.mySpace.dt)

def integrate(self, a):
        return 0.0

def totalPower(self):
        "Integral of the absolute value squared of the amplitude"
        DT = self.mySpace.dt #/ 10.0
        tValues = np.arange(self.minTime(), self.maxTime(), DT)
        Values = np.abs(self.valueAtTime(tValues))**2.0
        return scipy.integrate.simps(Values, dx = DT)

def totalIntegral(self):
        tValues = self.natural_time
        DT = tValues[1] - tValues[0]
        Values = np.abs(self.myFunction(self.natural_time))
        return scipy.integrate.simps(Values, dx = DT)

def totalPower(self):
        tValues = self.natural_time
        DT = tValues[1] - tValues[0]
        Values = np.abs(self.myFunction(self.natural_time))**2
        return scipy.integrate.simps(Values, dx = DT)

def xIntegration(self, Amplitude):
        "Integrates a given amplitude in x-space"
        toBeIntegrated = Amplitude
        for i in range(self.nuclearDimensionality):
            toBeIntegrated = scipy.integrate.simps(toBeIntegrated, dx=self.Dx_values[i])
        return toBeIntegrated

def kIntegration(self, Amplitude):
        "Integrates a given amplitude in k-space"
        toBeIntegrated = Amplitude
        for i in range(self.nuclearDimensionality):
            toBeIntegrated = scipy.integrate.simps(toBeIntegrated, dx=self.Dk_values[i])
        return toBeIntegrated

def tIntegration(self, timeAmplitude):
        "Integrates a given amplitude in time along the first axis"
        toBeIntegrated = timeAmplitude
        toBeIntegrated = scipy.integrate.simps(toBeIntegrated, dx=self.dt)
        return toBeIntegrated



    #@staticmethod

def timeOverlapWithOtherBraEWFOfPolarizationAndOverlapWithElectricField(self, braEWF, dipoleTuple, electricFieldTuple):
        """This is the workhorse function which treats self as the ket, and calculates
        the overlap with the supplied wavefunction after applying the dipole operator
        then takes the dot product with the supplied electric field and integrates"""
        output = []
        for i, EWF in enumerate(self):
            xVal = EWF.overlap(dipoleTuple[0] * braEWF[i]) * electricFieldTuple[0].valueAtTime(self.timeSeries[i])
            yVal = EWF.overlap(dipoleTuple[1] * braEWF[i]) * electricFieldTuple[1].valueAtTime(self.timeSeries[i])
            zVal = EWF.overlap(dipoleTuple[2] * braEWF[i]) * electricFieldTuple[2].valueAtTime(self.timeSeries[i])
            output.append(xVal + yVal + zVal)
        return scipy.integrate.simps(output, dx = self.mySpace.dt)

def propagate_ODE(psi_0,time,E_rot,E_0_squared_max,sigma,V0,V1,V2,abstol=1e-8,reltol=1e-8):
     ''' Same as propagate, except it uses an ODE solver instead
     of the matrix method.

     psi_0: initial wave function
     time: array of timesteps to integrate. psi_0 is given at time[0].
           The time step must be constant!
     E_rot: array of rotational energies for each J.
     E_0_squared_max: Max value of the square of the E field during pulse
     sigma: temporal variance of the gaussian
     V0,V1,V2: the three bands of the symmetric 5 diagonal interaction matrix,
               V0 being the diagonal.
     abstol, reltol: Error tolerances used for the ODE solver.

     '''
     if (numpy.any(numpy.abs(numpy.diff(numpy.diff(time)))>1e-20)):
         raise RuntimeError("Pulse time steps must be equidistant.");
     dt = time[1]-time[0];

     psi_t = numpy.empty((len(time),len(psi_0)), dtype=numpy.complex)
     psi_t[0,:] = psi_0;
     try:
         res = libpropagation.propagate_field_ODE(len(time),len(psi_0),time[0],dt,E_0_squared_max,sigma,psi_t,V0,V1,V2,E_rot, abstol, reltol);
     except:
         raise RuntimeError("For ODE propagation, you need the libpropagation C library, compiled with GSL support.");

     if (res != 0):
         raise RuntimeError("Basis size too small");

     return psi_t;

def norm(self):
        # Numerically integrate the pdf
        return 1./self.integrate()

def _cache(self, name=None):
        if name in [None,'extension','ellipticity','truncate']:
            self._norm = 1./self.integrate() * 1./self.jacobian
        else:
            return

def _cache(self, name=None):
        if name in ['extension','ellipticity','truncate']:
            self._norm = 1./self.integrate()
        else:
            return

def propagate(psi_0,time,E_rot,E_0_squared_max,sigma,eig,vec):
     ''' Propagate the wave function psi_0 (given as coefficients)
     through a gaussian laser pulse centered at time t = 0.

     psi_0: initial wave function
     time: array of timesteps to integrate. psi_0 is given at time[0].
           The time step must be constant!
     E_rot: array of rotational energies for each J.
     E_0_squared_max: Max value of the square of the E field during pulse
     sigma: temporal variance of the gaussian
     eig, vec: diagonalized interaction matrix V = vec * diag(eig) * vec.T

     '''
     if (numpy.any(numpy.abs(numpy.diff(numpy.diff(time)))>1e-20)):
         raise RuntimeError("Pulse time steps must be equidistant.");

     # Within the maximum time step, a Gaussian is approximately constant.
     if (len(E_rot) > 3):
        dt_max = min(2.4*sigma/150,1/(E_rot[-1]-E_rot[-3]));
     else:
        dt_max = min(2.4*sigma/150,1/E_rot[-1]);
     dt = time[1]-time[0];
     # If our given time step is larger than this, use a smaller time step
     if (dt > dt_max):
         scale = int(numpy.ceil(dt/dt_max));
         dt = dt/scale;
     else:
         scale = 1;

     expRot = numpy.exp(-1j*dt*E_rot);
     expRot2 = numpy.exp(-1j*dt*E_rot/2);

     psi_t = numpy.empty((len(time),len(psi_0)), dtype=numpy.complex)
     psi_t[0,:] = psi_0;
     vecT = numpy.ascontiguousarray(vec.T);
     if (libpropagation):
         res = libpropagation.propagate_field(len(time),scale,len(psi_0),time[0],dt,E_0_squared_max,sigma,psi_t,eig,vec,vecT,expRot,expRot2);
         if (res != 0):
             raise RuntimeError("Basis size too small");
     else:
         i = 1;
         for t in time[:-1]: # use split step
            psi_t[i,:] = expRot2*psi_t[i-1,:];
            for k in range(scale):
                # I(t), gaussian beam
                if (k > 0): # Double half step:
                    psi_t[i,:] = expRot*psi_t[i,:];
                tp = t+k*dt;
                E_0_squared = E_0_squared_max * numpy.exp(-(tp**2)/(2*sigma**2)); 
                psi_t[i,:] = ((numpy.exp(-dt*E_0_squared*1j*eig))*vec).dot(vecT.dot(psi_t[i,:]));
                if (numpy.max(numpy.abs(psi_t[i,-2:])**2) > 1e-5):
                    raise RuntimeError("Basis size too small");
            psi_t[i,:] = expRot2*psi_t[i,:];
            i = i + 1;

     return psi_t;