diff --git a/main/astrodatagui/AstrodataGUI.py b/main/astrodatagui/AstrodataGUI.py index 2038677..f079be2 100644 --- a/main/astrodatagui/AstrodataGUI.py +++ b/main/astrodatagui/AstrodataGUI.py @@ -62,6 +62,7 @@ class AstrodataGUI(QtWidgets.QMainWindow): self.btOpenSavedFlareCountFile.clicked.connect(self.btOpenSavedFlareCountFileClicked) self.flaredetectorPreview.periodsCalculated.connect(self.updatePeriods) + self.flaredetectorPreview.epochCalculated.connect(self.updateEpochPeriod) self.setupCustomSimbadQueries() self.show() @@ -205,6 +206,9 @@ class AstrodataGUI(QtWidgets.QMainWindow): def updatePeriods(self, periods: list): self.edPlotFoldPeriod.setText(str(periods[0].value)) + def updateEpochPeriod(self, epoch: float): + self.edPlotFoldEpochTime.setText(str(epoch)) + def btCombineSeqClicked(self): #self.gbPlotOptionsNormalize.setEnabled(False) mainName = self.listDownloaded.currentItem().text() diff --git a/main/astrodatagui/CalcAllFlaresThread.py b/main/astrodatagui/CalcAllFlaresThread.py index 2ee84cd..edc1171 100644 --- a/main/astrodatagui/CalcAllFlaresThread.py +++ b/main/astrodatagui/CalcAllFlaresThread.py @@ -4,32 +4,54 @@ import concurrent.futures from ..util.MinimalLightCurve import read import pandas as pd from ..flaredetector.flaredetector import calculateFlareFitsForLightcurve -from ..flaredetector.util import getTotalValidDataInSeconds +from ..flaredetector.util import * def getFlareCount(filesDict): lc = read(filesDict["FilePath"]) lc.flux = lc["sap_flux"] lc.flux_err = lc["sap_flux_err"] - sapPeaks, sapFits = calculateFlareFitsForLightcurve(lc.flatten(), minimalLC=True) + sapPeaks, sapFits = calculateFlareFitsForLightcurve(lc.flatten()) sapValSec, sapTds = getTotalValidDataInSeconds(lc, "sap_flux") + sapPeriodogram = lc.to_periodogram() + sapPeakPeriod = sapPeriodogram.period[findMaxIndices(sapPeriodogram, num=4, distance=100, sortByHighest=True)[0]] + sapEpochTime = getEpochTime(lc) + sapFoldedLC = lc.fold(period=sapPeakPeriod, epoch_time=sapEpochTime) + sapPhase, sapSineFit, sapFitType = getFoldedBestFit(sapFoldedLC) + sapMinima, sapMaxima = getFoldedFitPeakValley(sapSineFit) + sapminPhasesBoundsIndices, sapmaxPhasesBoundsIndices = getPhaseRangesNearPeak((sapMinima, sapMaxima), sapPhase) lc.flux = lc["pdcsap_flux"] lc.flux_err = lc["pdcsap_flux_err"] - pdcsapPeaks, pdcsapFits = calculateFlareFitsForLightcurve(lc.flatten(), minimalLC=True) + pdcsapPeaks, pdcsapFits = calculateFlareFitsForLightcurve(lc.flatten()) pdcsapValSec, pdcsapTds = getTotalValidDataInSeconds(lc, "pdcsap_flux") + pdcsapPeriodogram = lc.to_periodogram() + pdcsapPeakPeriod = pdcsapPeriodogram.period[findMaxIndices(pdcsapPeriodogram, num=4, distance=100, sortByHighest=True)[0]] + pdcsapEpochTime = getEpochTime(lc) + pdcsapFoldedLC = lc.fold(period=pdcsapPeakPeriod, epoch_time=pdcsapEpochTime) + pdcsapPhase, pdcsapSineFit, pdcsapFitType = getFoldedBestFit(pdcsapFoldedLC) + pdcsapMinima, pdcsapMaxima = getFoldedFitPeakValley(pdcsapSineFit) + pdcsapminPhasesBoundsIndices, pdcsapmaxPhasesBoundsIndices = getPhaseRangesNearPeak((pdcsapMinima, pdcsapMaxima), pdcsapPhase) filesDict["sapPeaks"] = sapPeaks filesDict["sapPeaksCount"] = len(sapPeaks) filesDict["sapFits"] = sapFits filesDict["sapValidTimespans"] = sapTds filesDict["sapValidSeconds"] = sapValSec + filesDict["sapPeriod"] = sapPeakPeriod.value + filesDict["sapPeriodMinima"] = sapMinima + filesDict["sapPeriodMinimaBoundaries"] = sapminPhasesBoundsIndices + filesDict["sapPeriodMaxima"] = sapMaxima + filesDict["sapPeriodMaximaBoundaries"] = sapmaxPhasesBoundsIndices filesDict["pdcsapPeaks"] = pdcsapPeaks filesDict["pdcsapPeaksCount"] = len(pdcsapPeaks) filesDict["pdcsapFits"] = pdcsapFits filesDict["pdcsapValidTimespans"] = pdcsapTds filesDict["pdcsapValidSeconds"] = pdcsapValSec + filesDict["pdcsapPeriod"] = pdcsapPeakPeriod.value + filesDict["pdcsapPeriodMinima"] = pdcsapMinima + filesDict["pdcsapPeriodMinimaBoundaries"] = pdcsapminPhasesBoundsIndices + filesDict["pdcsapPeriodMaxima"] = pdcsapMaxima + filesDict["pdcsapPeriodMaximaBoundaries"] = pdcsapmaxPhasesBoundsIndices del lc - if(filesDict["StarName"] == "CD-51 13128" and filesDict["Sequence"] == 1): - print(filesDict) return filesDict class CalcAllFlaresThread(QThread): diff --git a/main/astrodatagui/FlareSummaryPlotGUI.py b/main/astrodatagui/FlareSummaryPlotGUI.py index 211f5a1..e01880d 100644 --- a/main/astrodatagui/FlareSummaryPlotGUI.py +++ b/main/astrodatagui/FlareSummaryPlotGUI.py @@ -43,6 +43,8 @@ class FlareSummaryPlotGUI(QWidget): self.btShowFlaresPerStar.clicked.connect(self.btShowFlaresPerStarClicked) self.btShowFlaresPerStarNormalized = QPushButton("Show Flares/Star Normalized") self.btShowFlaresPerStarNormalized.clicked.connect(self.btShowFlaresPerStarNormalizedClicked) + self.btShowPeriods = QPushButton("Show Mean Periods") + self.btShowPeriods.clicked.connect(self.btShowPeriodsClicked) self.cbKepler = QCheckBox("Kepler") self.cbKepler.setChecked(True) @@ -68,6 +70,7 @@ class FlareSummaryPlotGUI(QWidget): self.buttonGridLayout.addWidget(self.btShowFlaresPerFile, 0, 1) self.buttonGridLayout.addWidget(self.btShowFlaresPerStar, 0, 2) self.buttonGridLayout.addWidget(self.btShowFlaresPerStarNormalized, 0, 3) + self.buttonGridLayout.addWidget(self.btShowPeriods, 0, 4) self.buttonGridLayout.addWidget(QLabel("Sources: "), 1, 0) self.buttonGridLayout.addWidget(self.cbKepler, 1, 1) @@ -169,7 +172,18 @@ class FlareSummaryPlotGUI(QWidget): 'pdcsapFits', 'pdcsapPeaks', 'sapFits', - 'sapPeaks']) + 'sapPeaks', + 'sapPeriod', + 'sapPeriodMinima', + 'sapPeriodMinimaBoundaries', + 'sapPeriodMaxima', + 'sapPeriodMaximaBoundaries', + 'pdcsapValidTimespans', + 'pdcsapPeriod', + 'pdcsapPeriodMinima', + 'pdcsapPeriodMinimaBoundaries', + 'pdcsapPeriodMaxima', + 'pdcsapPeriodMaximaBoundaries']) showSourceFilter = np.full(len(data), False) @@ -248,7 +262,19 @@ class FlareSummaryPlotGUI(QWidget): 'pdcsapFits', 'pdcsapPeaks', 'sapFits', - 'sapPeaks']) + 'sapValidTimespans', + 'sapPeaks', + 'sapPeriod', + 'sapPeriodMinima', + 'sapPeriodMinimaBoundaries', + 'sapPeriodMaxima', + 'sapPeriodMaximaBoundaries', + 'pdcsapValidTimespans', + 'pdcsapPeriod', + 'pdcsapPeriodMinima', + 'pdcsapPeriodMinimaBoundaries', + 'pdcsapPeriodMaxima', + 'pdcsapPeriodMaximaBoundaries']) showSourceFilter = np.full(len(data), False) @@ -329,4 +355,121 @@ class FlareSummaryPlotGUI(QWidget): self.figureAxis.set_ylabel("Flare count per Week") self.figureAxis.set_xlabel("Star Number") self.figureAxis.legend() - self.figure.canvas.draw_idle() \ No newline at end of file + self.figure.canvas.draw_idle() + + def btShowPeriodsClicked(self): + data = self.starFLareDictList.drop(columns=['Distance', + 'DistanceUnit', + 'FilePath', + 'RotVel', + 'RotVelUnit', + 'Sequence', + 'pdcsapFits', + 'pdcsapPeaks', + 'sapFits', + 'sapValidTimespans', + 'sapPeaks', + 'sapPeriodMinima', + 'sapPeriodMinimaBoundaries', + 'sapPeriodMaxima', + 'sapPeriodMaximaBoundaries', + 'pdcsapValidTimespans', + 'pdcsapPeriodMinima', + 'pdcsapPeriodMinimaBoundaries', + 'pdcsapPeriodMaxima', + 'pdcsapPeriodMaximaBoundaries', + 'sapPeaksCount', + 'pdcsapPeaksCount', + 'sapValidSeconds', + 'pdcsapValidSeconds']) + + showSourceFilter = np.full(len(data), False) + + if(self.cbKepler.isChecked()): + showKepler = data["Source"] == "Kepler" + showSourceFilter |= showKepler + if(self.cbK2.isChecked()): + showK2 = data["Source"] == "K2" + showSourceFilter |= showK2 + if(self.cbTESS.isChecked()): + showTESS = data["Source"] == "TESS" + showSourceFilter |= showTESS + + print(data["sapPeriod"]) + + data = data[showSourceFilter].groupby(["StarName", "SpType"], + as_index=False).agg({"sapPeriod": "mean", + "pdcsapPeriod": "mean"}) + x = np.arange(start=1, stop=len(data)+1, step=1, dtype=int) + + if(self.cbSpTypeL.isChecked()): + Lfilter = data["SpType"].str.startswith("L") + Lfilter &= showSourceFilter + if(self.cbSpTypeM.isChecked()): + Mfilter = data["SpType"].str.startswith("M") + Mfilter &= showSourceFilter + if(self.cbSpTypeK.isChecked()): + Kfilter = data["SpType"].str.startswith("K") + Kfilter &= showSourceFilter + if(self.cbSpTypeG.isChecked()): + Gfilter = data["SpType"].str.startswith("G") + Gfilter &= showSourceFilter + if(self.cbSpTypeF.isChecked()): + Ffilter = data["SpType"].str.startswith("F") + Ffilter &= showSourceFilter + if(self.cbSpTypeUnknown.isChecked()): + Unknownfilter = data["SpType"].str.startswith("-") + Unknownfilter &= showSourceFilter + + self.figureAxis.clear() + if(self.cbSpTypeL.isChecked()): + if(np.any(Lfilter)): + self.figureAxis.scatter(x[Lfilter], data[Lfilter]["sapPeriod"], + marker="o", color="brown", label="L SAP Period") + self.figureAxis.scatter(x[Lfilter], data[Lfilter]["pdcsapPeriod"], + marker="x", color="brown", label="L PDCSAP Period") + if(self.cbSpTypeM.isChecked()): + if(np.any(Mfilter)): + self.figureAxis.scatter(x[Mfilter], data[Mfilter]["sapPeriod"], + marker="o", color="red", label="M SAP Period") + self.figureAxis.scatter(x[Mfilter], data[Mfilter]["pdcsapPeriod"], + marker="x", color="red", label="M PDCSAP Period") + if(self.cbSpTypeK.isChecked()): + if(np.any(Kfilter)): + self.figureAxis.scatter(x[Kfilter], data[Kfilter]["sapPeriod"], + marker="o", color="orange", label="K SAP Period") + self.figureAxis.scatter(x[Kfilter], data[Kfilter]["pdcsapPeriod"], + marker="x", color="orange", label="K PDCSAP Period") + if(self.cbSpTypeG.isChecked()): + if(np.any(Gfilter)): + self.figureAxis.scatter(x[Gfilter], data[Gfilter]["sapPeriod"], + marker="o", color="yellow", label="G SAP Period") + self.figureAxis.scatter(x[Gfilter], data[Gfilter]["pdcsapPeriod"], + marker="x", color="yellow", label="G PDCSAP Period") + if(self.cbSpTypeF.isChecked()): + if(np.any(Ffilter)): + self.figureAxis.scatter(x[Ffilter], data[Ffilter]["sapPeriod"], + marker="o", color="greenyellow", label="F SAP Period") + self.figureAxis.scatter(x[Ffilter], data[Ffilter]["pdcsapPeriod"], + marker="x", color="greenyellow", label="F PDCSAP Period") + if(self.cbSpTypeUnknown.isChecked()): + if(np.any(Unknownfilter)): + self.figureAxis.scatter(x[Unknownfilter], data[Unknownfilter]["sapPeriod"], + marker="o", color="gray", label="Unknown SAP Period") + self.figureAxis.scatter(x[Unknownfilter], data[Unknownfilter]["pdcsapPeriod"], + marker="x", color="gray", label="Unknown PDCSAP Period") + + self.figureAxis.set_ylabel("Period / days") + self.figureAxis.set_xlabel("Star Number") + self.figureAxis.legend() + self.figure.canvas.draw_idle() + pass + + def btShowNumMinimaMaximaClicked(self): + pass + + def btShowFlaresInMinimaMaximaClicked(self): + pass + + def btShowFlaresInMinimaMaximaPerMinimaMaximaClicked(self): + pass diff --git a/main/astrodatagui/ui/FlaredetectorWidget.py b/main/astrodatagui/ui/FlaredetectorWidget.py index 292517e..24d9e47 100644 --- a/main/astrodatagui/ui/FlaredetectorWidget.py +++ b/main/astrodatagui/ui/FlaredetectorWidget.py @@ -254,7 +254,7 @@ class FlaredetectorWidget(QtWidgets.QWidget): cycle, foldedIndex = convertStarndardIndexToFoldedIndex(lc, peak["StandardIndex"]) self.figureAxis.plot(lc.phase[lc.cycle == cycle][foldedIndex].value, lc.flux[lc.cycle == cycle][foldedIndex], "x", color="red") - phase, sineFit = getFoldedBestFit(lc) + phase, sineFit, _ = getFoldedBestFit(lc) self.figureAxis.plot(phase, sineFit, color="red") minPhasesBoundsIndices, maxPhasesBoundsIndices = getPhaseRangesNearPeak(getFoldedFitPeakValley(sineFit), phase) plotPhaseRangesNearPeak((minPhasesBoundsIndices, maxPhasesBoundsIndices), phase, ax=self.figureAxis) diff --git a/main/flaredetector/flaredetector.py b/main/flaredetector/flaredetector.py index 754ef7f..f56fbca 100644 --- a/main/flaredetector/flaredetector.py +++ b/main/flaredetector/flaredetector.py @@ -37,12 +37,19 @@ def calculateFitForFlare(lightcurve, peakIndex): popt, pcov = curve_fit(lambda x, l, w: GaussExpo(x, Peak, l, w), flareFitDataPoints, flareData) pcovDiag = np.diag(pcov) fit = GaussExpo(flareFitDataPoints, Peak, *popt) + rss = compureRSS(fit, flareData) + aic = computeAIC(rss, 2, len(flareData)) + bic = computeBIC(rss, 2, len(flareData)) + r2 = 1 - rss/computeTotalSoS(flareData) + #print("Fit r2: ", r2) + if(r2 < 0.8): + raise Exception("Fit varies too much") flareFitDataPoints = flareFitDataPoints + peakIndex return np.array(Peak+1), flareFitDataPoints, np.array(fit+1) -def calculateFlareFitsForLightcurve(lightcurve, num=100, distance=1, height=0.05, minimalLC=False): +def calculateFlareFitsForLightcurve(lightcurve, num=100, distance=1, height=0.05): maxIndices = findMaxIndices(lightcurve, num, distance=distance, height=height, - sortByHighest=True, minimalLC=minimalLC) + sortByHighest=True) peaks = [] fits = [] for index in maxIndices: diff --git a/main/flaredetector/util.py b/main/flaredetector/util.py index babbcfc..7a207e8 100644 --- a/main/flaredetector/util.py +++ b/main/flaredetector/util.py @@ -1,21 +1,15 @@ import numpy as np from numpy.polynomial.polynomial import Polynomial -from scipy.signal import find_peaks +from scipy.signal import find_peaks, argrelextrema from scipy.optimize import curve_fit -def findMaxIndices(lc, num=3, distance=100, height=(None, None), sortByHighest=False, - minimalLC=False): - if(minimalLC): +def findMaxIndices(lc, num=3, distance=100, height=(None, None), sortByHighest=False): + if(hasattr(lc, "power")): + data = lc.power + elif(hasattr(lc, "flux")): data = lc.flux else: - from lightkurve.periodogram import Periodogram as OrigPG - from lightkurve.lightcurve import LightCurve as OrigLC - if(isinstance(lc, OrigPG)): - data = lc.power - elif(isinstance(lc, OrigLC)): - data = lc.flux - else: - return np.zeros(num)-1 + return np.zeros(num)-1 peak_indices, peak_dict = find_peaks(data, height=height, distance=distance) peak_heights = peak_dict["peak_heights"] @@ -88,16 +82,22 @@ def singleSine(t, A, w, p, c): 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) + popt, _ = curve_fit(singleSine, phase, flux, p0=initGuess, maxfev=300) return singleSine(phase, *popt) +def polynomial(phase, flux, degree): + return Polynomial.fit(phase, flux, degree).convert().coef + def fitPolynomial(phase, flux, degree): - return sum(p * phase**i for i, p in enumerate(Polynomial.fit(phase, flux, degree).convert().coef)) + return sum(p * phase**i for i, p in enumerate(polynomial(phase, flux, degree))) def compureRSS(fit, flux): residuals = flux - fit return np.sum(residuals**2) +def computeTotalSoS(flux): + return np.sum((flux - np.mean(flux))**2) + def computeAIC(rss, numParams, numDataPoints): return 2 * numParams + numDataPoints * np.log(rss/numDataPoints) @@ -113,30 +113,46 @@ def getFoldedBestFit(foldedLc): flux = foldedLc.normalize().flux phase = foldedLc.phase[filt].value flux = flux[filt] - #popt, pcov = curve_fit(sine, phase, flux[filt], maxfev=100000) - singleFit = fitSingleSine(phase, flux) - polyFit = fitPolynomial(phase, flux, 10) + polyDegree = 7 + fitThreshold = 0.0 + smoothed_flux = np.convolve(flux, np.ones(len(flux)//100)/(len(flux)//100), mode="valid") + peaks, _ = find_peaks(smoothed_flux, height=np.mean(smoothed_flux)) + #print("Num peaks: ", peaks) - singleRSS = compureRSS(singleFit, flux) - polyRSS = compureRSS(polyFit, flux) + try: + retFit = fitSingleSine(phase, flux) + r2 = 1 - (compureRSS(retFit, flux)/computeTotalSoS(flux)) + fitType = "sine" + #print("R2: ", r2) + if(r2 < fitThreshold): + raise Exception("Sine fit is suboptimal") + #print("Single sine preferred") + except: + retFit = fitPolynomial(phase, flux, polyDegree) + fitType = "poly" + #print("Polynomial preferred") - singleAIC = computeAIC(singleRSS, 4, len(flux)) - singleBIC = computeBIC(singleRSS, 4, len(flux)) - polyAIC = computeAIC(polyRSS, 10, len(flux)) - polyBIC = computeBIC(polyRSS, 10, len(flux)) - - if (singleAIC < polyAIC and singleBIC < polyBIC): - print("Single sine preferred") - sineFit = singleFit - else: - print("Polynomial preferred") - sineFit = polyFit - - sineFit *= multi - return phase, sineFit + retFit *= multi + return phase, retFit, fitType def getFoldedFitPeakValley(sineFit): - return np.argmin(sineFit), np.argmax(sineFit) + minima = argrelextrema(sineFit, np.less)[0] + maxima = argrelextrema(sineFit, np.greater)[0] + #print(minima) + #print(maxima) + if(len(minima) > 0 and len(maxima) > 0): + if(minima[0] < maxima[0] and minima[-1] > maxima[-1]): + if(sineFit[minima[0]] > sineFit[minima[-1]]): + minima = np.delete(minima, 0) + else: + minima = np.delete(minima, -1) + elif(maxima[0] < minima[0] and maxima[-1] > minima[-1]): + if(sineFit[maxima[0]] > sineFit[maxima[-1]]): + maxima = np.delete(maxima, -1) + else: + maxima = np.delete(maxima, 0) + + return minima, maxima def findNearestIndexOfValue(array, value): array = np.asarray(array) @@ -144,57 +160,76 @@ def findNearestIndexOfValue(array, value): return idx def getPhaseRangesNearPeak(maxArgs, phase): - minPhase = phase[maxArgs[0]] - maxPhase = phase[maxArgs[1]] + minPhases = phase[maxArgs[0]] + maxPhases = phase[maxArgs[1]] totalPhase = abs(phase[0]) + abs(phase[-1]) - phasePart = totalPhase * 0.3 / 2 - minPhaseBounds = [] - if(minPhase - phasePart < phase[0]): - minPhaseBounds.append((phase[0], minPhase + phasePart)) - minPhaseBounds.append((phase[-1] + (minPhase - phasePart - phase[0]), phase[-1])) - elif(minPhase + phasePart > phase[-1]): - minPhaseBounds.append((minPhase - phasePart, phase[-1])) - minPhaseBounds.append((phase[0], phase[0] + (minPhase + phasePart - phase[-1]))) - else: - minPhaseBounds.append((minPhase - phasePart, minPhase + phasePart)) + phasePart = totalPhase * 0.3 / (len(minPhases) + len(maxPhases)) + minPhasesBounds = [] + for minPhase in minPhases: + minPhaseBounds = [] + if(minPhase - phasePart < phase[0]): + minPhaseBounds.append((phase[0], minPhase + phasePart)) + minPhaseBounds.append((phase[-1] + (minPhase - phasePart - phase[0]), phase[-1])) + elif(minPhase + phasePart > phase[-1]): + minPhaseBounds.append((minPhase - phasePart, phase[-1])) + minPhaseBounds.append((phase[0], phase[0] + (minPhase + phasePart - phase[-1]))) + else: + minPhaseBounds.append((minPhase - phasePart, minPhase + phasePart)) + minPhasesBounds.append(minPhaseBounds) - maxPhaseBounds = [] - if(maxPhase - phasePart < phase[0]): - maxPhaseBounds.append((phase[0], maxPhase + phasePart)) - maxPhaseBounds.append((phase[-1] + (maxPhase - phasePart - phase[0]), phase[-1])) - elif(maxPhase + phasePart > phase[-1]): - maxPhaseBounds.append((maxPhase - phasePart, phase[-1])) - maxPhaseBounds.append((phase[0], phase[0] + (maxPhase + phasePart - phase[-1]))) - else: - maxPhaseBounds.append((maxPhase - phasePart, maxPhase + phasePart)) + maxPhasesBounds = [] + for maxPhase in maxPhases: + maxPhaseBounds = [] + if(maxPhase - phasePart < phase[0]): + maxPhaseBounds.append((phase[0], maxPhase + phasePart)) + maxPhaseBounds.append((phase[-1] + (maxPhase - phasePart - phase[0]), phase[-1])) + elif(maxPhase + phasePart > phase[-1]): + maxPhaseBounds.append((maxPhase - phasePart, phase[-1])) + maxPhaseBounds.append((phase[0], phase[0] + (maxPhase + phasePart - phase[-1]))) + else: + maxPhaseBounds.append((maxPhase - phasePart, maxPhase + phasePart)) + maxPhasesBounds.append(maxPhaseBounds) - minPhaseBoundsIndices = [] - for l in minPhaseBounds: - minPhaseBoundsIndices.append([findNearestIndexOfValue(phase, lv) for lv in l]) + minPhasesBoundsIndices = [] + for minPhaseBounds in minPhasesBounds: + minPhaseBoundsIndices = [] + for l in minPhaseBounds: + minPhaseBoundsIndices.append([findNearestIndexOfValue(phase, lv) for lv in l]) + minPhasesBoundsIndices.append(minPhaseBoundsIndices) - maxPhaseBoundsIndices = [] - for l in maxPhaseBounds: - maxPhaseBoundsIndices.append([findNearestIndexOfValue(phase, lv) for lv in l]) + maxPhasesBoundsIndices = [] + for maxPhaseBounds in maxPhasesBounds: + maxPhaseBoundsIndices = [] + for l in maxPhaseBounds: + maxPhaseBoundsIndices.append([findNearestIndexOfValue(phase, lv) for lv in l]) + maxPhasesBoundsIndices.append(maxPhaseBoundsIndices) - return minPhaseBoundsIndices, maxPhaseBoundsIndices + return minPhasesBoundsIndices, maxPhasesBoundsIndices def plotPhaseRangesNearPeak(indices, phase, ax=None): - minPhaseBoundsIndices = indices[0] - maxPhaseBoundsIndices = indices[1] + minPhasesBoundsIndices = indices[0] + maxPhasesBoundsIndices = indices[1] if(ax is None): import matplotlib.pyplot as plt - for pair in minPhaseBoundsIndices: - for l in pair: - plt.axes.axvline(phase[l], color="violet") - for pair in maxPhaseBoundsIndices: - for l in pair: - plt.axes.axvline(phase[l], color="orange") + for minPhaseBoundsIndices in minPhasesBoundsIndices: + for pair in minPhaseBoundsIndices: + for l in pair: + plt.axes.axvline(phase[l], color="violet") + for maxPhaseBoundsIndices in maxPhasesBoundsIndices: + for pair in maxPhaseBoundsIndices: + for l in pair: + plt.axes.axvline(phase[l], color="orange") else: - for pair in minPhaseBoundsIndices: - for l in pair: - ax.axvline(phase[l], color="violet") - for pair in maxPhaseBoundsIndices: - for l in pair: - ax.axvline(phase[l], color="orange") \ No newline at end of file + for minPhaseBoundsIndices in minPhasesBoundsIndices: + for pair in minPhaseBoundsIndices: + for l in pair: + ax.axvline(phase[l], color="violet") + for maxPhaseBoundsIndices in maxPhasesBoundsIndices: + for pair in maxPhaseBoundsIndices: + for l in pair: + ax.axvline(phase[l], color="orange") + +def getEpochTime(lc): + return lc.time.value[argrelextrema(np.asarray(lc.flux.value), np.less, order=500)[0]][0] \ No newline at end of file diff --git a/main/util/MinimalLightCurve.py b/main/util/MinimalLightCurve.py index 9df8ff5..1a6702c 100644 --- a/main/util/MinimalLightCurve.py +++ b/main/util/MinimalLightCurve.py @@ -17,6 +17,10 @@ from astropy.time.formats import TimeFromEpoch log = logging.getLogger(__name__) +class LightkurveWarning(Warning): + """Class for all Lightkurve warnings.""" + pass + class TimeBKJD(TimeFromEpoch): name = 'bkjd' unit = 1.0 @@ -402,6 +406,14 @@ class TessQualityFlags(QualityFlags): 32768: "Insufficient Targets for Error Correction Exclude", } +def validate_method(method, supported_methods): + method = method.lower() + if method in supported_methods: + return method + raise ValueError( + "method '{}' is not supported; " + "must be one of {}".format(method, supported_methods) + ) def _is_dict_like(data1): return hasattr(data1, "keys") and callable(getattr(data1, "keys")) @@ -879,6 +891,195 @@ class LightCurve(TimeSeries): return flatten_lc, trend_lc return flatten_lc + def to_periodogram(self, method="lombscargle", **kwargs): + supported_methods = ["ls", "bls", "lombscargle", "boxleastsquares"] + method = validate_method(method.replace(" ", ""), supported_methods) + if method in ["bls", "boxleastsquares"]: + from .MinimalPeriodogram import BoxLeastSquaresPeriodogram + + return BoxLeastSquaresPeriodogram.from_lightcurve(lc=self, **kwargs) + else: + from .MinimalPeriodogram import LombScarglePeriodogram + + return LombScarglePeriodogram.from_lightcurve(lc=self, **kwargs) + + def fold( + self, + period=None, + epoch_time=None, + epoch_phase=0, + wrap_phase=None, + normalize_phase=False, + ): + # Lightkurve v1.x assumed that `period` was given in days if no unit + # was specified. We maintain this behavior for backwards-compatibility. + if period is not None and not isinstance(period, Quantity): + period *= u.day + if epoch_time is not None and not isinstance(epoch_time, Time): + epoch_time = Time( + epoch_time, format=self.time.format, scale=self.time.scale + ) + if ( + epoch_phase is not None + and not isinstance(epoch_phase, Quantity) + and not normalize_phase + ): + epoch_phase *= u.day + if wrap_phase is not None and not isinstance(wrap_phase, Quantity): + wrap_phase *= u.day + + # Warn if `epoch_time` appears to use the wrong format + if epoch_time is not None and epoch_time.value > 2450000: + if self.time.format == "bkjd": + warnings.warn( + "`epoch_time` appears to be given in JD, " + "however the light curve time uses BKJD " + "(i.e. JD - 2454833).", + LightkurveWarning, + ) + elif self.time.format == "btjd": + warnings.warn( + "`epoch_time` appears to be given in JD, " + "however the light curve time uses BTJD " + "(i.e. JD - 2457000).", + LightkurveWarning, + ) + + ts = super().fold( + period=period, + epoch_time=epoch_time, + epoch_phase=epoch_phase, + wrap_phase=wrap_phase, + normalize_phase=normalize_phase, + ) + + # The folded time would pass the `TimeSeries` validation check if + # `normalize_phase=True`, so creating a `FoldedLightCurve` object + # requires the following three-step workaround: + # 1. Give the folded light curve a valid time column again + with ts._delay_required_column_checks(): + folded_time = ts.time.copy() + ts.remove_column("time") + ts.add_column(self.time, name="time", index=0) + # 2. Create the folded object + lc = FoldedLightCurve(data=ts) + # 3. Restore the folded time + with lc._delay_required_column_checks(): + lc.remove_column("time") + lc.add_column(folded_time, name="time", index=0) + + # Add extra column and meta data specific to FoldedLightCurve + lc.add_column( + self.time.copy(), name="time_original", index=len(self._required_columns) + ) + lc.meta["PERIOD"] = period + lc.meta["EPOCH_TIME"] = epoch_time + lc.meta["EPOCH_PHASE"] = epoch_phase + lc.meta["WRAP_PHASE"] = wrap_phase + lc.meta["NORMALIZE_PHASE"] = normalize_phase + lc.sort("time") + + return lc + + def normalize(self, unit="unscaled"): + validate_method(unit, ["unscaled", "percent", "ppt", "ppm"]) + median_flux = np.nanmedian(self.flux) + std_flux = np.nanstd(self.flux) + + # If the median flux is within half a standard deviation from zero, the + # light curve is likely zero-centered and normalization makes no sense. + if (median_flux == 0) or ( + np.isfinite(std_flux) and (np.abs(median_flux) < 0.5 * std_flux) + ): + warnings.warn( + "The light curve appears to be zero-centered " + "(median={:.2e} +/- {:.2e}); `normalize()` will divide " + "the light curve by a value close to zero, which is " + "probably not what you want." + "".format(median_flux, std_flux), + LightkurveWarning, + ) + # If the median flux is negative, normalization will invert the light + # curve and makes no sense. + if median_flux < 0: + warnings.warn( + "The light curve has a negative median flux ({:.2e});" + " `normalize()` will therefore divide by a negative " + "number and invert the light curve, which is probably" + "not what you want".format(median_flux), + LightkurveWarning, + ) + + # Create a new light curve instance and normalize its values + lc = self.copy() + lc.flux = lc.flux / median_flux + lc.flux_err = lc.flux_err / median_flux + if not lc.flux.unit: + lc.flux *= u.dimensionless_unscaled + if not lc.flux_err.unit: + lc.flux_err *= u.dimensionless_unscaled + + # Set the desired relative (dimensionless) units + if unit == "percent": + lc.flux = lc.flux.to(u.percent) + lc.flux_err = lc.flux_err.to(u.percent) + elif unit in ("ppt", "ppm"): + lc.flux = lc.flux.to(unit) + lc.flux_err = lc.flux_err.to(unit) + + lc.meta["NORMALIZED"] = True + return lc + + def remove_nans(self, column: str = "flux"): + return self[~np.isnan(self[column])] # This will return a sliced copy + +class FoldedLightCurve(LightCurve): + + @property + def phase(self): + """Alias for `LightCurve.time`.""" + return self.time + + @property + def cycle(self): + cycle_epoch_start = self.epoch_time - self.period / 2 + result = np.asarray(np.floor(((self.time_original - cycle_epoch_start) / self.period).value), dtype=int) + result = result - result.min() + return result + + @property + def odd_mask(self): + return self.cycle % 2 == 1 + + @property + def even_mask(self): + return ~self.odd_mask + + def _set_xlabel(self, kwargs): + if "xlabel" not in kwargs: + kwargs["xlabel"] = "Phase" + if isinstance(self.time, TimeDelta): + kwargs["xlabel"] += f" [{self.time.format.upper()}]" + return kwargs + + def plot(self, **kwargs): + kwargs = self._set_xlabel(kwargs) + return super(FoldedLightCurve, self).plot(**kwargs) + + def scatter(self, **kwargs): + kwargs = self._set_xlabel(kwargs) + return super(FoldedLightCurve, self).scatter(**kwargs) + + def errorbar(self, **kwargs): + kwargs = self._set_xlabel(kwargs) + return super(FoldedLightCurve, self).errorbar(**kwargs) + + def plot_river(self, **kwargs): + ax = super(FoldedLightCurve, self).plot_river( + period=self.period, epoch_time=self.epoch_time, **kwargs + ) + return ax + class KeplerLightCurve(LightCurve): """Subclass of :class:`LightCurve ` to represent data from NASA's Kepler and K2 mission.""" @@ -1418,3 +1619,5 @@ try: except registry.IORegistryError as exc: print(exc) pass # necessary to enable autoreload during debugging + + diff --git a/main/util/MinimalPeriodogram.py b/main/util/MinimalPeriodogram.py new file mode 100644 index 0000000..19bfa93 --- /dev/null +++ b/main/util/MinimalPeriodogram.py @@ -0,0 +1,1103 @@ +"""Defines the Periodogram class and associated tools.""" +from __future__ import division, print_function + +import copy +import logging +import math +import re +import warnings + +import numpy as np +from matplotlib import pyplot as plt + +import astropy +from astropy.table import Table +from astropy import units as u +from astropy.units import cds +from astropy.convolution import convolve, Box1DKernel +from astropy.time import Time + +from astropy.timeseries import LombScargle +from astropy.timeseries.periodograms.lombscargle import implementations # for .main._is_regular + +from .MinimalLightCurve import LightkurveWarning, validate_method +from .MinimalLightCurve import LightCurve + +log = logging.getLogger(__name__) + +__all__ = ["Periodogram", "LombScarglePeriodogram", "BoxLeastSquaresPeriodogram"] + + +class Periodogram(object): + """Generic class to represent a power spectrum (frequency vs power data). + + The Periodogram class represents a power spectrum, with values of + frequency on the x-axis (in any frequency units) and values of power on the + y-axis (in units of flux^2 / [frequency units]). + + Attributes + ---------- + frequency : `~astropy.units.Quantity` + Array of frequencies as an AstroPy Quantity object. + power : `~astropy.units.Quantity` + Array of power-spectral-densities. The Quantity must have units of + `flux^2 / freq_unit`, where freq_unit is the unit of the frequency + attribute. + nyquist : float + The Nyquist frequency of the lightcurve. In units of freq_unit, where + freq_unit is the unit of the frequency attribute. + label : str + Human-friendly object label, e.g. "KIC 123456789". + targetid : str + Identifier of the target. + default_view : "frequency" or "period" + Should plots be shown in frequency space or period space by default? + meta : dict + Free-form metadata associated with the Periodogram. + """ + + frequency = None + """The array of frequency values.""" + + power = None + """The array of power values.""" + + def __init__( + self, + frequency, + power, + nyquist=None, + label=None, + targetid=None, + default_view="frequency", + meta={}, + ): + # Input validation + if not isinstance(frequency, u.quantity.Quantity): + raise ValueError("frequency must be an `astropy.units.Quantity` object.") + if not isinstance(power, u.quantity.Quantity): + raise ValueError("power must be an `astropy.units.Quantity` object.") + # Frequency must have frequency units + try: + frequency.to(u.Hz) + except u.UnitConversionError: + raise ValueError("Frequency must be in units of 1/time.") + # Frequency and power must have sensible shapes + if frequency.shape[0] <= 1: + raise ValueError("frequency and power must have a length greater than 1.") + if frequency.shape != power.shape: + raise ValueError("frequency and power must have the same length.") + + self.frequency = frequency + self.power = power + self.nyquist = nyquist + self.label = label + self.targetid = targetid + self.default_view = self._validate_view(default_view) + self.meta = meta + + def _validate_view(self, view): + """Verifies whether `view` is is one of {"frequency", "period"} and + raises a helpful `ValueError` if not. + """ + if view is None and hasattr(self, "default_view"): + view = self.default_view + return validate_method(view, ["frequency", "period"]) + + def _is_evenly_spaced(self): + """Returns true if the values in ``frequency`` are evenly spaced. + + This helper method exists because some features, such as ``smooth()``, + ``estimate_numax()``, and ``estimate_deltanu()``, require a grid of + evenly-spaced frequencies. + """ + # verify that the first differences are all equal + freqdiff = np.diff(self.frequency.value) + if np.allclose(freqdiff[0], freqdiff): + return True + return False + + @property + def period(self): + """The array of periods, i.e. 1/frequency.""" + return 1.0 / self.frequency + + @property + def max_power(self): + """Power of the highest peak in the periodogram.""" + return np.nanmax(self.power) + + @property + def frequency_at_max_power(self): + """Frequency value corresponding to the highest peak in the periodogram.""" + return self.frequency[np.nanargmax(self.power)] + + @property + def period_at_max_power(self): + """Period value corresponding to the highest peak in the periodogram.""" + return 1.0 / self.frequency_at_max_power + + def bin(self, binsize=10, method="mean"): + """Bins the power spectrum. + + Parameters + ---------- + binsize : int + The factor by which to bin the power spectrum, in the sense that + the power spectrum will be smoothed by taking the mean in bins + of size N / binsize, where N is the length of the original + frequency array. Defaults to 10. + method : str, one of 'mean' or 'median' + Method to use for binning. Default is 'mean'. + + Returns + ------- + binned_periodogram : a `Periodogram` object + Returns a new `Periodogram` object which has been binned. + """ + # Input validation + if binsize < 1: + raise ValueError("binsize must be larger than or equal to 1") + method = validate_method(method, ["mean", "median"]) + + m = int(len(self.power) / binsize) # length of the binned arrays + if method == "mean": + binned_freq = self.frequency[: m * binsize].reshape((m, binsize)).mean(1) + binned_power = self.power[: m * binsize].reshape((m, binsize)).mean(1) + elif method == "median": + binned_freq = np.nanmedian( + self.frequency[: m * binsize].reshape((m, binsize)), axis=1 + ) + binned_power = np.nanmedian( + self.power[: m * binsize].reshape((m, binsize)), axis=1 + ) + + binned_pg = self.copy() + binned_pg.frequency = binned_freq + binned_pg.power = binned_power + return binned_pg + + def smooth(self, method="boxkernel", filter_width=0.1): + """Smooths the power spectrum using the 'boxkernel' or 'logmedian' method. + + If `method` is set to 'boxkernel', this method will smooth the power + spectrum by convolving with a numpy Box1DKernel with a width of + `filter_width`, where `filter width` is in units of frequency. + This is best for filtering out noise while maintaining seismic mode + peaks. This method requires the Periodogram to have an evenly spaced + grid of frequencies. A `ValueError` exception will be raised if this is + not the case. + + If `method` is set to 'logmedian', it smooths the power spectrum using + a moving median which moves across the power spectrum in a steps of + + log10(x0) + 0.5 * filter_width + + where `filter width` is in log10(frequency) space. This is best for + estimating the noise background, as it filters over the seismic peaks. + + Periodograms that are unsmoothed have multiplicative noise that is + distributed as chi squared 2 degrees of freedom. This noise + distribution has a well defined mean and median but the two are not + equivalent. The mean of a chi squared 2 dof distribution is 2, but the + median is 2(8/9)**3. + (see https://en.wikipedia.org/wiki/Chi-squared_distribution) + In order to maintain consistency between 'boxkernel' and 'logmedian' a + correction factor of (8/9)**3 is applied to (i.e., the median is divided + by the factor) to the median values. + + In addition to consistency with the 'boxkernel' method, the correction + of the median values is useful when applying the periodogram flatten + method. The flatten method divides the periodgram by the smoothed + periodogram using the 'logmedian' method. By appyling the correction + factor we follow asteroseismic convention that the signal-to-noise + power has a mean value of unity. (note the signal-to-noise power is + really the signal plus noise divided by the noise and hence should be + unity in the absence of any signal) + + Parameters + ---------- + method : str, one of 'boxkernel' or 'logmedian' + The smoothing method to use. Defaults to 'boxkernel'. + filter_width : float + If `method` = 'boxkernel', this is the width of the smoothing filter + in units of frequency. + If method = `logmedian`, this is the width of the smoothing filter + in log10(frequency) space. + + Returns + ------- + smoothed_pg : `Periodogram` object + Returns a new `Periodogram` object in which the power spectrum + has been smoothed. + """ + method = validate_method(method, ["boxkernel", "logmedian"]) + + if method == "boxkernel": + if filter_width <= 0.0: + raise ValueError( + "the `filter_width` parameter must be " + "larger than 0 for the 'boxkernel' method." + ) + try: + filter_width = u.Quantity(filter_width, self.frequency.unit) + except u.UnitConversionError: + raise ValueError( + "the `filter_width` parameter must have " "frequency units." + ) + + # Check to see if we have a grid of evenly spaced periods instead. + if not self._is_evenly_spaced(): + raise ValueError( + "the 'boxkernel' method requires the periodogram " + "to have a grid of evenly spaced frequencies." + ) + + fs = np.mean(np.diff(self.frequency)) + box_kernel = Box1DKernel(math.ceil((filter_width / fs).value)) + smooth_power = convolve(self.power.value, box_kernel) + smooth_pg = self.copy() + smooth_pg.power = u.Quantity(smooth_power, self.power.unit) + return smooth_pg + + if method == "logmedian": + if isinstance(filter_width, astropy.units.quantity.Quantity): + raise ValueError( + "the 'logmedian' method requires a dimensionless " + "value for `filter_width` in log10(frequency) space." + ) + count = np.zeros(len(self.frequency.value), dtype=int) + bkg = np.zeros_like(self.frequency.value) + x0 = np.log10(self.frequency[0].value) + corr_factor = (8.0 / 9.0) ** 3 + while x0 < np.log10(self.frequency[-1].value): + m = np.abs(np.log10(self.frequency.value) - x0) < filter_width + if len(bkg[m] > 0): + bkg[m] += np.nanmedian(self.power[m].value) / corr_factor + count[m] += 1 + x0 += 0.5 * filter_width + bkg /= count + smooth_pg = self.copy() + smooth_pg.power = u.Quantity(bkg, self.power.unit) + return smooth_pg + + def flatten(self, method="logmedian", filter_width=0.01, return_trend=False): + """Estimates the Signal-To-Noise (SNR) spectrum by dividing out an + estimate of the noise background. + + This method divides the power spectrum by a background estimated + using a moving filter in log10 space by default. For details on the + `method` and `filter_width` parameters, see `Periodogram.smooth()` + + Dividing the power through by the noise background produces a spectrum + with no units of power. Since the signal is divided through by a measure + of the noise, we refer to this as a `Signal-To-Noise` spectrum. + + Parameters + ---------- + method : str, one of 'boxkernel' or 'logmedian' + Background estimation method passed on to `Periodogram.smooth()`. + Defaults to 'logmedian'. + filter_width : float + If `method` = 'boxkernel', this is the width of the smoothing filter + in units of frequency. + If method = `logmedian`, this is the width of the smoothing filter + in log10(frequency) space. + return_trend : bool + If True, then the background estimate, alongside the SNR spectrum, + will be returned. + + Returns + ------- + snr_spectrum : `Periodogram` object + Returns a periodogram object where the power is an estimate of the + signal-to-noise of the spectrum, creating by dividing the powers + with a simple estimate of the noise background using a smoothing filter. + bkg : `Periodogram` object + The estimated power spectrum of the background noise. This is only + returned if `return_trend = True`. + """ + bkg = self.smooth(method=method, filter_width=filter_width) + snr_pg = self / bkg.power + snr = SNRPeriodogram( + snr_pg.frequency, + snr_pg.power, + nyquist=self.nyquist, + targetid=self.targetid, + label=self.label, + meta=self.meta, + ) + if return_trend: + return snr, bkg + return snr + + def to_table(self): + """Exports the Periodogram as an Astropy Table. + + Returns + ------- + table : `~astropy.table.Table` object + An AstroPy Table with columns 'frequency', 'period', and 'power'. + """ + return Table( + data=(self.frequency, self.period, self.power), + names=("frequency", "period", "power"), + meta=self.meta, + ) + + def copy(self): + """Returns a copy of the Periodogram object. + + This method uses the `copy.deepcopy` function to ensure that all + objects stored within the Periodogram are copied. + + Returns + ------- + pg_copy : Periodogram + A new `Periodogram` object which is a copy of the original. + """ + return copy.deepcopy(self) + + def __repr__(self): + return "Periodogram(ID: {})".format(self.label) + + def __getitem__(self, key): + copy_self = self.copy() + copy_self.frequency = self.frequency[key] + copy_self.power = self.power[key] + return copy_self + + def __add__(self, other): + copy_self = self.copy() + copy_self.power = copy_self.power + u.Quantity(other, self.power.unit) + return copy_self + + def __radd__(self, other): + return self.__add__(other) + + def __sub__(self, other): + return self.__add__(-other) + + def __rsub__(self, other): + copy_self = self.copy() + copy_self.power = other - copy_self.power + return copy_self + + def __mul__(self, other): + copy_self = self.copy() + copy_self.power = other * copy_self.power + return copy_self + + def __rmul__(self, other): + return self.__mul__(other) + + def __truediv__(self, other): + return self.__mul__(1.0 / other) + + def __rtruediv__(self, other): + copy_self = self.copy() + copy_self.power = other / copy_self.power + return copy_self + + def __div__(self, other): + return self.__truediv__(other) + + def __rdiv__(self, other): + return self.__rtruediv__(other) + + def show_properties(self): + """Prints a summary of the non-callable attributes of the Periodogram object. + + Prints in order of type (ints, strings, lists, arrays and others). + Prints in alphabetical order. + """ + attrs = {} + for attr in dir(self): + if not attr.startswith("_"): + res = getattr(self, attr) + if callable(res): + continue + + if isinstance(res, astropy.units.quantity.Quantity): + unit = res.unit + res = res.value + attrs[attr] = {"res": res} + attrs[attr]["unit"] = unit.to_string() + else: + attrs[attr] = {"res": res} + attrs[attr]["unit"] = "" + + if attr == "hdu": + attrs[attr] = {"res": res, "type": "list"} + for idx, r in enumerate(res): + if idx == 0: + attrs[attr]["print"] = "{}".format(r.header["EXTNAME"]) + else: + attrs[attr]["print"] = "{}, {}".format( + attrs[attr]["print"], "{}".format(r.header["EXTNAME"]) + ) + continue + + if isinstance(res, int): + attrs[attr]["print"] = "{}".format(res) + attrs[attr]["type"] = "int" + elif isinstance(res, float): + attrs[attr]["print"] = "{}".format(np.round(res, 4)) + attrs[attr]["type"] = "float" + elif isinstance(res, np.ndarray): + attrs[attr]["print"] = "array {}".format(res.shape) + attrs[attr]["type"] = "array" + elif isinstance(res, list): + attrs[attr]["print"] = "list length {}".format(len(res)) + attrs[attr]["type"] = "list" + elif isinstance(res, str): + if res == "": + attrs[attr]["print"] = "{}".format("None") + else: + attrs[attr]["print"] = "{}".format(res) + attrs[attr]["type"] = "str" + elif attr == "wcs": + attrs[attr]["print"] = "astropy.wcs.wcs.WCS" + attrs[attr]["type"] = "other" + else: + attrs[attr]["print"] = "{}".format(type(res)) + attrs[attr]["type"] = "other" + + output = Table( + names=["Attribute", "Description", "Units"], dtype=[object, object, object] + ) + idx = 0 + types = ["int", "str", "float", "list", "array", "other"] + for typ in types: + for attr, dic in attrs.items(): + if dic["type"] == typ: + output.add_row([attr, dic["print"], dic["unit"]]) + idx += 1 + print("lightkurve.Periodogram properties:") + output.pprint(max_lines=-1, max_width=-1) + + def to_seismology(self, **kwargs): + """Returns a `~lightkurve.seismology.Seismology` object to analyze the periodogram. + + Returns + ------- + seismology : `~lightkurve.seismology.Seismology` + Helper object to run asteroseismology methods. + """ + from .seismology import Seismology + + return Seismology(self) + + +class SNRPeriodogram(Periodogram): + """Defines a Signal-to-Noise Ratio (SNR) Periodogram class. + + This class is nearly identical to the standard :class:`Periodogram` class, + but has different plotting defaults. + """ + + def __init__(self, *args, **kwargs): + super(SNRPeriodogram, self).__init__(*args, **kwargs) + + def __repr__(self): + return "SNRPeriodogram(ID: {})".format(self.label) + + def plot(self, **kwargs): + """Plot the SNR spectrum using matplotlib's `plot` method. + See `Periodogram.plot` for details on the accepted arguments. + + Parameters + ---------- + kwargs : dict + Dictionary of arguments ot be passed to `Periodogram.plot`. + + Returns + ------- + ax : `~matplotlib.axes.Axes` + The matplotlib axes object. + """ + ax = super(SNRPeriodogram, self).plot(**kwargs) + if "ylabel" not in kwargs: + ax.set_ylabel("Signal to Noise Ratio (SNR)") + return ax + + +class LombScarglePeriodogram(Periodogram): + """Subclass of :class:`Periodogram ` + representing a power spectrum generated using the Lomb Scargle method. + """ + + def __init__(self, *args, **kwargs): + self._LS_object = kwargs.pop("ls_obj", None) + self.nterms = kwargs.pop("nterms", 1) + self.ls_method = kwargs.pop("ls_method", "fastchi2") + super(LombScarglePeriodogram, self).__init__(*args, **kwargs) + + def __repr__(self): + return "LombScarglePeriodogram(ID: {})".format(self.label) + + @staticmethod + def from_lightcurve( + lc, + minimum_frequency=None, + maximum_frequency=None, + minimum_period=None, + maximum_period=None, + frequency=None, + period=None, + nterms=1, + nyquist_factor=1, + oversample_factor=None, + freq_unit=None, + normalization="amplitude", + ls_method="fast", + **kwargs + ): + # Input validation + normalization = validate_method(normalization, ["psd", "amplitude"]) + if np.isnan(lc.flux).any() or (hasattr(lc.flux, 'unmasked') and np.isnan(lc.flux.unmasked).any()): + lc = lc.remove_nans() + log.debug( + "Lightcurve contains NaN values." + "These are removed before creating the periodogram." + ) + + # Setting default frequency units + if freq_unit is None: + freq_unit = 1 / u.day if normalization == "amplitude" else u.microhertz + + # Default oversample factor + if oversample_factor is None: + oversample_factor = 5.0 if normalization == "amplitude" else 1.0 + + if "min_period" in kwargs: + warnings.warn( + "`min_period` keyword is deprecated, " + "please use `minimum_period` instead.", + LightkurveWarning, + ) + minimum_period = kwargs.pop("min_period", None) + if "max_period" in kwargs: + warnings.warn( + "`max_period` keyword is deprecated, " + "please use `maximum_period` instead.", + LightkurveWarning, + ) + maximum_period = kwargs.pop("max_period", None) + if "min_frequency" in kwargs: + warnings.warn( + "`min_frequency` keyword is deprecated, " + "please use `minimum_frequency` instead.", + LightkurveWarning, + ) + minimum_frequency = kwargs.pop("min_frequency", None) + if "max_frequency" in kwargs: + warnings.warn( + "`max_frequency` keyword is deprecated, " + "please use `maximum_frequency` instead.", + LightkurveWarning, + ) + maximum_frequency = kwargs.pop("max_frequency", None) + + # Check if any values of period have been passed and set format accordingly + if not all(b is None for b in [period, minimum_period, maximum_period]): + default_view = "period" + else: + default_view = "frequency" + + # If period and frequency keywords have both been set, throw an error + if (not all(b is None for b in [period, minimum_period, maximum_period])) & ( + not all( + b is None for b in [frequency, minimum_frequency, maximum_frequency] + ) + ): + raise ValueError( + "You have input keyword arguments for both frequency and period. " + "Please only use one." + ) + + time = lc.time.copy() + + # Approximate Nyquist Frequency and frequency bin width in terms of days + nyquist = 0.5 * (1.0 / (np.median(np.diff(time.value)))) * (1 / cds.d) + fs = (1.0 / (time[-1] - time[0])) / oversample_factor + + # Convert these values to requested frequency unit + nyquist = nyquist.to(freq_unit) + fs = fs.to(freq_unit) + + # Warn if there is confusing input + if (frequency is not None) & ( + any([a is not None for a in [minimum_frequency, maximum_frequency]]) + ): + log.warning( + "You have passed both a grid of frequencies " + "and min_frequency/maximum_frequency arguments; " + "the latter will be ignored." + ) + if (period is not None) & ( + any([a is not None for a in [minimum_period, maximum_period]]) + ): + log.warning( + "You have passed a grid of periods " + "and minimum_period/maximum_period arguments; " + "the latter will be ignored." + ) + + # Tidy up the period stuff... + if maximum_period is not None: + # minimum_frequency MUST be none by this point. + minimum_frequency = 1.0 / maximum_period + if minimum_period is not None: + # maximum_frequency MUST be none by this point. + maximum_frequency = 1.0 / minimum_period + # If the user specified a period, copy it into the frequency. + if period is not None: + frequency = 1.0 / period + + # Do unit conversions if user input min/max frequency or period + if frequency is None: + if minimum_frequency is not None: + minimum_frequency = u.Quantity(minimum_frequency, freq_unit) + if maximum_frequency is not None: + maximum_frequency = u.Quantity(maximum_frequency, freq_unit) + if (minimum_frequency is not None) & (maximum_frequency is not None): + if minimum_frequency > maximum_frequency: + if default_view == "frequency": + raise ValueError( + "minimum_frequency cannot be larger than maximum_frequency" + ) + if default_view == "period": + raise ValueError( + "minimum_period cannot be larger than maximum_period" + ) + # If nothing has been passed in, set them to the defaults + if minimum_frequency is None: + minimum_frequency = fs + if maximum_frequency is None: + maximum_frequency = nyquist * nyquist_factor + + # Create frequency grid evenly spaced in frequency + frequency = np.arange( + minimum_frequency.value, maximum_frequency.value, fs.value + ) + + # Convert to desired units + frequency = u.Quantity(frequency, freq_unit) + + # Change to compatible ls method if sampling not even in frequency + if not implementations.main._is_regular(frequency) and ls_method in [ + "fastchi2", + "fast", + ]: + oldmethod = ls_method + ls_method = {"fastchi2": "chi2", "fast": "slow"}[ls_method] + log.warning( + "The requested periodogram is not evenly sampled in frequency.\n" + "Method has been changed from '{}' to '{}' to allow for this.".format( + oldmethod, ls_method + ) + ) + + if (nterms > 1) and (ls_method not in ["fastchi2", "chi2"]): + warnings.warn( + "Building a Lomb Scargle Periodogram using the `slow` method. " + "`nterms` has been set to >1, however this is not supported under the `{}` method. " + "To run with higher nterms, set `ls_method` to either 'fastchi2', or 'chi2'. " + "Please refer to the `astropy.timeseries.periodogram.LombScargle` documentation.".format( + ls_method + ), + LightkurveWarning, + ) + nterms = 1 + + if float(astropy.__version__[0]) >= 3: + LS = LombScargle( + time, lc.flux, nterms=nterms, normalization="psd", **kwargs + ) + power = LS.power(frequency, method=ls_method) + else: + LS = LombScargle(time, lc.flux, nterms=nterms, **kwargs) + power = LS.power(frequency, method=ls_method, normalization="psd") + + if normalization == "psd": # Power spectral density + # Rescale from the unnormalized power output by Astropy's + # Lomb-Scargle function to units of flux_variance / [frequency unit] + # that may be of more interest for asteroseismology. + power *= 2.0 / (len(time) * oversample_factor * fs) + elif normalization == "amplitude": + power = np.sqrt(power) * np.sqrt(4.0 / len(lc.time)) + + # Periodogram needs properties + return LombScarglePeriodogram( + frequency=frequency, + power=power, + nyquist=nyquist, + targetid=lc.meta.get("TARGETID"), + label=lc.meta.get("LABEL"), + default_view=default_view, + ls_obj=LS, + nterms=nterms, + ls_method=ls_method, + meta=lc.meta, + ) + + def model(self, time, frequency=None): + """Obtain the flux model for a given frequency and time + + Parameters + ---------- + time : np.ndarray + Time points to evaluate model. + frequency : frequency to evaluate model. Default is the frequency at + max power. + + Returns + ------- + result : lightkurve.LightCurve + Model object with the time and flux model + """ + if self._LS_object is None: + raise ValueError("No `astropy` Lomb Scargle object exists.") + if frequency is None: + frequency = self.frequency_at_max_power + f = self._LS_object.model(time, frequency) + lc = LightCurve( + time=time, + flux=f, + meta={"FREQUENCY": frequency}, + label="LS Model", + targetid="{} LS Model".format(self.targetid), + ) + return lc.normalize() + + +class BoxLeastSquaresPeriodogram(Periodogram): + """Subclass of :class:`Periodogram ` + representing a power spectrum generated using the Box Least Squares (BLS) method. + """ + + def __init__(self, *args, **kwargs): + self.duration = kwargs.pop("duration", None) + self.depth = kwargs.pop("depth", None) + self.snr = kwargs.pop("snr", None) + self._BLS_result = kwargs.pop("bls_result", None) + self._BLS_object = kwargs.pop("bls_obj", None) + + self.transit_time = kwargs.pop("transit_time", None) + self.time = kwargs.pop("time", None) + self.flux = kwargs.pop("flux", None) + self.time_unit = kwargs.pop("time_unit", None) + super(BoxLeastSquaresPeriodogram, self).__init__(*args, **kwargs) + + def __repr__(self): + return "BoxLeastSquaresPeriodogram(ID: {})".format(self.label) + + @staticmethod + def from_lightcurve(lc, **kwargs): + """Creates a `Periodogram` from a LightCurve using the Box Least Squares (BLS) method. + + Parameters + ---------- + lc : `LightCurve` object + The LightCurve from which to compute the Periodogram. + duration : float, array_like, or `~astropy.units.Quantity`, optional + The set of durations that will be considered. + Default to `[0.05, 0.10, 0.15, 0.20, 0.25, 0.33]` if not specified. + period : array_like or `~astropy.units.Quantity`, optional + The periods where the Periodogram should be computed. + If not provided, a default will be created using + `BoxLeastSquares.autoperiod() `. + minimum_period, maximum_period : float or `~astropy.units.Quantity`, optional + If ``period`` is not provided, the minimum/maximum periods to search. + The defaults will be computed as described in the notes below. + frequency_factor : float, optional + If ``period`` is not provided, a factor to control the frequency spacing of periods + to be considered. + kwargs : dict + Keyword arguments passed to + `BoxLeastSquares.power() ` + + Returns + ------- + Periodogram : `Periodogram` object + Returns a Periodogram object extracted from the lightcurve. + + Notes + ----- + If ``period`` is not provided, the default minimum period is computed from maximum duration and + the median observation time gap as + + .. code-block:: python + + minimum_period = max(median(diff(lc.time)) * 4, + max(duration) + median(diff(lc.time))) + + The default maximum period is computed as + + .. code-block:: python + + maximum_period = (max(lc.time) - min(lc.time)) / 3 + + ensuring that any systems with at least 3 transits are within the range of searched periods. + + """ + # BoxLeastSquares was added to `astropy.stats` in AstroPy v3.1 and then + # moved to `astropy.timeseries` in v3.2, which makes the import below + # somewhat complicated. + try: + from astropy.timeseries import BoxLeastSquares + except ImportError: + try: + from astropy.stats import BoxLeastSquares + except ImportError: + raise ImportError("BLS requires AstroPy v3.1 or later") + + # Validate user input for `lc` + # (BoxLeastSquares will not work if flux or flux_err contain NaNs) + lc = lc.remove_nans() + if np.isfinite(lc.flux_err).all(): + dy = lc.flux_err + else: + dy = None + + # Validate user input for `duration` + duration = kwargs.pop("duration", [0.05, 0.10, 0.15, 0.20, 0.25, 0.33]) + if duration is not None and ~np.all(np.isfinite(duration)): + raise ValueError( + "`duration` parameter contains illegal nan or inf value(s)" + ) + + # Validate user input for `period` + period = kwargs.pop("period", None) + minimum_period = kwargs.pop("minimum_period", None) + maximum_period = kwargs.pop("maximum_period", None) + if period is not None and ~np.all(np.isfinite(period)): + raise ValueError("`period` parameter contains illegal nan or inf value(s)") + if minimum_period is None: + if period is None: + minimum_period = np.max( + [ + np.median(np.diff(lc.time.value)) * 4, + np.max(duration) + np.median(np.diff(lc.time.value)), + ] + ) + else: + minimum_period = np.min(period) + if maximum_period is None: + if period is None: + maximum_period = (np.max(lc.time.value) - np.min(lc.time.value)) / 3.0 + else: + maximum_period = np.max(period) + + # Validate user input for `time_unit` + time_unit = kwargs.pop("time_unit", "day") + if time_unit not in dir(u): + raise ValueError( + "{} is not a valid value for `time_unit`".format(time_unit) + ) + + # Validate user input for `frequency_factor` + frequency_factor = kwargs.pop("frequency_factor", 10) + df = ( + frequency_factor + * np.min(duration) + / (np.max(lc.time.value) - np.min(lc.time.value)) ** 2 + ) + npoints = int(((1 / minimum_period) - (1 / maximum_period)) / df) + if npoints > 1e7: + raise ValueError( + "`period` contains {} points." + "Periodogram is too large to evaluate. " + "Consider setting `frequency_factor` to a higher value." + "".format(np.round(npoints, 4)) + ) + elif npoints > 1e5: + log.warning( + "`period` contains {} points." + "Periodogram is likely to be large, and slow to evaluate. " + "Consider setting `frequency_factor` to a higher value." + "".format(np.round(npoints, 4)) + ) + + # Create BLS object and run the BLS search + bls = BoxLeastSquares(lc.time, lc.flux, dy) + if period is None: + period = bls.autoperiod( + duration, + minimum_period=minimum_period, + maximum_period=maximum_period, + frequency_factor=frequency_factor, + ) + result = bls.power(period, duration, **kwargs) + if not isinstance(result.period, u.quantity.Quantity): + result.period = u.Quantity(result.period, time_unit) + if not isinstance(result.power, u.quantity.Quantity): + result.power = result.power * u.dimensionless_unscaled + if not isinstance(result.duration, u.quantity.Quantity): + result.duration = u.Quantity(result.duration, time_unit) + + return BoxLeastSquaresPeriodogram( + frequency=1.0 / result.period, + power=result.power, + default_view="period", + label=lc.meta.get("LABEL"), + targetid=lc.meta.get("TARGETID"), + transit_time=result.transit_time, + duration=result.duration, + depth=result.depth, + bls_result=result, + snr=result.depth_snr, + bls_obj=bls, + time=lc.time, + flux=lc.flux, + time_unit=time_unit, + ) + + def compute_stats(self, period=None, duration=None, transit_time=None): + """Computes commonly used vetting statistics for a transit model. + + See `~astropy.timeseries.BoxLeastSquares` docs for further details. + + Parameters + ---------- + period : float or Quantity + Period of the transits. Default is `period_at_max_power` + duration : float or Quantity + Duration of the transits. Default is `duration_at_max_power` + transit_time : float or Quantity + Transit midpoint of the transits. Default is `transit_time_at_max_power` + + Returns + ------- + stats : dict + Dictionary of vetting statistics + """ + if period is None: + period = self.period_at_max_power + log.warning("No period specified. Using period at max power") + if duration is None: + duration = self.duration_at_max_power + log.warning("No duration specified. Using duration at max power") + if transit_time is None: + transit_time = self.transit_time_at_max_power + log.warning("No transit time specified. Using transit time at max power") + if not isinstance(transit_time, Time): + transit_time = Time( + transit_time, format=self.time.format, scale=self.time.scale + ) + + return self._BLS_object.compute_stats( + u.Quantity(period, "d").value, u.Quantity(duration, "d").value, transit_time + ) + + def get_transit_model(self, period=None, duration=None, transit_time=None): + """Computes the transit model using the BLS, returns a lightkurve.LightCurve + + See `~astropy.timeseries.BoxLeastSquares` docs for further details. + + Parameters + ---------- + period : float or Quantity + Period of the transits. Default is `period_at_max_power` + duration : float or Quantity + Duration of the transits. Default is `duration_at_max_power` + transit_time : float or Quantity + Transit midpoint of the transits. Default is `transit_time_at_max_power` + + Returns + ------- + model : lightkurve.LightCurve + Model of transit + """ + from .lightcurve import LightCurve + + if period is None: + period = self.period_at_max_power + log.warning("No period specified. Using period at max power") + if duration is None: + duration = self.duration_at_max_power + log.warning("No duration specified. Using duration at max power") + if transit_time is None: + transit_time = self.transit_time_at_max_power + log.warning("No transit time specified. Using transit time at max power") + if not isinstance(transit_time, Time): + transit_time = Time( + transit_time, format=self.time.format, scale=self.time.scale + ) + + model_flux = self._BLS_object.model( + self.time, + u.Quantity(period, "d").value, + u.Quantity(duration, "d").value, + transit_time, + ) + model = LightCurve(time=self.time, flux=model_flux, label="Transit Model Flux") + return model + + def get_transit_mask(self, period=None, duration=None, transit_time=None): + """Returns a boolean array that is ``True`` during transits and + ``False`` elsewhere. + + Parameters + ---------- + period : float or Quantity + Period of the transits. Default is `period_at_max_power` + duration : float or Quantity + Duration of the transits. Default is `duration_at_max_power` + transit_time : float or Quantity + Transit midpoint of the transits. Default is `transit_time_at_max_power` + + Returns + ------- + transit_mask : np.array of bool + Mask that flags transits. Mask is ``True`` where there are transits. + """ + model = self.get_transit_model( + period=period, duration=duration, transit_time=transit_time + ) + return model.flux != np.median(model.flux) + + @property + def transit_time_at_max_power(self): + """Returns the transit time corresponding to the highest peak in the periodogram.""" + return self.transit_time[np.nanargmax(self.power)] + + @property + def duration_at_max_power(self): + """Returns the duration corresponding to the highest peak in the periodogram.""" + return self.duration[np.nanargmax(self.power)] + + @property + def depth_at_max_power(self): + """Returns the depth corresponding to the highest peak in the periodogram.""" + return self.depth[np.nanargmax(self.power)] + + def plot(self, **kwargs): + """Plot the BoxLeastSquaresPeriodogram spectrum using matplotlib's `plot` method. + See `Periodogram.plot` for details on the accepted arguments. + + Parameters + ---------- + kwargs : dict + Dictionary of arguments ot be passed to `Periodogram.plot`. + + Returns + ------- + ax : `~matplotlib.axes.Axes` + The matplotlib axes object. + """ + ax = super(BoxLeastSquaresPeriodogram, self).plot(**kwargs) + if "ylabel" not in kwargs: + ax.set_ylabel("BLS Power") + return ax + + def flatten(self, **kwargs): + raise NotImplementedError( + "`flatten` is not implemented for `BoxLeastSquaresPeriodogram`." + ) + + def smooth(self, **kwargs): + raise NotImplementedError( + "`smooth` is not implemented for `BoxLeastSquaresPeriodogram`. " + )