def r2(self, nostab=True, nosquare=False):
"""
Evaluates the spline, compares it with the data points and returns a weighted sum of residuals r2.
If nostab = False, stab points are included
This is precisely the same r2 as is used by splrep for the fit, and thus the same value as
returned by optc !
This method can set lastr2nostab, so be sure to end any optimization with it.
If nostab = True, we don't count the stab points
"""
if nostab == True :
splinemags = self.eval(nostab = True, jds = None)
errs = self.datapoints.mags[self.datapoints.mask] - splinemags
werrs = errs/self.datapoints.magerrs[self.datapoints.mask]
if nosquare:
r2 = np.sum(np.fabs(werrs))
else:
r2 = np.sum(werrs * werrs)
self.lastr2nostab = r2
else :
splinemags = self.eval(nostab = False, jds = None)
errs = self.datapoints.mags - splinemags
werrs = errs/self.datapoints.magerrs
if nosquare:
r2 = np.sum(np.fabs(werrs))
else:
r2 = np.sum(werrs * werrs)
self.lastr2stab = r2
return r2
#if red:
# return chi2/len(self.datapoints.jds)
python类splrep()的实例源码
def _interpolate(self, xa, xb, data):
from scipy import interpolate
f = interpolate.splrep(xa, data)
return interpolate.splev(xb, f)
def turning_points(mu, vv, Gv, Bv, dv=0.1):
DD = np.sqrt(const.h/(8*mu*const.m_u*const.c*100))*1.0e10/np.pi
# Gv spline
gsp = splrep(vv, Gv, s=0)
# Bv spline
bsp = splrep(vv, Bv, s=0)
# vibrational QN at which to evaluate turning points
V = np.arange(dv-1/2, vv[-1], dv)
# add a point close to v=-0.5, the bottom of the well
V = np.append([-1/2+0.0001], V)
Rmin = []; Rmax = []; E = []
# compute turning points using RKR method
print(u"RKR: v Rmin(A) Rmax(A) E(cm-1)")
for vib in V:
E.append(G(vib, gsp)) # energy of vibrational level
ff = fg_integral(vib, gsp, bsp, lambda x, y: 1)
gg = fg_integral(vib, gsp, bsp, B)
fg = np.sqrt(ff/gg + ff**2)
Rmin.append((fg - ff)*DD) # turning points
Rmax.append((fg + ff)*DD)
frac, integ = np.modf(vib)
if frac > 0 and frac < dv:
print(u" {:d} {:6.4f} {:6.4f} {:6.2f}".
format(int(vib), Rmin[-1], Rmax[-1], np.float(E[-1])))
return Rmin, Rmax, E
def formPEC(R, Rmin, Rmax, E, De, limb):
evcm = const.e/(const.c*const.h*100) # converts cm-1 to eV
# combine Rmin with Rmax to form PEC
Re = (Rmin[0] + Rmax[0])/2
print(u"RKR: Re = {:g}".format(Re))
RTP = np.append(Rmin[::-1], Rmax, 0) # radial positions of turning-points
PTP = np.append(E[::-1], E, 0) # potential energy at turning point
# interpolate
psp = splrep(RTP, PTP, s=0)
# Interpolate RKR curve to this grid
PEC = splev(R, psp, der=0)
# extrapolate using analytical function
inner_limb_Morse(R, PEC, RTP, PTP, Re, De)
if limb=='L':
outer_limb_LeRoy(R, PEC, RTP, PTP, De)
else:
outer_limb_Morse(R, PEC, RTP, PTP, De)
PTP /= evcm
PEC /= evcm # convert to eV
return PEC, RTP, PTP
# analytical functions
def p2eRatio(self):
"""
Calculates the proton density to electron density ratio using Eq. 7 of [1]_.
Notes
------
Uses the abundance and ionization equilibrium.
References
----------
.. [1] Young, P. R. et al., 2003, ApJS, `144, 135 <http://adsabs.harvard.edu/abs/2003ApJS..144..135Y>`_
"""
if hasattr(self, 'Temperature'):
temperature = self.Temperature
else:
temperature = self.IoneqAll['ioneqTemperature']
if not hasattr(self, 'AbundanceName'):
AbundanceName = self.Defaults['abundfile']
else:
AbundanceName = self.AbundanceName
tmp_abundance = io.abundanceRead(abundancename=AbundanceName)
abundance = tmp_abundance['abundance'][tmp_abundance['abundance']>0]
denominator = np.zeros(len(self.IoneqAll['ioneqTemperature']))
for i in range(len(abundance)):
for z in range(1,i+2):
denominator += z*self.IoneqAll['ioneqAll'][i,z,:]*abundance[i]
p2eratio = abundance[0]*self.IoneqAll['ioneqAll'][0,1,:]/denominator
nots = interpolate.splrep(np.log10(self.IoneqAll['ioneqTemperature']),p2eratio,s=0)
self.ProtonDensityRatio = interpolate.splev(np.log10(temperature),nots,der=0,ext=1)
def ioneqOne(self):
'''
Provide the ionization equilibrium for the selected ion as a function of temperature.
returned in self.IoneqOne
'''
#
if hasattr(self, 'Temperature'):
temperature = self.Temperature
else:
return
#
if hasattr(self, 'IoneqAll'):
ioneqAll = self.IoneqAll
else:
ioneqAll = io.ioneqRead(ioneqname = self.Defaults['ioneqfile'])
self.ioneqAll = self.IoneqAll
#
ioneqTemperature = ioneqAll['ioneqTemperature']
Z = self.Z
Ion = self.Ion
Dielectronic = self.Dielectronic
ioneqOne = np.zeros_like(temperature)
#
thisIoneq = ioneqAll['ioneqAll'][Z-1,Ion-1 + Dielectronic].squeeze()
gioneq = thisIoneq > 0.
goodt1 = self.Temperature >= ioneqTemperature[gioneq].min()
goodt2 = self.Temperature <= ioneqTemperature[gioneq].max()
goodt = np.logical_and(goodt1,goodt2)
y2 = interpolate.splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0)
#
if goodt.sum() > 0:
if self.Temperature.size > 1:
gIoneq = interpolate.splev(np.log(self.Temperature[goodt]),y2) #,der=0)
ioneqOne[goodt] = np.exp(gIoneq)
else:
gIoneq = interpolate.splev(np.log(self.Temperature),y2)
ioneqOne = np.exp(gIoneq)*np.ones(self.NTempDen, 'float64')
self.IoneqOne = ioneqOne
else:
self.IoneqOne = np.zeros_like(self.Temperature)
def calculate_free_free_loss(self, **kwargs):
"""
Calculate the free-free energy loss rate of an ion. The result is returned to the
`free_free_loss` attribute.
The free-free radiative loss rate is given by Eq. 5.15a of [1]_. Writing the numerical
constant in terms of the fine structure constant :math:`\\alpha`,
.. math::
\\frac{dW}{dtdV} = \\frac{4\\alpha^3h^2}{3\pi^2m_e}\left(\\frac{2\pi k_B}{3m_e}\\right)^{1/2}Z^2T^{1/2}\\bar{g}_B
where where :math:`Z` is the nuclear charge, :math:`T` is the electron temperature, and
:math:`\\bar{g}_{B}` is the wavelength-averaged and velocity-averaged Gaunt factor. The
Gaunt factor is calculated using the methods of [2]_. Note that this expression for the
loss rate is just the integral over wavelength of Eq. 5.14a of [1]_, the free-free emission, and
is expressed in units of erg :math:`\mathrm{cm}^3\,\mathrm{s}^{-1}`.
References
----------
.. [1] Rybicki and Lightman, 1979, Radiative Processes in Astrophysics,
`(Wiley-VCH) <http://adsabs.harvard.edu/abs/1986rpa..book.....R>`_
.. [2] Karzas and Latter, 1961, ApJSS, `6, 167
<http://adsabs.harvard.edu/abs/1961ApJS....6..167K>`_
"""
# interpolate wavelength-averaged K&L gaunt factors
gf_kl_info = ch_io.gffintRead()
gamma_squared = self.ionization_potential/ch_const.boltzmann/self.Temperature
for i, atemp in enumerate(self.Temperature):
print('%s T: %10.2e gamma_squared %10.2e'%(self.ion_string, atemp, gamma_squared[i]))
gaunt_factor = splev(np.log(gamma_squared),
splrep(gf_kl_info['g2'],gf_kl_info['gffint']), ext=3)
# calculate numerical constant
prefactor = (4.*(ch_const.fine**3)*(ch_const.planck**2)/3./(np.pi**2)/ch_const.emass
* np.sqrt(2.*np.pi*ch_const.boltzmann/3./ch_const.emass))
# include abundance and ionization equilibrium
prefactor *= self.abundance*self.ioneq_one(self.stage, **kwargs)
self.free_free_loss = prefactor*(self.Z**2)*np.sqrt(self.Temperature)*gaunt_factor
def freeFreeLoss(self, **kwargs):
"""
Calculate the free-free energy loss rate of an ion. The result is returned to the
`free_free_loss` attribute.
The free-free radiative loss rate is given by Eq. 5.15a of [1]_. Writing the numerical
constant in terms of the fine structure constant :math:`\\alpha`,
.. math::
\\frac{dW}{dtdV} = \\frac{4\\alpha^3h^2}{3\pi^2m_e}\left(\\frac{2\pi k_B}{3m_e}\\right)^{1/2}Z^2T^{1/2}\\bar{g}_B
where where :math:`Z` is the nuclear charge, :math:`T` is the electron temperature, and
:math:`\\bar{g}_{B}` is the wavelength-averaged and velocity-averaged Gaunt factor. The
Gaunt factor is calculated using the methods of [2]_. Note that this expression for the
loss rate is just the integral over wavelength of Eq. 5.14a of [1]_, the free-free emission, and
is expressed in units of erg :math:`\mathrm{cm}^3\,\mathrm{s}^{-1}`.
References
----------
.. [1] Rybicki and Lightman, 1979, Radiative Processes in Astrophysics,
`(Wiley-VCH) <http://adsabs.harvard.edu/abs/1986rpa..book.....R>`_
.. [2] Karzas and Latter, 1961, ApJSS, `6, 167
<http://adsabs.harvard.edu/abs/1961ApJS....6..167K>`_
"""
# interpolate wavelength-averaged K&L gaunt factors
gf_kl_info = ch_io.gffintRead()
gamma_squared = self.IprErg/ch_const.boltzmann/self.Temperature
# for i, atemp in enumerate(self.Temperature):
# print('%s T: %10.2e gamma_squared %10.2e'%(self.ion_string, atemp, gamma_squared[i]))
gaunt_factor = splev(np.log(gamma_squared),
splrep(gf_kl_info['g2'],gf_kl_info['gffint']), ext=3)
# calculate numerical constant
prefactor = (4.*(ch_const.fine**3)*(ch_const.planck**2)/3./(np.pi**2)/ch_const.emass
* np.sqrt(2.*np.pi*ch_const.boltzmann/3./ch_const.emass))
# include abundance and ionization equilibrium
prefactor *= self.abundance*self.IoneqOne
self.FreeFreeLoss = {'rate':prefactor*(self.Z**2)*np.sqrt(self.Temperature)*gaunt_factor}
def ioneqOne(self):
'''
Provide the ionization equilibrium for the selected ion as a function of temperature.
Similar to but not identical to ion.ioneqOne()
returned in self.IoneqOne
'''
#
if hasattr(self, 'Temperature'):
temperature = self.Temperature
else:
return
#
if hasattr(self, 'IoneqAll'):
ioneqAll = self.IoneqAll
else:
self.IoneqAll = ch_io.ioneqRead(ioneqname = self.Defaults['ioneqfile'])
ioneqAll = self.IoneqAll
#
ioneqTemperature = ioneqAll['ioneqTemperature']
Z = self.Z
stage = self.stage
ioneqOne = np.zeros_like(temperature)
#
thisIoneq = ioneqAll['ioneqAll'][Z-1,stage-1].squeeze()
gioneq = thisIoneq > 0.
goodt1 = self.Temperature >= ioneqTemperature[gioneq].min()
goodt2 = self.Temperature <= ioneqTemperature[gioneq].max()
goodt = np.logical_and(goodt1,goodt2)
y2 = splrep(np.log(ioneqTemperature[gioneq]),np.log(thisIoneq[gioneq]),s=0)
#
if goodt.sum() > 0:
if self.Temperature.size > 1:
gIoneq = splev(np.log(self.Temperature[goodt]),y2) #,der=0)
ioneqOne[goodt] = np.exp(gIoneq)
else:
gIoneq = splev(np.log(self.Temperature),y2)
ioneqOne = np.exp(gIoneq)
ioneqOne = np.atleast_1d(ioneqOne)
self.IoneqOne = ioneqOne
else:
self.IoneqOne = np.zeros_like(self.Temperature)
def __init__(self, t, x, y, omega, origin):
self.t = t
self.x = x
self.y = y
self.omega = omega
self.origin = origin
self.n = len(t)
self.x_spl = interpolate.splrep(self.t, self.x)
self.y_spl = interpolate.splrep(self.t, self.y)
self.omega_spl = interpolate.splrep(self.t, self.omega)
preview_dataset.py 文件源码
项目：ego-lane-analysis-system
作者: rodrigoberriel
项目源码
文件源码
阅读 18
收藏 0
点赞 0
评论 0
def draw_spline(image, lane, ys, color=[0, 0, 255]):
pts = [[x, ys[i]]for i, x in enumerate(lane) if not math.isnan(x)]
if len(pts)-1 > 0:
spline = interpolate.splrep([pt[1] for pt in pts], [pt[0] for pt in pts], k=len(pts)-1)
for i in range(min([pt[1] for pt in pts]), max([pt[1] for pt in pts])):
image[i, int(interpolate.splev(i, spline)), :] = color
return image
def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3):
self.A = A
self.h = h
self.Asp = 2*self.h**2/self.A
self.rho = rho
# Input lift and drag data
self.n = len(alpha)
self.alpha = alpha
self.CL = CL
self.CD = CD
self.k_spline = k_spline
# -------- Spline interpolation -----------------------------
if len(self.alpha.shape) == 1:
self.interpolationMethod = 'spline'
self.nrControlParameters = 1
self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline)
self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline)
elif len(self.alpha.shape) == 2:
self.interpolationMethod = 'griddata'
self.nrControlParameters = 2
self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing)
self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)
def setCalmWaterResistanceData(self, U, Fr, Re, CR, Cf):
self.Fr = Fr
self.U = U
self.Re = Re
self.CR = CR
self.Cf = Cf
self.CRspl = interpolate.splrep(self.Fr, self.CR, s=1e-9)
self.Cfspl = interpolate.splrep(self.Fr, self.Cf, s=1e-9)
def necessaryDriftAngle(self, U, Fy, scale=1):
alpha_test = np.linspace(0, np.max(self.alpha), 5)
Fx_test, Fy_test, Mz_test = self.driftInducedForces(U, alpha_test, scale=scale)
alpha_spline = interpolate.splrep(Fy_test, alpha_test)
return interpolate.splev(Fy, alpha_spline)
def reqRevolutions(self, T_req, U):
J_min = 0.05
n_min = U/(self.J_max*self.D)
n_max = U/(J_min*self.D)
n_test = 10
n = np.linspace(n_min, n_max, n_test)
T = self.Thrust(U, n)
n_spl = interpolate.splrep(T, n)
return float(interpolate.splev(T_req, n_spl))
def figure(self,mode,data,name,**options):
'''
Generate a figure to view the data.
Parameters
----------
mode : 'L','P'
'L' for lines and 'P' for pcolor.
data : ndarray
The data to be viewed.
name : str
The name of the figure.
options : dict
Other options.
'''
assert mode in ('L','P')
plt.title(os.path.basename(name))
if mode=='L':
if options.get('interpolate',False):
plt.plot(data[:,0],data[:,1],'r.')
X=np.linspace(data[:,0].min(),data[:,0].max(),10*data.shape[0])
for i in xrange(1,data.shape[1]):
tck=interpolate.splrep(data[:,0],data[:,i],k=3)
Y=interpolate.splev(X,tck,der=0)
plt.plot(X,Y,label=options['legend'][i-1] if 'legend' in options else None)
if 'legend' in options:
leg=plt.legend(fancybox=True,loc=options.get('legendloc',None))
leg.get_frame().set_alpha(0.5)
else:
plt.plot(data[:,0],data[:,1:])
if 'legend' in options:
leg=plt.legend(options['legend'],fancybox=True,loc=options.get('legendloc',None))
leg.get_frame().set_alpha(0.5)
elif mode=='P':
plt.colorbar(plt.pcolormesh(data[:,:,0],data[:,:,1],data[:,:,2]))
if 'axis' in options: plt.axis(options.get('axis','equal'))
if self.show and self.suspend: plt.show()
if self.show and not self.suspend: plt.pause(App.SUSPEND_TIME)
if self.savefig: plt.savefig('%s.png'%name)
plt.close()
def interpolate_to_mesh_indices(self, loop):
"""
Return interpolated loop indices to the temperature and density meshes defined for
the atomic data. For use with `~scipy.ndimage.map_coordinates`.
"""
nots_itemperature = splrep(self.temperature.value, np.arange(self.temperature.shape[0]))
nots_idensity = splrep(self.density.value, np.arange(self.density.shape[0]))
itemperature = splev(np.ravel(loop.electron_temperature.value), nots_itemperature)
idensity = splev(np.ravel(loop.density.value), nots_idensity)
return itemperature, idensity
def process(self, spectrogram):
mult = np.multiply(self.weights, np.ones((1, spectrogram.shape[1])))
square_mag = np.power(spectrogram.magnitude(),2)
square_mag = np.vstack((square_mag[1:,:],np.flipud(square_mag[1:,:])))
square_mag_sum = np.sum(square_mag, 0)
energy_threshold = np.percentile(square_mag_sum, 25)
res = np.fft.rfft(square_mag.T, 2 * (self.n_bins - 1)).T
resN = np.abs(res)
resP = np.angle(res)
yin = np.zeros((spectrogram.shape))
yin[0,:] = 1
tmp = np.zeros((spectrogram.shape[1],))
for tau in range(1, self.n_bins):
yin[tau,:] = square_mag_sum - resN[tau,:] * np.cos(resP[tau,:])
tmp += yin[tau,:]
yin[tau,:] *= tau/tmp
tau = self.tau_min + np.argmin(yin[self.tau_min:self.tau_max,:],0)
y_min = np.min(yin[self.tau_min:self.tau_max,:],0)
tau = tau[:,np.newaxis]
if self.interp:
ranges = np.hstack((tau - 5, tau - 4, tau - 3, tau-2,tau-1,tau,
tau+1,tau + 2, tau + 3, tau + 4, tau + 5))
new_tau = np.empty_like(tau)
for frame in range(spectrogram.shape[1]):
r = ranges[frame]
r[0] = max(r[0], 0)
r[1] = min(r[1], self.n_bins-1)
val = yin[r.astype(int), frame]
tck = interpolate.splrep(r, val, s=0, k = 2)
xnew = np.arange(r[0],r[-1], (r[-1] - r[0]) /10)
ynew = interpolate.splev(xnew, tck, der=0)
new_tau[frame] = xnew[np.argmin(ynew)]
y_min[frame] = np.min(ynew)
tau = new_tau
tau = tau[:,0]
pitch_confidence = 1- y_min
pitch = np.zeros((spectrogram.shape[1],))
pitch[tau!=0] = np.nan_to_num(self.sample_rate * 1.0 / (tau[tau!=0]))
pitch[tau==0] = 0
pitch_confidence[tau==0] = 0
pitch[square_mag_sum < energy_threshold] = 0
pitch_confidence[square_mag_sum < energy_threshold] = 0
return (pitch, pitch_confidence)
def interp_helper(all_data, trend_data, time_from):
'performs lf spline + hf fft interpolation of radial distance'
all_times, all_values = zip(*all_data)
trend_times, trend_values = zip(*trend_data)
split_time = int(time_to_index(time_from, all_times[0]))
trend_indices = array([time_to_index(item, all_times[0]) for item in trend_times])
spline = splrep(trend_indices, array(trend_values))
all_indices = array([time_to_index(item, all_times[0]) for item in all_times])
trend = splev(all_indices, spline)
detrended = array(all_values) - trend
trend_add = splev(arange(split_time, all_indices[-1]+1), spline)
dense_samples = detrended[:split_time]
sparse_samples = detrended[split_time:]
sparse_indices = (all_indices[split_time:]-split_time).astype(int)
amp = log(absolute(rfft(dense_samples)))
dense_freq = rfftfreq(dense_samples.size, 5)
periods = (3000.0, 300.0)
ind_from = int(round(1/(periods[0]*dense_freq[1])))
ind_to = int(round(1/(periods[1]*dense_freq[1])))
slope, _ = polyfit(log(dense_freq[ind_from:ind_to]), amp[ind_from:ind_to], 1)
params = {
't_max': periods[0],
'slope': slope,
'n_harm': 9,
'scale': [20, 4, 2*pi]
}
series_func, residual_func = make_residual_func(sparse_samples, sparse_indices, **params)
x0 = array([0.5]*(params["n_harm"]+2))
bounds = [(0, 1)]*(params["n_harm"]+2)
result = minimize(residual_func, x0, method="L-BFGS-B", bounds=bounds, options={'eps':1e-2})
interp_values = [trend + high_freq for trend, high_freq in
zip(trend_add, series_func(result.x)[:sparse_indices[-1]+1])]
#make_qc_plot(arange(sparse_indices[-1]+1), interp_values,
# sparse_indices, array(all_values[split_time:]))
interp_times = [index_to_time(ind, time_from) for ind in range(sparse_indices[-1]+1)]
return list(zip(interp_times, interp_values))
def set_values(self, values, points = None, dimension = None):
"""Calcualte the bspline parameter for the data points y
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])
dimension (int, list or None): the dimension(s) at which to change the curve, if None change dimension to values.shape[0]
"""
values = np.array(values, dtype = float);
if values.ndim == 1:
values = values[:,np.newaxis];
vdims = range(values.shape[1]);
# determine the dimensions at which to change curve
if dimension is None:
pdims = range(values.shape[1]);
self.ndim = values.shape[1];
else:
pdims = np.array(dimension, dtype = int);
if len(vdims) != len(pdims) or len(pdims) != values.shape[1] or max(pdims) > self.ndim:
raise RuntimeError('inconsistent number of dimensions %d, values %d and parameter %d and curve %d' % (values.shape[1], len(vdims), len(pdims), self.ndim));
#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[:,pdims] = self._projection_inverse.dot(values);
#self._values[:,pdims] = values;
else:
#tck,u = splprep(values, u = self.points, t = self.knots, task = -1, k = self.degree, s = 0); # splprep crashes due to erros in fitpack
#for d in range(self.ndim):
# self.parameter[d] = tck[1][d];
# self.values[d] = self.basis.dot(self.parameter[d]);
for v,p in zip(vdims, pdims):
tck = splrep(self.points, values[:,v], t = self.knots, task = -1, k = self.degree);
self.parameter[:,p] = tck[1][:self.nparameter];
# values will change
self._values = None;
#self._values = values; #fast but imprecise as values of spline due approximation might differ!
