Python numpy 模块，roots() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.roots()。
def get_Z(self, T, p):
#A = self.a * p / self.R**2 / T**2.5
#B = self.b * p / self.R / T
A = self.get_A(T, p)
B = self.get_B(T, p)
# Solve the cubic equation for compressibility factor z
# Z^3 - Z^2 + (A-B-B**2)*Z - A*B = 0
coeffs = [1, -1, A-B-B**2, -A*B]
roots = np.roots(coeffs)
real_roots = roots[np.isreal(roots)].real
valid_roots = real_roots[real_roots > p*self.b/self.R/T]
return valid_roots
# dZ/dT at const. p
def phormants(x, Fs):
N = len(x)
w = numpy.hamming(N)
# Apply window and high pass filter.
x1 = x * w
x1 = lfilter([1], [1., 0.63], x1)
# Get LPC.
ncoeff = 2 + Fs / 1000
A, e, k = lpc(x1, ncoeff)
#A, e, k = lpc(x1, 8)
# Get roots.
rts = numpy.roots(A)
rts = [r for r in rts if numpy.imag(r) >= 0]
# Get angles.
angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))
# Get frequencies.
frqs = sorted(angz * (Fs / (2 * math.pi)))
return frqs
def lpf(x, cutoff, fs, order=5):
"""
low pass filters signal with Butterworth digital
filter according to cutoff frequency
filter uses Gustafsson’s method to make sure
forward-backward filt == backward-forward filt
Note that edge effects are expected
Args:
x (array): signal data (numpy array)
cutoff (float): cutoff frequency (Hz)
fs (int): sample rate (Hz)
order (int): order of filter (default 5)
Returns:
filtered (array): low pass filtered data
"""
nyquist = fs / 2
b, a = butter(order, cutoff / nyquist)
if not np.all(np.abs(np.roots(a)) < 1):
raise PsolaError('Filter with cutoff at {} Hz is unstable given '
'sample frequency {} Hz'.format(cutoff, fs))
filtered = filtfilt(b, a, x, method='gust')
return filtered
def estimate_time_constant(y, p=2, sn=None, lags=5, fudge_factor=1.):
"""
Estimate AR model parameters through the autocovariance function
Parameters
----------
y : array, shape (T,)
One dimensional array containing the fluorescence intensities with
one entry per time-bin.
p : positive integer
order of AR system
sn : float
sn 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
Returns
-------
g : estimated coefficients of the AR process
"""
if sn is None:
sn = GetSn(y)
lags += p
xc = axcov(y, 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()
def find_closest(t, v, t0, v0):
""" Find the closest point on the curve f = a + b/x
to the given point (t,v)
"""
a = v0
b = v0*t0
# Solve for intersection points
eqn_coefs = [1/b, -t/b, 0, v-a, -b]
tis = np.roots(eqn_coefs)
tis = tis[abs(tis.imag/tis.real)<0.01].real # We care only real solutions
tis = tis[tis>0] # and positive ones
# Choose the shortest among solutions
ds = abs(tis-t)*np.sqrt(1 + np.power(tis,4)/(b*b)) # Distance from solutions to given point (t,v)
idx = np.argmin(ds)
ti = tis[idx]
vi = a + b/ti
return ti, vi
def get_Z(self, T, p):
kappa = 0.37464 + 1.54226*self.acentric - 0.26992*self.acentric**2
Tr = T/self.T_crit
alpha = (1 + kappa*(1 - Tr**0.5))**2
A = alpha * self.a * p / self.R**2 / T**2
B = self.b * p / self.R / T
# Solve the cubic equation for compressibility factor z
coeffs = [1, -(1-B), A-2*B-3*B**2, -(A*B-B**2-B**3)]
#print(coeffs)
roots = np.roots(coeffs)
real_roots = roots[np.isreal(roots)].real
valid_roots = real_roots[real_roots > p*self.b/self.R/T]
return valid_roots
def fit_cubic(y0, y1, g0, g1):
"""Fit cubic polynomial to function values and derivatives at x = 0, 1.
Returns position and function value of minimum if fit succeeds. Fit does
not succeeds if
1. polynomial doesn't have extrema or
2. maximum is from (0,1) or
3. maximum is closer to 0.5 than minimum
"""
a = 2*(y0-y1)+g0+g1
b = -3*(y0-y1)-2*g0-g1
p = np.array([a, b, g0, y0])
r = np.roots(np.polyder(p))
if not np.isreal(r).all():
return None, None
r = sorted(x.real for x in r)
if p[0] > 0:
maxim, minim = r
else:
minim, maxim = r
if 0 < maxim < 1 and abs(minim-0.5) > abs(maxim-0.5):
return None, None
return minim, np.polyval(p, minim)
def get_mu_tensor(self):
const_fact = self._dist_to_opt_avg**2 * self._h_min**2 / 2 / self._grad_var
coef = tf.Variable([-1.0, 3.0, 0.0, 1.0], dtype=tf.float32, name="cubic_solver_coef")
coef = tf.scatter_update(coef, tf.constant(2), -(3 + const_fact) )
roots = tf.py_func(np.roots, [coef], Tout=tf.complex64, stateful=False)
# filter out the correct root
root_idx = tf.logical_and(tf.logical_and(tf.greater(tf.real(roots), tf.constant(0.0) ),
tf.less(tf.real(roots), tf.constant(1.0) ) ), tf.less(tf.abs(tf.imag(roots) ), 1e-5) )
# in case there are two duplicated roots satisfying the above condition
root = tf.reshape(tf.gather(tf.gather(roots, tf.where(root_idx) ), tf.constant(0) ), shape=[] )
tf.assert_equal(tf.size(root), tf.constant(1) )
dr = self._h_max / self._h_min
mu = tf.maximum(tf.real(root)**2, ( (tf.sqrt(dr) - 1)/(tf.sqrt(dr) + 1) )**2)
return mu
def polyroots(p, realroots=False, condition=lambda r: True):
"""
Returns the roots of a polynomial with coefficients given in p.
p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
INPUT:
p - Rank-1 array-like object of polynomial coefficients.
realroots - a boolean. If true, only real roots will be returned and the
condition function can be written assuming all roots are real.
condition - a boolean-valued function. Only roots satisfying this will be
returned. If realroots==True, these conditions should assume the roots
are real.
OUTPUT:
A list containing the roots of the polynomial.
NOTE: This uses np.isclose and np.roots"""
roots = np.roots(p)
if realroots:
roots = [r.real for r in roots if isclose(r.imag, 0)]
roots = [r for r in roots if condition(r)]
duplicates = []
for idx, (r1, r2) in enumerate(combinations(roots, 2)):
if isclose(r1, r2):
duplicates.append(idx)
return [r for idx, r in enumerate(roots) if idx not in duplicates]
def polyroots(p, realroots=False, condition=lambda r: True):
"""
Returns the roots of a polynomial with coefficients given in p.
p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
INPUT:
p - Rank-1 array-like object of polynomial coefficients.
realroots - a boolean. If true, only real roots will be returned and the
condition function can be written assuming all roots are real.
condition - a boolean-valued function. Only roots satisfying this will be
returned. If realroots==True, these conditions should assume the roots
are real.
OUTPUT:
A list containing the roots of the polynomial.
NOTE: This uses np.isclose and np.roots"""
roots = np.roots(p)
if realroots:
roots = [r.real for r in roots if isclose(r.imag, 0)]
roots = [r for r in roots if condition(r)]
duplicates = []
for idx, (r1, r2) in enumerate(combinations(roots, 2)):
if isclose(r1, r2):
duplicates.append(idx)
return [r for idx, r in enumerate(roots) if idx not in duplicates]
def calc_cp(alpha,beta,k):
gamma = []
for i in range(k+1):
num = factorial(alpha+k)*factorial(alpha+beta+k+i)
denom = factorial(alpha+i)*factorial(k-i)*factorial(i)
gamma.insert(i,num/denom)
poly = []
for i in range(k+1):
if i == 0:
poly.insert(i,gamma[i])
else:
prod = [1]
j=1
while j<=i:
prod=conv(prod,[1,-1])
j=j+1
while len(poly)<len(prod):
poly.insert(0,0)
prod = [gamma[i]*t for t in prod]
poly = [sum(pair) for pair in zip(poly,prod)]
cp = numpy.roots(poly)
return cp
def pitch_from_relden(relden, cf, sw):
"""
This function calculates the pitch of cuboct of a given relative density, chamfer factor, and strut width.
:param relden: float. Desired relative density
:param cf: float. Chamfer factor of voxel
:param sw: float. strut width of voxel
:return: lattice pitch
"""
chamheight = float(sw) / cf
l_2 = sw / 2.0 + chamheight
l_3 = l_2 + sw * np.cos(math.radians(45)) # horizontal position of points
l_4 = np.sqrt(2) * (l_3 - sw / 2.0)
tan_theta = ((l_3 - l_2) / ((l_4 / 2.0) - (np.sqrt(2) * chamheight / 2.0)))
v1 = l_2 * (sw * sw + 4 * sw * (l_3 - sw / 2.0) + 2 * (l_3 - sw / 2.0) * (l_3 - sw / 2.0))
h = (l_4 / 2.0) * tan_theta
hs = chamheight * tan_theta * np.sqrt(2) / 2.0
v2 = ((l_4 * l_4 * h) - (2 * (chamheight * chamheight * hs))) / 3.0
v3 = 4 * sw * (0.5 * (l_3 - l_2) * (l_3 - l_2) + (l_3 - l_2) * chamheight)
v4 = sw * sw * (l_3 - l_2)
node_volume = v1 + v2 + v3 + v4
c1 = relden
c2 = (-6) * np.sqrt(2)*sw *sw
c3 = -6*node_volume + 12*sw*sw*np.sqrt(2)*(l_2 + l_3)
return max(np.roots([c1, 0, c2, c3]))
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 test_roots(self):
assert_array_equal(np.roots([1, 0, 0]), [0, 0])
def GetRootsPolynomial(a, b, xE, yE, xC, yC, r):
'''
'''
# Define some stuff
r2 = r * r;
a2 = a * a;
b2 = b * b;
a2b2 = a2 / b2;
x0 = xE - xC;
y0 = yE - yC;
y2 = y0 * y0;
x2 = x0 * x0;
# Get the coefficients
A = a2b2 - 1.;
B = -2. * x0 * a2b2;
C = r2 - y2 - a2 + a2b2 * x2;
D = 4. * y2 * a2b2;
c4 = A * A;
c3 = 2. * A * B;
c2 = 2. * A * C + B * B + D;
c1 = 2. * B * C - 2. * D * x0;
c0 = C * C - (b2 - x2) * D;
# Get the real roots
roots = [r.real + xC for r in np.roots([c4, c3, c2, c1, c0])
if np.abs(r.imag) < tol]
return roots
def filter_is_stable(a):
"""
Check if filter coefficients of IIR filter are stable.
Parameters
----------
a: list or 1darray of number
Denominator filter coefficients a.
Returns
-------
is_stable: bool
Filter is stable or not.
Notes
----
Filter is stable if absolute value of all roots is smaller than 1,
see [1]_.
References
----------
.. [1] HYRY, "SciPy 'lfilter' returns only NaNs" StackOverflow,
http://stackoverflow.com/a/8812737/1469195
"""
assert a[0] == 1.0, (
"a[0] should normally be zero, did you accidentally supply b?\n"
"a: {:s}".format(str(a)))
# from http://stackoverflow.com/a/8812737/1469195
return np.all(np.abs(np.roots(a))<1)
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 get_Z(self, T, p):
"""
Returns the compressibility factor of a real gas.
:param T: temperature (K)
:type T: float
:param p: pressure (Pa)
:type p: float
:return: compressibility factor
:rtype: float
"""
# a = self.get_a(T)
# coeffs = [
# 1,
# (-self.R*T - self.b*p + self.c*p + self.d*p)/(self.R*T),
# (-self.R*T*self.d*p + a*p - self.b*self.d*p**2 + self.c*self.d*p**2 + self.e*p**2)/(self.R**2*T**2),
# (-self.R*T*self.e*p**2 - a*self.b*p**2 + a*self.c*p**2 - self.b*self.e*p**3 + self.c*self.e*p**3)/(self.R**3*T**3)
# ]
coeffs = [1, self.get_A(T, p), self.get_B(T, p), self.get_C(T, p)]
#print(coeffs)
roots = np.roots(coeffs)
real_roots = roots[np.isreal(roots)].real
if len(real_roots) == 1:
real_roots = real_roots[0]
return real_roots
# Partial derivative of Z wrt T at constant p
def get_Z(self, T, p):
A = self.get_A(T, p)
B = self.get_B(T, p)
# Solve the cubic equation for compressibility factor z
coeffs = [1, -1, A-B-B**2, -A*B]
roots = np.roots(coeffs)
real_roots = roots[np.isreal(roots)].real
valid_roots = real_roots[real_roots > p*self.b/self.R/T]
return valid_roots
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 test_roots(self):
assert_array_equal(np.roots([1, 0, 0]), [0, 0])
def fit_quartic(y0, y1, g0, g1):
"""Fit constrained quartic polynomial to function values and erivatives at x = 0,1.
Returns position and function value of minimum or None if fit fails or has
a maximum. Quartic polynomial is constrained such that it's 2nd derivative
is zero at just one point. This ensures that it has just one local
extremum. No such or two such quartic polynomials always exist. From the
two, the one with lower minimum is chosen.
"""
def g(y0, y1, g0, g1, c):
a = c+3*(y0-y1)+2*g0+g1
b = -2*c-4*(y0-y1)-3*g0-g1
return np.array([a, b, c, g0, y0])
def quart_min(p):
r = np.roots(np.polyder(p))
is_real = np.isreal(r)
if is_real.sum() == 1:
minim = r[is_real][0].real
else:
minim = r[(r == max(-abs(r))) | r == -max(-abs(r))][0].real
return minim, np.polyval(p, minim)
D = -(g0+g1)**2-2*g0*g1+6*(y1-y0)*(g0+g1)-6*(y1-y0)**2 # discriminant of d^2y/dx^2=0
if D < 1e-11:
return None, None
else:
m = -5*g0-g1-6*y0+6*y1
p1 = g(y0, y1, g0, g1, .5*(m+np.sqrt(2*D)))
p2 = g(y0, y1, g0, g1, .5*(m-np.sqrt(2*D)))
if p1[0] < 0 and p2[0] < 0:
return None, None
[minim1, minval1] = quart_min(p1)
[minim2, minval2] = quart_min(p2)
if minval1 < minval2:
return minim1, minval1
else:
return minim2, minval2
def signalMinimum(img, fitParams=None, n_std=3):
'''
intersection between signal and background peak
'''
if fitParams is None:
fitParams = FitHistogramPeaks(img).fitParams
assert len(fitParams) > 1, 'need 2 peaks so get minimum signal'
i = signalPeakIndex(fitParams)
signal = fitParams[i]
bg = getBackgroundPeak(fitParams)
smn = signal[1] - n_std * signal[2]
bmx = bg[1] + n_std * bg[2]
if smn > bmx:
return smn
# peaks are overlapping
# define signal min. as intersection between both Gaussians
def solve(p1, p2):
s1, m1, std1 = p1
s2, m2, std2 = p2
a = (1 / (2 * std1**2)) - (1 / (2 * std2**2))
b = (m2 / (std2**2)) - (m1 / (std1**2))
c = (m1**2 / (2 * std1**2)) - (m2**2 / (2 * std2**2)) - \
np.log(((std2 * s1) / (std1 * s2)))
return np.roots([a, b, c])
i = solve(bg, signal)
try:
return i[np.logical_and(i > bg[1], i < signal[1])][0]
except IndexError:
# this error shouldn't occur... well
return max(smn, bmx)
def getSignalMinimum(fitParams, n_std=3):
assert len(fitParams) > 0, 'need min. 1 peak so get minimum signal'
if len(fitParams) == 1:
signal = fitParams[0]
return signal[1] - n_std * signal[2]
i = signalPeakIndex(fitParams)
signal = fitParams[i]
bg = fitParams[i - 1]
#bg = getBackgroundPeak(fitParams)
smn = signal[1] - n_std * signal[2]
bmx = bg[1] + n_std * bg[2]
if smn > bmx:
return smn
# peaks are overlapping
# define signal min. as intersection between both Gaussians
def solve(p1, p2):
s1, m1, std1 = p1
s2, m2, std2 = p2
a = (1 / (2 * std1**2)) - (1 / (2 * std2**2))
b = (m2 / (std2**2)) - (m1 / (std1**2))
c = (m1**2 / (2 * std1**2)) - (m2**2 / (2 * std2**2)) - \
np.log(((std2 * s1) / (std1 * s2)))
return np.roots([a, b, c])
i = solve(bg, signal)
try:
return i[np.logical_and(i > bg[1], i < signal[1])][0]
except IndexError:
# something didnt work out - fallback
return smn
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 update_X(self, X, mu, k=20):
U, S, VT = svdp(X, k=k)
P = np.c_[np.ones((k, 1)), 1-S, 1./2./mu-S]
sigma_star = np.zeros(k)
for t in range(k):
p = P[t, :]
delta = p[1]**2 - 4 * p[0] * p[2]
if delta <= 0:
sigma_star[t] = 0.
else:
solution = np.roots(p)
solution = sorted(solution, key=abs)
solution = np.array(solution)
if solution[0] * solution[1] <= 0:
sigma_star[t] = solution[1]
elif solution[1] < 0:
sigma_star[t] = 0.
else:
f = np.log(1 + solution[1]) + mu * (solution[1] - s[t])**2
if f > mu * s[t]**2:
sigma_star[t] = 0.
else:
sigma_star[t] = solution[1]
sigma_star = np.diag(sigma_star)
sigma_star = np.dot(np.dot(U, sigma_star), VT)
return sigma_star
def update_X(X, mu, k=6):
#U, S, VT = svdp(X, k=k)
U, S, VT = svds(X, k=k, which='LM')
P = np.c_[np.ones((k, 1)), 1-S, 1./2./mu-S]
sigma_star = np.zeros(k)
for t in range(k):
p = P[t, :]
delta = p[1]**2 - 4 * p[0] * p[2]
if delta <= 0:
sigma_star[t] = 0.
else:
solution = np.roots(p)
solution = solution.tolist()
solution.sort(key=abs)
solution = np.array(solution)
if solution[0] * solution[1] <= 0:
sigma_star[t] = solution[1]
elif solution[1] < 0:
sigma_star[t] = 0.
else:
f = np.log(1 + solution[1]) + mu * (solution[1] - s[t])**2
if f > mu *s[t]**2:
sigma_star[t] = 0.
else:
sigma_star[t] = solution[1]
sigma_star = sp.csr_matrix(np.diag(sigma_star))
sigma_star = safe_sparse_dot(safe_sparse_dot(U, sigma_star), VT)
sigma_star[abs(sigma_star)<1e-10] = 0
return sp.lil_matrix(sigma_star)
def real_roots(coeffs):
"""Get real roots of a polynomial.
Args:
coeffs (List[Float]): List of polynomial coefficients.
Returns:
numpy.ndarray: The (sorted) real roots of the polynomial.
"""
all_roots = np.roots(coeffs)
filtered = all_roots[all_roots.imag == 0.0].real
return np.sort(filtered)
def roots(a):
return matlabarray(np.roots(np.asarray(a).ravel()))
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()
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)