def estimate_time_constant(fluor, p = 2, sn = None, lags = 5, fudge_factor = 1.):
"""
Estimate AR model parameters through the autocovariance function
Inputs
----------
fluor : nparray
One dimensional array containing the fluorescence intensities with
one entry per time-bin.
p : positive integer
order of AR system
sn : float
noise standard deviation, estimated if not provided.
lags : positive integer
number of additional lags where he autocovariance is computed
fudge_factor : float (0< fudge_factor <= 1)
shrinkage factor to reduce bias
Return
-----------
g : estimated coefficients of the AR process
"""
if sn is None:
sn = GetSn(fluor)
lags += p
xc = axcov(fluor,lags)
xc = xc[:,np.newaxis]
A = scipy.linalg.toeplitz(xc[lags+np.arange(lags)],xc[lags+np.arange(p)]) - sn**2*np.eye(lags,p)
g = np.linalg.lstsq(A,xc[lags+1:])[0]
gr = np.roots(np.concatenate([np.array([1]),-g.flatten()]))
gr = (gr+gr.conjugate())/2.
gr[gr>1] = 0.95 + np.random.normal(0,0.01,np.sum(gr>1))
gr[gr<0] = 0.15 + np.random.normal(0,0.01,np.sum(gr<0))
g = np.poly(fudge_factor*gr)
g = -g[1:]
return g.flatten()
python类roots()的实例源码
def _compute_minimum(self, a1, p1, dp1, a2, p2, dp2):
cubic_mtx = np.zeros((4, 4))
cubic_mtx[0, :] = [1., a1, a1 ** 2, a1 ** 3]
cubic_mtx[1, :] = [1., a2, a2 ** 2, a2 ** 3]
cubic_mtx[2, :] = [0., 1., 2 * a1, 3 * a1 ** 2]
cubic_mtx[3, :] = [0., 1., 2 * a2, 3 * a2 ** 2]
c0, c1, c2, c3 = np.linalg.solve(cubic_mtx, [p1, p2, dp1, dp2])
d0 = c1
d1 = 2 * c2
d2 = 3 * c3
r1, r2 = np.roots([d2, d1, d0])
a = None
p = max(p1, p2)
if (a1 <= r1 <= a2 or a2 <= r1 <= a1) and np.isreal(r1):
px = self._phi(r1)
if px < p:
a = r1
p = px
dp = self._dphi(r1)
if (a1 <= r2 <= a2 or a2 <= r2 <= a1) and np.isreal(r2):
px = self._phi(r2)
if px < p:
a = r2
p = px
dp = self._dphi(r2)
return a, p, dp
def roots_of_characteristic(self):
"""
This function calculates z_0 and the 2m roots of the characteristic equation
associated with the Euler equation (1.7)
Note:
------
numpy.poly1d(roots, True) defines a polynomial using its roots that can be
evaluated at any point. If x_1, x_2, ... , x_m are the roots then
p(x) = (x - x_1)(x - x_2)...(x - x_m)
"""
m = self.m
? = self.?
# Calculate the roots of the 2m-polynomial
roots = np.roots(?)
# sort the roots according to their length (in descending order)
roots_sorted = roots[np.argsort(abs(roots))[::-1]]
z_0 = ?.sum() / np.poly1d(roots, True)(1)
z_1_to_m = roots_sorted[:m] # we need only those outside the unit circle
? = 1 / z_1_to_m
return z_1_to_m, z_0, ?
def test_roots(self):
assert_array_equal(np.roots([1, 0, 0]), [0, 0])
def tune_everything(x0squared, C, T, gmin, gmax):
# First tune based on dynamic range
if C==0:
dr=gmax/gmin
mustar=((np.sqrt(dr)-1)/(np.sqrt(dr)+1))**2
alpha_star = (1+np.sqrt(mustar))**2/gmax
return alpha_star,mustar
dist_to_opt = x0squared
grad_var = C
max_curv = gmax
min_curv = gmin
const_fact = dist_to_opt * min_curv**2 / 2 / grad_var
coef = [-1, 3, -(3 + const_fact), 1]
roots = np.roots(coef)
roots = roots[np.real(roots) > 0]
roots = roots[np.real(roots) < 1]
root = roots[np.argmin(np.imag(roots) ) ]
assert root > 0 and root < 1 and np.absolute(root.imag) < 1e-6
dr = max_curv / min_curv
assert max_curv >= min_curv
mu = max( ( (np.sqrt(dr) - 1) / (np.sqrt(dr) + 1) )**2, root**2)
lr_min = (1 - np.sqrt(mu) )**2 / min_curv
lr_max = (1 + np.sqrt(mu) )**2 / max_curv
alpha_star = lr_min
mustar = mu
return alpha_star, mustar
def tfps(beta=0.6, ca=1., de2=0.000545, theta=0., kk=1.):
"""
beta: Total plasma beta,
ca: Alfven speed based on mean field,
de2: me/mi,
theta: Angle of propagation in degrees
kk: Wavenumber of interest in units of kdi
Output is frequencies of the roots and the phase speeds w/k
The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
"""
tht = theta*pi/180.
ct = np.cos(tht)
st = np.sin(tht)
tt = st/ct
cs=np.sqrt(beta/2.)*ca
di =1.
caksq=np.cos(tht)**2*ca**2
cmsq=ca**2 + cs**2
pcs=np.zeros(4)
D = 1 + kk**2*de2
# Find out the root of the quadratic dispersion relation
pcs[0] = 1.
pcs[1] = -(ca**2/D + cs**2 + caksq*(1+kk**2*di**2/D)/D)*kk**2
pcs[2] = caksq*kk**4*(ca**2/D + cs**2 + cs**2*(1+kk**2*di**2/D))/D
pcs[3] = -(cs**2*kk**6*caksq**2)/D**2
w2 = np.roots(pcs); w = np.sqrt(w2)
speeds= w/kk
return w,speeds
def tfps(beta=0.6, ca=1., de2=0.000545, theta=0., kk=1.):
"""
beta: Total plasma beta,
ca: Alfven speed based on mean field,
de2: me/mi,
theta: Angle of propagation as fraction of pi/2),
kk: Wavenumber of interest in units of kdi
Output is frequencies of the roots and the phase speeds w/k
The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
"""
# theta=float(raw_input("Angle of propagation as fraction of pi/2: "))
# kk=float(raw_input("What k value? "))
# ca = float(raw_input("Alfven Speed: "))
# beta=float(raw_input("Total plasma beta?: "))
# de2= float(raw_input("me/mi : "))
ct = cos(theta*pi/2)
st = sin(theta*pi/2)
tt = st/ct
cs=sqrt(beta/2.)*ca
di =1.
caksq=cos(theta*pi/2.)**2*ca**2
cmsq=ca**2 + cs**2
pcs=np.zeros(4)
D = 1 + kk**2*de2
# Find out the root of the quadratic dispersion relation
pcs[0] = 1.
pcs[1] = -(ca**2/D + cs**2 + caksq*(1+kk**2*di**2/D)/D)*kk**2
pcs[2] = caksq*kk**4*(ca**2/D + cs**2 + cs**2*(1+kk**2*di**2/D))/D
pcs[3] = -(cs**2*kk**6*caksq**2)/D**2
w2 = np.roots(pcs); w = sqrt(w2)
speeds= w/kk
return w,speeds
def get_mu(self):
coef = [-1.0, 3.0, 0.0, 1.0]
coef[2] = -(3 + self._dist_to_opt**2 * self._h_min**2 / 2 / self._grad_var)
roots = np.roots(coef)
root = roots[np.logical_and(np.logical_and(np.real(roots) > 0.0,
np.real(roots) < 1.0), np.imag(roots) < 1e-5) ]
assert root.size == 1
dr = self._h_max / self._h_min
self._mu_t = max(np.real(root)[0]**2, ( (np.sqrt(dr) - 1) / (np.sqrt(dr) + 1) )**2 )
return
def test_roots(self):
assert_array_equal(np.roots([1, 0, 0]), [0, 0])
def get_H(self, mass_upper):
"""Solve for H+, the concentration of hydrogen ions
(the (pH) of seawater).
:param mass_upper: Carbon mass in ocenas in GtC
:type mass_upper: float
:return: H
:rtype: float
"""
p0 = 1
p1 = (self.k_1 - mass_upper * self.k_1 / self.Alk)
p2 = (1 - 2 * mass_upper / self.Alk) * self.k_1 * self.k_2
return max(np.roots([p0, p1, p2]))
def get_H(self, mass_upper):
"""Solve for H+, the concentration of hydrogen ions
(the (pH) of seawater).
:param mass_upper: Carbon mass in ocenas in GtC
:type mass_upper: float
:return: H
:rtype: float
"""
p0 = 1
p1 = (self.k_1 - mass_upper * self.k_1 / self.Alk)
p2 = (1 - 2 * mass_upper / self.Alk) * self.k_1 * self.k_2
return max(np.roots([p0, p1, p2]))
def test_roots(self):
assert_array_equal(np.roots([1, 0, 0]), [0, 0])
def polyroots01(p):
"""Returns the real roots between 0 and 1 of the polynomial with
coefficients given in p,
p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
p can also be a np.poly1d object. See polyroots for more information."""
return polyroots(p, realroots=True, condition=lambda tval: 0 <= tval <= 1)
def polyroots01(p):
"""Returns the real roots between 0 and 1 of the polynomial with
coefficients given in p,
p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
p can also be a np.poly1d object. See polyroots for more information."""
return polyroots(p, realroots=True, condition=lambda tval: 0 <= tval <= 1)
def circleby2ptsradius(pt1Vec, pt2Vec, radius):
if radius < la.norm(pt1Vec-pt2Vec)/2:
centerVec1 = np.array(['NaN', 'NaN'])
centerVec2 = np.array(['NaN', 'NaN'])
else:
a = pt1Vec[0]**2-pt2Vec[0]**2
b = pt1Vec[1]**2-pt2Vec[1]**2
c = -2*(pt1Vec[0]-pt2Vec[0])
d = -2*(pt1Vec[1]-pt2Vec[1])
e = a+b
Coeff1 = 1+(d/c)**2
Coeff2 = ((2*d*e/c**2) +(pt2Vec[0]*2*d/c) -2*pt2Vec[1])
Coeff3 = ((e/c)**2+(2*pt2Vec[0]*e/c)+pt2Vec[0]**2+pt2Vec[1]**2-\
radius**2)
All_coeff = np.array([Coeff1, Coeff2, Coeff3])
Eq_root = np.roots(All_coeff)
# centersArray--center of the circle. It is a 2x2 matrix. First row
# represents first possible center (x1, y1) and second row is the
# second possible center.
centersArray = np.zeros((len(Eq_root),2))
for i in list(range(len(Eq_root))):
x = -(e+d*Eq_root[i])/c
centersArray[i,0] = x
centersArray[i,1] = Eq_root[i]
# Check if values in centerArray are real or unreal numbers
if la.norm(centersArray.imag) != 0:
print('Circle with the specified radius does not pass through',
'specified points', sep=" ")
centerVec1 = np.array(['NaN', 'NaN'])
centerVec2 = np.array(['NaN', 'NaN'])
else:
centerVec1 = centersArray[0, :]
centerVec2 = centersArray[1, :]
return centerVec1, centerVec2
def tuneEverything(self, x0squared, c, t, gmin, gmax):
# First tune based on dynamic range
if c == 0:
dr = gmax / gmin
mustar = ((np.sqrt(dr) - 1) / (np.sqrt(dr) + 1))**2
alpha_star = (1 + np.sqrt(mustar))**2/gmax
return alpha_star, mustar
dist_to_opt = x0squared
grad_var = c
max_curv = gmax
min_curv = gmin
const_fact = dist_to_opt * min_curv**2 / 2 / grad_var
coef = [-1, 3, -(3 + const_fact), 1]
roots = np.roots(coef)
roots = roots[np.real(roots) > 0]
roots = roots[np.real(roots) < 1]
root = roots[np.argmin(np.imag(roots))]
assert root > 0 and root < 1 and np.absolute(root.imag) < 1e-6
dr = max_curv / min_curv
assert max_curv >= min_curv
mu = max(((np.sqrt(dr) - 1) / (np.sqrt(dr) + 1))**2, root**2)
lr_min = (1 - np.sqrt(mu))**2 / min_curv
alpha_star = lr_min
mustar = mu
return alpha_star, mustar
def check_mcs_step_septum(septum, proton, step_in):
"""
From K2 scattering routine. It is basically treatment of edge effect inside the septum jaws
@param septum: septum object
@param proton: Proton object
@param step_in: Random step chosen for this iteration
"""
theta_0 = 13.6e-3 / (proton.mom_p) / np.sqrt(septum.rad_length) * np.sqrt(step_in) * (1. + 0.038 * np.log(step_in / septum.rad_length))
x_temp = proton.x - septum.aper_u if proton.x <= (septum.aper_u + septum.thickness / 2.) else proton.x - septum.aper_u - septum.thickness
a = 9 * theta_0 ** 2 / 3.
b = -1 * proton.x_p ** 2
c = -2 * proton.x_p * x_temp
d = -x_temp ** 2
# print 'a', a, 'b', b, 'c', c, 'd', d
# test = np.linspace(-10, 10, 100)
# test_y = a * test ** 3 + b * test ** 2 + test * c + d
# plt.plot(test, test_y, 'lime')
sol = np.roots([a, b, c, d])
# print sol
sol_re = [100000.]
for i in sol:
if np.imag(i) == 0:
sol_re.append(float(np.real(i)))
min_sol = min(sol_re)
if min_sol < step_in:
if min_sol < septum.rad_length * 0.1e-2:
return septum.rad_length * 0.1e-2
else:
return min(sol_re)
else:
return step_in
def lpc_to_lsf(all_lpc):
if len(all_lpc.shape) < 2:
all_lpc = all_lpc[None]
order = all_lpc.shape[1] - 1
all_lsf = np.zeros((len(all_lpc), order))
for i in range(len(all_lpc)):
lpc = all_lpc[i]
lpc1 = np.append(lpc, 0)
lpc2 = lpc1[::-1]
sum_filt = lpc1 + lpc2
diff_filt = lpc1 - lpc2
if order % 2 != 0:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
deconv_sum = sum_filt
else:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])
roots_diff = np.roots(deconv_diff)
roots_sum = np.roots(deconv_sum)
angle_diff = np.angle(roots_diff[::2])
angle_sum = np.angle(roots_sum[::2])
lsf = np.sort(np.hstack((angle_diff, angle_sum)))
if len(lsf) != 0:
all_lsf[i] = lsf
return np.squeeze(all_lsf)
def lpc_to_lsf(all_lpc):
if len(all_lpc.shape) < 2:
all_lpc = all_lpc[None]
order = all_lpc.shape[1] - 1
all_lsf = np.zeros((len(all_lpc), order))
for i in range(len(all_lpc)):
lpc = all_lpc[i]
lpc1 = np.append(lpc, 0)
lpc2 = lpc1[::-1]
sum_filt = lpc1 + lpc2
diff_filt = lpc1 - lpc2
if order % 2 != 0:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
deconv_sum = sum_filt
else:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])
roots_diff = np.roots(deconv_diff)
roots_sum = np.roots(deconv_sum)
angle_diff = np.angle(roots_diff[::2])
angle_sum = np.angle(roots_sum[::2])
lsf = np.sort(np.hstack((angle_diff, angle_sum)))
if len(lsf) != 0:
all_lsf[i] = lsf
return np.squeeze(all_lsf)
def lpc_to_lsf(all_lpc):
if len(all_lpc.shape) < 2:
all_lpc = all_lpc[None]
order = all_lpc.shape[1] - 1
all_lsf = np.zeros((len(all_lpc), order))
for i in range(len(all_lpc)):
lpc = all_lpc[i]
lpc1 = np.append(lpc, 0)
lpc2 = lpc1[::-1]
sum_filt = lpc1 + lpc2
diff_filt = lpc1 - lpc2
if order % 2 != 0:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, 0, -1])
deconv_sum = sum_filt
else:
deconv_diff, _ = sg.deconvolve(diff_filt, [1, -1])
deconv_sum, _ = sg.deconvolve(sum_filt, [1, 1])
roots_diff = np.roots(deconv_diff)
roots_sum = np.roots(deconv_sum)
angle_diff = np.angle(roots_diff[::2])
angle_sum = np.angle(roots_sum[::2])
lsf = np.sort(np.hstack((angle_diff, angle_sum)))
if len(lsf) != 0:
all_lsf[i] = lsf
return np.squeeze(all_lsf)
