Python UnivariateSpline.derivatives方法代码示例

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
def calc_conductance_curve(V_list,T,R_T,C_sigma):
    #test_voltages = arange(-v_max,v_max,v_step)
    test_currents = []
    for V in V_list:
        test_currents.append(calc_current(V,T,R_T,C_sigma))
        #print "V: %g, current %g"%(V,test_currents[-1])

    ## calc conductances manually
    #test_conductances = []
    #for idx,V in enumerate (test_currents[1:-2]):
    #    if idx==0:
    #        print idx
    #    test_conductances.append((test_currents[idx+2]-test_currents[idx])/(2.0*v_step))
    #
    #test_voltages_G = test_voltages[1:-2]

    #
    # SPLINE
    #
    spline = UnivariateSpline(V_list,test_currents,s=0)
    #print "test_conductances"
    #indices = [x for x, y in enumerate(col1) if (y >0.7 or y<-0.7)]
    test_conductances = []
    for v_iter in V_list:
        test_conductances.append(spline.derivatives(v_iter)[1])
    return test_conductances

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
def QuasiPWeight(ReSE_A):
	''' calculating the Fermi-liquid quasiparticle weight (residue) Z '''
	N = len(En_A)
	#M = int(1e-3/dE) if dE < 1e-3 else 1	# very fine grids lead to oscillations
	# replace 1 with M below to dilute the grid
	ReSE = UnivariateSpline(En_A[int(N/2-10):int(N/2+10):1],ReSE_A[int(N/2-10):int(N/2+10):1])
	dReSEdw = ReSE.derivatives(0.0)[1]
	Z = 1.0/(1.0-dReSEdw)
	return sp.array([Z,dReSEdw])

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
def calc_conductance_curve_full(sigma,N,V_list,R_T,C_sigma,T_p,island_volume):
    """
    Calculates full conductance curve taking into account self-heating
    """
    test_currents = []
    for V in V_list:
        test_currents.append(calc_current_full(sigma,N,V,R_T,C_sigma,T_p,island_volume))
    # SPLINE
    spline = UnivariateSpline(V_list,test_currents,s=0)
    test_conductances = []
    for v_iter in V_list:
        test_conductances.append(spline.derivatives(v_iter)[1])
    return test_conductances

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
def _find_peak(x, y, ye, k=3, s=None):
    spl = UnivariateSpline(x, y, ye ** -1, k=k, s=s)
    f = lambda k : -spl(k)
    fprime = np.vectorize(lambda k : - spl.derivatives(k)[1])
    xp_best = None
    yp_best = -np.inf
    bounds = [(x.min(), x.max())]
    for i in xrange(5):
        x0 = (x.ptp() * np.random.rand() + x.min(),)
        xp, nfeval, rc = fmin_tnc(f, x0, fprime=fprime, bounds=bounds,
                messages=tnc.MSG_NONE)
        xp = xp.item()
        yp = spl(xp)
        if yp >= yp_best:
            xp_best = xp
            yp_best = yp
    return xp_best, yp_best.item()

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
 def add_Nsquared(self, rhokey="rho", depthkey="z", N2key="N2", s=0.2):
     """ Calculate the squared buoyancy frequency, based on in-situ density.
     Uses a smoothing spline to compute derivatives.
     
     rhokey::string              Data key to use for in-situ density
     depthkey::string            Data key to use for depth
     N2key::string               Data key to use for N^2
     s::float                    Spline smoothing factor (smaller values
                                 give a noisier result)
     """
     if rhokey not in self.fields:
         raise FieldError("add_Nsquared requires in-situ density")
     msk = self.nanmask((rhokey, depthkey))
     rho = self[rhokey][~msk]
     z = self[depthkey][~msk]
     rhospl = UnivariateSpline(z, rho, s=s)
     drhodz = np.asarray([-rhospl.derivatives(_z)[1] for _z in z])
     N2 = np.empty(len(self), dtype=np.float64)
     N2[msk] = np.nan
     N2[~msk] = -G / rho * drhodz
     return self._addkeydata(N2key, N2)

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
spl = UnivariateSpline(mass, luminosity, k=4, s=0)

