flaredetector: util: remove double sine fit

it nnever fits better than a polynomial, and often takes a lot
longer to compute

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