def integrate(self, attrname, radio, x_box, y_box):
element = radio.labels[radio.active]
source_local = getattr(self, attrname+"_"+element+"_source")
lower_xlim = float(x_box.value)
lower_ylim = float(y_box.value)
x = np.array(source_local.data["x"])
y = np.array(source_local.data["y"])
x_change = x[x>lower_xlim]*1e-4
y_change = y[len(y)-len(x_change):]
integral = np.trapz(y_change, x = x_change)
comparsion = np.sum(y) * x[-1]*1e-4
print(integral, comparsion)
python类trapz()的实例源码
def band_area(spectrum, low_endmember=None, high_endmember=None):
"""
Compute the area under the curve between two endpoints where the
y-value <= 1.
"""
x = spectrum.index
y = spectrum
if not low_endmember:
low_endmember = x[0]
if not high_endmember:
high_endmember = x[-1]
ny = y[low_endmember:high_endmember]
return np.trapz(-ny[ny <= 1.0])
def add_filtered_planck(self, wavelength, trans):
vals = np.zeros(self.x.shape, 'float64')
nu = clight/(wavelength*1e-8)
nu = nu[::-1]
for i,logT in enumerate(self.x):
T = 10**logT
# Black body at this nu, T
Bnu = ((2.0 * hcgs * nu**3) / clight**2.0) / \
(np.exp(hcgs * nu / (kboltz * T)) - 1.0)
# transmission
f = Bnu * trans[::-1]
# integrate transmission over nu
vals[i] = np.trapz(f,nu)
# normalize by total transmission over filter
self.y = vals/trans.sum() #/np.trapz(trans[::-1],nu)
#self.y = np.clip(np.maximum(vals, self.y), 0.0, 1.0)
def convolve(self, wavelengths, densities):
# define short names for the involved wavelength grids
wa = wavelengths
wb = self._Wavelengths
# create a combined wavelength grid, restricted to the overlapping interval
w1 = wa[ (wa>=wb[0]) & (wa<=wb[-1]) ]
w2 = wb[ (wb>=wa[0]) & (wb<=wa[-1]) ]
w = np.unique(np.hstack((w1,w2)))
if len(w) < 2: return 0
# log-log interpolate SED and transmission on the combined wavelength grid
# (use scipy interpolation function for SED because np.interp does not support broadcasting)
F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w)))
T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.))
# perform the integration
if self._PhotonCounter:
return np.trapz(x=w, y=w*F*T) / self._IntegratedTransmission
else:
return np.trapz(x=w, y=F*T) / self._IntegratedTransmission
## This function calculates and returns the integrated value for a given spectral energy distribution over the
# filter's wavelength range,
def integrate(self, wavelengths, densities):
# define short names for the involved wavelength grids
wa = wavelengths
wb = self._Wavelengths
# create a combined wavelength grid, restricted to the overlapping interval
w1 = wa[(wa >= wb[0]) & (wa <= wb[-1])]
w2 = wb[(wb >= wa[0]) & (wb <= wa[-1])]
w = np.unique(np.hstack((w1, w2)))
if len(w) < 2: return 0
# log-log interpolate SED and transmission on the combined wavelength grid
# (use scipy interpolation function for SED because np.interp does not support broadcasting)
F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w)))
T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.))
# perform the integration
if self._PhotonCounter: return np.trapz(x=w, y=w * F * T)
else: return np.trapz(x=w, y=F * T)
## This private helper function returns the natural logarithm for positive values, and a large negative number
# (but not infinity) for zero or negative values.
def convolve(self, wavelengths, densities):
# define short names for the involved wavelength grids
wa = wavelengths
wb = self._Wavelengths
# create a combined wavelength grid, restricted to the overlapping interval
w1 = wa[ (wa>=wb[0]) & (wa<=wb[-1]) ]
w2 = wb[ (wb>=wa[0]) & (wb<=wa[-1]) ]
w = np.unique(np.hstack((w1,w2)))
if len(w) < 2: return 0
# log-log interpolate SED and transmission on the combined wavelength grid
# (use scipy interpolation function for SED because np.interp does not support broadcasting)
F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w)))
T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.))
# perform the integration
if self._PhotonCounter:
return np.trapz(x=w, y=w*F*T) / self._IntegratedTransmission
else:
return np.trapz(x=w, y=F*T) / self._IntegratedTransmission
## This function calculates and returns the integrated value for a given spectral energy distribution over the
# filter's wavelength range,
def integrate(self, wavelengths, densities):
# define short names for the involved wavelength grids
wa = wavelengths
wb = self._Wavelengths
# create a combined wavelength grid, restricted to the overlapping interval
w1 = wa[(wa >= wb[0]) & (wa <= wb[-1])]
w2 = wb[(wb >= wa[0]) & (wb <= wa[-1])]
w = np.unique(np.hstack((w1, w2)))
if len(w) < 2: return 0
# log-log interpolate SED and transmission on the combined wavelength grid
# (use scipy interpolation function for SED because np.interp does not support broadcasting)
F = np.exp(interp1d(np.log(wa), _log(densities), copy=False, bounds_error=False, fill_value=0.)(np.log(w)))
T = np.exp(np.interp(np.log(w), np.log(wb), _log(self._Transmission), left=0., right=0.))
# perform the integration
if self._PhotonCounter: return np.trapz(x=w, y=w * F * T)
else: return np.trapz(x=w, y=F * T)
## This private helper function returns the natural logarithm for positive values, and a large negative number
# (but not infinity) for zero or negative values.
def obj_func_disc_nofit(xx, e_g, e_p, g, rho, Q_g, Q_p, head_g, head_p, optimizing = True):
H_T = int(price_duration.Frequency.max()) # total duration (100%)
prc_g, prc_p, freq_g, freq_p = [],[],[],[]
for i,x in enumerate(price_duration.Frequency):
if x < xx: # Power Generation price and duration
prc_g.append(price_duration.Price[i]), freq_g.append(x)
if H_T - xx < x < H_T: # Pumping price and duration
prc_p.append(price_duration.Price[i]), freq_p.append(x)
prc_g = np.array(prc_g) # generation price
prc_p = np.array(prc_p) # pumping price
freq_g = np.array(freq_g) # generation duration
freq_p = np.array(freq_p) # pumping duration
# Use numerical integration to integrate (Trapezoidal rule)
Power_Revenue = np.trapz(prc_g, freq_g, dx=0.1, axis = -1)*e_g*rho*g*Q_g*head_g/(10**6)
Pumping_Cost = np.trapz(prc_p, freq_p, dx=0.1, axis = -1)/e_p*rho*g*Q_p*head_p/(10**6)
z = Power_Revenue - Pumping_Cost # profit
return z if optimizing else -z
# fit a curve
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 nena(locs, info, callback=None):
bin_centers, dnfl_ = next_frame_neighbor_distance_histogram(locs, callback)
def func(d, a, s, ac, dc, sc):
f = a * (d / s**2) * _np.exp(-0.5 * d**2 / s**2)
fc = ac * (d / sc**2) * _np.exp(-0.5 * (d**2 + dc**2) / sc**2) * _iv(0, d * dc / sc)
return f + fc
pdf_model = _lmfit.Model(func)
params = _lmfit.Parameters()
area = _np.trapz(dnfl_, bin_centers)
median_lp = _np.mean([_np.median(locs.lpx), _np.median(locs.lpy)])
params.add('a', value=area/2, min=0)
params.add('s', value=median_lp, min=0)
params.add('ac', value=area/2, min=0)
params.add('dc', value=2*median_lp, min=0)
params.add('sc', value=median_lp, min=0)
result = pdf_model.fit(dnfl_, params, d=bin_centers)
return result, result.best_values['s']
def test_gamma():
tsample = 0.005 / 1000
with pytest.raises(ValueError):
t, g = utils.gamma(0, 0.1, tsample)
with pytest.raises(ValueError):
t, g = utils.gamma(2, -0.1, tsample)
with pytest.raises(ValueError):
t, g = utils.gamma(2, 0.1, -tsample)
for tau in [0.001, 0.01, 0.1]:
for n in [1, 2, 5]:
t, g = utils.gamma(n, tau, tsample)
npt.assert_equal(np.arange(0, t[-1] + tsample / 2.0, tsample), t)
if n > 1:
npt.assert_equal(g[0], 0.0)
# Make sure area under the curve is normalized
npt.assert_almost_equal(np.trapz(np.abs(g), dx=tsample), 1.0,
decimal=2)
# Make sure peak sits correctly
npt.assert_almost_equal(g.argmax() * tsample, tau * (n - 1))
def mag2diam(ref_mag, ref_band, teff, nbr_pts=10):
"""
Returns the DIAMETER (not radius) of a black body of a given
magnitude (U,B,V,R,I,J,H,K)
"""
filt = ALL_FILTERS[ref_band]
if nbr_pts != 10:
lam = _gen_lam_for_filter(filt, nbr_pts)
else:
lam = filt['span_10pts']
# flux in watts from magnitude
res = filt['flux_wm2m']*10**(-0.39809*ref_mag)*filt['delta_wl']
# integrated flux
res /= np.trapz(blackbody_spectral_irr(teff, lam), lam)
return 2*RAD2MAS*np.sqrt(res)
def checkFacts_pdf_Mu(nu=0, tau=0, mu_grid=None, **kwargs):
cPrior = c_Prior(nu=nu, tau=tau)
if mu_grid is None:
mu_grid = make_mu_grid(**kwargs)
phi_grid = mu2phi(mu_grid)
logpdf_grid = tau * phi_grid - nu * c_Phi(phi_grid) - cPrior
logJacobian_grid = logJacobian_mu(mu_grid)
pdf_grid = np.exp(logpdf_grid + logJacobian_grid)
Integral = np.trapz(pdf_grid, mu_grid)
E_mu_numeric = np.trapz(pdf_grid * mu_grid, mu_grid)
E_mu_formula = tau/nu
E_phi_numeric = np.trapz(pdf_grid * phi_grid, mu_grid)
E_phi_formula = E_phi(nu,tau)
print "nu=%7.3f tau=%7.3f" % (nu, tau)
print " Integral=% 7.3f should be % 7.3f" % (Integral, 1.0)
print " E[mu]=% 7.3f should be % 7.3f" % (E_mu_numeric, E_mu_formula)
print " E[phi]=% 7.3f should be % 7.3f" % (E_phi_numeric, E_phi_formula)
def checkFacts_pdf_Phi(
nu=0, tau=0, phi_grid=None, **kwargs):
pdf_grid, phi_grid = make_pdfgrid_phi(nu, tau, phi_grid, **kwargs)
mu_grid = phi2mu(phi_grid, **kwargs)
IntegralVal = np.trapz(pdf_grid, phi_grid)
E_phi_numeric = np.trapz(pdf_grid * phi_grid, phi_grid)
E_phi_formula = E_phi(nu=nu, tau=tau, **kwargs)
E_mu_numeric = np.trapz(pdf_grid * mu_grid, phi_grid)
E_mu_formula = tau/nu
mode_phi_numeric = phi_grid[np.argmax(pdf_grid)]
mode_phi_formula = mu2phi(tau/nu, **kwargs)
print "nu=%7.3f tau=%7.3f" % (nu, tau)
print " Integral=% 7.3f should be % 7.3f" % (IntegralVal, 1.0)
print " E[mu]=% 7.3f should be % 7.3f" % (E_mu_numeric, E_mu_formula)
print " E[phi]=% 7.3f should be % 7.3f" % (E_phi_numeric, E_phi_formula)
print " mode[phi]=% 7.3f should be % 7.3f" % (
mode_phi_numeric, mode_phi_formula)
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
# Record max overlap value for each gt box
# Return vector of overlap values
gt_overlaps = np.zeros(0)
for i in xrange(self.num_images):
gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
gt_boxes = self.roidb[i]['boxes'][gt_inds, :]
boxes = candidate_boxes[i]
if boxes.shape[0] == 0:
continue
overlaps = bbox_overlaps(boxes.astype(np.float),
gt_boxes.astype(np.float))
# gt_overlaps = np.hstack((gt_overlaps, overlaps.max(axis=0)))
_gt_overlaps = np.zeros((gt_boxes.shape[0]))
for j in xrange(gt_boxes.shape[0]):
argmax_overlaps = overlaps.argmax(axis=0)
max_overlaps = overlaps.max(axis=0)
gt_ind = max_overlaps.argmax()
gt_ovr = max_overlaps.max()
assert(gt_ovr >= 0)
box_ind = argmax_overlaps[gt_ind]
_gt_overlaps[j] = overlaps[box_ind, gt_ind]
assert(_gt_overlaps[j] == gt_ovr)
overlaps[box_ind, :] = -1
overlaps[:, gt_ind] = -1
gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))
num_pos = gt_overlaps.size
gt_overlaps = np.sort(gt_overlaps)
step = 0.001
thresholds = np.minimum(np.arange(0.5, 1.0 + step, step), 1.0)
recalls = np.zeros_like(thresholds)
for i, t in enumerate(thresholds):
recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
ar = 2 * np.trapz(recalls, thresholds)
return ar, gt_overlaps, recalls, thresholds
def make_curve(x, y, xlabel, ylabel, filename):
"""Draw and save ROC or PR curves."""
auc = np.trapz(x, y, dx=0.001) # area under the curve
sns.plt.figure()
sns.plt.clf()
sns.plt.plot(x, y)
sns.plt.xlim([0, 1])
sns.plt.ylim([0, 1])
sns.plt.xlabel(xlabel)
sns.plt.ylabel(ylabel)
sns.plt.title("AUC: {}".format(auc))
sns.plt.savefig(filename)
return auc
def plotROC(self):
for classifier in self.structures:
x, y = [], []
for threshold in np.arange(self.minScore, self.maxScore, 0.1):
tp, tn, fp, fn = 0, 0, 0, 0
for realS in self.structures:
for score in self.scores[(classifier, realS)]:
if score >= threshold: # Positive
if classifier == realS: # True
tp += 1
else: # False
fp += 1
else: # Negative
if classifier != realS: # True
tn += 1
else: # False
fn += 1
tpr = tp / (tp + fn)
fpr = fp / (tn + fp)
x.append(fpr)
y.append(tpr)
print("----- Receiver Operating Characteristic Curve -----")
print("Classifier:", classifier)
x.sort()
y.sort()
auc = np.trapz(y, x) # Area Under Curve
print("AUC:", auc)
pyplot.title("Classifier " + classifier)
pyplot.xlabel('False Positive Rate')
pyplot.ylabel('True Positive Rate')
pyplot.plot([0, 1], [0, 1])
pyplot.plot(x, y)
pyplot.show()
def test_limbdark(tol = 1e-4):
'''
Test the limb darkening normalization by comparing the total flux of
the star to what you get with the Stefan-Boltzmann law.
'''
# Instantiate the star
r = 0.1
teff = 3200
star = Star('A', m = 0.1, r = r, nz = 31, color = 'k',
limbdark = [u1], teff = teff)
# True luminosity
truth = 4 * np.pi * (r * RSUN) ** 2 * SBOLTZ * teff ** 4
# System
system = System(star, quiet = True)
system.distance = 12.2
# Compute the light curve
time = np.arange(0., 1., 10)
system.compute(time, lambda1 = 0.01, lambda2 = 1000, R = 3000)
# Luminosity
bol = np.trapz(system.flux[0], system.A.wavelength * 1e6)
lum = 4 * np.pi * bol * (12.2 * PARSEC) ** 2
# Check!
assert np.abs(lum - truth) / truth < tol, \
"Incorrect bolometric flux: %.10e != %.10e." % (lum, truth)
def autocorrelation_function(self, r):
"""compute the real space autocorrelation function for the Gaussian random field model
"""
# compute the cut-level parameter beta
beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1)
# the covariance of the GRF
acf_psi = (np.exp(-r/self.corr_length) * (1 + r / self.corr_length)
* np.sinc(2*r/self.repeat_distance))
# integral discretization. henning says: the resolution 1e-2 is ad hoc, test required,
# the integrand has a (integrable) singularity for t=1 and acf_psi = 1, so an adaptive
# discretization seems preferable -> TODO
dt = 1e-2
t = np.arange(0, 1, dt)
# the gridded integrand, via change of integration variable
# compared to the wp-2 docu, to enable array-based computation
t_gridded, acf_psi_gridded = np.meshgrid(t, acf_psi)
integrand_gridded = (acf_psi_gridded / np.sqrt(1 - (t_gridded * acf_psi_gridded)**2)
* np.exp( - beta**2 / (1 + t_gridded * acf_psi_gridded)))
acf = 1.0 / (2 * np.pi) * np.trapz(integrand_gridded, x=t_gridded)
return acf
# ft not known analytically: deligate
# def ft_autocorrelation_function(self, k):
def trapz2(f, x=None, y=None, dx=1.0, dy=1.0):
"""Double integrate."""
return numpy.trapz(numpy.trapz(f, x=y, dx=dy), x=x, dx=dx)
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
# Record max overlap value for each gt box
# Return vector of overlap values
gt_overlaps = np.zeros(0)
for i in xrange(self.num_images):
gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
gt_boxes = self.roidb[i]['boxes'][gt_inds, :]
boxes = candidate_boxes[i]
if boxes.shape[0] == 0:
continue
overlaps = bbox_overlaps(boxes.astype(np.float),
gt_boxes.astype(np.float))
# gt_overlaps = np.hstack((gt_overlaps, overlaps.max(axis=0)))
_gt_overlaps = np.zeros((gt_boxes.shape[0]))
for j in xrange(gt_boxes.shape[0]):
argmax_overlaps = overlaps.argmax(axis=0)
max_overlaps = overlaps.max(axis=0)
gt_ind = max_overlaps.argmax()
gt_ovr = max_overlaps.max()
assert(gt_ovr >= 0)
box_ind = argmax_overlaps[gt_ind]
_gt_overlaps[j] = overlaps[box_ind, gt_ind]
assert(_gt_overlaps[j] == gt_ovr)
overlaps[box_ind, :] = -1
overlaps[:, gt_ind] = -1
gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))
num_pos = gt_overlaps.size
gt_overlaps = np.sort(gt_overlaps)
step = 0.001
thresholds = np.minimum(np.arange(0.5, 1.0 + step, step), 1.0)
recalls = np.zeros_like(thresholds)
for i, t in enumerate(thresholds):
recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
ar = 2 * np.trapz(recalls, thresholds)
return ar, gt_overlaps, recalls, thresholds
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
gt_overlaps = np.zeros(0)
for i in xrange(self.num_images):
gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
gt_boxes = self.roidb[i]['boxes'][gt_inds, :]
boxes = candidate_boxes[i]
if boxes.shape[0] == 0:
continue
overlaps = bbox_overlaps(boxes.astype(np.float),
gt_boxes.astype(np.float))
_gt_overlaps = np.zeros((gt_boxes.shape[0]))
for j in xrange(gt_boxes.shape[0]):
argmax_overlaps = overlaps.argmax(axis=0)
max_overlaps = overlaps.max(axis=0)
gt_ind = max_overlaps.argmax()
gt_ovr = max_overlaps.max()
assert(gt_ovr >= 0)
box_ind = argmax_overlaps[gt_ind]
_gt_overlaps[j] = overlaps[box_ind, gt_ind]
assert(_gt_overlaps[j] == gt_ovr)
overlaps[box_ind, :] = -1
overlaps[:, gt_ind] = -1
gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))
num_pos = gt_overlaps.size
gt_overlaps = np.sort(gt_overlaps)
step = 0.05
thresholds = np.minimum(np.arange(0.2, 1.0 + step, step), 1.0)
recalls = np.zeros_like(thresholds)
for i, t in enumerate(thresholds):
recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
ar = 2 * np.trapz(recalls, thresholds)
return ar, gt_overlaps, recalls, thresholds
def compute_ess_transmission(beam_sigma, slits, dispersion):
"""Compute the transmission as a function of the momentum offset (in %) from a simple analytical model."""
n_steps = 10000
dx = 3.0/n_steps
sigma = beam_sigma/2.8
slits_at = slits/dispersion
error = np.arange(-1.5, 1.5, dx)
slits = np.zeros(n_steps)
slits[np.where((error < slits_at) & (error > -slits_at))] = 1.0
beam = np.exp(-(error/sigma)**2/2)
return np.roll(np.convolve(slits, beam, mode="same"), -1)/np.trapz(beam)
def bpm_fit(a, p):
if a is None:
return None
x = a['channel']['data'][p][:, 1]
y = a['channel']['data'][p][:, 2]
ar = np.trapz(y / np.sum(y) * len(y), x)
mean = np.mean(x * y / np.sum(y) * len(y))
rms = np.std(x * y / np.sum(y) * len(y))
popt, pcov = curve_fit(gaussian, x, y, p0=[ar, mean, rms])
return [popt, np.sqrt(pcov.diagonal())]
imdb2.py 文件源码
项目：Automatic_Group_Photography_Enhancement
作者: Yuliang-Zou
项目源码
文件源码
阅读 30
收藏 0
点赞 0
评论 0
def evaluate_recall(self, candidate_boxes, ar_thresh=0.5):
# Record max overlap value for each gt box
# Return vector of overlap values
gt_overlaps = np.zeros(0)
for i in xrange(self.num_images):
gt_inds = np.where(self.roidb[i]['gt_classes'] > 0)[0]
gt_boxes = self.roidb[i]['boxes'][gt_inds, :]
boxes = candidate_boxes[i]
if boxes.shape[0] == 0:
continue
overlaps = bbox_overlaps(boxes.astype(np.float),
gt_boxes.astype(np.float))
# gt_overlaps = np.hstack((gt_overlaps, overlaps.max(axis=0)))
_gt_overlaps = np.zeros((gt_boxes.shape[0]))
for j in xrange(gt_boxes.shape[0]):
argmax_overlaps = overlaps.argmax(axis=0)
max_overlaps = overlaps.max(axis=0)
gt_ind = max_overlaps.argmax()
gt_ovr = max_overlaps.max()
assert(gt_ovr >= 0)
box_ind = argmax_overlaps[gt_ind]
_gt_overlaps[j] = overlaps[box_ind, gt_ind]
assert(_gt_overlaps[j] == gt_ovr)
overlaps[box_ind, :] = -1
overlaps[:, gt_ind] = -1
gt_overlaps = np.hstack((gt_overlaps, _gt_overlaps))
num_pos = gt_overlaps.size
gt_overlaps = np.sort(gt_overlaps)
step = 0.001
thresholds = np.minimum(np.arange(0.5, 1.0 + step, step), 1.0)
recalls = np.zeros_like(thresholds)
for i, t in enumerate(thresholds):
recalls[i] = (gt_overlaps >= t).sum() / float(num_pos)
ar = 2 * np.trapz(recalls, thresholds)
return ar, gt_overlaps, recalls, thresholds
def avg_mole_per_sec(raw, start_point=0, end_point='eq', constV=False):
#print type(end_point)
T = raw['temperature']
x = raw['axis0']
rr = raw['net_reaction_rate']
V = raw['volume']
n_rxn = rr.shape[1] - 1
i_start = start_point
if type(end_point) is int:
i_end = end_point
if i_end == i_start:
mole_per_sec = rr[i_start, :] * float(V[i_start])
return np.transpose(mole_per_sec)
else:
if end_point is 'ign':
T_end = T[0] + 400
elif end_point is 'eq':
T_end = T[-1] - 50
i_end = np.argmin(abs(T - T_end))
avg_rr = np.ndarray([n_rxn+1,1])
for id_rxn in range(1, n_rxn + 1):
if constV:
mole_per_sec = rr[i_start:i_end,id_rxn] * V[0]
else:
mole_per_sec = np.multiply(rr[i_start:i_end, id_rxn], V[i_start:i_end])
int_rr = np.trapz(np.transpose(mole_per_sec), np.transpose(x[i_start:i_end]))
avg_rr[id_rxn] = float(int_rr) / (x[i_end] - x[i_start])
return avg_rr
def WRelease(self,kr,r0,kt,t0,kb,b0,T):
"""
Do the analytical "release" of guest to standard concentration
"""
### SETUP
# Example: print WRelease(5.0,24.00,100.0,180.0,100.0,180.0,298.15)
# kr, r0: distance restraint (r) force constant, target value
# kt, t0: angle restraint (theta) force constant, target value
# kb, b0: angle restraint (b) force constant, target value
# T: Temperature
R = 1.987204118e-3 # kcal/mol-K, a.k.a. boltzman constant
beta = 1/(T*R)
rlb,rub,rst = [0.0,100.0,0.0001] # r lower bound, upper bound, step size
tlb,tub,tst = [0.0,np.pi,0.00005] # theta ", ", "
blb,bub,bst = [0.0,np.pi,0.00005] # b ", ", "
def fr(val):
return (val**2)*np.exp(-beta*kr*(val-r0)**2)
def ft(val):
return np.sin(val)*np.exp(-beta*kt*(val-np.radians(t0))**2)
def fb(val):
return np.sin(val)*np.exp(-beta*kb*(val-np.radians(b0))**2)
### Integrate
rint,tint,bint = [0.0,0.0,0.0]
intrange = np.arange(rlb,rub,rst)
rint = np.trapz(fr(intrange),intrange)
intrange = np.arange(tlb,tub,tst)
tint = np.trapz(ft(intrange),intrange)
intrange = np.arange(blb,bub,bst)
bint = np.trapz(fb(intrange),intrange)
return R*T*np.log(np.pi*(1.0/1660.0)*rint*tint*bint)
def quad(fintegrand, xmin, xmax, n=1e4):
spacings = np.logspace(np.log10(xmin), np.log10(xmax), n)
integrand_arr = fintegrand(spacings)
val = np.trapz(integrand_arr, dx=np.diff(spacings))
return val
def a2t(at, Om0=0.27, Oml0=0.73, h=0.700):
integrand = lambda x: 1./(x*sqrt(Oml0+Om0*x**-3.0))
# current_time,err = si.quad(integrand,0.0,at,epsabs=1e-6,epsrel=1e-6)
current_time = quad(integrand, 1e-4, at)
# spacings = np.logspace(-5,np.log10(at),1e5)
# integrand_arr = integrand(spacings)
# current_time = np.trapz(integrand_arr,dx=np.diff(spacings))
current_time *= 9.779/h
return current_time
def integrate_inf(fcn, error=1e-3, initial_guess=10):
"""
Integrate a function *fcn* from zero to infinity, stopping when the answer
changes by less than *error*. Hopefully someday we can do something
better than this!
"""
xvals = np.logspace(0,np.log10(initial_guess), initial_guess+1)-.9
yvals = fcn(xvals)
xdiffs = xvals[1:] - xvals[:-1]
# Trapezoid rule, but with different dxes between values, so np.trapz
# will not work.
areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
area0 = np.sum(areas)
# Next guess.
next_guess = 10 * initial_guess
xvals = np.logspace(0,np.log10(next_guess), 2*initial_guess**2+1)-.99
yvals = fcn(xvals)
xdiffs = xvals[1:] - xvals[:-1]
# Trapezoid rule.
areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
area1 = np.sum(areas)
# Now we refine until the error is smaller than *error*.
diff = area1 - area0
area_last = area1
one_pow = 3
while diff > error:
next_guess *= 10
xvals = np.logspace(0,np.log10(next_guess), one_pow*initial_guess**one_pow+1) - (1 - 0.1**one_pow)
yvals = fcn(xvals)
xdiffs = xvals[1:] - xvals[:-1]
# Trapezoid rule.
areas = (yvals[1:] + yvals[:-1]) * xdiffs / 2.0
area_next = np.sum(areas)
diff = area_next - area_last
area_last = area_next
one_pow+=1
return area_last
