Python UnivariateSpline.get_knots方法代码示例

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
def fit_and_plot_mpl_bspline(x,y):
# Fit B-spline with matplotlib and plot it
    SPLINE_ORDER = 3
    spl = UnivariateSpline(x,y,k=SPLINE_ORDER)
    plt.plot(x,y,'r')
    plt.plot(x,spl(x),'g')

    knots  = spl.get_knots()
    coeffs = spl.get_coeffs()
    yr = plt.ylim()

    for knot in knots:
        plt.axvline(knot,color='g')
    #for coeff in coeffs:
        #plt.text(???,(yr[0]+yr[1])/2,"%f" % coeff)
    plt.show()

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
	def spline(self, k=3, s=0.5):
		""" Interpolates uneven observation data into daily data using scipy.interpolate.UnivariateSpline with parameters k and s given"""
		# requires filtered_data to exist, check by making sure it isn't empty
		if len(self.filtered_data) == 0:
			raise ValueError

		# first we need to redo the dates as days from start
		for patient in self.filtered_data:
			self.filtered_data[patient]['days'] = []
			for date in self.filtered_data[patient]['dates']:
				since_beginning = (date - self.filtered_data[patient]['dates'][0]).days
				self.filtered_data[patient]['days'].append(since_beginning)

		# now we spline
		error = 0
		splined_data = dict()
		for patient in self.filtered_data:
			day_indices = np.array(self.filtered_data[patient]['days'])
			splined_data[patient] = []
			for question in range(self.num_items):
				item_data = np.array([i[question] for i in self.filtered_data[patient]['data']])
				if len(item_data) == 0:
					print "Uh, no responses for question %d patient %s" % (question+1, patient)
					raise ValueError
				# Now we need to take out the data with None values... not certain how to handle it if the patient ends up not having enough data
				# because of this filtering, but for UPittSSRI data (primary focus as of 7/29/13) we don't need to worry about it
				good_indices = np.array([ind for ind, resp in enumerate(item_data) if resp != None])
				filtered_days = day_indices[good_indices]
				filtered_data = item_data[good_indices]
				full_space = np.linspace(0,self.keep_days,num=self.keep_days)
				spl = UnivariateSpline(filtered_days, filtered_data, k=k, s=s)
				self.residual = spl.get_residual()
				self.knots = spl.get_knots()
				splined = spl(full_space)
				splined_data[patient].append(splined)
				# Measure error
				ms = math.sqrt(sum((splined[filtered_days] - filtered_data)**2)/len(filtered_data))
				error += ms
			splined_data[patient] = np.array(splined_data[patient]).T
		error /= (self.num_items)*(len(splined_data))
		self.spline_err = error

		self.splined_data = splined_data
		self.data_splined = True
		self.data = splined_data
		return splined_data

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
def break_spline(x, y, **kwargs):
    '''
    Calculate the break in 2 linear trends using a spline.
    '''

    s = UnivariateSpline(x, y, **kwargs)
    knots = s.get_knots()

    if len(knots) > 3:
        print "Many knots"
        knot_expec = -0.8
        posn = np.argmin(knots - knot_expec)
        return knots[posn]  # Return the knot closest to the expected value

    elif len(knots) == 2:
        print "No knots"
        return -0.6  # Set the threshold at the expected

    else:
        return knots[1]

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
def fit_trace(trace, deg=5):
    spl = UnivariateSpline(trace.sample_f, trace.trace_f, k=deg)
    return [deg, spl.get_coeffs().tolist(), spl.get_knots().tolist()]

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
class DualSplineSmoother(object):
	'''
	claselfdocs
	'''


	def __init__(self, yp, workdir, scale, sm=200):
		'''
		Constructor
		'''
		yp = np.array(yp)
		self.l = len(yp)/2
		self.xPos = (self.l-2)/2 #fPos = (self.l-2)/2 + 2
		tnsc = 2/scale
		print tnsc
		plt.rcParams['font.size'] = 24
		plt.rcParams['lines.linewidth'] = 2.4
		self.workdir = workdir
		
		avProfilePoints = yp[:self.l]
		self.avx = np.append(np.append([0], np.sort(np.tanh(tnsc*avProfilePoints[:self.xPos]))),[1])
		self.av = avProfilePoints[self.xPos:]
		
		sigmaProfilePoints = yp[self.l:]
		self.sigmax = np.append(np.append([0], np.sort(np.tanh(tnsc*sigmaProfilePoints[:self.xPos]))),[1])
		self.sigma = sigmaProfilePoints[self.xPos:]
		
		self.m = UnivariateSpline(self.avx, self.av)
		print "Created spline with " + str(len(self.m.get_knots())) + " knots"

		self.s = UnivariateSpline(self.sigmax, self.sigma)
		print "Created spline with " + str(len(self.s.get_knots())) + " knots"

	def saveSpline(self, filename):
		tp = np.linspace(0, 1, 1000)
		with open(filename ,"w+") as f:
			for i in range(0, 1000):
				f.write( str(tp[i]) + " , " + str(self.m(tp[i])) )
				if i < 999:
					f.write("\n")
		f.close()
	
	def saveSigmaSpline(self, filename):
		tp = np.linspace(0, 1, 1000)
		with open(filename ,"w+") as f:
			for i in range(0, 1000):
				f.write( str(tp[i]) + " , " + str(self.s(tp[i])) )
				if i < 999:
					f.write("\n")
		f.close()
	
	def showSpline(self, order=0):
		plt.clf()
		print "Spline full information:"
		print self.m.get_knots()
		print self.m.get_coeffs()
		print self.m.get_residual()
		tp = np.linspace(0, 1, 1000)
		plt.subplot(211)
		plt.scatter(self.avx,self.av)
		plt.plot(tp,self.m(tp))
		plt.subplot(212)
		plt.scatter(self.sigmax,self.sigma)
		plt.plot(tp,self.s(tp))
		plt.savefig(self.workdir+"/splineFit.pdf")
		if order > 0:
			plt.clf()
			plt.subplot(211)
			plt.plot(tp,self.m(tp,1))
			plt.subplot(212)
			plt.plot(tp,self.s(tp,1))
			plt.savefig(self.workdir+"/splineDerivative.pdf")
	
	def plotSplineData(self, dataContainer, yscale):
		plt.clf()
		plt.xlim(0,1)
		plt.ylim(0,yscale)
		tp = np.linspace(0, 1, 100)
		plt.scatter(dataContainer.points[0],dataContainer.points[1]+dataContainer.background, c='b', marker='o', s=5)
		plt.plot(tp, self.m(tp)+dataContainer.background,'r', linewidth=2)
		plt.plot(tp, self.m(tp)+np.sqrt(self.s(tp))+dataContainer.background,'r--', linewidth=2)
		plt.plot(tp, self.m(tp)-np.sqrt(self.s(tp))+dataContainer.background,'r--', linewidth=2)
		plt.plot(tp, np.zeros(100) + dataContainer.background, '--', c='#BBBBBB', alpha=0.8)
		plt.savefig(self.workdir+"/splineVsData.pdf")
		
	def plotBinnedData(self, dataContainer):
		plt.clf()
		tp = np.linspace(0, 1, dataContainer.numBins)
		tpHD = np.linspace(0, 1, 500)
		plt.plot(self.m(tpHD), self.s(tpHD))
		plt.plot(dataContainer.avs, np.power(dataContainer.stds,2), 'o')
		plt.savefig(self.workdir+"/noiseVsBins.pdf")
		plt.clf()
		plt.plot(tpHD, self.m(tpHD),'r', linewidth=2)
		plt.plot(tp, dataContainer.avs, 'o')
		plt.savefig(self.workdir+"/splineVsBins.pdf")
		plt.clf()
		plt.plot(tpHD, self.s(tpHD),'r', linewidth=2)
		plt.plot(tp, np.power(dataContainer.stds,2), 'o')
		plt.savefig(self.workdir+"/spatialNoiseVsBins.pdf")
