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
python类polyval()的实例源码
def __call__(self, x, model_dispersion, model_normalized_flux, *parameters):
# Marginalize over all models.
assert len(parameters) == len(self.parameter_names)
# Continuum.
O = self.metadata["continuum_order"]
continuum = 1.0 if 0 > O \
else np.polyval(parameters[-(O + 1):][::-1], model_dispersion)
y = model_normalized_flux * continuum
# Smoothing?
try:
index = self.parameter_names.index("smoothing")
except ValueError:
None
else:
kernel = abs(parameters[index])
y = ndimage.gaussian_filter1d(y, kernel, axis=-1)
# Redshift?
try:
index = self.parameter_names.index("redshift")
except ValueError:
z = 0
else:
v = parameters[index]
z = v/299792.458 # km/s
return np.interp(x, model_dispersion * (1 + z), y)
def FIT(p, seq):
k = len(p)
predict = polyval(p, k+1)
return predict
def project_xy(xy_coords, pvec):
# get cubic polynomial coefficients given
#
# f(0) = 0, f'(0) = alpha
# f(1) = 0, f'(1) = beta
alpha, beta = tuple(pvec[CUBIC_IDX])
poly = np.array([
alpha + beta,
-2*alpha - beta,
alpha,
0])
xy_coords = xy_coords.reshape((-1, 2))
z_coords = np.polyval(poly, xy_coords[:, 0])
objpoints = np.hstack((xy_coords, z_coords.reshape((-1, 1))))
image_points, _ = cv2.projectPoints(objpoints,
pvec[RVEC_IDX],
pvec[TVEC_IDX],
K, np.zeros(5))
return image_points
def calcmlmags(self, lightcurve):
"""
Returns a "lc.mags"-like array made using the ml-parameters.
It has the same size as lc.mags, and contains the microlensing to be added to them.
The lightcurve object is not changed !
For normal use, call getmags() from the lightcurve.
Idea : think about only returning the seasons mags to speed it up ? Not sure if reasonable, as no seasons defined outside ?
"""
jds = lightcurve.jds[self.season.indices] # Is this already a copy ? It seems so. So no need for an explicit copy().
# We do not need to apply shifts (i.e. getjds()), as anyway we "center" the jds.
# Old method :
if self.mltype == "poly":
refjd = np.mean(jds)
jds -= refjd # This is apparently safe, it does not shifts the lightcurves jds.
allmags = np.zeros(len(lightcurve.jds))
allmags[self.season.indices] = np.polyval(self.params, jds) # probably faster then +=
return allmags
# Legendre polynomials :
if self.mltype == "leg":
rjd = (np.max(jds) - np.min(jds))/2.0
cjd = (np.max(jds) + np.min(jds))/2.0
jds = (jds - cjd)/rjd
allmags = np.zeros(len(lightcurve.jds))
for (n, p) in enumerate(self.params):
allmags[self.season.indices] += p * ss.legendre(n)(jds)
return allmags
def test_polyfit(self):
c = np.array([3., 2., 1.])
x = np.linspace(0, 2, 7)
y = np.polyval(c, x)
err = [1, -1, 1, -1, 1, -1, 1]
weights = np.arange(8, 1, -1)**2/7.0
# check 1D case
m, cov = np.polyfit(x, y+err, 2, cov=True)
est = [3.8571, 0.2857, 1.619]
assert_almost_equal(est, m, decimal=4)
val0 = [[2.9388, -5.8776, 1.6327],
[-5.8776, 12.7347, -4.2449],
[1.6327, -4.2449, 2.3220]]
assert_almost_equal(val0, cov, decimal=4)
m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
val = [[8.7929, -10.0103, 0.9756],
[-10.0103, 13.6134, -1.8178],
[0.9756, -1.8178, 0.6674]]
assert_almost_equal(val, cov2, decimal=4)
# check 2D (n,1) case
y = y[:, np.newaxis]
c = c[:, np.newaxis]
assert_almost_equal(c, np.polyfit(x, y, 2))
# check 2D (n,2) case
yy = np.concatenate((y, y), axis=1)
cc = np.concatenate((c, c), axis=1)
assert_almost_equal(cc, np.polyfit(x, yy, 2))
m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
assert_almost_equal(est, m[:, 0], decimal=4)
assert_almost_equal(est, m[:, 1], decimal=4)
assert_almost_equal(val0, cov[:, :, 0], decimal=4)
assert_almost_equal(val0, cov[:, :, 1], decimal=4)
def _ir_calibrate(self, data):
"""IR calibration
"""
cwl = self._header['block5']["central_wave_length"][0] * 1e-6
c__ = self._header['calibration']["speed_of_light"][0]
h__ = self._header['calibration']["planck_constant"][0]
k__ = self._header['calibration']["boltzmann_constant"][0]
a__ = (h__ * c__) / (k__ * cwl)
#b__ = ((2 * h__ * c__ ** 2) / (1.0e6 * cwl ** 5 * data.data)) + 1
data.data[:] *= 1.0e6 * cwl ** 5
data.data[:] **= -1
data.data[:] *= (2 * h__ * c__ ** 2)
data.data[:] += 1
#Te_ = a__ / np.log(b__)
data.data[:] = a__ / np.log(data.data)
c0_ = self._header['calibration']["c0_rad2tb_conversion"][0]
c1_ = self._header['calibration']["c1_rad2tb_conversion"][0]
c2_ = self._header['calibration']["c2_rad2tb_conversion"][0]
#data.data[:] = c0_ + c1_ * Te_ + c2_ * Te_ ** 2
data.data[:] = np.polyval([c2_, c1_, c0_], data.data)
data.mask[data.data < 0] = True
data.mask[np.isnan(data.data)] = True
def _discretize(self, observation):
if not self._is_discrete:
buckets = np.array(observation, dtype=np.int)
for i in range(len(observation)):
buckets[i] = np.digitize(observation[i], self._separators[i])
observation = np.polyval(buckets, self._bins)
return observation
def test_polyfit(self):
c = np.array([3., 2., 1.])
x = np.linspace(0, 2, 7)
y = np.polyval(c, x)
err = [1, -1, 1, -1, 1, -1, 1]
weights = np.arange(8, 1, -1)**2/7.0
# check 1D case
m, cov = np.polyfit(x, y+err, 2, cov=True)
est = [3.8571, 0.2857, 1.619]
assert_almost_equal(est, m, decimal=4)
val0 = [[2.9388, -5.8776, 1.6327],
[-5.8776, 12.7347, -4.2449],
[1.6327, -4.2449, 2.3220]]
assert_almost_equal(val0, cov, decimal=4)
m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
val = [[8.7929, -10.0103, 0.9756],
[-10.0103, 13.6134, -1.8178],
[0.9756, -1.8178, 0.6674]]
assert_almost_equal(val, cov2, decimal=4)
# check 2D (n,1) case
y = y[:, np.newaxis]
c = c[:, np.newaxis]
assert_almost_equal(c, np.polyfit(x, y, 2))
# check 2D (n,2) case
yy = np.concatenate((y, y), axis=1)
cc = np.concatenate((c, c), axis=1)
assert_almost_equal(cc, np.polyfit(x, yy, 2))
m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
assert_almost_equal(est, m[:, 0], decimal=4)
assert_almost_equal(est, m[:, 1], decimal=4)
assert_almost_equal(val0, cov[:, :, 0], decimal=4)
assert_almost_equal(val0, cov[:, :, 1], decimal=4)
def op(i, j, x, y):
p0, a, b, c = orthopy.line.recurrence_coefficients.jacobi(
i, 0, 0,
# standardization='monic'
standardization='p(1)=(n+alpha over n)'
)
val1 = orthopy.line.tools.evaluate_orthogonal_polynomial(
(x-y)/(x+y), p0, a, b, c
)
val1 = numpy.polyval(scipy.special.jacobi(i, 0, 0), (x-y)/(x+y))
# treat x==0, y==0 separately
if isinstance(val1, numpy.ndarray):
idx = numpy.where(numpy.logical_and(x == 0, y == 0))[0]
val1[idx] = numpy.polyval(scipy.special.jacobi(i, 0, 0), 0.0)
else:
if numpy.isnan(val1):
val1 = numpy.polyval(scipy.special.jacobi(i, 0, 0), 0.0)
p0, a, b, c = orthopy.line.recurrence_coefficients.jacobi(
j, 2*i+1, 0,
# standardization='monic'
standardization='p(1)=(n+alpha over n)'
)
val2 = orthopy.line.tools.evaluate_orthogonal_polynomial(
1-2*(x+y), p0, a, b, c
)
# val2 = numpy.polyval(scipy.special.jacobi(j, 2*i+1, 0), 1-2*(x+y))
flt = numpy.vectorize(float)
return flt(
numpy.sqrt(2*i + 1) * val1 * (x+y)**i
* numpy.sqrt(2*j + 2*i + 2) * val2
)
def test_eval(t, ref, tol=1.0e-14):
n = 5
p0, a, b, c = orthopy.line.recurrence_coefficients.legendre(
n, 'monic', symbolic=True
)
value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, c)
assert value == ref
# Evaluating the Legendre polynomial in this way is rather unstable, so
# don't go too far with n.
approx_ref = numpy.polyval(legendre(n, monic=True), t)
assert abs(value - approx_ref) < tol
return
def test_eval_vec(t, ref, tol=1.0e-14):
n = 5
p0, a, b, c = orthopy.line.recurrence_coefficients.legendre(
n, 'monic', symbolic=True
)
value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, c)
assert (value == ref).all()
# Evaluating the Legendre polynomial in this way is rather unstable, so
# don't go too far with n.
approx_ref = numpy.polyval(legendre(n, monic=True), t)
assert (abs(value - approx_ref) < tol).all()
return
def test_clenshaw(tol=1.0e-14):
n = 5
_, _, alpha, beta = \
orthopy.line.recurrence_coefficients.legendre(n, 'monic')
t = 1.0
a = numpy.ones(n+1)
value = orthopy.line.clenshaw(a, alpha, beta, t)
ref = math.fsum([
numpy.polyval(legendre(i, monic=True), t)
for i in range(n+1)])
assert abs(value - ref) < tol
return
def polyval(fit, points, der=0, avg=False):
"""Evaluate polynomial generated by ``polyfit()`` on `points`.
Parameters
----------
fit, points : see polyfit()
der : int, optional
Derivative order. Only for 1D, uses np.polyder().
avg : bool, optional
Internal hack, only used by ``avgpolyval()``.
Notes
-----
For 1D we provide "analytic" derivatives using np.polyder(). For ND, we
didn't implement an equivalent machinery. For 2D, you might get away with
fitting a bispline (see Interpol2D) and use it's derivs. For ND, try rbf.py's RBF
interpolator which has at least 1st derivatives for arbitrary dimensions.
See Also
--------
:class:`PolyFit`, :class:`PolyFit1D`, :func:`polyfit`
"""
pscale, pmin = fit['pscale'], fit['pmin']
vscale, vmin = fit['vscale'], fit['vmin']
if der > 0:
assert points.shape[1] == 1, "deriv only for 1d poly (ndim=1)"
# ::-1 b/c numpy stores poly coeffs in reversed order
dcoeffs = np.polyder(fit['coeffs'][::-1], m=der)
return np.polyval(dcoeffs, (points[:,0] - pmin[0,0]) / pscale[0,0]) / \
pscale[0,0]**der * vscale
else:
vand = vander((points - pmin) / pscale, fit['deg'])
if avg:
return np.dot(vand, fit['coeffs']) * vscale
else:
return np.dot(vand, fit['coeffs']) * vscale + vmin
def test_compare_numpy():
x = np.sort(np.random.rand(10))
y = np.random.rand(10)
yy1 = np.polyval(np.polyfit(x, y, 3), x)
for scale in [True,False]:
yy2 = num.PolyFit1D(x, y, 3, scale=scale)(x)
assert np.allclose(yy1, yy2)
def Legendre_matrix(N, ctheta):
r"""Matrix of weighted Legendre Polynominals.
Computes a matrix of weighted Legendre polynominals up to order N for
the given angles
.. math::
L_n(\theta) = \frac{2n+1}{4 \pi} P_n(\theta)
Parameters
----------
N : int
Maximum order.
ctheta : (Q,) array_like
Angles.
Returns
-------
Lmn : ((N+1), Q) numpy.ndarray
Matrix containing Legendre polynominals.
"""
if ctheta.ndim == 0:
M = 1
else:
M = len(ctheta)
Lmn = np.zeros([N+1, M], dtype=complex)
for n in range(N+1):
Lmn[n, :] = (2*n+1)/(4*np.pi) * np.polyval(special.legendre(n), ctheta)
return Lmn
def eval_fibmodel_poly(self, use_poly=False):
self.fibmodel = np.zeros((self.D, len(self.binx)))
for i in xrange(len(self.binx)):
if use_poly:
self.fibmodel[:,i] = np.polyval(self.fibmodel_polyvals[:,i],
1.* np.arange(self.D) / self.D)
else:
self.fibmodel[:,i] = np.interp(np.arange(self.D),
self.fibmodel_x,
self.fibmodel_y[:,i])
def eval_wave_poly(self):
self.wavelength = np.polyval(self.wave_polyvals,
1.* np.arange(self.D) / self.D)
