From 082145f139a296c5cbb620e4db0f0a3e6b453427 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 9 Dec 2024 14:40:09 +1100 Subject: [PATCH 01/28] Update identify.py --- astrosource/identify.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/identify.py b/astrosource/identify.py index bfc65d1..27b1f6a 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -331,7 +331,7 @@ def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=180 def process_file(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger): try: # Read the file - photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1) + photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1, invalid_raise=False) logger.info(fn) # Validate file shape and range From f8c59eaef33965bfe2ca26a442e3b57798135741 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sat, 21 Dec 2024 18:15:43 +1100 Subject: [PATCH 02/28] multiprocessing speedup and better variable search function using a sigmoid --- astrosource/analyse.py | 252 +++++++++++++++++++---- astrosource/comparison.py | 210 ++++++++++++++----- astrosource/identify.py | 410 +++++++++++++++++++++++--------------- astrosource/periodic.py | 331 +++++++++++++++++++++++++++--- 4 files changed, 924 insertions(+), 279 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index 670baef..a8db8df 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -13,12 +13,14 @@ from tqdm import tqdm #import traceback import logging -from multiprocessing import Pool +from multiprocessing import Pool, cpu_count #from astrosource.utils import photometry_files_to_array, AstrosourceException from astrosource.utils import AstrosourceException from astrosource.plots import plot_variability +from astropy.stats import sigma_clip +from scipy.optimize import curve_fit import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt @@ -100,12 +102,171 @@ def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCoun ) # Multiprocessing with Pool - with Pool() as pool: + with Pool(processes=max([cpu_count()-1,1])) as pool: results = pool.map(worker, targetFile) # Filter out None results return [res for res in results if res is not None] + +def fit_sigma_clipped_poly(x, y, order=3, sigma=2, parentPath=''): + """ + Fits a sigma-clipped polynomial to the data and identifies outliers. + + Parameters: + x (array-like): The x-values of the data points. + y (array-like): The y-values of the data points. + order (int): The order of the polynomial to fit (default: 3). + sigma (float): The number of standard deviations for sigma clipping (default: 2). + + Returns: + None (plots the data, fitted polynomial, and identified outliers). + """ + # Ensure inputs are numpy arrays + x = np.array(x) + y = np.array(y) + + # Sigma clipping + clipped_data = sigma_clip(y, sigma=sigma, maxiters=5) + mask = clipped_data.mask + + # Fit polynomial to non-masked data + poly_coeff = np.polyfit(x[~mask], y[~mask], order) + poly_func = np.poly1d(poly_coeff) + + # Calculate residuals and standard deviation + y_fit = poly_func(x) + residuals = y - y_fit + std_dev = np.std(residuals[~mask]) + + # Identify outliers + outliers = residuals > (2 * std_dev) + + # Plot results + plt.figure(figsize=(10, 6)) + plt.scatter(x, y, color='blue', label='Data') + plt.plot(x, y_fit, color='red', label=f'{order}-order Polynomial Fit') + plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5) + plt.xlabel('x') + plt.ylabel('y') + plt.title('Sigma-Clipped Polynomial Fit and Outliers') + plt.legend() + plt.grid(True) + #plt.show() + + plt.savefig(parentPath / "results/polifitdiagram.png") + + #breakpoint() + +def sigmoid_func(x, L, k, x0): + """Sigmoid function: y = L / (1 + exp(-k * (x - x0))) + 0.01""" + return L / (1 + np.exp(-k * (x - x0))) + 0.01 + +def fit_sigma_clipped_sigmoid(x, y, sigma=2, parentPath=''): + """ + Fits a sigma-clipped sigmoid function to the data with weighted fitting and identifies outliers. + + Parameters: + x (array-like): The x-values of the data points. + y (array-like): The y-values of the data points. + sigma (float): The number of standard deviations for sigma clipping (default: 2). + + Returns: + outlier_indices (array): Indices of the identified outliers. + """ + # Ensure inputs are numpy arrays + x = np.array(x) + y = np.array(y) + + # Sigma clipping + clipped_data = sigma_clip(y, sigma=sigma, maxiters=5) + mask = clipped_data.mask + + # Create weights that emphasize the middle of the x-range + weights = np.exp(-((x - np.median(x))**2) / (2 * (np.std(x) / 2)**2)) + + # Fit sigmoid function to non-masked data with weights + p0 = [max(y), 1, np.median(x)] # Initial guesses for L, k, x0 + popt, _ = curve_fit(sigmoid_func, x[~mask], y[~mask], p0=p0, sigma=1/weights[~mask]) + + # Calculate residuals and standard deviation + y_fit = sigmoid_func(x, *popt) + residuals = y - y_fit + std_dev = np.std(residuals[~mask]) + + # Identify outliers + outliers = residuals > (2 * std_dev) + outlier_indices = np.where(outliers)[0] + + # Plot results + plt.figure(figsize=(10, 6)) + plt.scatter(x, y, color='blue', label='Data') + plt.plot(x, y_fit, color='red', label='Sigmoid Fit') + plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5) + plt.xlabel('x') + plt.ylabel('y') + plt.title('Weighted Sigma-Clipped Sigmoid Fit and Outliers') + plt.legend() + plt.grid(True) + plt.show() + + plt.savefig(parentPath / "results/polifitdiagram.png") + + return outlier_indices + +def exponential_func(x, a, b, c): + """Exponential function: y = a * exp(b * x) + c""" + return a * np.exp(b * x) + c + +def fit_sigma_clipped_exp(x, y, sigma=2, parentPath=''): + """ + Fits a sigma-clipped exponential function to the data and identifies outliers. + + Parameters: + x (array-like): The x-values of the data points. + y (array-like): The y-values of the data points. + sigma (float): The number of standard deviations for sigma clipping (default: 2). + + Returns: + outlier_indices (array): Indices of the identified outliers. + """ + # Ensure inputs are numpy arrays + x = np.array(x) + y = np.array(y) + + # Sigma clipping + clipped_data = sigma_clip(y, sigma=sigma, maxiters=5) + mask = clipped_data.mask + + # Fit exponential function to non-masked data + popt, _ = curve_fit(exponential_func, x[~mask], y[~mask], p0=(1, 0.1, 1)) + + # Calculate residuals and standard deviation + y_fit = exponential_func(x, *popt) + residuals = y - y_fit + std_dev = np.std(residuals[~mask]) + + # Identify outliers + outliers = residuals > (2 * std_dev) + outlier_indices = np.where(outliers)[0] + + # Plot results + plt.figure(figsize=(10, 6)) + plt.scatter(x, y, color='blue', label='Data') + plt.plot(x, y_fit, color='red', label='Exponential Fit') + plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5) + plt.xlabel('x') + plt.ylabel('y') + plt.title('Sigma-Clipped Exponential Fit and Outliers') + plt.legend() + plt.grid(True) + plt.show() + + plt.savefig(parentPath / "results/polifitdiagram.png") + + return outlier_indices + + def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, varsearchglobalstdev=-99.9, varsearchthresh=10000, varsearchstdev=2.0, varsearchmagwidth=0.25, varsearchminimages=0.3, photCoords=None, photFileHolder=None, fileList=None): ''' Find stable comparison stars for the target photometry and remove variables @@ -218,10 +379,23 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, savetxt(parentPath / "results/starVariability.csv", outputVariableHolder, delimiter=",", fmt='%0.8f', header='RA,DEC,DiffMag,Variability,No_of_images_used') + + + + + + #breakpoint() + ## Routine that actually pops out potential variables. starVar = np.asarray(outputVariableHolder) + outliers=fit_sigma_clipped_sigmoid(starVar[:,2],starVar[:,3], parentPath=parentPath) + + #breakpoint() + + potentialVariables = np.delete(starVar, outliers, axis=0) + meanMags = starVar[:,2] variations = starVar[:,3] @@ -230,51 +404,51 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, xbins = np.arange(np.min(meanMags), np.max(meanMags), xStepSize) ybins = np.arange(np.min(variations), np.max(variations), yStepSize) - #split it into one array and identify variables in bins with centre - variationsByMag=[] - potentialVariables=[] - for xbinner in range(len (xbins)): - - starsWithin=[] - for q in range(len(meanMags)): - if meanMags[q] >= xbins[xbinner] and meanMags[q] < xbins[xbinner]+xStepSize: - starsWithin.append(variations[q]) - - meanStarsWithin= (np.mean(starsWithin)) - stdStarsWithin= (np.std(starsWithin)) - variationsByMag.append([xbins[xbinner]+0.5*xStepSize,meanStarsWithin,stdStarsWithin]) - - # At this point extract RA and Dec of stars that may be variable - for q in range(len(starVar[:,2])): - if starVar[q,2] >= xbins[xbinner] and starVar[q,2] < xbins[xbinner]+xStepSize: - if varsearchglobalstdev != -99.9: - if starVar[q,3] > varsearchglobalstdev : - potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) - elif starVar[q,3] > (meanStarsWithin + varsearchstdev*stdStarsWithin): - #print (starVar[q,3]) - potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) - - potentialVariables=np.array(potentialVariables) - logger.debug("Potential Variables Identified: " + str(potentialVariables.shape[0])) + # #split it into one array and identify variables in bins with centre + # variationsByMag=[] + # potentialVariables=[] + # for xbinner in range(len (xbins)): + + # starsWithin=[] + # for q in range(len(meanMags)): + # if meanMags[q] >= xbins[xbinner] and meanMags[q] < xbins[xbinner]+xStepSize: + # starsWithin.append(variations[q]) + + # meanStarsWithin= (np.mean(starsWithin)) + # stdStarsWithin= (np.std(starsWithin)) + # variationsByMag.append([xbins[xbinner]+0.5*xStepSize,meanStarsWithin,stdStarsWithin]) + + # # At this point extract RA and Dec of stars that may be variable + # for q in range(len(starVar[:,2])): + # if starVar[q,2] >= xbins[xbinner] and starVar[q,2] < xbins[xbinner]+xStepSize: + # if varsearchglobalstdev != -99.9: + # if starVar[q,3] > varsearchglobalstdev : + # potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) + # elif starVar[q,3] > (meanStarsWithin + varsearchstdev*stdStarsWithin): + # #print (starVar[q,3]) + # potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) + + # potentialVariables=np.array(potentialVariables) + # logger.debug("Potential Variables Identified: " + str(potentialVariables.shape[0])) if potentialVariables.shape[0] == 0: logger.info("No Potential Variables identified in this set of data using the parameters requested.") else: savetxt(parentPath / "results/potentialVariables.csv", potentialVariables , delimiter=",", fmt='%0.8f', header='RA,DEC,DiffMag,Variability') - plot_variability(outputVariableHolder, potentialVariables, parentPath, compFile) + plot_variability(outputVariableHolder, potentialVariables, parentPath, compFile) - plt.cla() - fig, ax = plt.subplots(figsize =(10, 7)) - plt.hist2d(meanMags, variations, bins =[xbins, ybins], norm=colors.LogNorm(), cmap = plt.cm.YlOrRd) - plt.colorbar() - plt.title("Variation Histogram") - ax.set_xlabel('Mean Differential Magnitude') - ax.set_ylabel('Variation (Standard Deviation)') - plt.plot(potentialVariables[:,2],potentialVariables[:,3],'bo') - plt.tight_layout() + plt.cla() + fig, ax = plt.subplots(figsize =(10, 7)) + plt.hist2d(meanMags, variations, bins =[xbins, ybins], norm=colors.LogNorm(), cmap = plt.cm.YlOrRd) + plt.colorbar() + plt.title("Variation Histogram") + ax.set_xlabel('Mean Differential Magnitude') + ax.set_ylabel('Variation (Standard Deviation)') + plt.plot(potentialVariables[:,2],potentialVariables[:,3],'bo') + plt.tight_layout() - plt.savefig(parentPath / "results/Variation2DHistogram.png") + plt.savefig(parentPath / "results/Variation2DHistogram.png") return outputVariableHolder diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 3684601..31167a9 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -13,7 +13,7 @@ import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt - +from multiprocessing import Pool, cpu_count from numpy import min, max, median, std, isnan, delete, genfromtxt, savetxt, load, \ asarray, add, append, log10, average, array, where, cos from astropy.units import degree, arcsecond @@ -282,77 +282,184 @@ def read_data_files(parentPath, fileList): def ensemble_comparisons(photFileArray, compFile, parentPath, photSkyCoord): - # get rid of dumb for loop - fileCount = [] - q=0 - for photFile in photFileArray: - allCounts = 0.0 + # # get rid of dumb for loop + # fileCount = [] + # q=0 + # for photFile in photFileArray: + # allCounts = 0.0 - if compFile.size ==2 and compFile.shape[0]==2: - matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) - else: - matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) + # if compFile.size ==2 and compFile.shape[0]==2: + # matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) + # else: + # matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) - idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) + # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - allCounts = add(allCounts,sum(photFile[idx,4])) + # allCounts = add(allCounts,sum(photFile[idx,4])) - fileCount.append(allCounts) - q=q+1 + # fileCount.append(allCounts) + # q=q+1 + + # logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(np.sum(np.array(fileCount)))) + # return fileCount + + #degree = np.deg2rad(1) # Degree to radians conversion + fileCount = np.zeros(len(photFileArray)) # Pre-allocate output array + + # Prepare the `matchCoord` SkyCoord object + if compFile.size == 2 and compFile.shape[0] == 2: + matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree) + else: + matchCoord = SkyCoord(ra=compFile[:,0] * degree, dec=compFile[:,1] * degree) - logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(np.sum(np.array(fileCount)))) + # Process all photSkyCoord entries in a vectorized way + for q, (photFile, skyCoord) in enumerate(zip(photFileArray, photSkyCoord)): + # Match coordinates + idx, _, _ = matchCoord.match_to_catalog_sky(skyCoord) + # Accumulate counts + fileCount[q] = np.sum(photFile[idx, 4]) + + total_counts = np.sum(fileCount) + logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(total_counts)) return fileCount + +# def process_phot_file_for_variation(args): +# """ +# Process a single photFileArray element. +# """ +# q, matchCoord, photSkyCoord, photFileArray, fileCount = args + +# photFile = photFileArray[q] +# idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) +# comp_diff_mag = 2.5 * np.log10(photFile[idx][4] / fileCount[q]) +# instr_mag = -2.5 * np.log10(photFile[idx][4]) + +# return comp_diff_mag, instr_mag + +def process_phot_file_for_variation(args): + q, matchCoord, photSkyCoord, photFileArray, fileCount = args + idx, _, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) + photFile = photFileArray[q] + compDiffMags = 2.5 * np.log10(photFile[idx, 4] / fileCount[q]) + instrMags = -2.5 * np.log10(photFile[idx, 4]) + return compDiffMags, instrMags + def calculate_comparison_variation(compFile, photFileArray, fileCount, parentPath, photSkyCoord): - stdCompStar=[] - sortStars=[] - logger.debug("Calculating Variation in Individual Comparisons") + # stdCompStar=[] + # sortStars=[] + # logger.debug("Calculating Variation in Individual Comparisons") - if compFile.size ==2 and compFile.shape[0]==2: - compDiffMags = [] + # if compFile.size ==2 and compFile.shape[0]==2: + # compDiffMags = [] - matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) - for q, photFile in enumerate(photFileArray): - idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - compDiffMags = append(compDiffMags,2.5 * log10(photFile[idx][4]/fileCount[q])) - instrMags = -2.5 * log10(photFile[idx][4]) + # # matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) + # # for q, photFile in enumerate(photFileArray): + # # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) + # # compDiffMags = append(compDiffMags,2.5 * log10(photFile[idx][4]/fileCount[q])) + # # instrMags = -2.5 * log10(photFile[idx][4]) + + # matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree) + # args = [(q, matchCoord, photSkyCoord, photFileArray, fileCount) for q in range(len(photFileArray))] + + # with Pool() as pool: + # results = pool.map(process_phot_file_for_variation, args) + + # # Unpack results + # compDiffMags = np.array([result[0] for result in results]) + # instrMags = np.array([result[1] for result in results]) + + # stdCompDiffMags=std(compDiffMags) + # medCompDiffMags=np.nanmedian(compDiffMags) + # medInstrMags= np.nanmedian(instrMags) + + # if np.isnan(stdCompDiffMags) : + # logger.error("Star Variability non rejected") + # stdCompDiffMags=99 + # stdCompStar.append(stdCompDiffMags) + # sortStars.append([compFile[0],compFile[1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) + + + + # else: + # sortStars=[] + # compDiffMags = [] + # instrMags=[] + # matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) + # for q, photFile in enumerate(photFileArray): + # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) + # compDiffMags.append(2.5 * log10(photFile[idx,4]/fileCount[q])) + # instrMags.append(-2.5 * log10(photFile[idx,4])) + + # compDiffMags=array(compDiffMags) + # instrMags=array(instrMags) + + # sortStars=[] + # for z in range(len(compDiffMags[0])): + # stdCompDiffMags=std(compDiffMags[:,z]) + # medCompDiffMags=np.nanmedian(compDiffMags[:,z]) + # medInstrMags=np.nanmedian(instrMags[:,z]) + # if np.isnan(stdCompDiffMags) : + # logger.error("Star Variability non rejected") + # stdCompDiffMags=99 + # stdCompStar.append(stdCompDiffMags) + # sortStars.append([compFile[z,0],compFile[z,1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) + + # return stdCompStar, sortStars + + stdCompStar = [] + sortStars = [] + logger.debug("Calculating Variation in Individual Comparisons") - stdCompDiffMags=std(compDiffMags) - medCompDiffMags=np.nanmedian(compDiffMags) - medInstrMags= np.nanmedian(instrMags) + if compFile.size == 2 and compFile.shape[0] == 2: + # Case for single comparison star + matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree) - if np.isnan(stdCompDiffMags) : - logger.error("Star Variability non rejected") - stdCompDiffMags=99 - stdCompStar.append(stdCompDiffMags) - sortStars.append([compFile[0],compFile[1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) + args = [(q, matchCoord, photSkyCoord, photFileArray, fileCount) for q in range(len(photFileArray))] + + with Pool(processes=max([cpu_count()-1,1])) as pool: + results = pool.map(process_phot_file_for_variation, args) + + compDiffMags = np.array([result[0] for result in results]) + instrMags = np.array([result[1] for result in results]) + + stdCompDiffMags = np.std(compDiffMags) + medCompDiffMags = np.nanmedian(compDiffMags) + medInstrMags = np.nanmedian(instrMags) + if np.isnan(stdCompDiffMags): + logger.error("Star Variability non rejected") + stdCompDiffMags = 99 + stdCompStar.append(stdCompDiffMags) + sortStars.append([compFile[0], compFile[1], stdCompDiffMags, medCompDiffMags, medInstrMags] + [0.0] * 8) else: - sortStars=[] + # Case for multiple comparison stars + matchCoord = SkyCoord(ra=compFile[:, 0] * degree, dec=compFile[:, 1] * degree) compDiffMags = [] - instrMags=[] - matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) + instrMags = [] + for q, photFile in enumerate(photFileArray): - idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - compDiffMags.append(2.5 * log10(photFile[idx,4]/fileCount[q])) - instrMags.append(-2.5 * log10(photFile[idx,4])) - - compDiffMags=array(compDiffMags) - instrMags=array(instrMags) - - sortStars=[] - for z in range(len(compDiffMags[0])): - stdCompDiffMags=std(compDiffMags[:,z]) - medCompDiffMags=np.nanmedian(compDiffMags[:,z]) - medInstrMags=np.nanmedian(instrMags[:,z]) - if np.isnan(stdCompDiffMags) : + idx, _, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) + compDiffMags.append(2.5 * np.log10(photFile[idx, 4] / fileCount[q])) + instrMags.append(-2.5 * np.log10(photFile[idx, 4])) + + compDiffMags = np.array(compDiffMags) + instrMags = np.array(instrMags) + + for z in range(compDiffMags.shape[1]): + stdCompDiffMags = np.std(compDiffMags[:, z]) + medCompDiffMags = np.nanmedian(compDiffMags[:, z]) + medInstrMags = np.nanmedian(instrMags[:, z]) + + if np.isnan(stdCompDiffMags): logger.error("Star Variability non rejected") - stdCompDiffMags=99 + stdCompDiffMags = 99 + stdCompStar.append(stdCompDiffMags) - sortStars.append([compFile[z,0],compFile[z,1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) + sortStars.append([compFile[z, 0], compFile[z, 1], stdCompDiffMags, medCompDiffMags, medInstrMags] + [0.0] * 8) return stdCompStar, sortStars @@ -1029,6 +1136,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n logger.debug(tabl) except: print ("something fishy in the calib mags table") + import traceback; logger.error(traceback.print_exc()) # Get the set of least variable stars to use as a comparison to calibrate the files (to eventually get the *ACTUAL* standards if asarray(calibStands).shape[0] == 0: diff --git a/astrosource/identify.py b/astrosource/identify.py index 27b1f6a..16ea204 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -4,8 +4,11 @@ import os import logging import pickle +from concurrent.futures import ProcessPoolExecutor -from numpy import genfromtxt, delete, asarray, save, savetxt, load, transpose, isnan, zeros, max, min, nan, where, average, cos, hstack, array, column_stack, copy, c_, sqrt, ptp +#import copy + +from numpy import genfromtxt, delete, argmax, asarray, save, savetxt, load, transpose, isnan, zeros, max, min, nan, where, average, cos, hstack, array, column_stack, copy, c_, sqrt, ptp from astropy import units as u from astropy.units import degree, arcsecond from astropy import wcs @@ -15,7 +18,7 @@ from barycorrpy import utc_tdb from tqdm import tqdm from prettytable import PrettyTable -from multiprocessing import Pool +from multiprocessing import Pool,cpu_count from astrosource.utils import AstrosourceException from astrosource.comparison import catalogue_call @@ -128,7 +131,7 @@ def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, de def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut): # Wrap arguments for multiprocessing args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut) - with Pool() as pool: + with Pool(processes=max([cpu_count()-1,1])) as pool: results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(filelist))) # Filter valid results @@ -229,6 +232,80 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr return outfile, photFile + +def process_convert_photometry_file(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger): + try: + # Read the file + photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1, invalid_raise=False) + logger.info(fn) + + # Validate file shape and range + if photFile.size <= 16 or photFile.shape[1] != 8: + logger.debug(f"REJECT {fn}: Invalid size or shape.") + return None + + if max(photFile[:, 0]) >= 360 or max(photFile[:, 1]) >= 90: + logger.debug(f"REJECT {fn}: Invalid RA/Dec range.") + return None + + # Reject NaN entries + photFile = photFile[~isnan(photFile).any(axis=1)] + + # Apply radial cut or image edge trimming + if racut != -99.9 and deccut != -99.9 and radiuscut != -99.9: + distanceArray = sqrt((photFile[:, 0] - racut)**2 + (photFile[:, 1] - deccut)**2) + photFile = delete(photFile, where(distanceArray > radiuscut / 60), axis=0) + else: + # Handle RA crossover + if (max(photFile[:, 0]) - min(photFile[:, 0])) > 180: + photFile[:, 0] = where(photFile[:, 0] > 180, photFile[:, 0] - 360, photFile[:, 0] + 360) + + # Clip edges + raRange, decRange = ptp(photFile[:, 0]), ptp(photFile[:, 1]) + raClip, decClip = raRange * ignoreedgefraction, decRange * ignoreedgefraction + photFile = photFile[ + (photFile[:, 0] > min(photFile[:, 0]) + raClip) & + (photFile[:, 0] < max(photFile[:, 0]) - raClip) & + (photFile[:, 1] > min(photFile[:, 1]) + decClip) & + (photFile[:, 1] < max(photFile[:, 1]) - decClip) + ] + + # Remove zero and low-count entries + photFile = photFile[(photFile[:, 4] > lowestcounts) & (photFile[:, 0] != 0) & (photFile[:, 1] != 0)] + + # Prepare output + photFile = c_[photFile, zeros((len(photFile), 4))] + photFile[:, 8:10] = nan + photFile[:, 11] = photFile[:, 4].copy() + + # Return processed data + if photFile.size > 16: + filepath = Path(fn).with_suffix('.npy') + photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree) + return filepath.name, photFile, photSkyCoord + + logger.debug(f"REJECT {fn}: Insufficient valid entries.") + return None + except Exception as e: + logger.error(f"Error processing file {fn}: {e}") + return None + + +def process_convert_files_multiprocessing(filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger): + # Pool for multiprocessing + with Pool(processes=max([cpu_count()-1,1])) as pool: + results = pool.starmap( + process_convert_photometry_file, + [(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger) for fn in filelist] + ) + + # Filter results + results = [res for res in results if res is not None] + new_files, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], []) + return list(new_files), list(photFileHolder), list(photSkyCoord) + + + def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=1800, racut=-99.9, deccut=-99.9, radiuscut=-99.9): new_files = [] @@ -328,83 +405,13 @@ def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=180 # logger.debug(fn) - def process_file(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger): - try: - # Read the file - photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1, invalid_raise=False) - logger.info(fn) - - # Validate file shape and range - if photFile.size <= 16 or photFile.shape[1] != 8: - logger.debug(f"REJECT {fn}: Invalid size or shape.") - return None - - if max(photFile[:, 0]) >= 360 or max(photFile[:, 1]) >= 90: - logger.debug(f"REJECT {fn}: Invalid RA/Dec range.") - return None - - # Reject NaN entries - photFile = photFile[~isnan(photFile).any(axis=1)] - - # Apply radial cut or image edge trimming - if racut != -99.9 and deccut != -99.9 and radiuscut != -99.9: - distanceArray = sqrt((photFile[:, 0] - racut)**2 + (photFile[:, 1] - deccut)**2) - photFile = delete(photFile, where(distanceArray > radiuscut / 60), axis=0) - else: - # Handle RA crossover - if (max(photFile[:, 0]) - min(photFile[:, 0])) > 180: - photFile[:, 0] = where(photFile[:, 0] > 180, photFile[:, 0] - 360, photFile[:, 0] + 360) - - # Clip edges - raRange, decRange = ptp(photFile[:, 0]), ptp(photFile[:, 1]) - raClip, decClip = raRange * ignoreedgefraction, decRange * ignoreedgefraction - photFile = photFile[ - (photFile[:, 0] > min(photFile[:, 0]) + raClip) & - (photFile[:, 0] < max(photFile[:, 0]) - raClip) & - (photFile[:, 1] > min(photFile[:, 1]) + decClip) & - (photFile[:, 1] < max(photFile[:, 1]) - decClip) - ] - - # Remove zero and low-count entries - photFile = photFile[(photFile[:, 4] > lowestcounts) & (photFile[:, 0] != 0) & (photFile[:, 1] != 0)] - - # Prepare output - photFile = c_[photFile, zeros((len(photFile), 4))] - photFile[:, 8:10] = nan - photFile[:, 11] = photFile[:, 4].copy() - - # Return processed data - if photFile.size > 16: - filepath = Path(fn).with_suffix('.npy') - photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree) - return filepath.name, photFile, photSkyCoord - - logger.debug(f"REJECT {fn}: Insufficient valid entries.") - return None - except Exception as e: - logger.error(f"Error processing file {fn}: {e}") - return None - - - def process_files_multiprocessing(filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger): - # Pool for multiprocessing - with Pool() as pool: - results = pool.starmap( - process_file, - [(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger) for fn in filelist] - ) - - # Filter results - results = [res for res in results if res is not None] - new_files, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], []) - return list(new_files), list(photFileHolder), list(photSkyCoord) + # Example usage - if __name__ == "__main__": - new_files, photFileHolder, photSkyCoord = process_files_multiprocessing( - filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger - ) + new_files, photFileHolder, photSkyCoord = process_convert_files_multiprocessing( + filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger + ) if new_files ==[] or photFileHolder==[] : @@ -466,6 +473,90 @@ def gather_files(paths, filelist=None, filetype="fz", bjd=False, ignoreedgefract return phot_list, filterCode, photFileHolder, photSkyCoord +def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): + try: + imgRejFlag = 0 + photReject = [] + + # Check minimum stars requirement + if referenceFrame.shape[0] < mincompstars: + return None, referenceFrame, photReject + + # Validate photometry file + if photFile.size > imgsize and photFile.size > 7: + photRAandDec = photCoords[index] + rejectStars = [] + + if (max(photFile[:, 0]) <= 360) and (photFile[0][0] != 'null') and (photFile[0][0] != 0.0): + testStars = SkyCoord(ra=referenceFrame[:, 0] * u.degree, dec=referenceFrame[:, 1] * u.degree) + idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec) + rejectStars = where(d2d.arcsecond > acceptDistance)[0] + else: + logger.debug(f"Image {index} rejected due to problematic entries.") + imgRejFlag = 1 + photReject.append(index) + + # Remove rejected stars from the reference frame + if rejectStars.size > 0: + if not ((len(rejectStars) / referenceFrame.shape[0]) > starreject and index > rejectStart): + referenceFrame = delete(referenceFrame, rejectStars, axis=0) + logger.debug(f"Image {index}: Removed {len(rejectStars)} stars.") + else: + logger.debug(f"Image {index} rejected due to too high a fraction of rejected stars.") + imgRejFlag = 1 + photReject.append(index) + + if imgRejFlag == 0: + logger.debug(f"Image {index}: All stars present.") + elif photFile.size <= 7: + logger.error(f"Image {index}: WCS coordinates broken or too few stars.") + photReject.append(index) + else: + logger.debug(f"Image {index}: Contains too few stars.") + photReject.append(index) + + return imgRejFlag, referenceFrame, photReject + except Exception as e: + logger.error(f"Error processing image {index}: {e}") + return None, referenceFrame, [index] + + +def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): + with Pool(processes=max([cpu_count()-1,1])) as pool: + # Prepare arguments + args = [ + (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger) + for index in range(len(photFileHolder)) + ] + + # Process files in parallel + results = pool.starmap(process_phot_file, args) + + # Collect results + updated_referenceFrame = referenceFrame + photReject = [] + + for imgRejFlag, refFrameUpdate, rejects in results: + if refFrameUpdate is not None: + updated_referenceFrame = refFrameUpdate + photReject.extend(rejects) + + return updated_referenceFrame, photReject + + +def evaluate_file_size(file_index, file_list, phot_file_holder, phot_coords): + """Evaluate a file to determine its size and suitability.""" + file = file_list[file_index] + phot_file = phot_file_holder[file_index] + file_size = phot_file.size + file_ra_dec = phot_coords[file_index] + return file_size, file_index, file_ra_dec + + +def compute_size(phot_file): + """Compute the size for a specific file.""" + return phot_file.size + def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskymapper=False,closerejectd=5.0, photCoords=None, photFileHolder=None, mincompstars=0.1, mincompstarstotal=-99, starreject=0.1 , acceptDistance=1.0, lowcounts=2000, hicounts=3000000, imageFracReject=0.0, rejectStart=3, maxcandidatestars=10000, restrictcompcolourcentre=-99.0, restrictcompcolourrange=-99.0, filterCode=None, restrictmagbrightest=-99.0, restrictmagdimmest=99.0, minfractionimages=0.5): """ Finds stars useful for photometry in each photometry/data file @@ -510,21 +601,89 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym if os.path.exists(paths['parent'] / "photFileHolder"): os.remove(paths['parent'] / "photFileHolder") - fileSizer=0 - logger.info("Finding image with most stars detected and reject ones with bad WCS") + # fileSizer=0 + # logger.info("Finding image with most stars detected and reject ones with bad WCS") + # referenceFrame = None + + # counter=0 + # for file in fileList: + # photFile = photFileHolder[counter] + + # # Sort through and find the largest file and use that as the reference file + # if photFile.size > fileSizer: + # referenceFrame = photFile + # fileSizer = photFile.size + # fileRaDec = photCoords[counter] + # logger.debug("{} - {}".format(photFile.size, file)) + # counter=counter+1 + + logger.info("Finding image with most stars detected and rejecting ones with bad WCS") referenceFrame = None - - counter=0 - for file in fileList: - photFile = photFileHolder[counter] - - # Sort through and find the largest file and use that as the reference file - if photFile.size > fileSizer: - referenceFrame = photFile - fileSizer = photFile.size - fileRaDec = photCoords[counter] - logger.debug("{} - {}".format(photFile.size, file)) - counter=counter+1 + file_sizer = 0 + fileRaDec = None + + + + + # list_of_sizes=[] + + + # with Pool() as pool: + # # Use pool.map to distribute the computation across multiple processes + # list_of_sizes = pool.starmap(compute_size, [(i, photFileHolder) for i in range(len(photFileHolder))]) + + + # def compute_size(phot_file): + # """Compute the size for a specific file.""" + # return phot_file.size + # import time + + + # with Pool() as pool: + # # Parallelize the size computation + # list_of_sizes = pool.map(compute_size, photFileHolder) + + + import time + # googtime=time.time() + + + + list_of_sizes=[] + for filething in photFileHolder: + # zingtime=time.time() + list_of_sizes.append(filething.size) + # print ("Z: " + str(time.time()-zingtime)) + + # print ("D: " + str(time.time()-googtime)) + + largest_file_index= argmax(list_of_sizes) + + referenceFrame=photFileHolder[largest_file_index] + fileRaDec = photCoords[largest_file_index] + fileSizer=referenceFrame.size + + # breakpoint() + + # for file in len(photFileHolder): + # list_of_sizes.append(photFileHolder.size) + + # with ProcessPoolExecutor() as executor: + # results = list(executor.map( + # evaluate_file_size, + # range(len(fileList)), + # [fileList] * len(fileList), + # [photFileHolder] * len(fileList), + # [photCoords] * len(fileList) + # )) + + # # Process the results + # for file_size, index, file_ra_dec in results: + # if file_size > file_sizer: + # referenceFrame = photFileHolder[index] + # file_sizer = file_size + # fileRaDec = file_ra_dec + # logger.debug("{} - {}".format(file_size, fileList[index])) if not referenceFrame.size: raise AstrosourceException("No suitable reference files found") @@ -670,75 +829,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym # q=q+1 - def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): - try: - imgRejFlag = 0 - photReject = [] - - # Check minimum stars requirement - if referenceFrame.shape[0] < mincompstars: - return None, referenceFrame, photReject - - # Validate photometry file - if photFile.size > imgsize and photFile.size > 7: - photRAandDec = photCoords[index] - rejectStars = [] - - if (max(photFile[:, 0]) <= 360) and (photFile[0][0] != 'null') and (photFile[0][0] != 0.0): - testStars = SkyCoord(ra=referenceFrame[:, 0] * u.degree, dec=referenceFrame[:, 1] * u.degree) - idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec) - rejectStars = where(d2d.arcsecond > acceptDistance)[0] - else: - logger.debug(f"Image {index} rejected due to problematic entries.") - imgRejFlag = 1 - photReject.append(index) - - # Remove rejected stars from the reference frame - if rejectStars.size > 0: - if not ((len(rejectStars) / referenceFrame.shape[0]) > starreject and index > rejectStart): - referenceFrame = delete(referenceFrame, rejectStars, axis=0) - logger.debug(f"Image {index}: Removed {len(rejectStars)} stars.") - else: - logger.debug(f"Image {index} rejected due to too high a fraction of rejected stars.") - imgRejFlag = 1 - photReject.append(index) - - if imgRejFlag == 0: - logger.debug(f"Image {index}: All stars present.") - elif photFile.size <= 7: - logger.error(f"Image {index}: WCS coordinates broken or too few stars.") - photReject.append(index) - else: - logger.debug(f"Image {index}: Contains too few stars.") - photReject.append(index) - - return imgRejFlag, referenceFrame, photReject - except Exception as e: - logger.error(f"Error processing image {index}: {e}") - return None, referenceFrame, [index] - - - def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): - with Pool() as pool: - # Prepare arguments - args = [ - (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger) - for index in range(len(photFileHolder)) - ] - - # Process files in parallel - results = pool.starmap(process_phot_file, args) - - # Collect results - updated_referenceFrame = referenceFrame - photReject = [] - - for imgRejFlag, refFrameUpdate, rejects in results: - if refFrameUpdate is not None: - updated_referenceFrame = refFrameUpdate - photReject.extend(rejects) - - return updated_referenceFrame, photReject # Example usage diff --git a/astrosource/periodic.py b/astrosource/periodic.py index 10263c4..96e4801 100644 --- a/astrosource/periodic.py +++ b/astrosource/periodic.py @@ -5,10 +5,12 @@ import sys import os +import traceback + import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt -from multiprocessing import Pool +from multiprocessing import Pool, cpu_count from functools import partial import logging @@ -36,47 +38,316 @@ # SOFTWARE. def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, paths, filterCode, logger): + # try: + # logger.debug(file) + # variableName = file.stem.split('_')[0] + # logger.debug("Variable Name: {}".format(variableName)) + + # # Load data + # varData = genfromtxt(file, dtype=float, delimiter=',') + # calibFile = file.parent / "{}{}".format(file.stem.replace('diff', 'calib'), file.suffix) + # logger.debug(calibFile) + + # calibData = None + # if calibFile.exists(): + # calibData = genfromtxt(calibFile, dtype=float, delimiter=',') + + # # Determine which dataset to use for analysis + # if calibFile.exists() and calibData is not None and calibData.size > 3: + # dataToUse = calibData + # logger.info("Using calibrated data") + # else: + # dataToUse = varData + # logger.info("Using differential data") + + # # Phase dispersion minimization + # pdm = phase_dispersion_minimization(dataToUse, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) + + # # Save plots and results + # plt.figure(figsize=(15, 5)) + # plt.plot(pdm["periodguess_array"], pdm["distance_results"]) + # plt.gca().invert_yaxis() + # plt.title("Range {0} d Steps: {1}".format(maxperiod - minperiod, periodsteps)) + # plt.xlabel(r"Trial Period") + # plt.ylabel(r"Likelihood of Period") + # plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") + # plt.clf() + + # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + # f.write(f"Variable: {variableName}\n") + # f.write(f"Distance Method Estimate (days): {pdm['distance_minperiod']}\n") + # f.write(f"Distance method error: {pdm['distance_error']}\n") + + # except Exception as e: + # logger.error(f"Error processing file {file}: {e}") + try: - logger.debug(file) - variableName = file.stem.split('_')[0] + trialRange=[minperiod, maxperiod] + + # logger.debug(file) + variableName=file.stem.split('_')[0] + #logger.debug(str(outcatPath).replace('//','')) logger.debug("Variable Name: {}".format(variableName)) - - # Load data varData = genfromtxt(file, dtype=float, delimiter=',') - calibFile = file.parent / "{}{}".format(file.stem.replace('diff', 'calib'), file.suffix) + calibFile = file.parent / "{}{}".format(file.stem.replace('diff','calib'), file.suffix) logger.debug(calibFile) - - calibData = None if calibFile.exists(): - calibData = genfromtxt(calibFile, dtype=float, delimiter=',') + calibData=genfromtxt(calibFile, dtype=float, delimiter=',') - # Determine which dataset to use for analysis - if calibFile.exists() and calibData is not None and calibData.size > 3: - dataToUse = calibData - logger.info("Using calibrated data") + #print (calibData.size) + if calibFile.exists(): + if (calibData.size > 3): + pdm=phase_dispersion_minimization(calibData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) + else: + logger.info("Calibration File not large enough to run period methods on Calibrated data, running on differential data") + pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) else: - dataToUse = varData - logger.info("Using differential data") - - # Phase dispersion minimization - pdm = phase_dispersion_minimization(dataToUse, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) + pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) - # Save plots and results plt.figure(figsize=(15, 5)) + + logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])) + logger.debug("Distance method error: " + str(pdm["distance_error"])) + + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("Variable : "+str(variableName) +"\n") + f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n") + f.write("Distance method error: " + str(pdm["distance_error"])+"\n") + plt.plot(pdm["periodguess_array"], pdm["distance_results"]) plt.gca().invert_yaxis() - plt.title("Range {0} d Steps: {1}".format(maxperiod - minperiod, periodsteps)) + plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) plt.xlabel(r"Trial Period") plt.ylabel(r"Likelihood of Period") plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") plt.clf() + + if (varData.size > 3): + phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1 + + plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') + plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(f"Differential {filterCode} Magnitude") + plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png") + plt.clf() - with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - f.write(f"Variable: {variableName}\n") - f.write(f"Distance Method Estimate (days): {pdm['distance_minperiod']}\n") - f.write(f"Distance method error: {pdm['distance_error']}\n") + if calibFile.exists(): + if (calibData.size > 3): + phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1 + plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') + plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(f"Calibrated {filterCode} Magnitude") + plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png") + plt.clf() + + tempPeriodCatOut=[] + for g in range(len(calibData[:,0])): + tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + if (varData.size > 3): + + tempPeriodCatOut=[] + for g in range(len(phaseTest)): + tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + tempPeriodCatOut=[] + for g in range(len(varData[:,0])): + tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4): + logger.info("No PDM results due to lack of datapoint coverage") + else: + logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])) + phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1 + logger.debug("PDM method error: " + str(pdm["stdev_error"])) + + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n") + f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n") + + + plt.plot(pdm["periodguess_array"], pdm["stdev_results"]) + plt.gca().invert_yaxis() + plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png") + plt.clf() + + if (varData.size > 3): + plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') + plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(r"Differential " + str(filterCode) + " Magnitude") + plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png") + plt.clf() + + if calibFile.exists(): + if (calibData.size > 3): + phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1 + plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') + plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude") + plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png") + plt.clf() + + tempPeriodCatOut=[] + for g in range(len(calibData[:,0])): + tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + if (varData.size > 3): + + tempPeriodCatOut=[] + for g in range(len(phaseTest)): + tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + tempPeriodCatOut=[] + for g in range(len(varData[:,0])): + tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + # Plot publication plots + + plt.figure(figsize=(5, 3)) + + plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5) + plt.gca().invert_yaxis() + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300) + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps") + + plt.clf() + + plt.figure(figsize=(5, 3)) + + plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5) + plt.gca().invert_yaxis() + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) + plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300) + plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps") + + plt.clf() + + + # ANOVA + + if calibFile.exists(): + if (calibData.size > 3): + if len(calibData[:,0]) < 75: + binsize=0.1 + else: + binsize=0.05 + minperbin=int((len(calibData[:,0])/10)) + elif (varData.size > 3): + if len(varData[:,0]) < 75: + binsize=0.1 + else: + binsize=0.05 + minperbin=int((len(varData[:,0])/10)) + + else: + if (varData.size > 3): + if len(varData[:,0]) < 75: + binsize=0.1 + else: + binsize=0.05 + minperbin=int((len(varData[:,0])/10)) + + if 'minperbin' in locals(): + if minperbin > 10: + minperbin=10 + else: + minperbin = 3 + + + # Theta Anova Method off for the moment until I put in a command-line option + + # if calibFile.exists(): + # if (calibData.size > 3): + # aovoutput=aov_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + # else: + # if (varData.size > 3): + # aovoutput=aov_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + + + # logger.debug("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])) + #with open(paths['parent'] / "periodEstimates.txt", "a+") as f: + # f.write("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])+"\n") + + if calibFile.exists(): + if (calibData.size > 3): + aovhmoutput=aovhm_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + elif (varData.size > 3): + aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + + else: + if (varData.size > 3): + aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + + logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])) + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n") + + # LOMB SCARGLE + for nts in range(2): + if calibFile.exists(): + if (calibData.size > 3): + lscargoutput = LombScargleMultiterm('periodifile', (calibData[:, 0]), (calibData[:, 1]), (calibData[:, 2]), + nterms=nts+1, + periodlower=minperiod, periodupper=maxperiod, samples=20, + disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + + logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)) + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n") + else: + if (varData.size > 3): + lscargoutput = LombScargleMultiterm('periodifile', (varData[:, 0]), (varData[:, 1]), (varData[:, 2]), + nterms=nts+1, + periodlower=minperiod, periodupper=maxperiod, samples=20, + disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) + + logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)) + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n") + except Exception as e: + print(traceback.format_exc()) logger.error(f"Error processing file {file}: {e}") def aov_theta(times, mags, errs, frequency, @@ -1819,8 +2090,10 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1 logger=logger ) + #breakpoint() + # Use multiprocessing - with Pool(processes=os.cpu_count()) as pool: + with Pool(processes=max([cpu_count()-1,1])) as pool: pool.map(worker, fileList) # Load in the files @@ -2087,7 +2360,7 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1 # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: # f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n") - if 'pdm' in locals(): - return pdm["distance_minperiod"] - else: - return 0.0 + # if 'pdm' in locals(): + # return pdm["distance_minperiod"] + # else: + return 0.0 From bb0c2d50f865cc8d14b18ffd49c92352085f1d55 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 22 Dec 2024 08:16:08 +1100 Subject: [PATCH 03/28] missing try except --- astrosource/comparison.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 31167a9..e876aab 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -22,7 +22,7 @@ from astroquery.vizier import Vizier from tqdm import tqdm from prettytable import PrettyTable - +import traceback import requests import http import urllib3 @@ -1136,7 +1136,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n logger.debug(tabl) except: print ("something fishy in the calib mags table") - import traceback; logger.error(traceback.print_exc()) + logger.error(traceback.print_exc()) # Get the set of least variable stars to use as a comparison to calibrate the files (to eventually get the *ACTUAL* standards if asarray(calibStands).shape[0] == 0: @@ -1218,7 +1218,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n colTemp=delete(colTemp, calibStandsReject, axis=0) # Pre-colour Reference plot - + try: plt.cla() fig = plt.gcf() outplotx=colTemp[:,4] @@ -1249,7 +1249,9 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n plt.savefig(parentPath / str("results/CalibrationSanityPlot_PreColour_Reference.eps")) logger.info('Estimated Colour Slope in Reference Frame: ' + str(m)) - + except: + print ("Couldn't make colour plot") + logger.error(traceback.print_exc()) # MAKE ALL PRE-COLOUR PLOTS if colourdetect == True and colourTerm == 0.0: From 6506490dd23817b26a35e7829d067760b201ebb4 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 07:47:37 +1100 Subject: [PATCH 04/28] smaller glitches in multiprocessing development --- astrosource/analyse.py | 13 ++++++++++--- astrosource/comparison.py | 14 ++++++++++++-- astrosource/plots.py | 19 ++++++++++++------- 3 files changed, 34 insertions(+), 12 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index a8db8df..a317a02 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -14,7 +14,7 @@ #import traceback import logging from multiprocessing import Pool, cpu_count - +import traceback #from astrosource.utils import photometry_files_to_array, AstrosourceException from astrosource.utils import AstrosourceException from astrosource.plots import plot_variability @@ -394,7 +394,9 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, #breakpoint() - potentialVariables = np.delete(starVar, outliers, axis=0) + #potentialVariables = np.delete(starVar, outliers, axis=0) + + potentialVariables=starVar[outliers] meanMags = starVar[:,2] variations = starVar[:,3] @@ -436,7 +438,12 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, else: savetxt(parentPath / "results/potentialVariables.csv", potentialVariables , delimiter=",", fmt='%0.8f', header='RA,DEC,DiffMag,Variability') - plot_variability(outputVariableHolder, potentialVariables, parentPath, compFile) + try: + plot_variability(outputVariableHolder, potentialVariables, parentPath, compFile) + except: + print ("MTF hunting this bug") + logger.error(traceback.print_exc()) + breakpoint() plt.cla() fig, ax = plt.subplots(figsize =(10, 7)) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index e876aab..813e9f7 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -924,8 +924,18 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n elif cat_name == 'SkyMapper' and noskymapper==True: logger.info("Skipping Skymapper") else: - - coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd) + try: + coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd) + except Exception as e: + + if 'Could not find RA' in str(e): + raise AstrosourceException("Empty catalogue produced from catalogue call") + #pass + else: + print (e) + print ("MTF hunting this bug") + print(traceback.print_exc()) + breakpoint() # If no results try next catalogue if len(coords.ra) == 0: coords=[] diff --git a/astrosource/plots.py b/astrosource/plots.py index 686c494..016da42 100644 --- a/astrosource/plots.py +++ b/astrosource/plots.py @@ -9,6 +9,7 @@ from astropy.units import degree, arcsecond import logging import numpy as np +import traceback #from astrosource.utils import photometry_files_to_array, AstrosourceException from astrosource.utils import AstrosourceException @@ -94,13 +95,17 @@ def plot_variability(output, variableID, parentPath, compFile): if (parentPath / 'results/calibCompsUsed.csv').exists(): logger.debug("Calibrated") calibCompFile=genfromtxt(parentPath / 'results/calibCompsUsed.csv', dtype=float, delimiter=',') - - if len(calibCompFile) == 5: - calibCompSkyCoord = SkyCoord(calibCompFile[0],calibCompFile[1], frame='icrs', unit=degree) - calibnumber=1 - else: - calibCompSkyCoord = SkyCoord(calibCompFile[:,0],calibCompFile[:,1], frame='icrs', unit=degree) - calibnumber=len(calibCompSkyCoord) + try: + if len(calibCompFile) == 5 and calibCompFile.size == 5: + calibCompSkyCoord = SkyCoord(calibCompFile[0],calibCompFile[1], frame='icrs', unit=degree) + calibnumber=1 + else: + calibCompSkyCoord = SkyCoord(calibCompFile[:,0],calibCompFile[:,1], frame='icrs', unit=degree) + calibnumber=len(calibCompSkyCoord) + except: + print ("MTF FIXING THIS ERROR") + logger.error(traceback.print_exc()) + breakpoint() calibCompExist=True calibCompStarPlot = [] From 3c501c026e622c4a2c361837b9bac9f853f2ec50 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 08:33:55 +1100 Subject: [PATCH 05/28] better memory management for huge datasets --- astrosource/analyse.py | 33 ++++++++++++++++++++++++++++----- astrosource/comparison.py | 7 ++++--- astrosource/periodic.py | 2 +- 3 files changed, 33 insertions(+), 9 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index a317a02..a848b78 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -13,7 +13,7 @@ from tqdm import tqdm #import traceback import logging -from multiprocessing import Pool, cpu_count +from multiprocessing import Pool, cpu_count, shared_memory import traceback #from astrosource.utils import photometry_files_to_array, AstrosourceException from astrosource.utils import AstrosourceException @@ -25,7 +25,7 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.colors as colors - +from functools import partial logger = logging.getLogger('astrosource') @@ -58,7 +58,14 @@ def get_total_counts(photFileArray, compFile, loopLength, photCoords): #logger.debug(allCountsArray) return allCountsArray -def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): +#def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): +def process_varsearch_target(target, photFileArray_shape, photFileArray_dtype, shm_name, allCountsArray, matchRadius, minimumNoOfObs): + + # Attach to the shared memory for photFileArray + existing_shm = shared_memory.SharedMemory(name=shm_name) + photFileArray = np.ndarray(photFileArray_shape, dtype=photFileArray_dtype, buffer=existing_shm.buf) + + diffMagHolder = [] for allcountscount, photFile in enumerate(photFileArray): @@ -84,6 +91,8 @@ def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, if diffMagHolder.size == sizeBefore: break + existing_shm.close() + # Append to output if sufficient observations are available if diffMagHolder.size > minimumNoOfObs: return [target[0], target[1], np.median(diffMagHolder), np.std(diffMagHolder), diffMagHolder.size] @@ -91,11 +100,21 @@ def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): + # Create shared memory for photFileArray + # As this can lead to gigantic RAM use for large datasets + shm = shared_memory.SharedMemory(create=True, size=photFileArray.nbytes) + shared_photFileArray = np.ndarray(photFileArray.shape, dtype=photFileArray.dtype, buffer=shm.buf) + shared_photFileArray[:] = photFileArray[:] + + # Partial function to pass shared arguments - from functools import partial + worker = partial( process_varsearch_target, - photFileArray=photFileArray, + #photFileArray=photFileArray, + photFileArray_shape=photFileArray.shape, + photFileArray_dtype=photFileArray.dtype, + shm_name=shm.name, allCountsArray=allCountsArray, matchRadius=matchRadius, minimumNoOfObs=minimumNoOfObs, @@ -105,6 +124,10 @@ def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCoun with Pool(processes=max([cpu_count()-1,1])) as pool: results = pool.map(worker, targetFile) + # Clean up shared memory + shm.close() + shm.unlink() + # Filter out None results return [res for res in results if res is not None] diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 813e9f7..a0ba96f 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -338,9 +338,10 @@ def ensemble_comparisons(photFileArray, compFile, parentPath, photSkyCoord): # return comp_diff_mag, instr_mag def process_phot_file_for_variation(args): - q, matchCoord, photSkyCoord, photFileArray, fileCount = args + #q, matchCoord, photSkyCoord, photFileArray, fileCount = args + q, matchCoord, photSkyCoord, photFile, fileCount = args idx, _, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - photFile = photFileArray[q] + #photFile = photFileArray[q] compDiffMags = 2.5 * np.log10(photFile[idx, 4] / fileCount[q]) instrMags = -2.5 * np.log10(photFile[idx, 4]) return compDiffMags, instrMags @@ -416,7 +417,7 @@ def calculate_comparison_variation(compFile, photFileArray, fileCount, parentPat # Case for single comparison star matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree) - args = [(q, matchCoord, photSkyCoord, photFileArray, fileCount) for q in range(len(photFileArray))] + args = [(q, matchCoord, photSkyCoord, photFileArray[q], fileCount) for q in range(len(photFileArray))] with Pool(processes=max([cpu_count()-1,1])) as pool: results = pool.map(process_phot_file_for_variation, args) diff --git a/astrosource/periodic.py b/astrosource/periodic.py index 96e4801..54268c4 100644 --- a/astrosource/periodic.py +++ b/astrosource/periodic.py @@ -19,7 +19,7 @@ from astropy.timeseries import LombScargle logger = logging.getLogger('astrosource') -NCPUS = 1 +#NCPUS = 1 # Note that the functions that calculate the ANOVA periodograms have been adapted from the astrobase codeset # These are aov_theta, resort_by_time, get_frequency_grid, sigclip_magseries, phase_magseries, aov_periodfind, phase_magseries_with_errs, aovhm_theta, aovhm_periodfind From 1abf758b408ee1686a6e74d56e4ea27d2d4afea4 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 09:32:45 +1100 Subject: [PATCH 06/28] windows needs __main__ to multiprocess --- astrosource/analyse.py | 6 +++++- astrosource/comparison.py | 6 +++++- astrosource/identify.py | 6 +++++- astrosource/periodic.py | 4 ++++ astrosource/plots.py | 4 ++++ 5 files changed, 23 insertions(+), 3 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index a848b78..a792447 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -810,4 +810,8 @@ def calibrated_photometry(paths, photometrydata, colourterm, colourerror, colour logger.info("then consider providing an appropriate colour for this filter using the --targetcolour option") logger.info("as well as an appropriate colour term for this filter (using --colourdetect or --colourterm).") - return pdata \ No newline at end of file + return pdata + +# Needed for windows to multiprocess appropriately +if __name__ == "__main__": + pass \ No newline at end of file diff --git a/astrosource/comparison.py b/astrosource/comparison.py index a0ba96f..5c572c8 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -1941,4 +1941,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n savetxt(parentPath / "results/calibCompsUsed.csv", compFile, delimiter=",", fmt='%0.8f',header='RA,Dec,Variability') sys.stdout.write('\n') - return colourTerm, colourError, compFile \ No newline at end of file + return colourTerm, colourError, compFile + +# Needed for windows to multiprocess appropriately +if __name__ == "__main__": + pass \ No newline at end of file diff --git a/astrosource/identify.py b/astrosource/identify.py index 16ea204..ebe94f9 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -1185,4 +1185,8 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym f.write(str(filename) +"\n") - return fileList, outputComps, photFileHolder, photCoords \ No newline at end of file + return fileList, outputComps, photFileHolder, photCoords + +# Needed for windows to multiprocess appropriately +if __name__ == "__main__": + pass \ No newline at end of file diff --git a/astrosource/periodic.py b/astrosource/periodic.py index 54268c4..c03a568 100644 --- a/astrosource/periodic.py +++ b/astrosource/periodic.py @@ -2364,3 +2364,7 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1 # return pdm["distance_minperiod"] # else: return 0.0 + +# Needed for windows to multiprocess appropriately +if __name__ == "__main__": + pass \ No newline at end of file diff --git a/astrosource/plots.py b/astrosource/plots.py index 016da42..444f1b5 100644 --- a/astrosource/plots.py +++ b/astrosource/plots.py @@ -373,3 +373,7 @@ def phased_plots(paths, filterCode, targets, period, phaseShift): return + +# Needed for windows to multiprocess appropriately +if __name__ == "__main__": + pass From 734dc2ef82e1ab9f66a0a71e1363b76936a20ee3 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 10:32:51 +1100 Subject: [PATCH 07/28] actual bug fix. Also windows multiprocessing workaround. --- astrosource/analyse.py | 71 ++++++++++++++++++++++++++++++++++++++---- astrosource/utils.py | 2 +- 2 files changed, 66 insertions(+), 7 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index a792447..084e842 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -13,6 +13,7 @@ from tqdm import tqdm #import traceback import logging +import multiprocessing as mp from multiprocessing import Pool, cpu_count, shared_memory import traceback #from astrosource.utils import photometry_files_to_array, AstrosourceException @@ -28,6 +29,16 @@ from functools import partial logger = logging.getLogger('astrosource') +import platform +# import platform +# # Check the operating system +# if platform.system() == "Windows": +# # Use 'forkserver' for Windows +# mp.set_start_method("forkserver", force=True) +# else: +# # Use the default method for other OS +# mp.set_start_method("fork", force=True) + def get_total_counts(photFileArray, compFile, loopLength, photCoords): @@ -98,6 +109,39 @@ def process_varsearch_target(target, photFileArray_shape, photFileArray_dtype, s return [target[0], target[1], np.median(diffMagHolder), np.std(diffMagHolder), diffMagHolder.size] return None +def process_varsearch_target_singleprocess(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): + + diffMagHolder = [] + + for allcountscount, photFile in enumerate(photFileArray): + # Efficiently calculate the closest match using numpy + idx = (np.abs(photFile[:, 0] - target[0]) + np.abs(photFile[:, 1] - target[1])).argmin() + d2d = np.sqrt((photFile[idx, 0] - target[0])**2 + (photFile[idx, 1] - target[1])**2) * 3600 + multTemp = -2.5 * np.log10(photFile[idx, 11] / allCountsArray[allcountscount][0]) + + if d2d < matchRadius and not np.isinf(multTemp): + diffMagHolder.append(multTemp) + + # Remove major outliers + diffMagHolder = np.array(diffMagHolder) + while True: + stdVar = np.std(diffMagHolder) + avgVar = np.mean(diffMagHolder) + sizeBefore = diffMagHolder.size + + # Mask outliers + mask = (diffMagHolder <= avgVar + 4 * stdVar) & (diffMagHolder >= avgVar - 4 * stdVar) + diffMagHolder = diffMagHolder[mask] + + if diffMagHolder.size == sizeBefore: + break + + # Append to output if sufficient observations are available + if diffMagHolder.size > minimumNoOfObs: + return [target[0], target[1], np.median(diffMagHolder), np.std(diffMagHolder), diffMagHolder.size] + + return None + def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): # Create shared memory for photFileArray @@ -394,9 +438,24 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, # outputVariableHolder.append( [target[0],target[1],median(diffMagHolder), std(diffMagHolder), diffMagHolder.shape[0]]) # process_varsearch_target - outputVariableHolder = process_varsearch_targets_multiprocessing( - targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs - ) + + + # Hack to get windows to not multiprocess until I figure out how to do it. + if platform.system() == "Windows": + + + outputVariableHolder=[] + for target in targetFile: + result= process_varsearch_target_singleprocess(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs) + if result == None: + pass + else: + outputVariableHolder.append(result) + else: + + outputVariableHolder = process_varsearch_targets_multiprocessing( + targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs + ) print ("Star Variability done in " + str(time.time()-taketime)) @@ -528,15 +587,15 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file starDistanceRejCount=0 logger.debug("****************************") logger.debug("Processing Variable {}".format(q+1)) - if int(len(targets)) == 4 and targets.size==4: + if int(len(targets)) == 5 and targets.size==5: logger.debug("RA {}".format(targets[0])) else: logger.debug("RA {}".format(targets[q][0])) - if int(len(targets)) == 4 and targets.size==4: + if int(len(targets)) == 5 and targets.size==5: logger.debug("Dec {}".format(targets[1])) else: logger.debug("Dec {}".format(targets[q][1])) - if int(len(targets)) == 4 and targets.size==4: + if int(len(targets)) == 5 and targets.size==5: #varCoord = SkyCoord(targets[0],(targets[1]), frame='icrs', unit=degree) # Need to remove target stars from consideration singleCoordRA=targets[0] singleCoordDEC=targets[1] diff --git a/astrosource/utils.py b/astrosource/utils.py index 33b1852..3cf5463 100644 --- a/astrosource/utils.py +++ b/astrosource/utils.py @@ -95,7 +95,7 @@ def get_targets(targetfile): # Remove any nan rows from targets targetRejecter=[] - if not (targets.shape[0] == 4 and targets.size == 4): + if not (targets.shape[0] == 5 and targets.size == 5): for z in range(targets.shape[0]): if isnan(targets[z][0]): targetRejecter.append(z) From 9cd81b67bbbbbb0939397633e1792429ff543cdb Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 14:46:34 +1100 Subject: [PATCH 08/28] Update analyse.py --- astrosource/analyse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index 084e842..c8a55ab 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -782,7 +782,7 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file except ValueError: #raise AstrosourceException("No target stars were detected in your dataset. Check your input target(s) RA/Dec") - import traceback; logger.error(traceback.print_exc()) + logger.error(traceback.print_exc()) logger.error("This target star was not detected in your dataset. Check your input target(s) RA/Dec") #logger.info("Rejected Stdev Measurements: : {}".format(stdevReject)) #logger.error("Rejected Error Measurements: : {}".format(starErrorRejCount)) From 5ce35e981d9c7ab8b44e49a670ba2c03f85fbc2d Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Mon, 23 Dec 2024 16:13:30 +1100 Subject: [PATCH 09/28] Some better error handling --- astrosource/periodic.py | 301 +++++++++++++++++++++------------------- 1 file changed, 162 insertions(+), 139 deletions(-) diff --git a/astrosource/periodic.py b/astrosource/periodic.py index c03a568..713d212 100644 --- a/astrosource/periodic.py +++ b/astrosource/periodic.py @@ -106,162 +106,173 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p plt.figure(figsize=(15, 5)) - logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])) - logger.debug("Distance method error: " + str(pdm["distance_error"])) - + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: f.write("Variable : "+str(variableName) +"\n") - f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n") - f.write("Distance method error: " + str(pdm["distance_error"])+"\n") - - plt.plot(pdm["periodguess_array"], pdm["distance_results"]) - plt.gca().invert_yaxis() - plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) - plt.xlabel(r"Trial Period") - plt.ylabel(r"Likelihood of Period") - plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") - plt.clf() - - if (varData.size > 3): - phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1 - - plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') - plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') - plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') - plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') + + + try: + logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])) + logger.debug("Distance method error: " + str(pdm["distance_error"])) + pdmfailed=False + except: + logger.debug("Distance Method Failed") + pdmfailed=True + + if not pdmfailed: + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("Variable : "+str(variableName) +"\n") + f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n") + f.write("Distance method error: " + str(pdm["distance_error"])+"\n") + + plt.plot(pdm["periodguess_array"], pdm["distance_results"]) plt.gca().invert_yaxis() - plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) - plt.xlabel(r"Phase ($\phi$)") - plt.ylabel(f"Differential {filterCode} Magnitude") - plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png") + plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") plt.clf() - - if calibFile.exists(): - if (calibData.size > 3): - phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1 - plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') - plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') - plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') - plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + + if (varData.size > 3): + phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1 + + plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') + plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') plt.gca().invert_yaxis() plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) plt.xlabel(r"Phase ($\phi$)") - plt.ylabel(f"Calibrated {filterCode} Magnitude") - plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png") + plt.ylabel(f"Differential {filterCode} Magnitude") + plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png") plt.clf() + if calibFile.exists(): + if (calibData.size > 3): + phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1 + plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') + plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(f"Calibrated {filterCode} Magnitude") + plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png") + plt.clf() + + tempPeriodCatOut=[] + for g in range(len(calibData[:,0])): + tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + if (varData.size > 3): + tempPeriodCatOut=[] - for g in range(len(calibData[:,0])): - tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]]) + for g in range(len(phaseTest)): + tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - if (varData.size > 3): - - tempPeriodCatOut=[] - for g in range(len(phaseTest)): - tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) - tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + tempPeriodCatOut=[] + for g in range(len(varData[:,0])): + tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - tempPeriodCatOut=[] - for g in range(len(varData[:,0])): - tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]]) - tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4): - logger.info("No PDM results due to lack of datapoint coverage") - else: - logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])) - phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1 - logger.debug("PDM method error: " + str(pdm["stdev_error"])) - - with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n") - f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n") - - - plt.plot(pdm["periodguess_array"], pdm["stdev_results"]) - plt.gca().invert_yaxis() - plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) - plt.xlabel(r"Trial Period") - plt.ylabel(r"Likelihood of Period") - plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png") - - plt.clf() - - if (varData.size > 3): - plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') - plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') - plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') - plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') + if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4): + logger.info("No PDM results due to lack of datapoint coverage") + else: + logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])) + phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1 + logger.debug("PDM method error: " + str(pdm["stdev_error"])) + + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n") + f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n") + + + plt.plot(pdm["periodguess_array"], pdm["stdev_results"]) plt.gca().invert_yaxis() - plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) - plt.xlabel(r"Phase ($\phi$)") - plt.ylabel(r"Differential " + str(filterCode) + " Magnitude") - plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png") + plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png") + plt.clf() - - if calibFile.exists(): - if (calibData.size > 3): - phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1 - plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') - plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') - plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') - plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + + if (varData.size > 3): + plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') + plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') + plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') plt.gca().invert_yaxis() plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) plt.xlabel(r"Phase ($\phi$)") - plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude") - plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png") + plt.ylabel(r"Differential " + str(filterCode) + " Magnitude") + plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png") plt.clf() + if calibFile.exists(): + if (calibData.size > 3): + phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1 + plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') + plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') + plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') + plt.gca().invert_yaxis() + plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) + plt.xlabel(r"Phase ($\phi$)") + plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude") + plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png") + plt.clf() + + tempPeriodCatOut=[] + for g in range(len(calibData[:,0])): + tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + if (varData.size > 3): + tempPeriodCatOut=[] - for g in range(len(calibData[:,0])): - tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]]) + for g in range(len(phaseTest)): + tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - if (varData.size > 3): - - tempPeriodCatOut=[] - for g in range(len(phaseTest)): - tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) - tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') + + tempPeriodCatOut=[] + for g in range(len(varData[:,0])): + tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]]) + tempPeriodCatOut=asarray(tempPeriodCatOut) + savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - tempPeriodCatOut=[] - for g in range(len(varData[:,0])): - tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]]) - tempPeriodCatOut=asarray(tempPeriodCatOut) - savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # Plot publication plots - - plt.figure(figsize=(5, 3)) - - plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5) - plt.gca().invert_yaxis() - plt.xlabel(r"Trial Period") - plt.ylabel(r"Likelihood of Period") - plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) - plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300) - plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps") - - plt.clf() - - plt.figure(figsize=(5, 3)) - - plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5) - plt.gca().invert_yaxis() - plt.xlabel(r"Trial Period") - plt.ylabel(r"Likelihood of Period") - plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) - plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300) - plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps") - - plt.clf() + # Plot publication plots + + plt.figure(figsize=(5, 3)) + + plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5) + plt.gca().invert_yaxis() + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300) + plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps") + + plt.clf() + + plt.figure(figsize=(5, 3)) + + plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5) + plt.gca().invert_yaxis() + plt.xlabel(r"Trial Period") + plt.ylabel(r"Likelihood of Period") + plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) + plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300) + plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps") + + plt.clf() # ANOVA @@ -319,9 +330,12 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p if (varData.size > 3): aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])) - with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n") + try: + logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])) + with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: + f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n") + except: + logger.debug("Harmonic Anova Method Estimate Failed") # LOMB SCARGLE for nts in range(2): @@ -1586,8 +1600,11 @@ def aovhm_periodfind(times, frequencies = nparange(startf, endf, stepsize) # map to parallel workers - if (not nworkers) or (nworkers > NCPUS): - nworkers = NCPUS + + + # if (not nworkers) or (nworkers > NCPUS): + # nworkers = NCPUS + #nworkers=1 # renormalize the working mags to zero and scale them so that the # variance = 1 for use with our LSP functions @@ -1849,7 +1866,13 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu distance_results = [] stdev_results = [] + + if len(varData) < 5: + return + + #try: (julian_dates, fluxes) = (varData[:,0],varData[:,1]) + normalizedFluxes = normalize(fluxes) for r in range(periodsteps): From 31d9e32e9fef2c6e26f1da2022f425bb928c91cf Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Tue, 1 Apr 2025 19:18:37 +1100 Subject: [PATCH 10/28] udpate to skymapper dr4 --- astrosource/analyse.py | 3 +++ astrosource/comparison.py | 2 +- astrosource/identify.py | 30 ++++++++++++++++++++++-------- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index c8a55ab..12c6721 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -576,6 +576,8 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file allcountscount=0 + #breakpoint() + if len(targets)== 4 and targets.size == 4: loopLength=1 else: @@ -784,6 +786,7 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file #raise AstrosourceException("No target stars were detected in your dataset. Check your input target(s) RA/Dec") logger.error(traceback.print_exc()) logger.error("This target star was not detected in your dataset. Check your input target(s) RA/Dec") + #breakpoint() #logger.info("Rejected Stdev Measurements: : {}".format(stdevReject)) #logger.error("Rejected Error Measurements: : {}".format(starErrorRejCount)) #logger.error("Rejected Distance Measurements: : {}".format(starDistanceRejCount)) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 5c572c8..86d1ec1 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -625,7 +625,7 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): TABLES = {'APASS':'II/336/apass9', 'SDSS' :'V/154/sdss16', 'PanSTARRS' : 'II/349/ps1', - 'SkyMapper' : 'II/358/smss' + 'SkyMapper' : 'II/379/smssdr4' } tbname = TABLES.get(cat_name, None) diff --git a/astrosource/identify.py b/astrosource/identify.py index ebe94f9..c54e417 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -20,6 +20,7 @@ from prettytable import PrettyTable from multiprocessing import Pool,cpu_count +import platform from astrosource.utils import AstrosourceException from astrosource.comparison import catalogue_call @@ -470,6 +471,8 @@ def gather_files(paths, filelist=None, filetype="fz", bjd=False, ignoreedgefract file1=open(paths['parent'] / "filterCode","wb") pickle.dump(filterCode, file1) file1.close + + #breakpoint() return phot_list, filterCode, photFileHolder, photSkyCoord @@ -478,6 +481,7 @@ def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistanc imgRejFlag = 0 photReject = [] + # Check minimum stars requirement if referenceFrame.shape[0] < mincompstars: return None, referenceFrame, photReject @@ -522,15 +526,24 @@ def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistanc def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): - with Pool(processes=max([cpu_count()-1,1])) as pool: - # Prepare arguments - args = [ - (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger) - for index in range(len(photFileHolder)) - ] - # Process files in parallel - results = pool.starmap(process_phot_file, args) + + # Hack to get windows to not multiprocess until I figure out how to do it. + if platform.system() == "Windows": + results=[] + for index in range(len(photFileHolder)): + results.append(process_phot_file(index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger)) + + else: + with Pool(processes=max([cpu_count()-1,1])) as pool: + # Prepare arguments + args = [ + (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger) + for index in range(len(photFileHolder)) + ] + + # Process files in parallel + results = pool.starmap(process_phot_file, args) # Collect results updated_referenceFrame = referenceFrame @@ -541,6 +554,7 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen updated_referenceFrame = refFrameUpdate photReject.extend(rejects) + return updated_referenceFrame, photReject From e27e8826115b78a6bc26bd3ac97d04c64a350e62 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Fri, 4 Apr 2025 11:35:55 +1100 Subject: [PATCH 11/28] swap sdss and panstarss. Also WCS fail rejection --- astrosource/identify.py | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/astrosource/identify.py b/astrosource/identify.py index c54e417..2aca2dd 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -213,14 +213,24 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr for entry in range(len(photFile[:,0])): if photFile[entry,0] < 0: photFile[entry,0] = photFile[entry,0] + 360 + + + #remove odd zero entries photFile[:,0][photFile[:,0] == 0.0 ] = nan photFile[:,0][photFile[:,0] == 0.0 ] = nan photFile[:,1][photFile[:,1] == 0.0 ] = nan photFile[:,1][photFile[:,1] == 0.0 ] = nan + num_rows_with_nans = sum(isnan(photFile).any(axis=1)) + + if len(photFile) == num_rows_with_nans: + print(f"REJECT {infile}: Very likely no wcs fit - RA and Dec columns are zero.") + photFile=photFile[~isnan(photFile).any(axis=1)] + + #remove lowcounts rejectStars=where(photFile[:,4] < lowestcounts)[0] photFile=delete(photFile, rejectStars, axis=0) @@ -272,6 +282,12 @@ def process_convert_photometry_file(fn, racut, deccut, radiuscut, ignoreedgefrac ] # Remove zero and low-count entries + + #num_rows_with_nans = sum(isnan(photFile).any(axis=1)) + #print (len((photFile[:, 0] != 0) & (photFile[:, 1] != 0))) + if len((photFile[:, 0] != 0) & (photFile[:, 1] != 0)) == 0: + print(f"REJECTED {fn}: Very likely no wcs fit - RA and Dec columns are zero.") + photFile = photFile[(photFile[:, 4] > lowestcounts) & (photFile[:, 0] != 0) & (photFile[:, 1] != 0)] # Prepare output @@ -481,6 +497,7 @@ def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistanc imgRejFlag = 0 photReject = [] + #breakpoint() # Check minimum stars requirement if referenceFrame.shape[0] < mincompstars: @@ -526,7 +543,7 @@ def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistanc def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): - + # Hack to get windows to not multiprocess until I figure out how to do it. if platform.system() == "Windows": @@ -544,7 +561,7 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen # Process files in parallel results = pool.starmap(process_phot_file, args) - + # Collect results updated_referenceFrame = referenceFrame photReject = [] @@ -553,7 +570,8 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen if refFrameUpdate is not None: updated_referenceFrame = refFrameUpdate photReject.extend(rejects) - + + #breakpoint() return updated_referenceFrame, photReject @@ -850,6 +868,8 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger ) + #breakpoint() + # Remove files and Hold the photSkyCoords in memory photCoords=asarray(photCoords, dtype=object) @@ -1094,9 +1114,9 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym 'SkyMapper' : {'filter' : 'iPSF', 'error' : 'e_iPSF', 'colmatch' : 'rPSF', 'colerr' : 'e_rPSF', 'colname' : 'r-i', 'colrev' : '1'}, 'PanSTARRS': {'filter' : 'imag', 'error' : 'e_imag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-i', 'colrev' : '1'}, 'APASS' : {'filter' : 'i_mag', 'error' : 'e_i_mag', 'colmatch' : 'r_mag', 'colerr' : 'e_r_mag', 'colname' : 'r-i', 'colrev' : '1'}}, - 'zs' : {'PanSTARRS': {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}, + 'zs' : {'SDSS' : {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}, 'SkyMapper' : {'filter' : 'zPSF', 'error' : 'e_zPSF', 'colmatch' : 'rPSF', 'colerr' : 'e_rPSF', 'colname' : 'r-zs', 'colrev' : '1'}, - 'SDSS' : {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}}, + 'PanSTARRS': {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}}, } try: From ef18b3d5cf1bf3b3a2707f197d2e11550a523c57 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Tue, 27 May 2025 18:24:58 +1000 Subject: [PATCH 12/28] Vizier Server hotfix #1 --- astrosource/comparison.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 86d1ec1..011223d 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -521,7 +521,7 @@ def remove_stars_targets(parentPath, compFile, acceptDistance, targetFile, remov tableFound=False # Check VSX for any known variable stars and remove them from the list logger.info("Searching for known variable stars in VSX......") - vServers = ['vizier.cfa.harvard.edu','vizier.u-strasbg.fr'] + vServers = ['vizier.cds.unistra.fr','vizier.cfa.harvard.edu','vizier.u-strasbg.fr'] vS=0 try: while vS < len(vServers): @@ -1945,4 +1945,4 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n # Needed for windows to multiprocess appropriately if __name__ == "__main__": - pass \ No newline at end of file + pass From 28eef1142a991af710134e82b1acb3137e9d2002 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Tue, 27 May 2025 18:31:08 +1000 Subject: [PATCH 13/28] small patches --- astrosource/analyse.py | 6 +++++- astrosource/main.py | 7 +++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index 12c6721..ee1f74d 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -254,7 +254,11 @@ def fit_sigma_clipped_sigmoid(x, y, sigma=2, parentPath=''): # Fit sigmoid function to non-masked data with weights p0 = [max(y), 1, np.median(x)] # Initial guesses for L, k, x0 - popt, _ = curve_fit(sigmoid_func, x[~mask], y[~mask], p0=p0, sigma=1/weights[~mask]) + try: + popt, _ = curve_fit(sigmoid_func, x[~mask], y[~mask], p0=p0, sigma=1/weights[~mask]) + except: + logger.info("Sigmoid curve fit failed") + return [] # Calculate residuals and standard deviation y_fit = sigmoid_func(x, *popt) diff --git a/astrosource/main.py b/astrosource/main.py index 0b48688..e1c3913 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -2,7 +2,7 @@ import click import sys import logging - +import os from numpy import array, all from astropy.utils.exceptions import AstropyWarning, AstropyDeprecationWarning @@ -217,7 +217,10 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree ts.find_variables() if variablehunt == True: - targets = get_targets(parentPath / 'results/potentialVariables.csv') + if os.path.exists(parentPath / 'results/potentialVariables.csv'): + targets = get_targets(parentPath / 'results/potentialVariables.csv') + else: + targets=None From a5a3e4c1c42322d604e7c8cb72a0ccdc473ab447 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Wed, 28 May 2025 06:35:37 +1000 Subject: [PATCH 14/28] vizier hot fix #2 --- astrosource/comparison.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 011223d..60286d9 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -630,11 +630,11 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): tbname = TABLES.get(cat_name, None) kwargs = {'radius': str(1.5*radius)+' deg'} - kwargs['catalog'] = cat_name + kwargs['catalog'] = cat_name - vServers = ['vizier.u-strasbg.fr', + vServers = ['vizier.cds.unistra.fr','vizier.u-strasbg.fr']#, #'vizier.nao.ac.jp', - 'vizier.cfa.harvard.edu'] + #'vizier.cfa.harvard.edu'] #'vizier.iucaa.in', #'vizier.china-vo.org', #'vizier.inasan.ru', @@ -701,10 +701,15 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): # tableFound=True cycler=0 + vS=0 try: while vS < len(vServers): + v.VIZIER_SERVER=vServers[vS] - query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs) + try: + query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs) + except: + logger.info("Failed vizier query)") if str(query)=="Empty TableList": vS=vS+1 else: @@ -1143,7 +1148,10 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n logger.debug('Calibration Stars Identified below') tabl = PrettyTable() tabl.field_names = ["RA","Dec","MagDiff","Mag"] - tabl.add_rows(calibStands[:,0:4]) + #abl.add_rows(calibStands[:,0:4]) + rows = calibStands[:,0:4].tolist() + if rows: + tabl.add_rows(rows) logger.debug(tabl) except: print ("something fishy in the calib mags table") @@ -1945,4 +1953,4 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n # Needed for windows to multiprocess appropriately if __name__ == "__main__": - pass + pass \ No newline at end of file From 4243a89150de752bde7fb8e593f9b8198adea903 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Wed, 28 May 2025 06:46:09 +1000 Subject: [PATCH 15/28] Update README.md --- README.md | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index bb1b3b4..4d5567a 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,17 @@ pip install astrosource ### Install development version -If you need to install the development branch, download from [GitHub](https://github.com/zemogle/astrosource) and from the root of the repo, run: +The development version tends to be more cutting edge with hotfixes applied to latest problems and with new features that are new enough that we are still monitoring for bugs. Having said that, it is probably the better choice to use, just be aware that there may be bugs and if you find one, let us know! + +If you need to install the development branch, there are two methods: + +(1) Use pip to install the development branch + +```bash +pip install -U git+https://github.com/zemogle/astrosource@dev +``` + +(2) download from [GitHub](https://github.com/zemogle/astrosource) and from the root of the repo, run: ```bash cd astrosource From 330588ac2e48abd82a3652c1c68c4c630e24b761 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Wed, 28 May 2025 06:47:09 +1000 Subject: [PATCH 16/28] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 4d5567a..93f3b37 100644 --- a/README.md +++ b/README.md @@ -32,7 +32,7 @@ If you need to install the development branch, there are two methods: pip install -U git+https://github.com/zemogle/astrosource@dev ``` -(2) download from [GitHub](https://github.com/zemogle/astrosource) and from the root of the repo, run: +(2) download from [GitHub](https://github.com/zemogle/astrosource/tree/dev) and from the root of the repo, run: ```bash cd astrosource From 1eb8545c52ffcd607a62c551a3e52673bdd86027 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 1 Jun 2025 07:23:56 +1000 Subject: [PATCH 17/28] significant tidy --- astrosource/analyse.py | 125 +----------- astrosource/comparison.py | 317 ++--------------------------- astrosource/detrend.py | 49 +---- astrosource/identify.py | 407 ++++++-------------------------------- astrosource/main.py | 4 - astrosource/periodic.py | 390 +----------------------------------- astrosource/plots.py | 5 - 7 files changed, 88 insertions(+), 1209 deletions(-) diff --git a/astrosource/analyse.py b/astrosource/analyse.py index ee1f74d..a4e2700 100644 --- a/astrosource/analyse.py +++ b/astrosource/analyse.py @@ -30,15 +30,6 @@ logger = logging.getLogger('astrosource') import platform -# import platform -# # Check the operating system -# if platform.system() == "Windows": -# # Use 'forkserver' for Windows -# mp.set_start_method("forkserver", force=True) -# else: -# # Use the default method for other OS -# mp.set_start_method("fork", force=True) - def get_total_counts(photFileArray, compFile, loopLength, photCoords): @@ -50,7 +41,6 @@ def get_total_counts(photFileArray, compFile, loopLength, photCoords): for photFile in photFileArray: allCounts = 0.0 allCountsErr = 0.0 - #fileRaDec = SkyCoord(ra=photFile[:, 0]*degree, dec=photFile[:, 1]*degree) fileRaDec = photCoords[counter] counter=counter+1 #Array of comp measurements @@ -66,7 +56,6 @@ def get_total_counts(photFileArray, compFile, loopLength, photCoords): if (compFile.shape[0] == 5 and compFile.size == 5) or (compFile.shape[0] == 3 and compFile.size == 3): break allCountsArray.append([allCounts, allCountsErr]) - #logger.debug(allCountsArray) return allCountsArray #def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs): @@ -74,9 +63,7 @@ def process_varsearch_target(target, photFileArray_shape, photFileArray_dtype, s # Attach to the shared memory for photFileArray existing_shm = shared_memory.SharedMemory(name=shm_name) - photFileArray = np.ndarray(photFileArray_shape, dtype=photFileArray_dtype, buffer=existing_shm.buf) - - + photFileArray = np.ndarray(photFileArray_shape, dtype=photFileArray_dtype, buffer=existing_shm.buf) diffMagHolder = [] for allcountscount, photFile in enumerate(photFileArray): @@ -155,7 +142,6 @@ def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCoun worker = partial( process_varsearch_target, - #photFileArray=photFileArray, photFileArray_shape=photFileArray.shape, photFileArray_dtype=photFileArray.dtype, shm_name=shm.name, @@ -223,8 +209,6 @@ def fit_sigma_clipped_poly(x, y, order=3, sigma=2, parentPath=''): plt.savefig(parentPath / "results/polifitdiagram.png") - #breakpoint() - def sigmoid_func(x, L, k, x0): """Sigmoid function: y = L / (1 + exp(-k * (x - x0))) + 0.01""" return L / (1 + np.exp(-k * (x - x0))) + 0.01 @@ -358,7 +342,7 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, photFileArray=photFileHolder - photFileCoords=photCoords + #photFileCoords=photCoords # allocate minimum images to detect minimumNoOfObs=int(varsearchminimages*len(fileList)) @@ -381,7 +365,7 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, compFile = genfromtxt(parentPath / "results/compsUsed.csv", dtype=float, delimiter=',') logger.debug("Stable Comparison Candidates below variability threshold") - outputPhot = [] + #outputPhot = [] # Get total counts for each file allCountsArray = get_total_counts(photFileArray, compFile, loopLength=compFile.shape[0], photCoords=photCoords) @@ -400,53 +384,15 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, ## NEED TO REMOVE COMPARISON STARS FROM TARGETLIST - allcountscount=0 + #allcountscount=0 # For each variable calculate the variability outputVariableHolder=[] logger.info("Measuring variability of stars...... ") taketime=time.time() - # for target in targetFile: - # diffMagHolder=[] - # allcountscount=0 - - - # for photFile in photFileArray: - # # A bit rougher than using SkyCoord, but way faster - # # The amount of calculations is too slow for SkyCoord - # idx=(np.abs(photFile[:,0] - target[0]) + np.abs(photFile[:,1] - target[1])).argmin() - # d2d=pow(pow(photFile[idx,0] - target[0],2) + pow(photFile[idx,1] - target[1],2),0.5) * 3600 - # multTemp=(multiply(-2.5,log10(divide(photFile[idx][11],allCountsArray[allcountscount][0])))) - # if less(d2d, matchRadius) and (multTemp != inf) : - # diffMagHolder=append(diffMagHolder,multTemp) - # allcountscount=add(allcountscount,1) - - # ## REMOVE MAJOR OUTLIERS FROM CONSIDERATION - # diffMagHolder=np.array(diffMagHolder) - # while True: - # stdVar=std(diffMagHolder) - # avgVar=average(diffMagHolder) - - # sizeBefore=diffMagHolder.shape[0] - # #print (sizeBefore) - - # diffMagHolder[diffMagHolder > avgVar+(4*stdVar) ] = np.nan - # diffMagHolder[diffMagHolder < avgVar-(4*stdVar) ] = np.nan - # diffMagHolder=diffMagHolder[~np.isnan(diffMagHolder)] - - # if diffMagHolder.shape[0] == sizeBefore: - # break - - # if (diffMagHolder.shape[0] > minimumNoOfObs): - # outputVariableHolder.append( [target[0],target[1],median(diffMagHolder), std(diffMagHolder), diffMagHolder.shape[0]]) - -# process_varsearch_target - - - # Hack to get windows to not multiprocess until I figure out how to do it. - if platform.system() == "Windows": - + # Hack to get windows to not multiprocess until a good solution is found. + if platform.system() == "Windows": outputVariableHolder=[] for target in targetFile: @@ -465,23 +411,12 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, savetxt(parentPath / "results/starVariability.csv", outputVariableHolder, delimiter=",", fmt='%0.8f', header='RA,DEC,DiffMag,Variability,No_of_images_used') - - - - - - #breakpoint() - ## Routine that actually pops out potential variables. starVar = np.asarray(outputVariableHolder) outliers=fit_sigma_clipped_sigmoid(starVar[:,2],starVar[:,3], parentPath=parentPath) - #breakpoint() - - #potentialVariables = np.delete(starVar, outliers, axis=0) - potentialVariables=starVar[outliers] meanMags = starVar[:,2] @@ -492,33 +427,6 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, xbins = np.arange(np.min(meanMags), np.max(meanMags), xStepSize) ybins = np.arange(np.min(variations), np.max(variations), yStepSize) - # #split it into one array and identify variables in bins with centre - # variationsByMag=[] - # potentialVariables=[] - # for xbinner in range(len (xbins)): - - # starsWithin=[] - # for q in range(len(meanMags)): - # if meanMags[q] >= xbins[xbinner] and meanMags[q] < xbins[xbinner]+xStepSize: - # starsWithin.append(variations[q]) - - # meanStarsWithin= (np.mean(starsWithin)) - # stdStarsWithin= (np.std(starsWithin)) - # variationsByMag.append([xbins[xbinner]+0.5*xStepSize,meanStarsWithin,stdStarsWithin]) - - # # At this point extract RA and Dec of stars that may be variable - # for q in range(len(starVar[:,2])): - # if starVar[q,2] >= xbins[xbinner] and starVar[q,2] < xbins[xbinner]+xStepSize: - # if varsearchglobalstdev != -99.9: - # if starVar[q,3] > varsearchglobalstdev : - # potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) - # elif starVar[q,3] > (meanStarsWithin + varsearchstdev*stdStarsWithin): - # #print (starVar[q,3]) - # potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]]) - - # potentialVariables=np.array(potentialVariables) - # logger.debug("Potential Variables Identified: " + str(potentialVariables.shape[0])) - if potentialVariables.shape[0] == 0: logger.info("No Potential Variables identified in this set of data using the parameters requested.") else: @@ -580,8 +488,6 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file allcountscount=0 - #breakpoint() - if len(targets)== 4 and targets.size == 4: loopLength=1 else: @@ -616,20 +522,11 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file allcountscount=0 for imgs, photFile in enumerate(tqdm(photFileArray)): - #fileRaDec = photCoordsFile[imgs] - #idx, d2d, _ = varCoord.match_to_catalog_sky(fileRaDec) - #breakpoint() idx=(np.abs(photFile[:,0] - singleCoordRA) + np.abs(photFile[:,1] - singleCoordDEC)).argmin() d2d=pow(pow(photFile[idx,0] - singleCoordRA,2) + pow(photFile[idx,1] - singleCoordDEC,2),0.5) * 3600 - #print (d2d) - starRejected=0 - if (less(d2d, targetRadius)): - #magErrVar = 1.0857 * (photFile[idx][5]/photFile[idx][4]) - - - + if (less(d2d, targetRadius)): # If the file hasn't been calibrated, then it still contains the countrate in it. # So convert these to mags, otherwise use the calibrated error. if photFile[idx][4] > 50: @@ -725,8 +622,6 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file fileCount.append(allCountsArray[allcountscount][0]) allcountscount=allcountscount+1 - #breakpoint() - # Check for dud images imageReject=[] for j in range(asarray(outputPhot).shape[0]): @@ -757,7 +652,6 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file break # Reject by outsized error elimination - while True: errorsArray=[] for j in range(asarray(outputPhot).shape[0]): @@ -787,13 +681,8 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file outputPhot=delete(outputPhot, starReject, axis=0) except ValueError: - #raise AstrosourceException("No target stars were detected in your dataset. Check your input target(s) RA/Dec") logger.error(traceback.print_exc()) logger.error("This target star was not detected in your dataset. Check your input target(s) RA/Dec") - #breakpoint() - #logger.info("Rejected Stdev Measurements: : {}".format(stdevReject)) - #logger.error("Rejected Error Measurements: : {}".format(starErrorRejCount)) - #logger.error("Rejected Distance Measurements: : {}".format(starDistanceRejCount)) # Add calibration columns outputPhot= np.c_[outputPhot, np.ones(outputPhot.shape[0]),np.ones(outputPhot.shape[0])] diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 60286d9..b6a094c 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -33,8 +33,6 @@ logger = logging.getLogger('astrosource') - - def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None, photSkyCoord=None, matchRadius=1.45): ''' Find stable comparison stars for the target photometry @@ -54,7 +52,6 @@ def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None, if fileList==None: fileList= genfromtxt(parentPath / 'results/usedImages.txt', dtype=str) - #print (fileList) usedImages=[] q=0 imageRemove=[] @@ -73,7 +70,6 @@ def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None, logger.debug('The Following file does not contain one or more of the provided comparison stars. It has been removed') logger.debug(file) logger.debug('**********************') - #fileList.remove(file) imageRemove.append(q) else: usedImages.append(file) @@ -282,28 +278,6 @@ def read_data_files(parentPath, fileList): def ensemble_comparisons(photFileArray, compFile, parentPath, photSkyCoord): - # # get rid of dumb for loop - # fileCount = [] - # q=0 - # for photFile in photFileArray: - # allCounts = 0.0 - - # if compFile.size ==2 and compFile.shape[0]==2: - # matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) - # else: - # matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) - - # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - - # allCounts = add(allCounts,sum(photFile[idx,4])) - - # fileCount.append(allCounts) - # q=q+1 - - # logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(np.sum(np.array(fileCount)))) - # return fileCount - - #degree = np.deg2rad(1) # Degree to radians conversion fileCount = np.zeros(len(photFileArray)) # Pre-allocate output array # Prepare the `matchCoord` SkyCoord object @@ -323,20 +297,6 @@ def ensemble_comparisons(photFileArray, compFile, parentPath, photSkyCoord): logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(total_counts)) return fileCount - -# def process_phot_file_for_variation(args): -# """ -# Process a single photFileArray element. -# """ -# q, matchCoord, photSkyCoord, photFileArray, fileCount = args - -# photFile = photFileArray[q] -# idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) -# comp_diff_mag = 2.5 * np.log10(photFile[idx][4] / fileCount[q]) -# instr_mag = -2.5 * np.log10(photFile[idx][4]) - -# return comp_diff_mag, instr_mag - def process_phot_file_for_variation(args): #q, matchCoord, photSkyCoord, photFileArray, fileCount = args q, matchCoord, photSkyCoord, photFile, fileCount = args @@ -347,68 +307,7 @@ def process_phot_file_for_variation(args): return compDiffMags, instrMags def calculate_comparison_variation(compFile, photFileArray, fileCount, parentPath, photSkyCoord): - # stdCompStar=[] - # sortStars=[] - # logger.debug("Calculating Variation in Individual Comparisons") - - # if compFile.size ==2 and compFile.shape[0]==2: - # compDiffMags = [] - - # # matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) - # # for q, photFile in enumerate(photFileArray): - # # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - # # compDiffMags = append(compDiffMags,2.5 * log10(photFile[idx][4]/fileCount[q])) - # # instrMags = -2.5 * log10(photFile[idx][4]) - - # matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree) - # args = [(q, matchCoord, photSkyCoord, photFileArray, fileCount) for q in range(len(photFileArray))] - - # with Pool() as pool: - # results = pool.map(process_phot_file_for_variation, args) - - # # Unpack results - # compDiffMags = np.array([result[0] for result in results]) - # instrMags = np.array([result[1] for result in results]) - - # stdCompDiffMags=std(compDiffMags) - # medCompDiffMags=np.nanmedian(compDiffMags) - # medInstrMags= np.nanmedian(instrMags) - - # if np.isnan(stdCompDiffMags) : - # logger.error("Star Variability non rejected") - # stdCompDiffMags=99 - # stdCompStar.append(stdCompDiffMags) - # sortStars.append([compFile[0],compFile[1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) - - - - # else: - # sortStars=[] - # compDiffMags = [] - # instrMags=[] - # matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) - # for q, photFile in enumerate(photFileArray): - # idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q]) - # compDiffMags.append(2.5 * log10(photFile[idx,4]/fileCount[q])) - # instrMags.append(-2.5 * log10(photFile[idx,4])) - - # compDiffMags=array(compDiffMags) - # instrMags=array(instrMags) - - # sortStars=[] - # for z in range(len(compDiffMags[0])): - # stdCompDiffMags=std(compDiffMags[:,z]) - # medCompDiffMags=np.nanmedian(compDiffMags[:,z]) - # medInstrMags=np.nanmedian(instrMags[:,z]) - # if np.isnan(stdCompDiffMags) : - # logger.error("Star Variability non rejected") - # stdCompDiffMags=99 - # stdCompStar.append(stdCompDiffMags) - # sortStars.append([compFile[z,0],compFile[z,1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]) - - # return stdCompStar, sortStars - stdCompStar = [] sortStars = [] logger.debug("Calculating Variation in Individual Comparisons") @@ -469,10 +368,10 @@ def remove_stars_targets(parentPath, compFile, acceptDistance, targetFile, remov logger.info("Removing Target Stars from potential Comparisons") tableFound=False - if not (compFile.shape[0] == 2 and compFile.size ==2): - fileRaDec = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) - else: - fileRaDec = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) + # if not (compFile.shape[0] == 2 and compFile.size ==2): + # fileRaDec = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree) + # else: + # fileRaDec = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree) # Remove any nan rows from targetFile if targetFile is not None: @@ -527,14 +426,7 @@ def remove_stars_targets(parentPath, compFile, acceptDistance, targetFile, remov while vS < len(vServers): v=Vizier(columns=['RAJ2000', 'DEJ2000']) # Skymapper by default does not report the error columns v.ROW_LIMIT=-1 - - #'vizier.nao.ac.jp', - #] - #'vizier.iucaa.in', - #'vizier.china-vo.org', - #'vizier.inasan.ru', - #'vizier.idia.ac.za'] - + v.VIZIER_SERVER=vServers[vS] vS=vS + 1 variableResult=v.query_region(avgCoord, radius=1.5*radius*degree, catalog='VSX') @@ -654,17 +546,13 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): radecname = {'ra' :'raj2000', 'dec': 'dej2000'} if cat_name in ['APASS']: - searchColumns=[] - #searchColumns=[radecname['ra'], radecname['dec'], 'i_mag', opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] - #print (searchColumns) + searchColumns=[] else: searchColumns=[radecname['ra'], radecname['dec'], opt['filter'], opt['error'], opt['colmatch'], opt['colerr']] - #print (searchColumns) - + v=Vizier(columns=searchColumns) # Skymapper by default does not report the error columns v.ROW_LIMIT=-1 - #queryConstraint= # Filter out bad data from catalogues if cat_name == 'PanSTARRS': queryConstraint={'Qual' : '52 || 60 || 61'} @@ -675,31 +563,6 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): elif cat_name == 'APASS': queryConstraint={} - - - # vServers = ['vizier.cfa.harvard.edu','vizier.u-strasbg.fr'] - # vS=0 - # try: - # while vS < len(vServers): - # v=Vizier(columns=['RAJ2000', 'DEJ2000']) # Skymapper by default does not report the error columns - # v.ROW_LIMIT=-1 - - # #'vizier.nao.ac.jp', - # #] - # #'vizier.iucaa.in', - # #'vizier.china-vo.org', - # #'vizier.inasan.ru', - # #'vizier.idia.ac.za'] - - # v.VIZIER_SERVER=vServers[vS] - # vS=vS + 1 - # variableResult=v.query_region(avgCoord, radius=1.5*radius*degree, catalog='VSX') - # if str(variableResult)=="Empty TableList": - # logger.info("VSX Returned an Empty Table from " + str(vServers[vS]) + ".") - # else: - # variableResult=variableResult['B/vsx/vsx'] - # tableFound=True - cycler=0 vS=0 try: @@ -768,13 +631,6 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): if len(resp) != 0: catReject=[] - # v=Vizier() # Skymapper by default does not report the error columns - # v.ROW_LIMIT=-1 - # query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs) - # resp = query[tbname] - - #breakpoint() - for q in range(len(resp)): if np.asarray(resp[opt['filter']][q]) == 0.0 or np.asarray(resp[opt['error']][q]) == 0.0 or np.isnan(resp[opt['filter']][q]) or np.isnan(resp[opt['error']][q]): catReject.append(q) @@ -854,8 +710,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n if not calibPath.exists(): os.makedirs(calibPath) - #Vizier.ROW_LIMIT = -1 - # Get List of Files Used fileList=[] for line in (parentPath / "results/usedImages.txt").read_text().strip().split('\n'): @@ -952,111 +806,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n else: max_sep=1.5 * arcsecond - # The following commented out code is temporary code to manually access Skymapper DR3 and APASS10 prior to them becoming publically available through vizier. - # It will be removed in later versions. - - # SKYMAPPER OVERRIDE - # max_sep=1.5 * arcsecond - # smapdr3=genfromtxt(parentPath / 'SkymapperDR3.csv', dtype=float, delimiter=',') - # coords = namedtuple(typename='data',field_names=['ra','dec','mag','emag','cat_name', 'colmatch', 'colerr']) - # coords.ra=smapdr3[:,0] - # coords.dec=smapdr3[:,1] - # print (opt['filter']) - # if opt['filter'] == 'uPSF': - # coords.mag = smapdr3[:,2] - # coords.emag = smapdr3[:,3] - # coords.colmatch = smapdr3[:,4] - # coords.colerr = smapdr3[:,5] - # coords.colrev = 0 - - # if opt['filter'] == 'gPSF': - # coords.mag = smapdr3[:,4] - # coords.emag = smapdr3[:,5] - # coords.colmatch = smapdr3[:,6] - # coords.colerr = smapdr3[:,7] - # coords.colrev = 0 - - # if opt['filter'] == 'rPSF': - # coords.mag = smapdr3[:,6] - # coords.emag = smapdr3[:,7] - # coords.colmatch = smapdr3[:,8] - # coords.colerr = smapdr3[:,9] - # coords.colrev = 0 - - # if opt['filter'] == 'iPSF': - # coords.mag = smapdr3[:,8] - # coords.emag = smapdr3[:,9] - # coords.colmatch = smapdr3[:,6] - # coords.colerr = smapdr3[:,7] - # coords.colrev = 1 - - # if opt['filter'] == 'zPSF': - # coords.mag = smapdr3[:,10] - # coords.emag = smapdr3[:,11] - # coords.colmatch = smapdr3[:,8] - # coords.colerr = smapdr3[:,9] - # coords.colrev = 1 - - # APASS OVERRIDE - # max_sep=2.5 * arcsecond - # smapdr3=genfromtxt(parentPath / 'ap10.csv', dtype=float, delimiter=',') - # coords = namedtuple(typename='data',field_names=['ra','dec','mag','emag','cat_name', 'colmatch', 'colerr']) - # coords.ra=smapdr3[:,0] - # coords.dec=smapdr3[:,1] - # print (opt['filter']) - - # if opt['filter'] == 'Bmag': - # coords.mag = smapdr3[:,2] - # coords.emag = smapdr3[:,3] - # coords.colmatch = smapdr3[:,4] - # coords.colerr = smapdr3[:,5] - # coords.colrev = 0 - - # if opt['filter'] == 'Vmag': - # coords.mag = smapdr3[:,4] - # coords.emag = smapdr3[:,5] - # coords.colmatch = smapdr3[:,2] - # coords.colerr = smapdr3[:,3] - # coords.colrev = 1 - - # if opt['filter'] == 'umag': - # coords.mag = smapdr3[:,6] - # coords.emag = smapdr3[:,7] - # coords.colmatch = smapdr3[:,8] - # coords.colerr = smapdr3[:,9] - # coords.colrev = 0 - - # if opt['filter'] == 'gmag': - # coords.mag = smapdr3[:,8] - # coords.emag = smapdr3[:,9] - # coords.colmatch = smapdr3[:,10] - # coords.colerr = smapdr3[:,11] - # coords.colrev = 0 - - # if opt['filter'] == 'rmag': - # coords.mag = smapdr3[:,10] - # coords.emag = smapdr3[:,11] - # coords.colmatch = smapdr3[:,12] - # coords.colerr = smapdr3[:,13] - # coords.colrev = 0 - - # if opt['filter'] == 'imag': - # coords.mag = smapdr3[:,12] - # coords.emag = smapdr3[:,13] - # coords.colmatch = smapdr3[:,10] - # coords.colerr = smapdr3[:,11] - # coords.colrev = 1 - - # if opt['filter'] == 'zmag': - # coords.mag = smapdr3[:,14] - # coords.emag = smapdr3[:,15] - # coords.colmatch = smapdr3[:,12] - # coords.colerr = smapdr3[:,13] - # coords.colrev = 1 - # Save catalogue search catalogueOut=np.transpose(np.vstack(np.array([[coords.ra],[coords.dec],[coords.mag],[coords.emag],[coords.colmatch],[coords.colerr]]))) - #catalogueOut=np.array([[coords.ra],[coords.dec],[coords.mag],[coords.emag],[coords.colmatch],[coords.colerr]]) savetxt(parentPath / "results/catalogueSearch.csv", np.asarray(catalogueOut), newline='\n' , delimiter=",", fmt='%0.8f',header='RA,Dec,Mag,err_Mag,colourmatch_mag,colourmatch_mag_err') del catalogueOut @@ -1113,13 +864,11 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n ### If looking for colour, remove those without matching colour information - calibStandsReject=[] if (asarray(calibStands).shape[0] != 9 and asarray(calibStands).size !=9) and len(calibStands) != 0: for q in range(len(asarray(calibStands)[:,0])): reject=0 - #if colourdetect == True: # There needs to be a colour in here to make various plots. if np.isnan(calibStands[q][6]): # if no matching colour reject=1 @@ -1192,8 +941,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n #Colour Term Estimation routine - - # use a temporary calibStands array for colour terms arrayCalibStands=np.asarray(calibStands) @@ -1254,9 +1001,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n x, residuals, rank, s = sqsol plt.plot(outplotx,m*outplotx+c,'r') except: - logger.info("Couldn't fit a line to pre-colour reference plot") - - + logger.info("Couldn't fit a line to pre-colour reference plot") plt.ylim(max(outploty)+0.05,min(outploty)-0.05) plt.xlim(min(outplotx)-0.05,max(outplotx)+0.05) @@ -1271,8 +1016,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n except: print ("Couldn't make colour plot") logger.error(traceback.print_exc()) + # MAKE ALL PRE-COLOUR PLOTS - if colourdetect == True and colourTerm == 0.0: logger.info('Estimating Colour Slope from all Frames...') @@ -1290,14 +1035,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n os.makedirs(colourPath) for qert in range(len(photFileHolder)): - - #breakpoint() - #photFrame = load(parentPath / filen) photFrame = photFileHolder[z] z=z+1 - #photFrame[:,5] = 1.0857 * (photFrame[:,5]/photFrame[:,4]) - #photFrame[:,4]=-2.5 * np.log10(photFrame[:,4]) - photCoords=SkyCoord(ra=photFrame[:,0]*degree, dec=photFrame[:,1]*degree) colTemp=[] for q in range(len(arrayCalibStands[:,0])): @@ -1473,14 +1212,9 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n # Get colour information into photFile # adding in colour columns to photfile - idx,d2d,_=catCoords.match_to_catalog_sky(photCoords) - - - + idx,d2d,_=catCoords.match_to_catalog_sky(photCoords) rejectStars=where(d2d > max_sep)[0] - idx=delete(idx, rejectStars, axis=0) - - + idx=delete(idx, rejectStars, axis=0) for q in range(len(idx)): if colrev == 1: @@ -1490,16 +1224,10 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n photFile[idx[q],9]=pow(pow(coords.colerr[q],2)+pow(coords.emag[q],2),0.5) photFile[idx[q],10]=1 - #breakpoint() - #Replace undetected colours with average colour of the field photFile[:,8]=np.nan_to_num(photFile[:,8],nan=np.nanmedian(photFile[:,8])) photFile[:,9]=np.nan_to_num(photFile[:,9],nan=np.nanmedian(photFile[:,9])) - # If the colour detection routine was run, this was calculated already - #if not colourdetect: - # photFile[:,5]=1.0857 * (photFile[:,5]/photFile[:,4]) - - + notNanTemp=photFile[np.isnan(photFile).any(axis=1)] #rows where any value is nan if colourdetect: notNanTemp[:,4]=notNanTemp[:,4] -(colourTerm*notNanTemp[:,8]) @@ -1508,9 +1236,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n emptyCalibFlag=0 NanTemp=photFile[~np.isnan(photFile).any(axis=1)] #rows where no value is nan - if len(NanTemp) != 0: - #if not colourdetect: - # NanTemp[:,4]=-2.5*log10(photFile[:,4]) + if len(NanTemp) != 0: NanTemp[:,10]=2 # 2 means that there was no colour to use to embed the colour. else: emptyCalibFlag=1 @@ -1527,9 +1253,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n else: calibCoord=SkyCoord(ra=calibStands[:,0]*degree,dec=calibStands[:,1]*degree) - - - if calibStands.size == 13 and calibStands.shape[0]== 13: for q in range(len(calibStands[:,0])): idx,d2d,_=calibCoord.match_to_catalog_sky(photCoords) @@ -1579,10 +1302,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n weights=1/(calibOut[:,1]) linA = np.vstack([outplotx,np.ones(len(outplotx))]).T * np.sqrt(weights[:,np.newaxis]) linB = outploty * np.sqrt(weights) - - #breakpoint() - - try: sqsol = np.linalg.lstsq(linA,linB, rcond=None) @@ -1600,8 +1319,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n plt.subplots_adjust(left=0.15, right=0.98, top=0.98, bottom=0.17, wspace=0.3, hspace=0.4) fig.set_size_inches(6,3) plt.savefig(colourPath / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.png")) - plt.savefig(colourPath / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.eps")) - + plt.savefig(colourPath / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.eps")) except: logger.info("\nFailed linear fit on a postcolour plot: ") @@ -1626,7 +1344,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n # PRINT FINAL REFERENCE IMAGE COLOUR PLOT - #referenceFrame = genfromtxt(parentPath / 'results/referenceFrame.csv', dtype=float, delimiter=',') POSTreferenceFrame=np.copy(referenceFrame) # add the colour error in here POSTreferenceFrame[:,5] = pow(pow(referenceFrame[:,5],2) + pow(colourError,2),0.5) @@ -1759,7 +1476,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n fig = plt.gcf() outplotx=calibOverlord[:,0] outploty=calibOverlord[:,5] - #breakpoint() sqsol = np.linalg.lstsq(np.vstack([calibOverlord[:,0],np.ones(len(calibOverlord[:,0]))]).T,calibOverlord[:,5], rcond=None) m, c = sqsol[0] x, residuals, rank, s = sqsol @@ -1797,7 +1513,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n except: logger.info("Could not create Difference versus Magnitude calibration plot") - # Non-linearity correction routine if linearise == False: logger.info("Skipping Non-linearity correction routine") @@ -1870,14 +1585,10 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n logger.debug("CORRECTING EACH FILE FOR NONLINEARITY") correctFileList = calibPath.glob("*.calibrated.csv") for file in correctFileList: - # NEED TO FIX UP COMPARED! - photFile = genfromtxt(file, dtype=float, delimiter=',') - for q in range(len(photFile[:,0])): photFile[q,4]=photFile[q,4] - nonlinearSlope*photFile[q,4]+nonlinearZero - savetxt(file, photFile, delimiter=",", fmt='%0.8f') # Difference vs time calibration plot diff --git a/astrosource/detrend.py b/astrosource/detrend.py index 0407aef..dc6acda 100644 --- a/astrosource/detrend.py +++ b/astrosource/detrend.py @@ -29,74 +29,35 @@ def detrend_data(paths, filterCode, detrendfraction=0.1): fileList = paths['outcatPath'].glob('*diffExcel*csv') r=0 - #logger.debug(fileList) for file in fileList: - #photFile = load(paths['parent'] / file) - #print (file) photFile = genfromtxt(file, dtype=float, delimiter=',') - #print (photFile) exists=os.path.isfile(str(file).replace('diff','calib')) if exists: - calibFile = genfromtxt(str(file).replace('diff','calib'), dtype=float, delimiter=',') - #logger.debug("Calibration difference") - #logger.debug(-(photFile[:,1]-calibFile[:,1])[0]) + calibFile = genfromtxt(str(file).replace('diff','calib'), dtype=float, delimiter=',') calibDiff=-((photFile[:,1]-calibFile[:,1])[0]) - #logger.debug(photFile[:,1]) - #logger.debug(photFile[:,0]) - - #print (calibFile) - #logger.debug(file) - #logger.debug(photFile[:,1]) - + baseSubDate=min(photFile[:,0]) - #logger.debug(baseSubDate) - #logger.debug(math.floor(baseSubDate)) - - photFile[:,0]=photFile[:,0]-baseSubDate - - - - - #leftMost = click.prompt("Enter left side most valid date:") - #leftFlat = click.prompt("Enter left side end of flat region:") - - #rightFlat = click.prompt("Enter right side start of flat region:") - #rightMost = click.prompt("Enter right side most valid date:") - - + photFile[:,0]=photFile[:,0]-baseSubDate detrendPercent=0.1 leftMost=min(photFile[:,0]) - leftFlat=min(photFile[:,0])+ (detrendPercent * max(photFile[:,0])-min(photFile[:,0])) - + leftFlat=min(photFile[:,0])+ (detrendPercent * max(photFile[:,0])-min(photFile[:,0])) rightFlat=max(photFile[:,0])- (detrendPercent * max(photFile[:,0])-min(photFile[:,0])) rightMost=max(photFile[:,0]) - - - #print (leftMost) - #print (leftFlat) - #print (rightFlat) - #print (rightMost) - # Clip off edges clipReject=[] for i in range(photFile.shape[0]): if photFile[i,0] < float(leftMost) or photFile[i,0] > float(rightMost): clipReject.append(i) - #logger.debug(photFile[i,1]) - #logger.debug("REJECT") logger.debug(photFile.shape[0]) photFile=delete(photFile, clipReject, axis=0) logger.debug(photFile.shape[0]) - # Get an array holding only the flat bits transitReject=[] flatFile=asarray(photFile) for i in range(flatFile.shape[0]): if not (flatFile[i,0] > float(leftMost) and flatFile[i,0] < float(leftFlat)) or (flatFile[i,0] > float(rightFlat) and flatFile[i,0] < float(rightMost)): transitReject.append(i) - #logger.debug(flatFile[i,0]) - #logger.debug("REJECT") logger.debug(flatFile.shape[0]) flatFile=delete(flatFile, transitReject, axis=0) logger.debug(flatFile.shape[0]) @@ -113,7 +74,6 @@ def detrend_data(paths, filterCode, detrendfraction=0.1): for i in range(flatFile.shape[0]): flatFile[i,1]=flatFile[i,1]-(polyFit[1]+(polyFit[0]*flatFile[i,0])) - #Remove trend from actual data if polyFitRequest == 2: for i in range(photFile.shape[0]): @@ -122,7 +82,6 @@ def detrend_data(paths, filterCode, detrendfraction=0.1): for i in range(photFile.shape[0]): photFile[i,1]=photFile[i,1]-(polyFit[1]+(polyFit[0]*photFile[i,0])) - #return basedate to the data photFile[:,0]=photFile[:,0]+baseSubDate diff --git a/astrosource/identify.py b/astrosource/identify.py index 2aca2dd..c34c893 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -26,6 +26,53 @@ logger = logging.getLogger('astrosource') +def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut): + try: + # Determine if the file is from S3 + s3 = False + try: + fitsobj = Path(f) + except TypeError: + fitsobj = f.open() + s3 = True + + # Extract photometry + filepath, photFile = extract_photometry( + fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, + lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut + ) + + if not filepath or not photFile: + return None # Skip if extraction failed + + # Close S3 file if applicable + if s3: + f.close() + filename = f.name + else: + filename = fitsobj.name + + # Return results if photFile is valid + if photFile.size > 100: + photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree) + return filepath, photFile, photSkyCoord + return None + except Exception as e: + print(f"Error processing file {f}: {e}") + return None + +def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut): + # Wrap arguments for multiprocessing + args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut) + with Pool(processes=max([cpu_count()-1,1])) as pool: + results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(filelist))) + + # Filter valid results + results = [res for res in results if res is not None] + phot_dict, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], []) + return list(phot_dict), list(photFileHolder), list(photSkyCoord) + + def rename_data_file(prihdr, bjd=False): prihdrkeys = prihdr.keys() @@ -51,8 +98,6 @@ def rename_data_file(prihdr, bjd=False): airMass=(str(prihdr['AIRMASS']).replace('.','a')) instruMe=(prihdr['INSTRUME']).replace(' ','').replace('/','').replace('-','') - - if (prihdr['MJD-OBS'] == 'UNKNOWN'): timeobs = 'UNKNOWN' elif bjd: @@ -67,86 +112,11 @@ def export_photometry_files(filelist, indir, filetype='csv', bjd=False, ignoreed phot_dict = [] new_files = [] photFileHolder=[] - photSkyCoord=[] - # filelist = list(filelist) - # for f in tqdm(filelist): - # s3 = False - # try: - # fitsobj = Path(f) - # except TypeError: - # fitsobj = f.open() - # s3 = True - - - # filepath, photFile = extract_photometry(fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut) - # if not filepath and not photFile : - # continue - # if s3: - # f.close() - # filename = f.name - # else: - # filename = fitsobj.name - - # if photFile.size > 100: - # phot_dict.append(filepath) - - # photFileHolder.append(photFile) - # photSkyCoord.append(SkyCoord(ra=photFile[:,0]*u.degree, dec=photFile[:,1]*u.degree)) - - def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut): - try: - # Determine if the file is from S3 - s3 = False - try: - fitsobj = Path(f) - except TypeError: - fitsobj = f.open() - s3 = True - - # Extract photometry - filepath, photFile = extract_photometry( - fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, - lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut - ) - - if not filepath or not photFile: - return None # Skip if extraction failed - - # Close S3 file if applicable - if s3: - f.close() - filename = f.name - else: - filename = fitsobj.name - - # Return results if photFile is valid - if photFile.size > 100: - photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree) - return filepath, photFile, photSkyCoord - return None - except Exception as e: - print(f"Error processing file {f}: {e}") - return None - - - def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut): - # Wrap arguments for multiprocessing - args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut) - with Pool(processes=max([cpu_count()-1,1])) as pool: - results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(filelist))) - - # Filter valid results - results = [res for res in results if res is not None] - phot_dict, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], []) - return list(phot_dict), list(photFileHolder), list(photSkyCoord) - - + photSkyCoord=[] # Example usage phot_dict, photFileHolder, photSkyCoord = process_files_multiprocessing( filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut - ) - - + ) return phot_dict, photFileHolder, photSkyCoord @@ -212,10 +182,8 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr photFile[:,0] = photFile[:,0] - 180.0 for entry in range(len(photFile[:,0])): if photFile[entry,0] < 0: - photFile[entry,0] = photFile[entry,0] + 360 + photFile[entry,0] = photFile[entry,0] + 360 - - #remove odd zero entries photFile[:,0][photFile[:,0] == 0.0 ] = nan photFile[:,0][photFile[:,0] == 0.0 ] = nan @@ -227,9 +195,7 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr if len(photFile) == num_rows_with_nans: print(f"REJECT {infile}: Very likely no wcs fit - RA and Dec columns are zero.") - photFile=photFile[~isnan(photFile).any(axis=1)] - - + photFile=photFile[~isnan(photFile).any(axis=1)] #remove lowcounts rejectStars=where(photFile[:,4] < lowestcounts)[0] @@ -281,10 +247,7 @@ def process_convert_photometry_file(fn, racut, deccut, radiuscut, ignoreedgefrac (photFile[:, 1] < max(photFile[:, 1]) - decClip) ] - # Remove zero and low-count entries - - #num_rows_with_nans = sum(isnan(photFile).any(axis=1)) - #print (len((photFile[:, 0] != 0) & (photFile[:, 1] != 0))) + # Remove zero and low-count entries if len((photFile[:, 0] != 0) & (photFile[:, 1] != 0)) == 0: print(f"REJECTED {fn}: Very likely no wcs fit - RA and Dec columns are zero.") @@ -328,108 +291,10 @@ def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=180 new_files = [] photFileHolder=[] photSkyCoord=[] - #breakpoint() - # for fn in filelist: - # photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1) - # logger.info(fn) - # # reject nan entries in file - # if photFile.size > 16 and photFile.shape[1] == 8: #ignore zero sized files and files with only one or two entries - # if max(photFile[:,0]) < 360 and max(photFile[:,1]) < 90: - - # if (photFile.size > 50): - # if not (( asarray(photFile[:,0]) > 360).sum() > 0) : - # if not(( asarray(photFile[:,1]) > 90).sum() > 0) : - - # photFile=photFile[~isnan(photFile).any(axis=1)] - - - - # # if radial cut do that otherwise chop off image edges - # if racut != -99.9 and deccut !=-99.9 and radiuscut !=-99.9: - # distanceArray=pow((pow(photFile[:,0]-float(racut),2)+pow(photFile[:,1]-float(deccut),2)),0.5) - # photFile=delete(photFile, where(distanceArray > float(radiuscut)/60), axis=0) - # else: - # #Routine to deal with RA 0 crossovers - # crossover=0 - # if (max(photFile[:,0])-min(photFile[:,0])) > 180: - # for entry in range(len(photFile[:,0])): - # if photFile[entry,0] > 180: - # photFile[entry,0] = photFile[entry,0] - 180 - # else: - # photFile[entry,0] = photFile[entry,0] + 180 - - # crossover=1 - - # # Remove edge detections - # raRange=(max(photFile[:,0])-min(photFile[:,0])) - # decRange=(max(photFile[:,1])-min(photFile[:,1])) - # raMin=min(photFile[:,0]) - # raMax=max(photFile[:,0]) - # decMin=min(photFile[:,1]) - # decMax=max(photFile[:,1]) - # raClip=raRange*ignoreedgefraction - # decClip=decRange*ignoreedgefraction - # photFile[:,0][photFile[:,0] < raMin + raClip ] = nan - # photFile[:,0][photFile[:,0] > raMax - raClip ] = nan - # photFile[:,1][photFile[:,1] > decMax - decClip ] = nan - # photFile[:,1][photFile[:,1] < decMin + decClip ] = nan - # if crossover == 1: - # photFile[:,0] = photFile[:,0] - 180.0 - # for entry in range(len(photFile[:,0])): - # if photFile[entry,0] < 0: - # photFile[entry,0] = photFile[entry,0] + 360 - # #remove odd zero entries - # photFile[:,0][photFile[:,0] == 0.0 ] = nan - # photFile[:,0][photFile[:,0] == 0.0 ] = nan - # photFile[:,1][photFile[:,1] == 0.0 ] = nan - # photFile[:,1][photFile[:,1] == 0.0 ] = nan - # photFile=photFile[~isnan(photFile).any(axis=1)] - - # #remove lowcounts - # rejectStars=where(photFile[:,4] < lowestcounts)[0] - # photFile=delete(photFile, rejectStars, axis=0) - - # filepath = Path(fn).with_suffix('.npy') - - # photFile=c_[photFile,zeros(len(photFile[:,0])),zeros(len(photFile[:,0])),zeros(len(photFile[:,0])),zeros(len(photFile[:,0]))] - # photFile[:,8]=nan - # photFile[:,9]=nan - # photFile[:,10]=0 - # photFile[:,11]=copy(photFile[:,4]) - - # if photFile.size > 16: - # new_files.append(filepath.name) - # photFileHolder.append(photFile) - # photSkyCoord.append(SkyCoord(ra=photFile[:,0]*u.degree, dec=photFile[:,1]*u.degree)) - - - - - # else: - # logger.debug("REJECT") - # logger.debug(fn) - # else: - # logger.debug("REJECT") - # logger.debug(fn) - # else: - # logger.debug("REJECT") - # logger.debug(fn) - # else: - # logger.debug("REJECT") - # logger.debug(fn) - # else: - # logger.debug("REJECT") - # logger.debug(fn) - - - - - # Example usage new_files, photFileHolder, photSkyCoord = process_convert_files_multiprocessing( filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger - ) - + ) if new_files ==[] or photFileHolder==[] : raise AstrosourceException("Either there are no files of this photometry type or radial Cut seems to be outside of the range of your images... check your racut, deccut and radiuscut") @@ -487,8 +352,6 @@ def gather_files(paths, filelist=None, filetype="fz", bjd=False, ignoreedgefract file1=open(paths['parent'] / "filterCode","wb") pickle.dump(filterCode, file1) file1.close - - #breakpoint() return phot_list, filterCode, photFileHolder, photSkyCoord @@ -543,8 +406,7 @@ def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistanc def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger): - - + # Hack to get windows to not multiprocess until I figure out how to do it. if platform.system() == "Windows": results=[] @@ -571,8 +433,6 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen updated_referenceFrame = refFrameUpdate photReject.extend(rejects) - #breakpoint() - return updated_referenceFrame, photReject @@ -633,61 +493,15 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym if os.path.exists(paths['parent'] / "photFileHolder"): os.remove(paths['parent'] / "photFileHolder") - # fileSizer=0 - # logger.info("Finding image with most stars detected and reject ones with bad WCS") - # referenceFrame = None - - # counter=0 - # for file in fileList: - # photFile = photFileHolder[counter] - - # # Sort through and find the largest file and use that as the reference file - # if photFile.size > fileSizer: - # referenceFrame = photFile - # fileSizer = photFile.size - # fileRaDec = photCoords[counter] - # logger.debug("{} - {}".format(photFile.size, file)) - # counter=counter+1 - + logger.info("Finding image with most stars detected and rejecting ones with bad WCS") referenceFrame = None - file_sizer = 0 + #file_sizer = 0 fileRaDec = None - - - - # list_of_sizes=[] - - - # with Pool() as pool: - # # Use pool.map to distribute the computation across multiple processes - # list_of_sizes = pool.starmap(compute_size, [(i, photFileHolder) for i in range(len(photFileHolder))]) - - - # def compute_size(phot_file): - # """Compute the size for a specific file.""" - # return phot_file.size - # import time - - - # with Pool() as pool: - # # Parallelize the size computation - # list_of_sizes = pool.map(compute_size, photFileHolder) - - - import time - # googtime=time.time() - - - list_of_sizes=[] for filething in photFileHolder: - # zingtime=time.time() list_of_sizes.append(filething.size) - # print ("Z: " + str(time.time()-zingtime)) - - # print ("D: " + str(time.time()-googtime)) largest_file_index= argmax(list_of_sizes) @@ -695,28 +509,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym fileRaDec = photCoords[largest_file_index] fileSizer=referenceFrame.size - # breakpoint() - - # for file in len(photFileHolder): - # list_of_sizes.append(photFileHolder.size) - - # with ProcessPoolExecutor() as executor: - # results = list(executor.map( - # evaluate_file_size, - # range(len(fileList)), - # [fileList] * len(fileList), - # [photFileHolder] * len(fileList), - # [photCoords] * len(fileList) - # )) - - # # Process the results - # for file_size, index, file_ra_dec in results: - # if file_size > file_sizer: - # referenceFrame = photFileHolder[index] - # file_sizer = file_size - # fileRaDec = file_ra_dec - # logger.debug("{} - {}".format(file_size, fileList[index])) - if not referenceFrame.size: raise AstrosourceException("No suitable reference files found") @@ -749,7 +541,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym lowcounts=lowcounts+(0.05*lowcounts) hicounts=hicounts-(0.05*hicounts) - logger.debug("Number of stars post") referenceFrame = delete(referenceFrame, rejectStars, axis=0) logger.debug("Final count range, Low: " +str(int(lowcounts))+ " High: "+ str(int(hicounts))) @@ -778,102 +569,18 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym while (compchecker < mincompstars): # Keep going until you get the minimum number of Comp Stars imgsize=imageFracReject * fileSizer # set threshold size rejStartCounter = 0 - #usedImages=[] # Set up used images array imgReject = 0 # Number of images rejected due to high rejection rate loFileReject = 0 # Number of images rejected due to too few stars in the photometry file wcsFileReject=0 imgOverride=0 - # q=0 - # photReject=[] - - # for Nholder in range(len(photFileHolder)): - - # if ( not referenceFrame.shape[0] < mincompstars): - # rejStartCounter = rejStartCounter +1 - # #photFile = load(paths['parent'] / file) - # photFile = photFileHolder[q] - # logger.debug('Image Number: ' + str(rejStartCounter)) - # logger.debug(file) - # logger.debug("Image threshold size: "+str(imgsize)) - # logger.debug("Image catalogue size: "+str(photFile.size)) - # imgRejFlag=0 - # if photFile.size > imgsize and photFile.size > 7 : - # phottmparr = asarray(photFile) - # photRAandDec = photCoords[q] - - # # Checking existance of stars in all photometry files - # rejectStars=[] # A list to hold what stars are to be rejected - - # if (( phottmparr[:,0] > 360).sum() == 0) and ( phottmparr[0][0] != 'null') and ( phottmparr[0][0] != 0.0) : - # # Faster way than loop - # testStars=SkyCoord(ra = referenceFrame[:,0]*u.degree, dec = referenceFrame[:,1]*u.degree) - # idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec) - # rejectStars=where(d2d.arcsecond > acceptDistance)[0] - - # else: - # logger.debug('**********************') - # logger.debug('Image Rejected due to problematic entries in Photometry File') - # logger.debug('**********************') - # imgReject=imgReject+1 - # photReject.append(q) - # imgRejFlag=1 - - # # if the rejectstar list is not empty, remove the stars from the reference List - # if len(rejectStars) != 0: - - # if not (((len(rejectStars) / referenceFrame.shape[0]) > starreject) and rejStartCounter > rejectStart): - # referenceFrame = delete(referenceFrame, rejectStars, axis=0) - # logger.debug('**********************') - # logger.debug('Stars Removed : ' +str(len(rejectStars))) - # logger.debug('Remaining Stars: ' +str(referenceFrame.shape[0])) - # logger.debug('**********************') - - # else: - # logger.debug('**********************') - # logger.debug('Image Rejected due to too high a fraction of rejected stars') - # logger.debug(len(rejectStars) / referenceFrame.shape[0]) - # logger.debug('**********************') - # imgReject=imgReject+1 - # photReject.append(q) - - # elif imgRejFlag==0: - # logger.debug('**********************') - # logger.debug('All Stars Present') - # logger.debug('**********************') - - # # If we have removed all stars, we have failed! - # if (referenceFrame.shape[0]==0): - # logger.error("Problem file - {}".format(file)) - # logger.error("Running Loop again") - - # elif photFile.size < 7: - # logger.error('**********************') - # logger.error("WCS Coordinates broken") - # logger.error('**********************') - # wcsFileReject=wcsFileReject+1 - # photReject.append(q) - # else: - # logger.debug('**********************') - # logger.debug("CONTAINS TOO FEW STARS") - # logger.debug('**********************') - # loFileReject=loFileReject+1 - # photReject.append(q) - - # q=q+1 - - - # Example usage updated_referenceFrame, photReject = process_photometry_files_multiprocessing( photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger ) - #breakpoint() # Remove files and Hold the photSkyCoords in memory - photCoords=asarray(photCoords, dtype=object) - photCoords=delete(photCoords, photReject, axis=0) photFileHolder=asarray(photFileHolder, dtype=object) photFileHolder=delete(photFileHolder, photReject, axis=0) @@ -962,7 +669,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym photCoords=originalphotSkyCoord photFileHolder=originalphotFileHolder - elif starreject == 0.15 and imageFracReject == 0.8 and mincompstars !=1: logger.error("Maximum number of Candidate Comparison Stars found this cycle: " + str(compchecker)) logger.error("Failed to find sufficient comparison candidates with the maximum restrictions, trying with a lower value for mincompstars") @@ -1137,7 +843,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd) # If no results try next catalogue - #print (len(coords.ra)) if len(coords.ra) == 0: coords=[] raise AstrosourceException("Empty catalogue produced from catalogue call") diff --git a/astrosource/main.py b/astrosource/main.py index e1c3913..096e929 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -126,12 +126,8 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree logger.info('All output files removed') return if not (ra and dec) and not target_file and not variablehunt: - #logger.error("Either RA and Dec or a targetfile must be specified") logger.error("No specified RA or Dec nor targetfile nor request for a variable hunt. It is assumed you have no target to analyse.") notarget=True - #return - - if ra and dec: ra, dec = convert_coords(ra, dec) diff --git a/astrosource/periodic.py b/astrosource/periodic.py index 713d212..410a496 100644 --- a/astrosource/periodic.py +++ b/astrosource/periodic.py @@ -38,55 +38,10 @@ # SOFTWARE. def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, paths, filterCode, logger): - # try: - # logger.debug(file) - # variableName = file.stem.split('_')[0] - # logger.debug("Variable Name: {}".format(variableName)) - - # # Load data - # varData = genfromtxt(file, dtype=float, delimiter=',') - # calibFile = file.parent / "{}{}".format(file.stem.replace('diff', 'calib'), file.suffix) - # logger.debug(calibFile) - - # calibData = None - # if calibFile.exists(): - # calibData = genfromtxt(calibFile, dtype=float, delimiter=',') - - # # Determine which dataset to use for analysis - # if calibFile.exists() and calibData is not None and calibData.size > 3: - # dataToUse = calibData - # logger.info("Using calibrated data") - # else: - # dataToUse = varData - # logger.info("Using differential data") - - # # Phase dispersion minimization - # pdm = phase_dispersion_minimization(dataToUse, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) - - # # Save plots and results - # plt.figure(figsize=(15, 5)) - # plt.plot(pdm["periodguess_array"], pdm["distance_results"]) - # plt.gca().invert_yaxis() - # plt.title("Range {0} d Steps: {1}".format(maxperiod - minperiod, periodsteps)) - # plt.xlabel(r"Trial Period") - # plt.ylabel(r"Likelihood of Period") - # plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") - # plt.clf() - - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write(f"Variable: {variableName}\n") - # f.write(f"Distance Method Estimate (days): {pdm['distance_minperiod']}\n") - # f.write(f"Distance method error: {pdm['distance_error']}\n") - - # except Exception as e: - # logger.error(f"Error processing file {file}: {e}") try: - trialRange=[minperiod, maxperiod] - - # logger.debug(file) + trialRange=[minperiod, maxperiod] variableName=file.stem.split('_')[0] - #logger.debug(str(outcatPath).replace('//','')) logger.debug("Variable Name: {}".format(variableName)) varData = genfromtxt(file, dtype=float, delimiter=',') calibFile = file.parent / "{}{}".format(file.stem.replace('diff','calib'), file.suffix) @@ -94,7 +49,6 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p if calibFile.exists(): calibData=genfromtxt(calibFile, dtype=float, delimiter=',') - #print (calibData.size) if calibFile.exists(): if (calibData.size > 3): pdm=phase_dispersion_minimization(calibData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) @@ -248,8 +202,7 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p tempPeriodCatOut=asarray(tempPeriodCatOut) savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - # Plot publication plots - + # Plot publication plots plt.figure(figsize=(5, 3)) plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5) @@ -274,9 +227,7 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p plt.clf() - - # ANOVA - + # ANOVA if calibFile.exists(): if (calibData.size > 3): if len(calibData[:,0]) < 75: @@ -306,8 +257,7 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p minperbin = 3 - # Theta Anova Method off for the moment until I put in a command-line option - + # Theta Anova Method off for the moment until I put in a command-line option # if calibFile.exists(): # if (calibData.size > 3): # aovoutput=aov_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) @@ -315,7 +265,6 @@ def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, p # if (varData.size > 3): # aovoutput=aov_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - # logger.debug("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])) #with open(paths['parent'] / "periodEstimates.txt", "a+") as f: # f.write("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])+"\n") @@ -587,7 +536,6 @@ def sigclip_magseries(times, mags, errs, # fake the errors if they don't exist # this is inconsequential to sigma-clipping # we don't return these dummy values if the input errs are None - #print (errs) if errs is None: # assume 0.1% errors if not given # this should work for mags and fluxes @@ -1068,16 +1016,6 @@ def aov_periodfind(times, # default end period is length of time series startf = 1.0/(stimes.max() - stimes.min()) - # # if we're not using autofreq, then use the provided frequencies - # if not autofreq: - # frequencies = nparange(startf, endf, stepsize) - - # else: - # # this gets an automatic grid of frequencies to use - # frequencies = get_frequency_grid(stimes, - # minfreq=startf, - # maxfreq=endf) - stepsize = (endf-startf)/periodsteps frequencies = nparange(startf, endf, stepsize) @@ -1099,8 +1037,6 @@ def aov_periodfind(times, periods = 1.0/frequencies plt.plot(periods, lsp) - #plt.gca().invert_yaxis() - #plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) plt.xlabel(r"Trial Period") plt.ylabel(r"Likelihood of Period") plt.savefig(periodPath / f"{variableName}_ANOVAThetaLikelihoodPlot.png") @@ -1167,10 +1103,6 @@ def aov_periodfind(times, perioddiff = abs(period - prevperiod) bestperiodsdiff = [abs(period - x) for x in nbestperiods] - # print('prevperiod = %s, thisperiod = %s, ' - # 'perioddiff = %s, peakcount = %s' % - # (prevperiod, period, perioddiff, peakcount)) - # this ensures that this period is different from the last # period and from all the other existing best periods by # periodepsilon to make sure we jump to an entire different peak @@ -1583,29 +1515,9 @@ def aovhm_periodfind(times, # default end period is length of time series startf = 1.0/(stimes.max() - stimes.min()) - - - - # # if we're not using autofreq, then use the provided frequencies - # if not autofreq: - # frequencies = nparange(startf, endf, stepsize) - - # else: - # # this gets an automatic grid of frequencies to use - # frequencies = get_frequency_grid(stimes, - # minfreq=startf, - # maxfreq=endf) - stepsize = (endf-startf)/periodsteps frequencies = nparange(startf, endf, stepsize) - # map to parallel workers - - - # if (not nworkers) or (nworkers > NCPUS): - # nworkers = NCPUS - #nworkers=1 - # renormalize the working mags to zero and scale them so that the # variance = 1 for use with our LSP functions if normalize: @@ -1619,9 +1531,6 @@ def aovhm_periodfind(times, magvariance_bot = (nmags.size - 1)*npsum(1.0/(serrs*serrs)) / nmags.size magvariance = magvariance_top/magvariance_bot - #tasks = [(stimes, nmags, serrs, x, nharmonics, magvariance) - # for x in frequencies] - lsp=[] for x in frequencies: lsp.append(aovhm_theta(times, mags, errs, x, nharmonics, magvariance)) @@ -1630,8 +1539,6 @@ def aovhm_periodfind(times, periods = 1.0/frequencies plt.plot(periods, lsp) - #plt.gca().invert_yaxis() - #plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) plt.xlabel(r"Trial Period") plt.ylabel(r"Likelihood of Period") plt.savefig(periodPath / f"{variableName}_ANOVAharmonic_LikelihoodPlot.png") @@ -1698,10 +1605,6 @@ def aovhm_periodfind(times, perioddiff = abs(period - prevperiod) bestperiodsdiff = [abs(period - x) for x in nbestperiods] - # print('prevperiod = %s, thisperiod = %s, ' - # 'perioddiff = %s, peakcount = %s' % - # (prevperiod, period, perioddiff, peakcount)) - # this ensures that this period is different from the last # period and from all the other existing best periods by # periodepsilon to make sure we jump to an entire different peak @@ -1722,8 +1625,6 @@ def aovhm_periodfind(times, plt.gca().invert_yaxis() plt.title("Period: " + str(finperiods[bestperiodind])) plt.xlim(-0.01,2.01) - #plt.gca().invert_yaxis() - #plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) plt.xlabel(r"Phase") plt.ylabel(r"Magnitude") plt.savefig(periodPath / f"{variableName}_ANOVAHarmonicLightcurve.png") @@ -1940,8 +1841,6 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu break stepper=stepper+1 - - #print ("Stdev method error: " + str((righthandP - lefthandP)/2)) pdm["stdev_error"] = (righthandP - lefthandP)/2 @@ -1983,10 +1882,7 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu def LombScargleMultiterm(infile, t, m, d, periodlower=0.2, periodupper=2.5, nterms=1, multisearch=False, samples=5, disablelightcurve=False, periodPath=False, variableName="NoName", periodsteps=10000): - #print( - # 'using ' + str(samples) + ' samples per peak, start P = ' + str(periodlower) + ', end P = ' + str(periodupper)) - # Calculate the Lomb-Scargle periodogram values - + ls = LombScargle(t, m, d, nterms=nterms, fit_mean=True) # create equally spaced test in time space, then swap to unequal spaces in frequency @@ -1994,10 +1890,7 @@ def LombScargleMultiterm(infile, t, m, d, periodlower=0.2, periodupper=2.5, nter freq = 1/freq - try: - - #freq, power = ls.autopower(samples_per_peak=samples, minimum_frequency=1 / periodupper, - # maximum_frequency=1 / periodlower) + try: power = ls.power(freq) except: print ("Lomb Scargle failed.") # Need to hunt down a very rare but existant memory problem @@ -2087,7 +1980,7 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1 logger.info("Maximum Period Tested : " +str(maxperiod)) logger.info("Number of Period Trials: " +str(periodsteps)) - trialRange=[minperiod, maxperiod] + #trialRange=[minperiod, maxperiod] # Get list of phot files periodPath = paths['periods'] @@ -2113,279 +2006,10 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1 logger=logger ) - #breakpoint() - # Use multiprocessing with Pool(processes=max([cpu_count()-1,1])) as pool: pool.map(worker, fileList) - - # Load in the files - # for file in fileList: - # logger.debug(file) - # variableName=file.stem.split('_')[0] - # #logger.debug(str(outcatPath).replace('//','')) - # logger.debug("Variable Name: {}".format(variableName)) - # varData = genfromtxt(file, dtype=float, delimiter=',') - # calibFile = file.parent / "{}{}".format(file.stem.replace('diff','calib'), file.suffix) - # logger.debug(calibFile) - # if calibFile.exists(): - # calibData=genfromtxt(calibFile, dtype=float, delimiter=',') - - # #print (calibData.size) - # if calibFile.exists(): - # if (calibData.size > 3): - # pdm=phase_dispersion_minimization(calibData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) - # else: - # logger.info("Calibration File not large enough to run period methods on Calibrated data, running on differential data") - # pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) - # else: - # pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName) - - # plt.figure(figsize=(15, 5)) - - # logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])) - # logger.debug("Distance method error: " + str(pdm["distance_error"])) - - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write("Variable : "+str(variableName) +"\n") - # f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n") - # f.write("Distance method error: " + str(pdm["distance_error"])+"\n") - - # plt.plot(pdm["periodguess_array"], pdm["distance_results"]) - # plt.gca().invert_yaxis() - # plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) - # plt.xlabel(r"Trial Period") - # plt.ylabel(r"Likelihood of Period") - # plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png") - # plt.clf() - - # if (varData.size > 3): - # phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1 - - # plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') - # plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') - # plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') - # plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') - # plt.gca().invert_yaxis() - # plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) - # plt.xlabel(r"Phase ($\phi$)") - # plt.ylabel(f"Differential {filterCode} Magnitude") - # plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png") - # plt.clf() - - # if calibFile.exists(): - # if (calibData.size > 3): - # phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1 - # plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') - # plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') - # plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') - # plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') - # plt.gca().invert_yaxis() - # plt.title("Period: {0} d Steps: {1}".format(pdm["distance_minperiod"], periodsteps)) - # plt.xlabel(r"Phase ($\phi$)") - # plt.ylabel(f"Calibrated {filterCode} Magnitude") - # plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png") - # plt.clf() - - # tempPeriodCatOut=[] - # for g in range(len(calibData[:,0])): - # tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # if (varData.size > 3): - - # tempPeriodCatOut=[] - # for g in range(len(phaseTest)): - # tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # tempPeriodCatOut=[] - # for g in range(len(varData[:,0])): - # tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4): - # logger.info("No PDM results due to lack of datapoint coverage") - # else: - # logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])) - # phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1 - # logger.debug("PDM method error: " + str(pdm["stdev_error"])) - - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n") - # f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n") - - - # plt.plot(pdm["periodguess_array"], pdm["stdev_results"]) - # plt.gca().invert_yaxis() - # plt.title("Range {0} d Steps: {1}".format(trialRange, periodsteps)) - # plt.xlabel(r"Trial Period") - # plt.ylabel(r"Likelihood of Period") - # plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png") - - # plt.clf() - - # if (varData.size > 3): - # plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None') - # plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None') - # plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None') - # plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None') - # plt.gca().invert_yaxis() - # plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) - # plt.xlabel(r"Phase ($\phi$)") - # plt.ylabel(r"Differential " + str(filterCode) + " Magnitude") - # plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png") - # plt.clf() - - # if calibFile.exists(): - # if (calibData.size > 3): - # phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1 - # plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None') - # plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None') - # plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None') - # plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None') - # plt.gca().invert_yaxis() - # plt.title("Period: {0} d Steps: {1}".format(pdm["stdev_minperiod"], periodsteps)) - # plt.xlabel(r"Phase ($\phi$)") - # plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude") - # plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png") - # plt.clf() - - # tempPeriodCatOut=[] - # for g in range(len(calibData[:,0])): - # tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # if (varData.size > 3): - - # tempPeriodCatOut=[] - # for g in range(len(phaseTest)): - # tempPeriodCatOut.append([phaseTest[g],varData[g,1]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # tempPeriodCatOut=[] - # for g in range(len(varData[:,0])): - # tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]]) - # tempPeriodCatOut=asarray(tempPeriodCatOut) - # savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f') - - # # Plot publication plots - - # plt.figure(figsize=(5, 3)) - - # plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5) - # plt.gca().invert_yaxis() - # plt.xlabel(r"Trial Period") - # plt.ylabel(r"Likelihood of Period") - # plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) - # plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300) - # plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps") - - # plt.clf() - - # plt.figure(figsize=(5, 3)) - - # plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5) - # plt.gca().invert_yaxis() - # plt.xlabel(r"Trial Period") - # plt.ylabel(r"Likelihood of Period") - # plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4) - # plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300) - # plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps") - - # plt.clf() - - - # # ANOVA - - # if calibFile.exists(): - # if (calibData.size > 3): - # if len(calibData[:,0]) < 75: - # binsize=0.1 - # else: - # binsize=0.05 - # minperbin=int((len(calibData[:,0])/10)) - # elif (varData.size > 3): - # if len(varData[:,0]) < 75: - # binsize=0.1 - # else: - # binsize=0.05 - # minperbin=int((len(varData[:,0])/10)) - - # else: - # if (varData.size > 3): - # if len(varData[:,0]) < 75: - # binsize=0.1 - # else: - # binsize=0.05 - # minperbin=int((len(varData[:,0])/10)) - - # if 'minperbin' in locals(): - # if minperbin > 10: - # minperbin=10 - # else: - # minperbin = 3 - - - # # Theta Anova Method off for the moment until I put in a command-line option - # # if calibFile.exists(): - # # if (calibData.size > 3): - # # aovoutput=aov_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - # # else: - # # if (varData.size > 3): - # # aovoutput=aov_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - - - # # logger.debug("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])) - # #with open(paths['parent'] / "periodEstimates.txt", "a+") as f: - # # f.write("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])+"\n") - - # if calibFile.exists(): - # if (calibData.size > 3): - # aovhmoutput=aovhm_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - # elif (varData.size > 3): - # aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - - # else: - # if (varData.size > 3): - # aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - - # logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])) - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n") - - # # LOMB SCARGLE - # for nts in range(2): - # if calibFile.exists(): - # if (calibData.size > 3): - # lscargoutput = LombScargleMultiterm('periodifile', (calibData[:, 0]), (calibData[:, 1]), (calibData[:, 2]), - # nterms=nts+1, - # periodlower=minperiod, periodupper=maxperiod, samples=20, - # disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - - # logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)) - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n") - # else: - # if (varData.size > 3): - # lscargoutput = LombScargleMultiterm('periodifile', (varData[:, 0]), (varData[:, 1]), (varData[:, 2]), - # nterms=nts+1, - # periodlower=minperiod, periodupper=maxperiod, samples=20, - # disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps) - - # logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)) - # with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f: - # f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n") - - # if 'pdm' in locals(): - # return pdm["distance_minperiod"] - # else: return 0.0 # Needed for windows to multiprocess appropriately diff --git a/astrosource/plots.py b/astrosource/plots.py index 444f1b5..f55011f 100644 --- a/astrosource/plots.py +++ b/astrosource/plots.py @@ -11,7 +11,6 @@ import numpy as np import traceback -#from astrosource.utils import photometry_files_to_array, AstrosourceException from astrosource.utils import AstrosourceException logger = logging.getLogger('astrosource') @@ -109,7 +108,6 @@ def plot_variability(output, variableID, parentPath, compFile): calibCompExist=True calibCompStarPlot = [] - #for q in range(len(calibCompSkyCoord)): if calibnumber==1: idx, d2d, _ = calibCompSkyCoord.match_to_catalog_sky(outputSkyCoord) calibCompStarPlot.append([output[idx][2],output[idx][3]]) @@ -140,7 +138,6 @@ def plot_variability(output, variableID, parentPath, compFile): plt.plot(np.asarray(compStarPlot)[:,0],np.asarray(compStarPlot)[:,1], 'yo') if calibCompExist: plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none') - # plt.plot(linex, liney) plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04) plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1) @@ -195,7 +192,6 @@ def plot_variability(output, variableID, parentPath, compFile): plt.plot(np.asarray(compStarPlot)[:,0],np.asarray(compStarPlot)[:,1], 'yo') if calibCompExist: plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none') - # plt.plot(linex, liney) plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04) plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1) plt.grid(True) @@ -215,7 +211,6 @@ def plot_variability(output, variableID, parentPath, compFile): if calibCompExist: plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none') fig.set_size_inches(16,9) - # plt.plot(linex, liney) plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04) plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1) plt.grid(True) From da84fb5473196fdf53048b59abfe6a1290285dc1 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 1 Jun 2025 08:16:24 +1000 Subject: [PATCH 18/28] better timeout handling and skymapper handling --- astrosource/comparison.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index b6a094c..d8c2d15 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -26,7 +26,7 @@ import requests import http import urllib3 - +from collections import OrderedDict from astrosource.utils import AstrosourceException import logging @@ -566,17 +566,25 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): cycler=0 vS=0 try: + vtimeout=180 while vS < len(vServers): v.VIZIER_SERVER=vServers[vS] + v.TIMEOUT = vtimeout try: query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs) - except: - logger.info("Failed vizier query)") - if str(query)=="Empty TableList": + if str(query)=="Empty TableList": + vS=vS+1 + else: + break + except requests.exceptions.ReadTimeout: + + logger.info("Read Timeout for Vizier, might be a very dense field or slow server today") vS=vS+1 - else: - break + except: + logger.info("Failed vizier query") + print(traceback.print_exc()) + except VOSError: raise AstrosourceException("Could not find RA {} Dec {} in {}".format(avgCoord.ra.value,avgCoord.dec.value, cat_name)) @@ -771,7 +779,15 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n raise AstrosourceException(f"{filterCode} is not accepted at present") # Look up in online catalogues and make sure there are sufficient comparison stars - + # If in the realm of skymapper and a skymapper filter, put skymapper at the front + if avgCoord.dec.deg < 15 and 'SkyMapper' in catalogues: + ordered = OrderedDict() + ordered['SkyMapper'] = catalogues['SkyMapper'] + for key, val in catalogues.items(): + if key != 'SkyMapper': + ordered[key] = val + catalogues = ordered + coords=[] for cat_name, opt in catalogues.items(): try: From d2886331b11e5a4fe95082e6f44e97235def9551 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 08:03:18 +1000 Subject: [PATCH 19/28] small bugs in the fz route --- astrosource/identify.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/astrosource/identify.py b/astrosource/identify.py index c34c893..dcd0571 100644 --- a/astrosource/identify.py +++ b/astrosource/identify.py @@ -42,7 +42,7 @@ def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, de lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut ) - if not filepath or not photFile: + if filepath is None or photFile is None or len(photFile) == 0: return None # Skip if extraction failed # Close S3 file if applicable @@ -65,7 +65,7 @@ def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowe # Wrap arguments for multiprocessing args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut) with Pool(processes=max([cpu_count()-1,1])) as pool: - results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(filelist))) + results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(list(filelist)))) # Filter valid results results = [res for res in results if res is not None] From b3eee081b231f2ce2c19bdadb2f49779fddaef61 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 11:53:02 +1000 Subject: [PATCH 20/28] booleans aren't floats, Michael --- astrosource/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/main.py b/astrosource/main.py index 096e929..d14f10b 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -55,7 +55,7 @@ -@click.option('--period', is_flag=True, type=float, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') +@click.option('--period', is_flag=True, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') @click.option('--periodlower', '-pl', type=float, default=-99.9, help='Shortest period to trial in days. Default is 0.05 days') @click.option('--periodupper', '-pu', type=float, default=-99.9, help='Longest period to trial in days. Default is one-third the observational baseline of your dataset.') @click.option('--periodtests', '-pt', type=int, default=10000, help='Number of different trial periods to run') From da6b1ec6e756d4559026fbdc08b113103645731c Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 11:58:51 +1000 Subject: [PATCH 21/28] Update main.py --- astrosource/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/main.py b/astrosource/main.py index d14f10b..e4f18df 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -55,7 +55,7 @@ -@click.option('--period', is_flag=True, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') +@click.option('--period', is_flag=True, default=True, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') @click.option('--periodlower', '-pl', type=float, default=-99.9, help='Shortest period to trial in days. Default is 0.05 days') @click.option('--periodupper', '-pu', type=float, default=-99.9, help='Longest period to trial in days. Default is one-third the observational baseline of your dataset.') @click.option('--periodtests', '-pt', type=int, default=10000, help='Number of different trial periods to run') From c2738de61ac3c10094d5854846fac56c285e4e07 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 13:00:31 +1000 Subject: [PATCH 22/28] Update main.py --- astrosource/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/main.py b/astrosource/main.py index e4f18df..483b2d5 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -55,7 +55,7 @@ -@click.option('--period', is_flag=True, default=True, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') +@click.option('--period', is_flag=True, default=False, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.') @click.option('--periodlower', '-pl', type=float, default=-99.9, help='Shortest period to trial in days. Default is 0.05 days') @click.option('--periodupper', '-pu', type=float, default=-99.9, help='Longest period to trial in days. Default is one-third the observational baseline of your dataset.') @click.option('--periodtests', '-pt', type=int, default=10000, help='Number of different trial periods to run') From c21a10af5b8ce685138f7c7d6f272f8a84bf0521 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 15:30:15 +1000 Subject: [PATCH 23/28] output targets list version 1 --- astrosource/main.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/astrosource/main.py b/astrosource/main.py index 483b2d5..c800592 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -217,15 +217,25 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree targets = get_targets(parentPath / 'results/potentialVariables.csv') else: targets=None - - - + if targets is not None: if full or phot: ts.photometry(filesave=True, targets=targets) if full or plot: ts.plot(detrend=detrend, period=period, eebls=eebls, filesave=True) - + + # Output the list of targets that was used in this run. + with open(csv_path, 'w', newline='') as csvfile: + writer = csv.writer(csvfile) + num_cols = len(targets[0]) + 1 + header = [f'col{i}' for i in range(1, num_cols + 1)] + writer.writerow(header) + + for idx, targ in enumerate(targets, start=1): + row = [f'V{idx}'] + list(targ) + writer.writerow(row) + + sys.stdout.write("✅ AstroSource analysis complete\n") except AstrosourceException as e: From db5f385461864377ad4541864bb9bc2cce5e0a5b Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 15:35:45 +1000 Subject: [PATCH 24/28] Update main.py --- astrosource/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/astrosource/main.py b/astrosource/main.py index c800592..35e0edb 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -225,7 +225,7 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree ts.plot(detrend=detrend, period=period, eebls=eebls, filesave=True) # Output the list of targets that was used in this run. - with open(csv_path, 'w', newline='') as csvfile: + with open((parentPath / 'results/TargetsUsed.csv', 'w', newline='') as csvfile: writer = csv.writer(csvfile) num_cols = len(targets[0]) + 1 header = [f'col{i}' for i in range(1, num_cols + 1)] From d93f7bd8e57e807f9bab74dc4b651b39a222c0da Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 15:41:44 +1000 Subject: [PATCH 25/28] output targets v3 --- astrosource/main.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/astrosource/main.py b/astrosource/main.py index 35e0edb..5687133 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -2,6 +2,7 @@ import click import sys import logging +import csv import os from numpy import array, all from astropy.utils.exceptions import AstropyWarning, AstropyDeprecationWarning @@ -225,15 +226,18 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree ts.plot(detrend=detrend, period=period, eebls=eebls, filesave=True) # Output the list of targets that was used in this run. - with open((parentPath / 'results/TargetsUsed.csv', 'w', newline='') as csvfile: + results_dir = parentPath / 'results' + #results_dir.mkdir(exist_ok=True) + + csv_path = results_dir / 'TargetsUsed.csv' + with open(csv_path, 'w', newline='') as csvfile: writer = csv.writer(csvfile) num_cols = len(targets[0]) + 1 header = [f'col{i}' for i in range(1, num_cols + 1)] writer.writerow(header) for idx, targ in enumerate(targets, start=1): - row = [f'V{idx}'] + list(targ) - writer.writerow(row) + writer.writerow([f'V{idx}'] + list(targ)) sys.stdout.write("✅ AstroSource analysis complete\n") From d42ebf17d5e601c39a1dcb64d382beec0fdeb397 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 16:11:13 +1000 Subject: [PATCH 26/28] update apass columns --- astrosource/comparison.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index d8c2d15..3bbf91a 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -546,7 +546,8 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): radecname = {'ra' :'raj2000', 'dec': 'dej2000'} if cat_name in ['APASS']: - searchColumns=[] + searchColumns=[searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] + ] else: searchColumns=[radecname['ra'], radecname['dec'], opt['filter'], opt['error'], opt['colmatch'], opt['colerr']] @@ -1680,4 +1681,4 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n # Needed for windows to multiprocess appropriately if __name__ == "__main__": - pass \ No newline at end of file + pass From 257f976cc54f4fd8fd277d9873624d9e6e8a2b55 Mon Sep 17 00:00:00 2001 From: Michael Fitzgerald Date: Sun, 15 Jun 2025 16:20:37 +1000 Subject: [PATCH 27/28] typo --- astrosource/comparison.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 3bbf91a..e0fc868 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -546,8 +546,7 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): radecname = {'ra' :'raj2000', 'dec': 'dej2000'} if cat_name in ['APASS']: - searchColumns=[searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] - ] + searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] else: searchColumns=[radecname['ra'], radecname['dec'], opt['filter'], opt['error'], opt['colmatch'], opt['colerr']] From 7d678cc8d8abd7b16873f8eca5034b8b553469cc Mon Sep 17 00:00:00 2001 From: mfitzasp Date: Sun, 15 Jun 2025 16:35:13 +1000 Subject: [PATCH 28/28] adjust targets output --- astrosource/astrosource.py | 17 +++++++++++++++++ astrosource/comparison.py | 4 ++-- astrosource/main.py | 14 -------------- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/astrosource/astrosource.py b/astrosource/astrosource.py index e4b0a83..80c97fe 100644 --- a/astrosource/astrosource.py +++ b/astrosource/astrosource.py @@ -4,6 +4,7 @@ import ssl import sys from pathlib import Path +import csv if (not os.environ.get('PYTHONHTTPSVERIFY', '') and getattr(ssl, '_create_unverified_context', None)): ssl._create_default_https_context = ssl._create_unverified_context @@ -240,6 +241,22 @@ def plot(self, detrend=False, eebls=False, period=False, phaseShift=0.0, filesav phased_plots(filterCode=self.filtercode, paths=self.paths, targets=self.targets, period=self.period, phaseShift=phaseShift) if eebls: plot_bls(paths=self.paths) + + # Output the list of targets that was used in this run. + #results_dir = parentPath / 'results' + #results_dir.mkdir(exist_ok=True) + + csv_path = self.paths['results'] / 'TargetsUsed.csv' + with open(csv_path, 'w', newline='') as csvfile: + writer = csv.writer(csvfile) + num_cols = len(self.targets[0]) + 1 + header = [f'col{i}' for i in range(1, num_cols + 1)] + writer.writerow(header) + + for idx, targ in enumerate(self.targets, start=1): + writer.writerow([f'V{idx}'] + list(targ)) + + def output(self, mode, data): output_files(paths=self.paths, photometrydata=data, mode=mode) diff --git a/astrosource/comparison.py b/astrosource/comparison.py index 3bbf91a..9962d4f 100644 --- a/astrosource/comparison.py +++ b/astrosource/comparison.py @@ -546,8 +546,8 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd): radecname = {'ra' :'raj2000', 'dec': 'dej2000'} if cat_name in ['APASS']: - searchColumns=[searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] - ] + searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')] + else: searchColumns=[radecname['ra'], radecname['dec'], opt['filter'], opt['error'], opt['colmatch'], opt['colerr']] diff --git a/astrosource/main.py b/astrosource/main.py index 5687133..18df835 100644 --- a/astrosource/main.py +++ b/astrosource/main.py @@ -224,20 +224,6 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree ts.photometry(filesave=True, targets=targets) if full or plot: ts.plot(detrend=detrend, period=period, eebls=eebls, filesave=True) - - # Output the list of targets that was used in this run. - results_dir = parentPath / 'results' - #results_dir.mkdir(exist_ok=True) - - csv_path = results_dir / 'TargetsUsed.csv' - with open(csv_path, 'w', newline='') as csvfile: - writer = csv.writer(csvfile) - num_cols = len(targets[0]) + 1 - header = [f'col{i}' for i in range(1, num_cols + 1)] - writer.writerow(header) - - for idx, targ in enumerate(targets, start=1): - writer.writerow([f'V{idx}'] + list(targ)) sys.stdout.write("✅ AstroSource analysis complete\n")