generate_plots_MKGF: generate plots individually per star too

generate_plots_MKGF.py

@@ -27,7 +27,7 @@ def mkdir_p(mypath):
             pass
         else: raise
 
-fileName = "datav3.cff"
+fileName = "datav4.cff"
 data = pd.read_pickle(fileName)
 
 binList = [10, 20, 30]
@@ -245,10 +245,10 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel
     for spTyp in spTypes:
         finalData = pd.DataFrame()
         
-        Mfilter = data["SpType"].str.startswith(spTyp)
+        typeFilter = data["SpType"].str.startswith(spTyp)
-        numStars = len(set(data[Mfilter]["StarName"]))
+        numStars = len(set(data[typeFilter]["StarName"]))
-        Mfilter &= showSourceFilter
+        typeFilter &= showSourceFilter
-        finalData = pd.concat([finalData, data[Mfilter]], ignore_index=True)
+        finalData = pd.concat([finalData, data[typeFilter]], ignore_index=True)
 
         pdcsapbinningData = []
         locFolder = f"{folderPath}/{spTyp}/"
@@ -342,7 +342,7 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel
                 figFlarepeakHist, ((axFlarepeakHist)) = plt.subplots(nrows=1, ncols=1)
                 axFlarepeakHist.set_ylabel("Flare peak")
                 axFlarepeakHist.set_xlabel("Phase")