#.........这里部分代码省略.........

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
    def fit_pspec(self, breaks=None, log_break=True, low_cut=None,
                  high_cut=None, fit_verbose=False, bootstrap=False,
                  **bootstrap_kwargs):
        '''
        Fit the 1D Power spectrum using a segmented linear model. Note that
        the current implementation allows for only 1 break point in the
        model. If the break point is estimated via a spline, the breaks are
        tested, starting from the largest, until the model finds a good fit.

        Parameters
        ----------
        breaks : float or None, optional
            Guesses for the break points. If given as a list, the length of
            the list sets the number of break points to be fit. If a choice is
            outside of the allowed range from the data, Lm_Seg will raise an
            error. If None, a spline is used to estimate the breaks.
        log_break : bool, optional
            Sets whether the provided break estimates are log-ed values.
        low_cut : `~astropy.units.Quantity`, optional
            Lowest frequency to consider in the fit.
        high_cut : `~astropy.units.Quantity`, optional
            Highest frequency to consider in the fit.
        fit_verbose : bool, optional
            Enables verbose mode in Lm_Seg.
        bootstrap : bool, optional
            Bootstrap using the model residuals to estimate the standard
            errors.
        bootstrap_kwargs : dict, optional
            Pass keyword arguments to `~turbustat.statistics.fitting_utils.residual_bootstrap`.
        '''

        # Make the data to fit to
        if low_cut is None:
            # Default to the largest frequency, since this is just 1 pixel
            # But cut out fitting the total power point.
            self.low_cut = 0.98 / (float(self.data.shape[0]) * u.pix)
        else:
            self.low_cut = \
                self._spectral_freq_unit_conversion(low_cut, u.pix**-1)

        if high_cut is None:
            # Set to something larger than
            self.high_cut = (self.freqs.max().value * 1.01) / u.pix
        else:
            self.high_cut = \
                self._spectral_freq_unit_conversion(high_cut, u.pix**-1)

        # We keep both the real and imag frequencies, but we only need the real
        # component when fitting. We'll cut off the total power frequency too
        # in case it causes an extra break point (which seems to happen).
        shape = self.freqs.size
        rfreqs = self.freqs[1:shape // 2].value
        ps1D = self.ps1D[1:shape // 2]

        y = np.log10(ps1D[clip_func(rfreqs, self.low_cut.value,
                                    self.high_cut.value)])
        x = np.log10(rfreqs[clip_func(rfreqs, self.low_cut.value,
                                      self.high_cut.value)])

        if breaks is None:
            from scipy.interpolate import UnivariateSpline

            # Need to order the points
            spline = UnivariateSpline(x, y, k=1, s=1)

            # The first and last are just the min and max of x
            breaks = spline.get_knots()[1:-1]

            if fit_verbose:
                print("Breaks found from spline are: " + str(breaks))

            # Take the number according to max_breaks starting at the
            # largest x.
            breaks = breaks[::-1]

            # Ensure a break doesn't fall at the max or min.
            if breaks.size > 0:
                if breaks[0] == x.max():
                    breaks = breaks[1:]
            if breaks.size > 0:
                if breaks[-1] == x.min():
                    breaks = breaks[:-1]

            if x.size <= 3 or y.size <= 3:
                raise Warning("There are no points to fit to. Try lowering "
                              "'lg_scale_cut'.")

            # If no breaks, set to half-way before last point
            if breaks.size == 0:
                x_copy = np.sort(x.copy())
                breaks = np.array([0.5 * (x_copy[-1] + x_copy[-3])])

            # Now try these breaks until a good fit including the break is
            # found. If none are found, it accepts that there wasn't a good
            # break and continues on.
            i = 0
            while True:
                self.fit = \
                    Lm_Seg(x, y, breaks[i])
                self.fit.fit_model(verbose=fit_verbose, cov_type='HC3')
#.........这里部分代码省略.........

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
class SplineSmoother(object):
	'''
	claselfdocs
	'''


	def __init__(self, yp, workdir, scale, sm=200):
		'''
		Constructor
		'''
		
		self.workdir = workdir
		yp = np.array(yp)
		self.l = len(yp)
		self.xPos = (self.l-2)/2 #fPos = (self.l-2)/2 + 2
		tnsc = 2/scale
		print tnsc
		plt.rcParams['font.size'] = 24
		plt.rcParams['lines.linewidth'] = 2.4
		self.workdir = workdir
		
		self.avx = np.append(np.append([0], np.sort(np.tanh(tnsc*yp[:self.xPos]))),[1])
		self.av = yp[self.xPos:]
		self.m = UnivariateSpline(self.avx, self.av)
		
		plt.rcParams['font.size'] = 24
		plt.rcParams['lines.linewidth'] = 2.4

		print "Created spline with " + str(len(self.m.get_knots())) + " knots"

	def saveSpline(self, filename):
		tp = np.linspace(0, 1, 1000)
		with open(filename ,"w+") as f:
			for i in range(0, 1000):
				f.write( str(tp[i]) + " , " + str(self.m(tp[i])) )
				if i < 999:
					f.write("\n")
		f.close()
	
	def showSpline(self, order=0):
		plt.clf()
		print "Spline full information:"
		print self.m.get_knots()
		print self.m.get_coeffs()
		print self.m.get_residual()
		tp = np.linspace(0, 1, 1000)
		plt.scatter(self.avx, self.av)
		plt.plot(tp,self.m(tp))
		plt.savefig(self.workdir+"/splineFit.pdf")
		if order > 0:
			plt.plot(tp,self.m(tp,1))
			if order > 1:
				plt.plot(tp,self.m(tp,2))
			plt.savefig(self.workdir+"/splineDerivative.pdf")
	
	def plotSplineData(self, dataContainer, s, p, se, yscale):
		plt.clf()
		plt.xlim(0,1)
		plt.ylim(0,yscale)
		tp = np.linspace(0, 1, 500)
		plt.scatter(dataContainer.points[0],dataContainer.points[1]+dataContainer.background, c='b', marker='o', s=5)
		plt.plot(tp, self.m(tp)+dataContainer.background,'r', linewidth=2)
		plt.plot(tp, self.m(tp)+np.sqrt(s*(p*self.m(tp)*self.m(tp)+self.m(tp)+se))+dataContainer.background,'r--', linewidth=2)
		plt.plot(tp, self.m(tp)-np.sqrt(s*(p*self.m(tp)*self.m(tp)+self.m(tp)+se))+dataContainer.background,'r--', linewidth=2)
		plt.plot(tp, np.zeros(500) + dataContainer.background, '--', c='#BBBBBB', alpha=0.8)
		plt.savefig(self.workdir+"/splineVsData.pdf")
		
	def plotBinnedData(self, dataContainer, s, p, se, xmin, xmax):
		plt.clf()
		sigma = lambda x: s * (p*x*x + x + se)
		t = np.linspace(xmin, xmax, 500)
		plt.xlim(xmin, xmax)
		plt.plot(t, sigma(t))
		plt.plot(dataContainer.avs, np.power(dataContainer.stds,2), 'o')
		plt.savefig(self.workdir+"/noiseVsBins.pdf")
		plt.clf()
		tp = np.linspace(0, 1, dataContainer.numBins)
		plt.plot(tp, self.m(tp),'r', linewidth=2)
		plt.plot(tp, dataContainer.avs, 'o')
		plt.savefig(self.workdir+"/splineVsBins.pdf")
	
	def plotFisherInfo(self, dataContainer, s, p, se, ymax, ymaxsq):
		plt.clf()
		t = np.linspace(0, 1, 1000)
		
		minf = lambda x: -1 * self.m(x)
		minx = fminbound(minf, 0, 1)
		fval = self.m(minx)
		s = s/fval
		se = se/fval
		p = p*fval
		print fval, s, p, se
		fi = lambda m, mp: 2*np.power(s * (p * np.power(m,2) + m),2)/(np.power(mp,2) * 
		(2 * s * (p * np.power(m,2) + m) + np.power(s * (2 * p * m + 1), 2)))
		fiapp = lambda m, mp: s * (p * np.power(m,2) + m) / np.power(mp,2)
		plt.xlim(0, 1)
		plt.ylim(0, ymaxsq)
		plt.plot(t, fi(self.m(t)/fval, self.m(t, 1)/fval))
		plt.plot(t, fiapp(self.m(t)/fval, self.m(t, 1)/fval), 'r')
		plt.savefig(self.workdir+"/variance.pdf")
#.........这里部分代码省略.........

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
    def fit_pspec(self, breaks=None, log_break=True, lg_scale_cut=2, verbose=False):
        """
        Fit the 1D Power spectrum using a segmented linear model. Note that
        the current implementation allows for only 1 break point in the
        model. If the break point is estimated via a spline, the breaks are
        tested, starting from the largest, until the model finds a good fit.

        Parameters
        ----------
        breaks : float or None, optional
            Guesses for the break points. If given as a list, the length of
            the list sets the number of break points to be fit. If a choice is
            outside of the allowed range from the data, Lm_Seg will raise an
            error. If None, a spline is used to estimate the breaks.
        log_break : bool, optional
            Sets whether the provided break estimates are log-ed values.
        lg_scale_cut : int, optional
            Cuts off largest scales, which deviate from the powerlaw.
        verbose : bool, optional
            Enables verbose mode in Lm_Seg.
        """

        if breaks is None:
            from scipy.interpolate import UnivariateSpline

            # Need to order the points
            shape = self.vel_freqs.size
            spline_y = np.log10(self.ps1D[1 : shape / 2])
            spline_x = np.log10(self.vel_freqs[1 : shape / 2])

            spline = UnivariateSpline(spline_x, spline_y, k=1, s=1)

            # The first and last are just the min and max of x
            breaks = spline.get_knots()[1:-1]

            if verbose:
                print "Breaks found from spline are: " + str(breaks)

            # Take the number according to max_breaks starting at the
            # largest x.
            breaks = breaks[::-1]

            x = np.log10(self.vel_freqs[lg_scale_cut + 1 : -lg_scale_cut])
            y = np.log10(self.ps1D[lg_scale_cut + 1 : -lg_scale_cut])

            # Now try these breaks until a good fit including the break is
            # found. If none are found, it accept that there wasn't a good
            # break and continues on.
            i = 0
            while True:
                self.fit = Lm_Seg(x, y, breaks[i])
                self.fit.fit_model(verbose=verbose)

                if self.fit.params.size == 5:
                    break
                i += 1
                if i >= breaks.shape:
                    warnings.warn(
                        "No good break point found. Returned fit\
                                   does not include a break!"
                    )
                    break

            return self

        if not log_break:
            breaks = np.log10(breaks)

        self.fit = Lm_Seg(
            np.log10(self.vel_freqs[lg_scale_cut + 1 : -lg_scale_cut]),
            np.log10(self.ps1D[lg_scale_cut + 1 : -lg_scale_cut]),
            breaks,
        )
        self.fit.fit_model(verbose=verbose)

        return self

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import get_knots [as 别名]
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
nItems = 500
x = np.linspace(-3, 3, nItems)
y = np.exp(-x**2) + 0.1 * np.random.randn(nItems)
plt.plot(x, y, 'ro', ms=5)

smoothFactor = float(nItems) / 50.0

spl = UnivariateSpline(x, y, k = 2, s = smoothFactor)
xs = np.linspace(-3, 3, 1000)
plt.plot(xs, spl(xs), 'g', lw=3)

coeffs = spl.get_coeffs()
knots  = spl.get_knots()

print 'coeffs', coeffs
print 'knots', knots

knotsY = spl(knots)
plt.plot(knots, knotsY, 'bo')

#spl.set_smoothing_factor(0.5)
#print 'coeffs', spl.get_coeffs()
#print 'knots', spl.get_knots()
#plt.plot(xs, spl(xs), 'b', lw=3)
plt.show()