diff --git a/main/flaredetector/util.py b/main/flaredetector/util.py index 5c82c5c..babbcfc 100644 --- a/main/flaredetector/util.py +++ b/main/flaredetector/util.py @@ -86,30 +86,11 @@ def convertStarndardIndexToFoldedIndex(foldedLc, standardIndex): def singleSine(t, A, w, p, c): return A * np.sin(w*t + p) + c -def doubleSine(t, A1, w1, p1, c1, A2, w2, p2, c2, split1, split2): - y = np.zeros_like(t) - m1 = t < split1 - m2 = (t >= split1) & (t < split2) - m3 = t >= split2 - y[m1] = A1 * np.sin(w1 * t[m1] + p1) + c1 - y[m2] = A2 * np.sin(w2 * t[m2] + p2) + (A1 * np.sin(w1 * split1 + p1) + c1 - A2 * np.sin(w2 * split1 + p2)) - y[m3] = A1 * np.sin(w1 * t[m3] + p1) + (A2 * np.sin(w2 * split2 + p2) + (A1 * np.sin(w1 * split1 + p1) + c1 - A2 * np.sin(w2 * split1 + p2)) - A1 * np.sin(w1 * split2 + p1)) - return y - def fitSingleSine(phase, flux): initGuess = [np.ptp(flux)/2, 2 * np.pi / (np.max(phase) - np.min(phase)), 0, np.mean(flux)] popt, _ = curve_fit(singleSine, phase, flux, p0=initGuess, maxfev=100000) return singleSine(phase, *popt) -def fitDoubleSine(phase, flux): - initSplit1 = abs(phase[0]) * 2 / 3 + phase[0] - initSplit2 = abs(phase[0]) * 4 / 3 + phase[0] - initGuess = [np.ptp(flux) / 2, 2 * np.pi / (np.max(phase) - np.min(phase)), 0, np.mean(flux), - np.ptp(flux) / 2, 2 * np.pi / (np.max(phase) - np.min(phase)), 0, np.mean(flux), - initSplit1, initSplit2] - popt, _ = curve_fit(doubleSine, phase, flux, p0=initGuess, maxfev=100000) - return doubleSine(phase, *popt) - def fitPolynomial(phase, flux, degree): return sum(p * phase**i for i, p in enumerate(Polynomial.fit(phase, flux, degree).convert().coef)) @@ -133,41 +114,24 @@ def getFoldedBestFit(foldedLc): phase = foldedLc.phase[filt].value flux = flux[filt] #popt, pcov = curve_fit(sine, phase, flux[filt], maxfev=100000) - import time - start = time.time() singleFit = fitSingleSine(phase, flux) - singleTime = time.time() - start - start = time.time() - doubleFit = fitDoubleSine(phase, flux) - doubleTime = time.time() - start - start = time.time() polyFit = fitPolynomial(phase, flux, 10) - polyTime = time.time() - start singleRSS = compureRSS(singleFit, flux) - doubleRSS = compureRSS(doubleFit, flux) polyRSS = compureRSS(polyFit, flux) singleAIC = computeAIC(singleRSS, 4, len(flux)) singleBIC = computeBIC(singleRSS, 4, len(flux)) - doubleAIC = computeAIC(doubleRSS, 10, len(flux)) - doubleBIC = computeBIC(doubleRSS, 10, len(flux)) polyAIC = computeAIC(polyRSS, 10, len(flux)) polyBIC = computeBIC(polyRSS, 10, len(flux)) - if (singleAIC < doubleAIC and singleAIC < polyAIC and singleBIC < doubleBIC and singleBIC < polyBIC): + if (singleAIC < polyAIC and singleBIC < polyBIC): print("Single sine preferred") sineFit = singleFit - elif (doubleAIC < singleAIC and doubleAIC < polyAIC and doubleBIC < singleBIC and doubleBIC < polyBIC): - print("Double sine preferred") - sineFit = doubleFit else: print("Polynomial preferred") sineFit = polyFit - print("Single sine time: ", singleTime) - print("Double sine time: ", doubleTime) - print("Polynomial time: ", polyTime) sineFit *= multi return phase, sineFit