-                H, xedges, yedges = np.histogram2d(xData, yData, bins=bins, range=[[0, 2], [min(yData), maxY]])
+                H, xedges, yedges = np.histogram2d(xData, yData, bins=bins, range=[[0, 2], [0.99, maxY]])
                 cmax = 11
                 H_clipped = np.clip(H, None, cmax)
                 im = axFlarepeakHist.imshow(H_clipped.T, origin='lower', interpolation='nearest',
@@ -371,3 +371,140 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel
         
             plt.savefig(f"{locFolder}/{''.join(spTyp)}-Flarepeaks_maxY-{maxY}.png")
             plt.close()
+
+starDB: StarDB = StarDB.getInstance("stars.db")
+
+for starName in starDB.getAllStars():
+    finalData = pd.DataFrame()
+    starNameR = starName.replace('*', '_star_')
+    nameFilter = data["StarName"] == starName
+    nameFilter &= showSourceFilter
+    finalData = pd.concat([finalData, data[nameFilter]], ignore_index=True)
+    color = ""
+
+    pdcsapbinningData = []
+    locFolder = f"{folderPath}/stars/{starNameR}/"
+    mkdir_p(f"{locFolder}/")
+    csvFile = open(f"{locFolder}/{starNameR}.csv", "a")
+    csvFile.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Normalized Phase of Peak,Peak in Period")
+    csvFile.write("\n")
+    for ind, row in finalData.reset_index().iterrows():
+        PDCSAPminOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][0]
+        PDCSAPmaxOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][1]
+        if(row["SpType"][0] == "M"):
+            color = "red"
+        elif(row["SpType"][0] == "K"):
+            color = "orange"
+        elif(row["SpType"][0] == "G"):
+            color = "yellow"
+        elif(row["SpType"][0] == "F"):
+            color = "greenyellow"
+        else:
+            color = "gray"
+
+        if(len(row["pdcsapFoldedPeaksPhasePair"]) > 0):
+                pdcsapVals = pd.DataFrame(row["pdcsapFoldedPeaksPhasePair"])
+                for td, peak, pv in zip(pdcsapVals["Phase"], pdcsapVals["Peak"], row["pdcsapPeaks"]):
+                    normPhase = normalizePhase(td.value, np.abs(PDCSAPminOrigPhase), np.abs(PDCSAPmaxOrigPhase))
+                    csvFile.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{normPhase},{peak['FlarePeak']}")
+                    csvFile.write("\n")
+                    pdcsapbinningData.append({"SpType": f'{row["SpType"][0:2] if len(row["SpType"]) > 1 else row["SpType"][0]}',
+                                         "PDCSAPNormPhase": normPhase,
+                                         "Peak": peak["FlarePeak"]})
+                    if(peak["FlarePeak"] > 100):
+                        print(row["StarName"], "has over 100 peak")
+    csvFile.close()
+    pdcsapbinningData = pd.DataFrame(pdcsapbinningData)
+    PDCSAPdataList = []
+    PDCSAPlabelList = []
+    PDCSAPcolorList = []
+    PDCSAPdataList2dhistPhase = []
+    PDCSAPdataList2dhistPeak = []
+
+    if(len(pdcsapbinningData) < 1):
+        continue
+
+    PDCSAPdataList2dhistPhase.append(pdcsapbinningData[:]["PDCSAPNormPhase"])
+    PDCSAPdataList2dhistPeak.append(pdcsapbinningData[:]["Peak"])
+
+    xData = pd.DataFrame()
+    yData = pd.DataFrame()
+    for aX, aY in zip(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak):
+        xData = pd.concat([xData, aX], ignore_index=True)
+        yData = pd.concat([yData, aY], ignore_index=True)
+
+    xData = np.asarray(xData.values)[:,0]
+    yData = np.asarray(yData.values)[:,0]
+    PDCSAPdataList.append(pdcsapbinningData[:]["PDCSAPNormPhase"])
+
+    PDCSAPlabelList.append(f"{starName}")
+    PDCSAPcolorList.append(color)
+
+    for bins in binList:
+        figHisto, ((axHisto)) = plt.subplots(nrows=1, ncols=1)
+        y, binEdges, _ = axHisto.hist(PDCSAPdataList, bins,
+                                      label=PDCSAPlabelList,
+                                      color=PDCSAPcolorList,
+                                      stacked=True,
+                                      range=[0, 2])
+        bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
+        if(isinstance(y[0], np.ndarray)):
+            y = y[-1]
+        n_i = y
+        m_i = bincenters * np.pi
+        N = np.sum(n_i)
+        mean = np.sum(n_i * m_i)/N
+        stdDev = np.sqrt(np.sum(((n_i - mean)**2)) / (N-1))
+        menStd = np.sqrt(y)
+        axHisto.bar(bincenters[y > 0], y[y > 0], width=0, color='r', yerr=stdDev)
+        if(~np.isnan(stdDev) & ~np.isinf(stdDev)):
+            axHisto.set_ylim(0, max(y[y > 0 & ~np.isnan(y) & ~np.isinf(y)] + stdDev))
+        axHisto.set_ylabel("Num. flares")
+        axHisto.set_xlabel("Phase")
+        axHisto.set_title(f"Flare count in phase of {starName} with {bins} bins")
+        axHisto.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
+        axHisto.xaxis.set_major_locator(MaxNLocator(5))
+        axHisto.legend()
+
+        axHistoPhase = axHisto.twinx()
+        secAxisXdata = np.linspace(0, 2, num=10000)
+        secAxisYdata = np.cos(secAxisXdata*np.pi) + 1
+        axHistoPhase.plot(secAxisXdata, secAxisYdata)
+        axHistoPhase.set_ylim(0, 7)
+
+        plt.savefig(f"{locFolder}/{starNameR}-Flarecount-{bins}_Bins.png")
+        plt.close()
+
+        for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
+            figFlarepeakHist, ((axFlarepeakHist)) = plt.subplots(nrows=1, ncols=1)
+            axFlarepeakHist.set_ylabel("Flare peak")
+            axFlarepeakHist.set_xlabel("Phase")
+            H, xedges, yedges = np.histogram2d(xData, yData, bins=bins, range=[[0, 2], [0.99, maxY]])
+            cmax = 11
+            H_clipped = np.clip(H, None, cmax)
+            im = axFlarepeakHist.imshow(H_clipped.T, origin='lower', interpolation='nearest',
+                            extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]],
+                            aspect='auto', cmap='viridis')
+            figFlarepeakHist.colorbar(im, label='Counts', ax=axFlarepeakHist)
+            axFlarepeakHist.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
+            axFlarepeakHist.xaxis.set_major_locator(MaxNLocator(5))
+            axFlarepeakHist.set_title(f"Flare peak per phase histogram of {starName} with {bins} bins")
+            plt.savefig(f"{locFolder}/{starNameR}-Flarepeaks-{bins}_Bins_maxY-{maxY}.png")
+            plt.close()
+
+    for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
+        figFlarePeaks, ((axFlarePeaks)) = plt.subplots(nrows=1, ncols=1)
+        axFlarePeaks.scatter(pdcsapbinningData[:]["PDCSAPNormPhase"],
+                             pdcsapbinningData[:]["Peak"],
+                             label=f"{starName}", color=color)
+        axFlarePeaks.set_xlim(0, 2)
+        axFlarePeaks.set_ylim(0.99, maxY)
+        axFlarePeaks.set_ylabel("Flare peak")
+        axFlarePeaks.set_xlabel("Phase")
+        axFlarePeaks.set_title(f"Flare peaks per phase of {starName}")
+        axFlarePeaks.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
+        axFlarePeaks.xaxis.set_major_locator(MaxNLocator(5))
+        axFlarePeaks.legend()
+
+        plt.savefig(f"{locFolder}/{starNameR}-Flarepeaks_maxY-{maxY}.png")
+        plt.close()