print mass, luminosity
lum_deriv = luminosity
#lum_deriv[0] = spl.derivatives(mass[0])[1]


xs = np.linspace(0.5, 100., 1000)
#xs = np.logspace(0.5, 100., 1000)

xs_deriv = np.zeros(len(mass))
#xs_deriv = np.zeros(len(xs))
print len(xs)
for i in range(len(mass)):
#    print i, xs[i]
    xs_deriv[i] = spl.derivatives(mass[i])[1]
#    xs_deriv[i] = spl.derivatives(xs[i])[1]
    print mass[i], luminosity[i], xs_deriv[i], xs_deriv[i] * mass[i] / luminosity[i]

#log_lum = np.zeros(len(mass))
#for i in range(len(mass)):#(0,1,2):#range(len(mass)):
#    lum_deriv[i] = spl.derivatives(mass[i])[1]
#    log_lum[i] = lum_deriv[i] / (mass[i]*mass[i])
#    print i, mass[i], luminosity[i], lum_deriv[i], log_lum[i]
#
#print lum_orig
#luminosity = lum_orig
#print mass, luminosity, lum_deriv, log_lum
#print len(xs), len(xs_deriv)

pl.clf()

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
	def estimate(self, observedLC):
		"""!
		Estimate intrinsicFlux, period, eccentricity, omega, tau, & a2sini 
		"""
		## intrinsicFluxEst
		maxPeriodFactor = 10.0
		model = LombScargleFast().fit(observedLC.t, observedLC.y, observedLC.yerr)
		periods, power = model.periodogram_auto(nyquist_factor = observedLC.numCadences)
		model.optimizer.period_range = (2.0*np.mean(observedLC.t[1:] - observedLC.t[:-1]), maxPeriodFactor*observedLC.T)
		periodEst = model.best_period
		numIntrinsicFlux = 100
		lowestFlux = np.min(observedLC.y[np.where(observedLC.mask == 1.0)])
		highestFlux = np.max(observedLC.y[np.where(observedLC.mask == 1.0)])
		intrinsicFlux = np.linspace(np.min(observedLC.y[np.where(observedLC.mask == 1.0)]), np.max(observedLC.y[np.where(observedLC.mask == 1.0)]), num = numIntrinsicFlux)
		intrinsicFluxList = list()
		totalIntegralList = list()
		for f in xrange(1, numIntrinsicFlux - 1):
			beamedLC = observedLC.copy()
			beamedLC.x = np.require(np.zeros(beamedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
			for i in xrange(beamedLC.numCadences):
				beamedLC.y[i] = observedLC.y[i]/intrinsicFlux[f]
				beamedLC.yerr[i] = observedLC.yerr[i]/intrinsicFlux[f]
			dopplerLC = beamedLC.copy()
			dopplerLC.x = np.require(np.zeros(dopplerLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
			for i in xrange(observedLC.numCadences):
				dopplerLC.y[i] = math.pow(beamedLC.y[i], 1.0/3.44)
				dopplerLC.yerr[i] = (1.0/3.44)*math.fabs(dopplerLC.y[i]*(beamedLC.yerr[i]/beamedLC.y[i]))
			dzdtLC = dopplerLC.copy()
			dzdtLC.x = np.require(np.zeros(dopplerLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
			for i in xrange(observedLC.numCadences):
				dzdtLC.y[i] = 1.0 - (1.0/dopplerLC.y[i])
				dzdtLC.yerr[i] = math.fabs((-1.0*dopplerLC.yerr[i])/math.pow(dopplerLC.y[i], 2.0))
			foldedLC = dzdtLC.fold(periodEst)
			foldedLC.x = np.require(np.zeros(foldedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
			integralSpline = UnivariateSpline(foldedLC.t[np.where(foldedLC.mask == 1.0)], foldedLC.y[np.where(foldedLC.mask == 1.0)], 1.0/foldedLC.yerr[np.where(foldedLC.mask == 1.0)], k = 3, s = None, check_finite = True)
			totalIntegral = math.fabs(integralSpline.integral(foldedLC.t[0], foldedLC.t[-1]))
			intrinsicFluxList.append(intrinsicFlux[f])
			totalIntegralList.append(totalIntegral)
		intrinsicFluxEst = intrinsicFluxList[np.where(np.array(totalIntegralList) == np.min(np.array(totalIntegralList)))[0][0]]

		## periodEst
		for i in xrange(beamedLC.numCadences):
			beamedLC.y[i] = observedLC.y[i]/intrinsicFluxEst
			beamedLC.yerr[i] = observedLC.yerr[i]/intrinsicFluxEst
			dopplerLC.y[i] = math.pow(beamedLC.y[i], 1.0/3.44)
			dopplerLC.yerr[i] = (1.0/3.44)*math.fabs(dopplerLC.y[i]*(beamedLC.yerr[i]/beamedLC.y[i]))
			dzdtLC.y[i] = 1.0 - (1.0/dopplerLC.y[i])
			dzdtLC.yerr[i] = math.fabs((-1.0*dopplerLC.yerr[i])/math.pow(dopplerLC.y[i], 2.0))
		model = LombScargleFast().fit(dzdtLC.t, dzdtLC.y, dzdtLC.yerr)
		periods, power = model.periodogram_auto(nyquist_factor = dzdtLC.numCadences)
		model.optimizer.period_range = (2.0*np.mean(dzdtLC.t[1:] - dzdtLC.t[:-1]), maxPeriodFactor*dzdtLC.T)
		periodEst = model.best_period

		## eccentricityEst & omega2Est
		# First find a full period going from rising to falling. 
		risingSpline = UnivariateSpline(dzdtLC.t[np.where(dzdtLC.mask == 1.0)], dzdtLC.y[np.where(dzdtLC.mask == 1.0)], 1.0/dzdtLC.yerr[np.where(dzdtLC.mask == 1.0)], k = 3, s = None, check_finite = True)
		risingSplineRoots = risingSpline.roots()
		firstRoot = risingSplineRoots[0]
		if risingSpline.derivatives(risingSplineRoots[0])[1] > 0.0:
			tRising = risingSplineRoots[0]
		else:
			tRising = risingSplineRoots[1]
		# Now fold the LC starting at tRising and going for a full period.
		foldedLC = dzdtLC.fold(periodEst, tStart = tRising)
		foldedLC.x = np.require(np.zeros(foldedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
		# Fit the folded LC with a spline to figure out alpha and beta
		fitLC = foldedLC.copy()
		foldedSpline = UnivariateSpline(foldedLC.t[np.where(foldedLC.mask == 1.0)], foldedLC.y[np.where(foldedLC.mask == 1.0)], 1.0/foldedLC.yerr[np.where(foldedLC.mask == 1.0)], k = 3, s = 2*foldedLC.numCadences, check_finite = True)
		for i in xrange(fitLC.numCadences):
			fitLC.x[i] = foldedSpline(fitLC.t[i])
		# Now get the roots and find the falling root
		tZeros = foldedSpline.roots()
		if tZeros.shape[0] == 1: # We have found just tFalling
			tFalling = tZeros[0]
			tRising = fitLC.t[0]
			startIndex = 0
			tFull = fitLC.t[-1]
			stopIndex = fitLC.numCadences
		elif tZeros.shape[0] == 2: # We have found tFalling and one of tRising or tFull
			if foldedSpline.derivatives(tZeros[0])[1] < 0.0:
				tFalling = tZeros[0]
				tFull = tZeros[1]
				stopIndex = np.where(fitLC.t < tFull)[0][-1]
				tRising = fitLC.t[0]
				startIndex = 0
			elif foldedSpline.derivatives(tZeros[0])[1] > 0.0:
				if foldedSpline.derivatives(tZeros[1])[1] < 0.0:
					tRising = tZeros[0]
					startIndex = np.where(fitLC.t > tRising)[0][0]
					tFalling = tZeros[1]
					tFull = fitLC.t[-1]
					stopIndex = fitLC.numCadences
				else:
					raise RuntimeError('Could not determine alpha & omega correctly because the first root is rising but the second root is not falling!')
		elif tZeros.shape[0] == 3:
			tRising = tZeros[0]
			startIndex = np.where(fitLC.t > tRising)[0][0]
			tFalling = tZeros[1]
			tFull = tZeros[2]
			stopIndex = np.where(fitLC.t < tFull)[0][-1]
#.........这里部分代码省略.........

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
        ####
        ydelta = Results_a[stage,1]
        threslist = Results_a[stage,0]
        ####
    
        #Intermittency profile

        Dv = np.empty((len(threslist),9),dtype=np.double)
    
        epsilon = 2**np.linspace(0,8,9)
        for ithres,thres in enumerate(threslist):
            sp = UnivariateSpline(np.log10(epsilon),
                                  np.log10(box_counting_vol[ithres,:9]),
                                  k=5)
            for j in range(len(epsilon)):
                Dv[ithres,j] = sp.derivatives(np.log10(epsilon[j]))[1]
    
        Ds = np.empty((len(threslist),9),dtype=np.double)
        epsilon = 2**np.linspace(0,8,9)
        for ithres,thres in enumerate(threslist):
            sp = UnivariateSpline(np.log10(epsilon),
                                  np.log10(box_counting_sur[ithres,:9]),
                                  k=5)
            for j in range(len(epsilon)):
                Ds[ithres,j] = sp.derivatives(np.log10(epsilon[j]))[1]
    
        gamma = np.array([r[0] for r in Results_a[stage,2:]])
        genusv = np.array([r[1] for r in Results_a[stage,2:]])
        genuss = np.array([r[3] for r in Results_a[stage,2:]])
        
        if stage==1:

# 需要导入模块: from scipy.interpolate import UnivariateSpline [as 别名]
# 或者: from scipy.interpolate.UnivariateSpline import derivatives [as 别名]
    def estimate(self, observedLC):
        """!
        Estimate intrinsicFlux, period, eccentricity, omega, tau, & a2sini
        """
        # fluxEst
        if observedLC.numCadences > 50:
            model = gatspy.periodic.LombScargleFast(optimizer_kwds={"quiet": True}).fit(observedLC.t,
                                                                                        observedLC.y,
                                                                                        observedLC.yerr)
        else:
            model = gatspy.periodic.LombScargle(optimizer_kwds={"quiet": True}).fit(observedLC.t,
                                                                                    observedLC.y,
                                                                                    observedLC.yerr)
        periods, power = model.periodogram_auto(nyquist_factor=observedLC.numCadences)
        model.optimizer.period_range = (
            2.0*np.mean(observedLC.t[1:] - observedLC.t[:-1]), observedLC.T)
        periodEst = model.best_period
        numIntrinsicFlux = 100
        lowestFlux = np.min(observedLC.y[np.where(observedLC.mask == 1.0)])
        highestFlux = np.max(observedLC.y[np.where(observedLC.mask == 1.0)])
        intrinsicFlux = np.linspace(np.min(observedLC.y[np.where(observedLC.mask == 1.0)]), np.max(
            observedLC.y[np.where(observedLC.mask == 1.0)]), num=numIntrinsicFlux)
        intrinsicFluxList = list()
        totalIntegralList = list()
        for f in xrange(1, numIntrinsicFlux - 1):
            beamedLC = observedLC.copy()
            beamedLC.x = np.require(np.zeros(beamedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
            for i in xrange(beamedLC.numCadences):
                beamedLC.y[i] = observedLC.y[i]/intrinsicFlux[f]
                beamedLC.yerr[i] = observedLC.yerr[i]/intrinsicFlux[f]
            dopplerLC = beamedLC.copy()
            dopplerLC.x = np.require(np.zeros(dopplerLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
            for i in xrange(observedLC.numCadences):
                dopplerLC.y[i] = math.pow(beamedLC.y[i], 1.0/3.44)
                dopplerLC.yerr[i] = (1.0/3.44)*math.fabs(dopplerLC.y[i]*(beamedLC.yerr[i]/beamedLC.y[i]))
            dzdtLC = dopplerLC.copy()
            dzdtLC.x = np.require(np.zeros(dopplerLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
            for i in xrange(observedLC.numCadences):
                dzdtLC.y[i] = 1.0 - (1.0/dopplerLC.y[i])
                dzdtLC.yerr[i] = math.fabs((-1.0*dopplerLC.yerr[i])/math.pow(dopplerLC.y[i], 2.0))
            foldedLC = dzdtLC.fold(periodEst)
            foldedLC.x = np.require(np.zeros(foldedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
            integralSpline = UnivariateSpline(
                foldedLC.t[np.where(foldedLC.mask == 1.0)], foldedLC.y[np.where(foldedLC.mask == 1.0)],
                1.0/foldedLC.yerr[np.where(foldedLC.mask == 1.0)], k=3, s=None, check_finite=True)
            totalIntegral = math.fabs(integralSpline.integral(foldedLC.t[0], foldedLC.t[-1]))
            intrinsicFluxList.append(intrinsicFlux[f])
            totalIntegralList.append(totalIntegral)
        fluxEst = intrinsicFluxList[
            np.where(np.array(totalIntegralList) == np.min(np.array(totalIntegralList)))[0][0]]

        # periodEst
        for i in xrange(beamedLC.numCadences):
            beamedLC.y[i] = observedLC.y[i]/fluxEst
            beamedLC.yerr[i] = observedLC.yerr[i]/fluxEst
            dopplerLC.y[i] = math.pow(beamedLC.y[i], 1.0/3.44)
            dopplerLC.yerr[i] = (1.0/3.44)*math.fabs(dopplerLC.y[i]*(beamedLC.yerr[i]/beamedLC.y[i]))
            dzdtLC.y[i] = 1.0 - (1.0/dopplerLC.y[i])
            dzdtLC.yerr[i] = math.fabs((-1.0*dopplerLC.yerr[i])/math.pow(dopplerLC.y[i], 2.0))
            if observedLC.numCadences > 50:
                model = gatspy.periodic.LombScargleFast(optimizer_kwds={"quiet": True}).fit(dzdtLC.t,
                                                                                            dzdtLC.y,
                                                                                            dzdtLC.yerr)
            else:
                model = gatspy.periodic.LombScargle(optimizer_kwds={"quiet": True}).fit(dzdtLC.t,
                                                                                        dzdtLC.y,
                                                                                        dzdtLC.yerr)
        periods, power = model.periodogram_auto(nyquist_factor=dzdtLC.numCadences)
        model.optimizer.period_range = (2.0*np.mean(dzdtLC.t[1:] - dzdtLC.t[:-1]), dzdtLC.T)
        periodEst = model.best_period

        # eccentricityEst & omega2Est
        # First find a full period going from rising to falling.
        risingSpline = UnivariateSpline(
            dzdtLC.t[np.where(dzdtLC.mask == 1.0)], dzdtLC.y[np.where(dzdtLC.mask == 1.0)],
            1.0/dzdtLC.yerr[np.where(dzdtLC.mask == 1.0)], k=3, s=None, check_finite=True)
        risingSplineRoots = risingSpline.roots()
        firstRoot = risingSplineRoots[0]
        if risingSpline.derivatives(risingSplineRoots[0])[1] > 0.0:
            tRising = risingSplineRoots[0]
        else:
            tRising = risingSplineRoots[1]
        # Now fold the LC starting at tRising and going for a full period.
        foldedLC = dzdtLC.fold(periodEst, tStart=tRising)
        foldedLC.x = np.require(np.zeros(foldedLC.numCadences), requirements=['F', 'A', 'W', 'O', 'E'])
        # Fit the folded LC with a spline to figure out alpha and beta
        fitLC = foldedLC.copy()
        foldedSpline = UnivariateSpline(
            foldedLC.t[np.where(foldedLC.mask == 1.0)], foldedLC.y[np.where(foldedLC.mask == 1.0)],
            1.0/foldedLC.yerr[np.where(foldedLC.mask == 1.0)], k=3, s=2*foldedLC.numCadences,
            check_finite=True)
        for i in xrange(fitLC.numCadences):
            fitLC.x[i] = foldedSpline(fitLC.t[i])
        # Now get the roots and find the falling root
        tZeros = foldedSpline.roots()

        # Find tRising, tFalling, tFull, startIndex, & stopIndex via DBSCAN #######################
        # Find the number of clusters
        '''dbsObj = DBSCAN(eps = periodEst/10.0, min_samples = 1)
        db = dbsObj.fit(tZeros.reshape(-1,1))
#.........这里部分代码省略.........