def compute_beta_splines(TT, minT, splinesoptions={}):
"""
This function computes the volumetric thermal expansion as a numerical
derivative of the volume as a function of temperature V(T) given in the
input array *minT*. This array can obtained
from the free energy minimization which should be done before.
This function uses splines for smoothing the derivative.
"""
betaT = np.zeros(len(TT))
x = np.array(TT)
y0 = np.array(minT)
if (splinesoptions == {}):
tck0 = interpolate.splrep(x, y0)
else:
tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])
ynew0 = interpolate.splev(x, tck0, der=0)
betaT = interpolate.splev(x, tck0, der=1)
return betaT/minT
python类splrep()的实例源码
def resample_photosphere(opacities, photosphere, opacity_index):
""" Resample photospheric quantities onto a new opacity scale. """
if opacity_index is None:
return photosphere
resampled_photosphere = np.zeros(photosphere.shape)
n_quantities = photosphere.shape[1]
for i in range(n_quantities):
if i == opacity_index: continue
# Create spline function.
tk = \
interpolate.splrep(photosphere[:, opacity_index], photosphere[:, i])
# Evaluate photospheric quantities at the new opacities
resampled_photosphere[:, i] = interpolate.splev(opacities.flatten(), tk)
# Update photosphere with new opacities
resampled_photosphere[:, opacity_index] = opacities
return resampled_photosphere
def resample_1d(data, npoints, smooth = 0, periodic = False, derivative = 0):
"""Resample 1d data using n equidistant points
Arguments:
data (array): data points
npoints (int): number of points in equidistant resampling
smooth (number): smoothness factor
periodic (bool): if True assumes the curve is a closed curve
Returns:
(array): resampled data points
"""
x = np.linspace(0, 1, data.shape[0]);
dinterp = splrep(x, data, s = smooth, per = periodic);
if npoints is all:
npoints = data.shape[0];
x2 = np.linspace(0, 1, npoints);
return splev(x2, dinterp, der = derivative)
def phi(self, points = None):
"""Returns a Spline representing the tangent angle along a 2d curve
Arguments:
points (int, array or None): sample points used to determine the phi
Returns:
Spline: spline of phi
"""
if self.ndim != 2:
raise RuntimeError('phi angle can only be computed for 2d curves');
points = self.get_points(points, error = 'cannot determine sample points needed for the calculation of phi');
#get the tangents and tangent angles
tgs = splev(points, self.tck(), der = 1);
tgs = np.vstack(tgs).T;
phi = np.arctan2(tgs[:,1], tgs[:,0]);
phi = np.mod(phi + np.pi, 2 * np.pi) - np.pi;
#return Spline(phi, points = self.points, knots = self.knots, degree = self.degree - 1);
tck = splrep(points, phi, s = 0.0, k = self.degree + 1);
return Spline(tck = tck, points = points);
def set_parameter_from_values(self, values, points = None):
"""Change the parameter to represent the values and resmaple on current sample points if necessary
Arguments:
values (array): values
points (array or None): sample points for the values
"""
if points is None:
self.set_values(values);
else:
if points is self._points or (self._points is not None and self._points.shape[0] == points.shape[0] and np.allclose(points, self._points)):
self.set_values(values);
else:
tck = splrep(points, values, t = self._knots, task = -1, k = self.degree);
self._parameter = tck[1][:self._nparameter];
# values changed
self._values = None;
def set_parameter_from_values(self, values, points = None):
"""Change the parameter to represent the values and resmaple on current sample points if necessary
Arguments:
values (array): values
points (array or None): sample points for the values
"""
if points is None:
self.set_values(values);
else:
if points is self._points or (self._points is not None and self._points.shape[0] == points.shape[0] and np.allclose(points, self._points)):
self.set_values(values);
else:
tck = splrep(points, values, t = self._knots, task = -1, k = self.degree);
self._parameter = tck[1][:self._nparameter];
# values changed
self._values = None;
def optc(self):
"""
Optimize the coeffs, don't touch the knots
This is the fast guy, one reason to use splines :-)
Returns the chi2 in case you want it (including stabilization points) !
Sets lastr2stab, but not lastr2nostab !
"""
out = si.splrep(self.datapoints.jds, self.datapoints.mags, w=1.0/self.datapoints.magerrs, xb=None, xe=None, k=self.k, task=-1, s=None, t=self.getintt(), full_output=1, per=0, quiet=1)
# We check if it worked :
if not out[2] <= 0:
raise RuntimeError("Problem with spline representation, message = %s" % (out[3]))
self.c = out[0][1] # save the coeffs
#import matplotlib.pyplot as plt
#plt.plot(self.datapoints.jds, self.datapoints.magerrs)
#plt.show()
self.lastr2stab = out[1]
return out[1]
def klgfbInterp(self, wvl, n, l):
'''A Python version of the CHIANTI IDL procedure karzas_xs.
Interpolates free-bound gaunt factor of Karzas and Latter, (1961, Astrophysical Journal
Supplement Series, 6, 167) as a function of wavelength (wvl).
'''
try:
klgfb = self.Klgfb
except:
self.Klgfb = ch_io.klgfbRead()
klgfb = self.Klgfb
# get log of photon energy relative to the ionization potential
sclE = np.log(self.Ip/(wvl*ch_const.ev2ang))
thisGf = klgfb['klgfb'][n-1, l]
spl = splrep(klgfb['pe'], thisGf)
gf = splev(sclE, spl)
return gf
def ioneq_one(self, stage, **kwargs):
"""
Calculate the equilibrium fractional ionization of the ion as a function of temperature.
Uses the `ChiantiPy.core.ioneq` module and does a first-order spline interpolation to the data. An
ionization equilibrium file can be passed as a keyword argument, `ioneqfile`. This can
be passed through as a keyword argument to any of the functions that uses the
ionization equilibrium.
Parameters
----------
stage : int
Ionization stage, e.g. 25 for Fe XXV
"""
tmp = ioneq(self.Z)
tmp.load(ioneqName=kwargs.get('ioneqfile', None))
ionization_equilibrium = splev(self.Temperature,
splrep(tmp.Temperature, tmp.Ioneq[stage-1,:], k=1), ext=1)
return np.where(ionization_equilibrium < 0., 0., ionization_equilibrium)
def setDriftData(self, alpha, CL, CDi, CM):
self.alpha = alpha
self.CL = CL
self.CDi = CDi
self.CM = CM
# Calculate coefficients for lift
fitParams, fitCovariances = curve_fit(liftFunc, self.alpha, self.CL)
self.CL_a1 = fitParams[0]
self.CL_a2 = fitParams[1]
# Calculate coefficients for induced drag
fitParams, fitCovariances = curve_fit(inducedDragFunc, self.alpha, self.CDi)
self.CDi_a2 = fitParams[0]
# Spline interpolation for moment
self.CMspl = interpolate.splrep(self.alpha, self.CM, s=1e-9)
def setStabilityData(self, heelAngles, GZ):
self.heelAngles = heelAngles
self.GZ = GZ
self.restoringMoment = self.GZ*self.Volume*self.rho*self.g
self.restoringMomentSpl = interpolate.splrep(self.heelAngles, self.restoringMoment)
# Calculate maximum heel angle
nTest = 100
heelAnglesTest = np.linspace(0, np.max(self.heelAngles), nTest)
restoringMomentTest = interpolate.splev(heelAnglesTest, self.restoringMomentSpl)
iMax = np.argmax(restoringMomentTest)
self.maxHeelAngle = heelAnglesTest[iMax]
def alt_sky_model(self):
fac = 10
mn = 1e9
mx = 0.
for fiber in self.good_fibers:
mn = np.min([mn, fiber.wavelength.min()])
mx = np.max([mx, fiber.wavelength.max()])
xs = np.linspace(0, 1, self.D*fac)
A = np.zeros((len(xs), len(self.good_fibers)))
for i, fiber in enumerate(self.good_fibers):
y = fiber.spectrum / fiber.fiber_to_fiber
xp = np.interp(fiber.wavelength, np.linspace(mn, mx, self.D*fac),
xs, left=0.0, right=0.0)
tck = splrep(xp, y)
A[:, i] = splev(xs, tck)
ys = biweight_location(A, axis=(1,))
self.masterwave = np.linspace(mn, mx, self.D*fac)
B, c = bspline_x0(self.masterwave, nknots=self.D)
sol = np.linalg.lstsq(c, ys)[0]
self.mastersky = np.dot(c, sol)
def derivatives(xs,ys,ders=(1,)):
'''
Calculate the numerical derivatives of `ys` at points `xs` by use of the third-order spline interpolation.
Parameters
----------
xs : 1d ndarray
The sample points of the argument.
ys: 1d ndarray
The corresponding sample points of the function.
ders : tuple of integer
The derivatives to calculate.
Returns
-------
2d ndarray
The derivatives at the sample points of the argument.
'''
assert len(xs)==len(ys)
xs,ys=np.asarray(xs),np.asarray(ys)
result=np.zeros((len(ders),len(xs)),dtype=ys.dtype)
tck=ip.splrep(xs,ys,s=0)
for i,der in enumerate(ders):
result[i]=ip.splev(xs,tck,der=der)
return result
def _setup_channels(self):
"""
Setup channel, specifically the wavelength or temperature response functions.
Notes
-----
This should be replaced once the response functions are available in SunPy. Probably should
configure wavelength response function interpolators also.
"""
aia_fn = pkg_resources.resource_filename('synthesizAR', 'instruments/data/sdo_aia.json')
with open(aia_fn, 'r') as f:
aia_info = json.load(f)
for channel in self.channels:
channel['name'] = str(channel['wavelength'].value).strip('.0')
channel['instrument_label'] = '{}_{}'.format(self.fits_template['detector'],
channel['telescope_number'])
channel['wavelength_range'] = None
x = aia_info[channel['name']]['temperature_response_x']
y = aia_info[channel['name']]['temperature_response_y']
channel['temperature_response_spline'] = splrep(x, y)
x = aia_info[channel['name']]['response_x']
y = aia_info[channel['name']]['response_y']
channel['wavelength_response_spline'] = splrep(x, y)
def lininterp2(x1, y1, x):
"""Linear interpolation at points x between numpy arrays (x1, y1).
Only y1 is allowed to be two-dimensional. The x1 values should be sorted
from low to high. Returns a numpy.array of y values corresponding to
points x.
"""
return splev(x, splrep(x1, y1, s=0, k=1))
def compute_alpha_splines(TT,minT,ibrav,splinesoptions):
"""
This function calculates the thermal expansions alphaT at different temperatures
as the previous function but using spline interpolation as implemented in
scipy.interpolate.
"""
alphaT = np.zeros(len(TT)*6)
alphaT.shape = (len(TT),6)
x = np.array(TT)
y0 = np.array(minT[:,0])
y1 = np.array(minT[:,1])
y2 = np.array(minT[:,2])
if (splinesoptions=={}):
tck0 = interpolate.splrep(x, y0)
tck1 = interpolate.splrep(x, y1)
tck2 = interpolate.splrep(x, y2)
else:
tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])
tck1 = interpolate.splrep(x, y1, k=splinesoptions['k1'], s=splinesoptions['s1'])
tck2 = interpolate.splrep(x, y2, k=splinesoptions['k2'], s=splinesoptions['s2'])
ynew0 = interpolate.splev(x, tck0, der=0)
alphaT[:,0] = interpolate.splev(x, tck0, der=1)
ynew1 = interpolate.splev(x, tck1, der=0)
alphaT[:,1] = interpolate.splev(x, tck1, der=1)
ynew2 = interpolate.splev(x, tck2, der=0)
alphaT[:,2] = interpolate.splev(x, tck2, der=1)
# now normalize the alphaTs properly. It must be different for different ibrav
# to avoid a divide by 0 error (minT is zero for lattice parameters not defined
# in the system)
if ibrav==4:
alphaT[:,0] = alphaT[:,0]/minT[:,0]
alphaT[:,2] = alphaT[:,2]/minT[:,2]
return alphaT
################################################################################
def compute_blackbody_response(self, Teffs=None):
"""
Computes blackbody intensities across the entire range of
effective temperatures.
@Teffs: an array of effective temperatures. If None, a default
array from ~300K to ~500000K with 97 steps is used. The default
array is uniform in log10 scale.
Returns: n/a
"""
if Teffs == None:
log10Teffs = np.linspace(2.5, 5.7, 97) # this corresponds to the 316K-501187K range.
Teffs = 10**log10Teffs
# Energy-weighted intensities:
log10ints_energy = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=False)) for Teff in Teffs])
self._bb_func_energy = interpolate.splrep(Teffs, log10ints_energy, s=0)
self._log10_Inorm_bb_energy = lambda Teff: interpolate.splev(Teff, self._bb_func_energy)
# Photon-weighted intensities:
log10ints_photon = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=True)) for Teff in Teffs])
self._bb_func_photon = interpolate.splrep(Teffs, log10ints_photon, s=0)
self._log10_Inorm_bb_photon = lambda Teff: interpolate.splev(Teff, self._bb_func_photon)
self.content.append('blackbody')
self.atmlist.append('blackbody')
def burgess_tully_descale(x, y, energy_ratio, c, scaling_type):
"""
Convert scaled Burgess-Tully parameters to physical quantities. For more details see
[1]_.
References
----------
.. [1] Burgess, A. and Tully, J. A., 1992, A&A, `254, 436 <http://adsabs.harvard.edu/abs/1992A%26A...254..436B>`_
"""
nots = splrep(x, y, s=0)
if scaling_type == 1:
x_new = 1.0 - np.log(c)/np.log(energy_ratio + c)
upsilon = splev(x_new, nots, der=0)*np.log(energy_ratio + np.e)
elif scaling_type == 2:
x_new = energy_ratio/(energy_ratio + c)
upsilon = splev(x_new, nots, der=0)
elif scaling_type == 3:
x_new = energy_ratio/(energy_ratio + c)
upsilon = splev(x_new, nots, der=0)/(energy_ratio + 1.0)
elif scaling_type == 4:
x_new = 1.0 - np.log(c)/np.log(energy_ratio + c)
upsilon = splev(x_new, nots, der=0)*np.log(energy_ratio + c)
elif scaling_type == 5:
# dielectronic
x_new = energy_ratio/(energy_ratio + c)
upsilon = splev(x_new, nots, der=0)/energy_ratio
elif scaling_type == 6:
# protons
x_new = energy_ratio/(energy_ratio + c)
upsilon = 10**splev(x_new, nots, der=0)
else:
raise ValueError('Unrecognized BT92 scaling option.')
return upsilon
def from_tck(self, tck, points = None):
"""Change spline parameter and knot structure using a tck object returned by splprep or splrep
Arguments:
tck (tuple): t,c,k tuple returned by splprep
points (int, array or None): optional sample points pecification
"""
t,c,k = tck;
self.degree = k;
# set knots
self._knots = t[k+1:-k-1];
self._knots_all = t;
self._nparameter = self._knots.shape[0] + k + 1;
#set parameter
if isinstance(c, list):
c = np.vstack(c);
elif isinstance(c, np.ndarray) and c.ndim == 1:
c = np.vstack([c]);
c = c[:,:self._nparameter].T;
self.ndim = c.shape[1];
self._parameter = c[:];
#set points
if points is not None:
self.set_points(points);
# values and projection matrix will need to be updated after change of the knots
self._values = None;
self._projection = None;
self._projection_inverse = None;
############################################################################
### Spline getter
def tck_1d(self, dimension = 0):
"""Returns tck tuple for use with 1d spline functions
Arguments:
dimension (int): dimension for which to return tck tuple
Returns:
tuple: tck couple as retyurned by splrep
"""
k = self.degree;
p = np.pad(self._parameter[:,dimension], (0,k+1), 'constant');
return (self._knots_all, p, k);
def set_values(self, values, points = None):
"""Calculate the bspline parameter for the given values and sample points
Arguments:
values (array): values of data points
points (array or None): sample points for the data values, if None use internal points or linspace(0,1,values.shape[0])
"""
values = np.array(values, dtype = float);
#set points
if points is None:
if self._points is None:
self.set_points(values.shape[0]);
else:
self.set_points(points);
if values.shape[0] != self._points.shape[0]:
raise ValueError('number of values %d mismatch number of points %d' % (values.shape[0], self._points.shape[0]));
#set parameter from values
if self.with_projection and self.projection_inverse is not False:
self._parameter = self._projection_inverse.dot(values);
self._values = values;
else:
tck = splrep(self._points, values, t = self._knots, task = -1, k = self.degree);
self._parameter = tck[1][:self._nparameter];
# values changed
self._values = None;
#self._values = values; #fast but imprecise as values on spline due approximation might differ!
def from_tck(self, tck, points = None):
"""Set spline parameter and knot structure using a tck object returned by splrep
Arguments:
tck (tuple): t,c,k tuple as returned by splrep
points (int, array or None): optional sample points pecification
"""
t,c,k = tck;
self.degree = k;
# set knots
self._knots = t[k+1:-k-1];
self._knots_all = t;
self._nparameter = self._knots.shape[0] + k + 1;
# set parameter
self._parameter = c[:self._nparameter];
#set points
if points is not None:
self.set_points(points);
# values and projection matrix will need to be updated after change of the knots
self._values = None;
self._projection = None;
self._projection_inverse = None;
############################################################################
### Spline getter
def tck(self):
"""Returns a tck tuple for use with spline functions
Note:
This returns a tck tuple as splrep
"""
k = self.degree;
p = np.pad(self._parameter, (0,k+1), 'constant');
return (self._knots_all, p, k);
def set_values(self, values, points = None):
"""Calculate the bspline parameter for the given values and sample points
Arguments:
values (array): values of data points
points (array or None): sample points for the data values, if None use internal points or linspace(0,1,values.shape[0])
"""
values = np.array(values, dtype = float);
#set points
if points is None:
if self._points is None:
self.set_points(values.shape[0]);
else:
self.set_points(points);
if values.shape[0] != self._points.shape[0]:
raise ValueError('number of values %d mismatch number of points %d' % (values.shape[0], self._points.shape[0]));
#set parameter from values
if self.with_projection and self.projection_inverse is not False:
self._parameter = self._projection_inverse.dot(values);
self._values = values;
else:
tck = splrep(self._points, values, t = self._knots, task = -1, k = self.degree);
self._parameter = tck[1][:self._nparameter];
# values changed
self._values = None;
#self._values = values; #fast but imprecise as values on spline due approximation might differ!
def from_tck(self, tck, points = None):
"""Set spline parameter and knot structure using a tck object returned by splrep
Arguments:
tck (tuple): t,c,k tuple as returned by splrep
points (int, array or None): optional sample points pecification
"""
t,c,k = tck;
self.degree = k;
# set knots
self._knots = t[k+1:-k-1];
self._knots_all = t;
self._nparameter = self._knots.shape[0] + k + 1;
# set parameter
self._parameter = c[:self._nparameter];
#set points
if points is not None:
self.set_points(points);
# values and projection matrix will need to be updated after change of the knots
self._values = None;
self._projection = None;
self._projection_inverse = None;
############################################################################
### Spline getter
def tck(self):
"""Returns a tck tuple for use with spline functions
Note:
This returns a tck tuple as splrep
"""
k = self.degree;
p = np.pad(self._parameter, (0,k+1), 'constant');
return (self._knots_all, p, k);
def convex_hull_removal(pixel, wvl):
"""
Remove the convex-hull of the signal by hull quotient.
Parameters:
pixel: `list`
1D HSI data (p), a pixel.
wvl: `list`
Wavelength of each band (p x 1).
Results: `list`
Data with convex hull removed (p).
Reference:
Clark, R.N. and T.L. Roush (1984) Reflectance Spectroscopy: Quantitative
Analysis Techniques for Remote Sensing Applications, J. Geophys. Res., 89,
6329-6340.
"""
points = list(zip(wvl, pixel))
# close the polygone
poly = [(points[0][0],0)]+points+[(points[-1][0],0)]
hull = _jarvis.convex_hull(poly)
# the last two points are on the x axis, remove it
hull = hull[:-2]
x_hull = [u for u,v in hull]
y_hull = [v for u,v in hull]
tck = interpolate.splrep(x_hull, y_hull, k=1)
iy_hull = interpolate.splev(wvl, tck, der=0)
norm = []
for ysig, yhull in zip(pixel, iy_hull):
if yhull != 0:
norm.append(ysig/yhull)
else:
norm.append(1)
return norm, x_hull, y_hull
def smooth_plot(figure,X,Y,color1,color2,xlabel="",ylabel="",filled=False,n_points=400,smooth_factor=1.0,spline_order=3,linewidth=3,alpha=1.0,label=""):
"""
"""
X_smooth = np.linspace(X.min(),X.max(),n_points)
tck = splrep(X,Y,s=smooth_factor,k=spline_order)
Y_smooth = splev(X_smooth,tck,der=0)
font = fm.FontProperties(family = 'Trebuchet', weight ='light')
#font = fm.FontProperties(family = 'CenturyGothic',fname = '/Library/Fonts/Microsoft/Century Gothic', weight ='light')
figure.patch.set_facecolor('white')
axes = figure.add_subplot(111)
axes.plot(X,Y,linewidth=1,color=tuple(color2),alpha=0.2)
if filled:
axes.fill_between(X_smooth,Y_smooth,0,color=color2,alpha=0.1)
for i in xrange(100):
color = tuple(color1*(1.0-i/100.0) + color2*(i/100.0))
if i == 0:
axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha,label=label)
else:
axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha)
axes.set_xlim(X.min(),X.max())
axes.set_xlabel(xlabel,fontproperties=font, size=10, style='italic')
axes.set_xticklabels(axes.get_xticks(),fontproperties=font, size=12)
if '%' in ylabel:
axes.set_ylim(0,np.minimum(2*Y.max(),100))
axes.set_ylabel(ylabel, fontproperties=font, size=10, style='italic')
axes.set_yticklabels(axes.get_yticks(),fontproperties=font, size=12)
def spline(self):
"""
I return a new pycs.gen.spl.Spline object corresponding to the source.
So this is a bit the inverse of the constructor.
You can then put this spline object as ML of a lightcurve, as source spline, or whatever.
Note that my output spline has LOTs of knots... it is an interpolating spline, not a regression spline !
..note:: This spline is a priori for display purposes only. To do an interpolation, it might be safer (and faster) to use
the above linear interpolation eval() function.
But making a plot, you'll see that this spline seems well accurate.
"""
x = self.ijds.copy()
y = self.imags.copy()
magerrs = np.zeros(len(x))
out = si.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=0.0, t=None, full_output=1, per=0, quiet=1)
# s = 0.0 means interpolating spline !
tck = out[0]
# From this we want to build a real Spline object.
datapoints = pycs.gen.spl.DataPoints(x, y, magerrs, splitup=False, sort=False, stab=False)
outspline = pycs.gen.spl.Spline(datapoints, t = tck[0], c = tck[1], k = tck[2], plotcolour=self.plotcolour)
# Conditions for the knots (no, everything is ok with them, no need to tweak anything).
#intt = outspline.getintt()
#outspline.setintt(intt)
#print tck[2]
#sys.exit()
outspline.knottype = "MadeBySource"
outspline.showknots = False # Usually we have a lot of them, slows down.
return outspline
def getintt(self):
"""
Returns the internal knots (i.e., not even the datapoints extrema)
This is what you need to feed into splrep !
There are nint - 1 such knots
"""
return self.t[(self.k+1):-(self.k+1)].copy() # We cut the outer knots.
