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
python类integrate()的实例源码
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
