From 7088011d01b8f0ac4144c963d3cc40dde18cb6b2 Mon Sep 17 00:00:00 2001 From: EdwardBerman <10emb29@gmail.com> Date: Tue, 22 Oct 2024 09:06:48 -0400 Subject: [PATCH] move --- shopt/analyticLBFGS.jl | 101 +++++ shopt/argparser.jl | 14 + shopt/chisq.jl | 12 + shopt/dataOutprocessing.jl | 426 ++++++++++++++++++ shopt/dataPreprocessing.jl | 414 +++++++++++++++++ shopt/dependencies.jl | 42 ++ shopt/fancyPrint.jl | 27 ++ shopt/interpolate.jl | 149 ++++++ shopt/kaisserSquires.jl | 49 ++ shopt/masks.jl | 79 ++++ shopt/outliers.jl | 46 ++ shopt/oversample.jl | 35 ++ shopt/pca.jl | 20 + shopt/pixelGridAutoencoder.jl | 22 + shopt/plot.jl | 455 +++++++++++++++++++ shopt/powerSpectrum.jl | 16 + shopt/radialProfiles.jl | 27 ++ shopt/reader.jl | 87 ++++ shopt/shopt.jl | 820 ++++++++++++++++++++++++++++++++++ 19 files changed, 2841 insertions(+) create mode 100644 shopt/analyticLBFGS.jl create mode 100644 shopt/argparser.jl create mode 100644 shopt/chisq.jl create mode 100644 shopt/dataOutprocessing.jl create mode 100644 shopt/dataPreprocessing.jl create mode 100644 shopt/dependencies.jl create mode 100644 shopt/fancyPrint.jl create mode 100644 shopt/interpolate.jl create mode 100644 shopt/kaisserSquires.jl create mode 100644 shopt/masks.jl create mode 100644 shopt/outliers.jl create mode 100644 shopt/oversample.jl create mode 100644 shopt/pca.jl create mode 100644 shopt/pixelGridAutoencoder.jl create mode 100644 shopt/plot.jl create mode 100644 shopt/powerSpectrum.jl create mode 100644 shopt/radialProfiles.jl create mode 100644 shopt/reader.jl create mode 100644 shopt/shopt.jl diff --git a/shopt/analyticLBFGS.jl b/shopt/analyticLBFGS.jl new file mode 100644 index 0000000..7c59bfc --- /dev/null +++ b/shopt/analyticLBFGS.jl @@ -0,0 +1,101 @@ +include("radialProfiles.jl") + +#= +Functions for Cost and Gradient used in the Optimize step with LBFGS +NB: Reparameterization for [s, g1, g2] via [σ, e1, e2] to constraint update steps inside R+ x B_2(r) +=# + +function cost(params; r = r, c= c, starL=starCatalog[iteration], radial=fGaussian, AnalyticStampSize=AnalyticStampSize, get_middle_nxn=get_middle_nxn) + Totalcost = 0 + σ = params[1] + s_guess = σ^2 + e1_guess = params[2] + e2_guess = params[3] + ellipticity = sqrt((e1_guess)^2 + (e2_guess)^2) + normG = sqrt(1 + 0.5*( (1/ellipticity^2) - sqrt( (4/ellipticity^2)+ (1/ellipticity^4) ) )) + ratio = ellipticity/normG + g1_guess = e1_guess/ratio + g2_guess = e2_guess/ratio + + starL = get_middle_nxn(starL, AnalyticStampSize) + r = AnalyticStampSize + c = AnalyticStampSize + + sum = 0 + for u in 1:r + for v in 1:c + try + sum += radial(u,v, g1_guess, g2_guess, s_guess, r/2,c/2) + catch + sum += 0 + end + end + end + A_guess = 1/sum + + for u in 1:r + for v in 1:c + if isnan(starL[u,v]) + Totalcost += 0 + else + Totalcost += 0.5*(A_guess*radial(u, v, g1_guess, g2_guess, s_guess, r/2, c/2) - starL[u, v])^2 + end + end + end + return Totalcost +end + + +function costD(params; r=r, c=c, starL=pixelGridFits[iteration], radial=fGaussian, AnalyticStampSize=AnalyticStampSize, get_middle_nxn=get_middle_nxn) + Totalcost = 0 + σ = params[1] + s_guess = σ^2 + e1_guess = params[2] + e2_guess = params[3] + ellipticity = sqrt((e1_guess)^2 + (e2_guess)^2) + normG = sqrt(1 + 0.5*( (1/ellipticity^2) - sqrt( (4/ellipticity^2)+ (1/ellipticity^4) ) )) + ratio = ellipticity/normG + g1_guess = e1_guess/ratio + g2_guess = e2_guess/ratio + + starL = get_middle_nxn(starL, AnalyticStampSize) + r = AnalyticStampSize + c = AnalyticStampSize + + sum = 0 + for u in 1:r + for v in 1:c + try + sum += radial(u,v, g1_guess, g2_guess, s_guess, r/2,c/2) + catch + sum += 0 + end + end + end + A_guess = 1/sum + + for u in 1:r + for v in 1:c + if isnan(starL[u,v]) + Totalcost += 0 + else + Totalcost += 0.5*(A_guess*radial(u, v, g1_guess, g2_guess, s_guess, r/2, c/2) - starL[u, v])^2 + end + end + end + return Totalcost +end + +function g!(storage, p) + grad_cost = ForwardDiff.gradient(cost, p) + storage[1] = grad_cost[1] + storage[2] = grad_cost[2] + storage[3] = grad_cost[3] +end + +function gD!(storage, p) + grad_cost = ForwardDiff.gradient(costD, p) + storage[1] = grad_cost[1] + storage[2] = grad_cost[2] + storage[3] = grad_cost[3] +end diff --git a/shopt/argparser.jl b/shopt/argparser.jl new file mode 100644 index 0000000..08d7dd1 --- /dev/null +++ b/shopt/argparser.jl @@ -0,0 +1,14 @@ +#= +Function to parse arguments from the command line +=# + +function process_arguments(args) + fancyPrint("Parsing Arguments") + configdir = args[1] + println("━ Config Directory: ", configdir) + outdir = args[2] + println("━ Output Directory: ", outdir) + catalog = args[3] + println("━ Catalog: ", catalog) +end + diff --git a/shopt/chisq.jl b/shopt/chisq.jl new file mode 100644 index 0000000..c22ee7b --- /dev/null +++ b/shopt/chisq.jl @@ -0,0 +1,12 @@ +function chisq_cost(params; starL=starCatalog[iteration], weight_map=errVignets[iteration]) + pixel_values = params + chisq = sum(x -> isfinite(x) ? x : 0, (pixel_values .- vec(starL)) .^ 2 ./ vec(weight_map)) + return chisq +end + + +function chisq_g!(storage, p) + chisq_grad_cost = ForwardDiff.gradient(chisq_cost, p) + storage[1:length(chisq_grad_cost)] = chisq_grad_cost +end + diff --git a/shopt/dataOutprocessing.jl b/shopt/dataOutprocessing.jl new file mode 100644 index 0000000..f758db4 --- /dev/null +++ b/shopt/dataOutprocessing.jl @@ -0,0 +1,426 @@ +## .shopt file +using PyCall +using Dates + +#= +Functions to be called to write output to a fits file and add more diagnostic plots +=# + +function writeData(size, shear1, shear2, sizeD, shear1D, shear2D) + df = DataFrame(star = 1:length(size), + s_model=size, + g1_model=shear1, + g2_model=shear2, + s_data=sizeD, + g1_data=shear1D, + g2_data=shear2D) + + CSV.write(joinpath("outdir", "df.shopt"), df) +end + + +function readData() + DataFrame(CSV.File(joinpath("outdir", "df.shopt"))) +end + +function writeFitsData(truncate_summary_file = truncate_summary_file, g2C = g2C, g1C = g1C, sC = sC, summary_name=summary_name, sampled_indices=sampled_indices, s_model=s_model, g1_model=g1_model, g2_model=g2_model, s_data=s_data, g1_data=g1_data, g2_data=g2_data, u_coordinates = u_coordinates, v_coordinates = v_coordinates, PolynomialMatrix = PolynomialMatrix, outdir = outdir, configdir=configdir, starCatalog = starCatalog, pixelGridFits=pixelGridFits, fancyPrint=fancyPrint, training_stars=training_stars, training_u_coords=training_u_coords, training_v_coords=training_v_coords, validation_stars=validation_stars, validation_u_coords=validation_u_coords, validation_v_coords=validation_v_coords, validation_star_catalog=validation_star_catalog, degree=degree, YAMLSAVE=YAMLSAVE, parametersHistogram=parametersHistogram, parametersScatterplot=parametersScatterplot, cairomakiePlots=cairomakiePlots, pythonPlots=pythonPlots, outlier_indices=outlier_indices , failedStars=failedStars, s_blacklist=s_blacklist) + + m, n = size(starCatalog[1]) + array_3d = zeros(m, n, length(starCatalog)) + for (i, array) in enumerate(starCatalog) + array_3d[:, :, i] = array + end + starCatalog = array_3d + + m, n = size(pixelGridFits[1]) + array_3d = zeros(m, n, length(pixelGridFits)) + for (i, array) in enumerate(pixelGridFits) + array_3d[:, :, i] = array + end + pixelgridfits = array_3d + + #= + m, n = size(errVignets[1]) + array_3d = zeros(m, n, length(errVignets)) + for (i, array) in enumerate(errVignets) + array_3d[:, :, i] = array + end + errVignets = array_3d + =# + summary_name = summary_name*"summary.shopt" + + if truncate_summary_file + py""" + from astropy.io import fits + import numpy as np + + PolynomialMatrix = np.array($PolynomialMatrix, dtype=np.float64) + deg_element = $degree + degree_array = np.array([deg_element], dtype=np.float64) + + c = fits.Column(name='POLYNOMIAL_DEGREE', array=degree_array, format='D') + + sC = np.array($sC, dtype=np.float64) + g1C = np.array($g1C, dtype=np.float64) + g2C = np.array($g2C, dtype=np.float64) + + summary_statistics_hdu = fits.BinTableHDU.from_columns([c]) + summary_statistics_hdu.header['EXTNAME'] = 'SUMMARY_STATISTICS' + + hdu1 = fits.PrimaryHDU(PolynomialMatrix) + hdu1.header['EXTNAME'] = 'POLYNOMIAL_MATRIX' + + s_col = fits.Column(name='s_MATRIX', array=sC, format='D') + g1_col = fits.Column(name='g1_MATRIX', array=g1C, format='D') + g2_col = fits.Column(name='g2_MATRIX', array=g2C, format='D') + + analytic_profile_hdu = fits.BinTableHDU.from_columns([s_col, g1_col, g2_col]) + analytic_profile_hdu.header['EXTNAME'] = 'ANALYTIC_PROFILE' + + hdul = fits.HDUList([hdu1, summary_statistics_hdu, analytic_profile_hdu]) + py_summary_name = $summary_name + hdul.writeto(py_summary_name, overwrite=True) + + """ + else + + py""" + from astropy.io import fits + import numpy as np + + sC = np.array($sC, dtype=np.float64) + g1C = np.array($g1C, dtype=np.float64) + g2C = np.array($g2C, dtype=np.float64) + s_model = np.array($s_model, dtype=np.float64) + g1_model = np.array($g1_model, dtype=np.float64) + g2_model = np.array($g2_model, dtype=np.float64) + s_data = np.array($s_data, dtype=np.float64) + g1_data = np.array($g1_data, dtype=np.float64) + g2_data = np.array($g2_data, dtype=np.float64) + u_coordinates = np.array($u_coordinates, dtype=np.float64) + v_coordinates = np.array($v_coordinates, dtype=np.float64) + PolynomialMatrix = np.array($PolynomialMatrix, dtype=np.float64) + starCatalog = np.array($starCatalog, dtype=np.float64) + pixelGridFits = np.array($pixelGridFits, dtype=np.float64) + #errVignets = np.array($errVignets, dtype=np.float64) + #meanRelativeError = np.array($meanRelativeError, dtype=np.float64) + training_u_coords = np.array($training_u_coords, dtype=np.float64) + training_v_coords = np.array($training_v_coords, dtype=np.float64) + validation_u_coords = np.array($validation_u_coords, dtype=np.float64) + validation_v_coords = np.array($validation_v_coords, dtype=np.float64) + training_stars = np.array($training_stars, dtype=np.float64) + validation_stars = np.array($validation_stars, dtype=np.float64) + deg_element = $degree + degree_array = np.array([deg_element], dtype=np.float64) + + hdu1 = fits.PrimaryHDU(PolynomialMatrix) + hdu1.header['EXTNAME'] = 'POLYNOMIAL_MATRIX' + + c00 = fits.Column(name='U_COORDINATES', array=u_coordinates, format='D') + c01 = fits.Column(name='V_COORDINATES', array=v_coordinates, format='D') + c02 = fits.Column(name='TRAINING_U_COORDS', array=training_u_coords, format='D') + c03 = fits.Column(name='TRAINING_V_COORDS', array=training_v_coords, format='D') + c04 = fits.Column(name='VALIDATION_U_COORDS', array=validation_u_coords, format='D') + c05 = fits.Column(name='VALIDATION_V_COORDS', array=validation_v_coords, format='D') + c1 = fits.Column(name='S_MODEL', array=s_model, format='D') + c2 = fits.Column(name='g1_MODEL', array=g1_model, format='D') + c3 = fits.Column(name='g2_MODEL', array=g2_model, format='D') + c4 = fits.Column(name='s_DATA', array=s_data, format='D') + c5 = fits.Column(name='g1_DATA', array=g1_data, format='D') + c6 = fits.Column(name='g2_DATA', array=g2_data, format='D') + #c7 = fits.Column(name='MEAN_RELATIVE_ERROR', array=meanRelativeError, format='D') + c7 = fits.Column(name='POLYNOMIAL_DEGREE', array=degree_array, format='D') + + + VIGNETS_hdu = fits.ImageHDU(starCatalog) + VIGNETS_hdu.header['EXTNAME'] = 'VIGNETS' + + #errVignets_hdu = fits.ImageHDU(errVignets) + #errVignets_hdu.header['EXTNAME'] = 'ERR_VIGNETS' + + pixelGridFits_hdu = fits.ImageHDU(pixelGridFits) + pixelGridFits_hdu.header['EXTNAME'] = 'PIXEL_GRID_FITS' + + validation_hdu = fits.ImageHDU(validation_stars) + validation_hdu.header['EXTNAME'] = 'VALIDATION_STARS' + + s_col = fits.Column(name='s_MATRIX', array=sC, format='D') + g1_col = fits.Column(name='g1_MATRIX', array=g1C, format='D') + g2_col = fits.Column(name='g2_MATRIX', array=g2C, format='D') + + analytic_profile_hdu = fits.BinTableHDU.from_columns([s_col, g1_col, g2_col]) + analytic_profile_hdu.header['EXTNAME'] = 'ANALYTIC_PROFILE' + + summary_statistics_hdu = fits.BinTableHDU.from_columns([c00, c01, c02, c03, c04, c05, c1, c2, c3, c4, c5, c6, c7]) + summary_statistics_hdu.header['EXTNAME'] = 'SUMMARY_STATISTICS' + + outlier_indices = np.array($outlier_indices, dtype=np.float64) + s_blacklist = np.array($s_blacklist, dtype=np.float64) + failedStars = np.array($failedStars, dtype=np.float64) + col00 = fits.Column(name='FLAG_SNR', array=outlier_indices, format='D') + col01 = fits.Column(name='FLAG_S_BLACKLIST_MODEL', array=s_blacklist, format='D') + col02 = fits.Column(name='FLAG_PIXEL_GRID_OR_S_BLACKLIST_LEARNED', array=failedStars, format='D') + flag_hdu = fits.BinTableHDU.from_columns([col00, col01, col02]) + flag_hdu.header['EXTNAME'] = 'FLAGS' + + hdul = fits.HDUList([hdu1, summary_statistics_hdu, VIGNETS_hdu, pixelGridFits_hdu, validation_hdu, flag_hdu, analytic_profile_hdu]) + py_summary_name = $summary_name + hdul.writeto(py_summary_name, overwrite=True) + """ + end + + command1 = `mv $summary_name $outdir` + run(command1) + + if YAMLSAVE + command2 = `cp $configdir/shopt.yml $outdir` + run(command2) + + current_time = now() + command3 = `mv $outdir/shopt.yml $outdir/$(replace(Dates.format(Dates.now(), "YY:mm:dd:HH:MM:SS"), ":" => "_")*"_shopt.yml")` + + #command3 = `mv $outdir/shopt.yml $outdir/$(Dates.format(Time(current_time), "HH:MM:SS")*"_shopt.yml")` + run(command3) + end + + #command4 = `python $outdir/diagnostics.py` + #run(command4) + py_outdir = outdir * "/" * replace(Dates.format(Dates.now(), "YY:mm:dd:HH:MM:SS"), ":" => "_") + #py_outdir = outdir*"/$(now())" + if pythonPlots + fancyPrint("Producing Additional Python Plots") + else + py""" + import os + py_outdir = $py_outdir + + if not os.path.exists(py_outdir): + os.makedirs(py_outdir) + + """ + end + println("Creating ", py_outdir, " directory...") + + deg = degree + + if pythonPlots + py""" + from astropy.io import fits + import numpy as np + import matplotlib.pyplot as plt + import matplotlib.colors as colors + import matplotlib.cbook as cbook + from matplotlib import cm + from matplotlib import animation + import imageio + import os + from mpl_toolkits.axes_grid1 import make_axes_locatable + + + plt.ion() + + #Generalize this + python_current_outdir = $outdir + + f = fits.open(python_current_outdir + '/summary.shopt') + #f = fits.open('/home/eddieberman/research/mcclearygroup/shopt/outdir/summary.shopt') + + polyMatrix = f[0].data + #print(polyMatrix[30,30,:]) + deg = $deg + def polynomial_interpolation_star(u,v, polynomialMatrix): + r,c = np.shape(f[0].data)[0], np.shape(f[0].data)[1] + star = np.zeros((r,c)) + for i in range(r): + for j in range(c): + def objective_function(p, x, y, degree): + num_coefficients = (degree + 1) * (degree + 2) // 2 + value = 0 + counter = 0 + for a in range(1, degree + 2): + for b in range(1, degree + 2): + if (a - 1) + (b - 1) <= degree: + value += p[counter] * x**(a - 1) * y**(b - 1) + counter += 1 + return value + star[i,j] = objective_function(polynomialMatrix[i,j], u, v, deg) + + star = star/np.sum(star) + return star + + python_sampled_indices = $sampled_indices + sample_one = python_sampled_indices[0] + sample_two = python_sampled_indices[1] + sample_three = python_sampled_indices[2] + + a = polynomial_interpolation_star(f[1].data['validation_u_coords'][sample_one], f[1].data['validation_v_coords'][sample_one] ,polyMatrix) + #print(np.shape(a)) + fig, axes = plt.subplots(2, 3) + + # Display the first image in the first subplot + axes[0, 0].imshow(f[5].data[sample_one, :, : ], norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[0, 0].set_title('Pixel Grid Fit') + + # Display the second image in the second subplot + axes[0, 1].imshow(a, norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[0, 1].set_title('Polynomial Interpolation') + + axes[0, 2].imshow(f[5].data[sample_one, :, : ] - a, norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[0, 2].set_title('Residuals') + + b = polynomial_interpolation_star(f[1].data['validation_u_coords'][sample_two], f[1].data['validation_v_coords'][sample_two], polyMatrix) + + axes[1, 0].imshow(f[5].data[sample_two, :, : ], norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[1, 0].set_title('Pixel Grid Fit') + + # Display the second image in the second subplot + axes[1, 1].imshow(b, norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[1, 1].set_title('Polynomial Interpolation') + + axes[1, 2].imshow(f[5].data[sample_two, :, : ] - b, norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes[1, 2].set_title('Residuals') + # Adjust the spacing between subplots + plt.tight_layout() + plt.show() + + + + vignets = f[2].data + a,b,c = np.shape(f[4].data) + pixelGrid = np.zeros((b,c,a)) + + for i in range(a): + pixelGrid[:,:,i] = f[4].data[i,:,:] + + #print(np.shape(pixelGrid)) + #print(np.shape(vignets)) + + fig2, axes2 = plt.subplots(1, 3) + axes2[0].imshow(vignets[:,:,0], norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes2[0].set_title('vignets') + axes2[1].imshow(pixelGrid[:,:,0], norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes2[1].set_title('pixel grid') + axes2[2].imshow(vignets[:,:,0] - pixelGrid[:,:,0], norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes2[2].set_title('Log Scale Residuals') + plt.tight_layout() + plt.show() + + + def meanRelativeError(vignets, pixelGrid): + meanRelativeError = np.zeros((vignets.shape[0], vignets.shape[1])) + for j in range(vignets.shape[0]): + for k in range(vignets.shape[1]): + RelativeError = [] + for i in range(vignets.shape[2]): + RelativeError.append((vignets[j,k,i] - pixelGrid[j,k,i]) / (vignets[j,k,i] + 1e-10)) + meanRelativeError[j,k] = np.nanmean(RelativeError) + return meanRelativeError + + fig3, axes3 = plt.subplots(1) + im = axes3.imshow(meanRelativeError(vignets, pixelGrid), cmap=plt.cm.bwr_r, norm=colors.SymLogNorm(linthresh=1*10**(-6))) + axes3.set_title('Mean Relative Error') + divider = make_axes_locatable(axes3) + cax = divider.append_axes("right", size="5%", pad=0.05) + fig3.colorbar(im, cax=cax) + + fft_image = np.fft.fft(vignets[:,:,0] - pixelGrid[:,:,0]) + fft_image = np.abs(fft_image) ** 2 + + pk = [] + for i in range(1, 11): + radius = np.linspace(1, max(vignets.shape[0]/2, vignets.shape[1]/2) - 1, num=10) + radiusPixels = [] + for u in range(fft_image.shape[0]): + for v in range(fft_image.shape[1]): + if round(np.sqrt((u - fft_image.shape[0]/2)**2 + (v - fft_image.shape[1]/2)**2) - radius[i-1]) == 0: + radiusPixels.append(fft_image[u, v]) + pk.append(np.nanmean(radiusPixels)) + + + fig4, axes4 = plt.subplots(1,2) + axes4[0].imshow(fft_image) + axes4[0].set_title('FFT of Residuals') + axes4[1].plot(radius, pk) + axes4[1].set_title('Power Spectra') + + def update(frame): + im.set_array(polyMatrix[:, :, frame]) # Update the plot data for each frame + return im, + + fig5 = plt.figure() + ax5 = fig5.add_subplot(111) + im = ax5.imshow(polyMatrix[:, :, 0], cmap=plt.cm.bwr_r, norm=colors.SymLogNorm(linthresh=1*10**(-6)), interpolation='nearest') + + frames = range(np.shape(polyMatrix)[2]) # Number of frames + animation = animation.FuncAnimation(fig5, update, frames=frames, interval=200, blit=True) + + def mreHist(vignets, pixelGrid): + RelativeError = [] + for j in range(vignets.shape[0]): + for k in range(vignets.shape[1]): + for i in range(vignets.shape[2]): + RelativeError.append((vignets[j,k,i] - pixelGrid[j,k,i]) / (vignets[j,k,i]) + 1e-10) + return RelativeError + + fig6, axes6 = plt.subplots(1, figsize=(30, 10)) + bins = np.linspace(-1, 1, 21) + axes6.hist(mreHist(vignets, pixelGrid), bins = bins, color = "lightblue", ec="red", lw=3) #bins = bins + axes6.set_xlabel('Relative Error', fontsize=20) + axes6.set_ylabel('Frequency',fontsize=20) + axes6.set_yscale('log') + axes6.set_title(f'Relative Error Histogram (logscale)\n $\mu = {np.nanmean(mreHist(vignets, pixelGrid))}$, $\sigma = {np.nanstd(mreHist(vignets, pixelGrid))/np.sqrt(len(mreHist(vignets, pixelGrid)) - np.sum(np.isnan(mreHist(vignets, pixelGrid))) )}$', fontsize=30) + + py_outdir = $py_outdir + + if not os.path.exists(py_outdir): + os.makedirs(py_outdir) + + fig.savefig(os.path.join(py_outdir,'pgfVsInterp.png')) + fig2.savefig(os.path.join(py_outdir, 'logV_P_R.png')) + fig3.savefig(os.path.join(py_outdir, 'MRE.png')) + fig4.savefig(os.path.join(py_outdir, 'fftPk.png')) + fig6.savefig(os.path.join(py_outdir, 'RelativeErrorHistogram.png')) + #fig6.savefig(os.path.join(py_outdir, 'RelativeErrorHistogram.pdf')) + + """ + end + + command4 = `mv $outdir/$(replace(Dates.format(Dates.now(), "YY:mm:dd:HH:MM:SS"), ":" => "_")*"_shopt.yml") $py_outdir` + run(command4) + + command5 = `cp $outdir/$summary_name $py_outdir` + run(command5) + + if cairomakiePlots + command6 = `mv s_uv.png $py_outdir` + run(command6) + + command7 = `mv g1_uv.png $py_outdir` + run(command7) + + command8 = `mv g2_uv.png $py_outdir` + run(command8) + end + + if parametersScatterplot + command9 = `mv parametersScatterplot.png $py_outdir` + run(command9) + + command10 = `mv parametersScatterplot.pdf $py_outdir` + run(command10) + end + + if parametersHistogram + command11 = `mv parametersHistogram.pdf $py_outdir` + run(command11) + + command12 = `mv parametersHistogram.png $py_outdir` + run(command12) + end + + # run on sampled indices, copy diagnostics.py to py""" """ here +end + +# function P(u,v) +# PSF at uv +# end diff --git a/shopt/dataPreprocessing.jl b/shopt/dataPreprocessing.jl new file mode 100644 index 0000000..f1f7614 --- /dev/null +++ b/shopt/dataPreprocessing.jl @@ -0,0 +1,414 @@ +#= +Load in Config shopt.yml +=# +nt = nthreads() +fancyPrintTwo("Julia is using $nt threads to run ShOpt.jl. To change this, on UNIX run export JULIA_NUM_THREADS=[Integer] and on Windows set JULIA_NUM_THREADS=[Integer].") +fancyPrintTwo("For example: export JULIA_NUM_THREADS=4 or set JULIA_NUM_THREADS=4") +config = YAML.load_file(joinpath(configdir, "shopt.yml")) +epochs = config["NNparams"]["epochs"] +degree = config["polynomialDegree"] +new_img_dim = config["stampSize"] +r, c = new_img_dim, new_img_dim +snrCutoff = config["dataProcessing"]["SnRPercentile"] +YAMLSAVE = config["saveYaml"] +minAnalyticGradientModel = config["AnalyticFitParams"]["minGradientAnalyticModel"] +minAnalyticGradientLearned = config["AnalyticFitParams"]["minGradientAnalyticLearned"] +AnalyticStampSize = config["AnalyticFitParams"]["analyticFitStampSize"] +minPixelGradient = config["NNparams"]["minGradientPixel"] +UnicodePlotsPrint = config["plots"]["unicodePlots"] +pythonPlots = config["plots"]["pythonPlots"] +parametersHistogram = config["plots"]["normalPlots"]["parametersHistogram"] +parametersScatterplot = config["plots"]["normalPlots"]["parametersScatterplot"] +cairomakiePlots = config["plots"]["cairomakiePlots"]["streamplots"] +sLowerBound = config["dataProcessing"]["sLowerBound"] +sUpperBound = config["dataProcessing"]["sUpperBound"] +comments = config["CommentsOnRun"] +unity_sum = config["sum_pixel_grid_and_inputs_to_unity"] +training_ratio = config["training_ratio"] +summary_name = config["summary_name"] +mode = config["mode"] # Options: auotoencoder, lanczos +PCAterms = config["PCAterms"] +polynomial_interpolation_stopping_gradient = config["polynomial_interpolation_stopping_gradient"] +lanczos = config["lanczos"] +truncate_summary_file = config["truncate_summary_file"] +iterationsPolyfit = config["iterations"] +alpha = config["alpha"] +chisq_stopping_gradient = config["chisq_stopping_gradient"] +median_constraint = config["median_constraint"] +outlier_sigma_clipping = config["outlier_sigma_clipping"] +#= +Log these config choices +=# + +println("Key Config Choices:") +if mode == "chisq" + println("━ Mode: χ2") +else + println("━ Mode: ", mode) +end +println("━ α:", alpha) +println("━ Chisq Stopping Gradient:", chisq_stopping_gradient) +println("━ Iterations: ", iterationsPolyfit) +println("━ Median Constraint: ", median_constraint) +println("━ Summary Name: ", summary_name) +println("━ PCA Terms: ", PCAterms) +println("━ Lanczos n: ", lanczos) +println("━ Polynomial Interpolation Stopping Gradient: ", polynomial_interpolation_stopping_gradient) +println("━ Max Epochs: ", epochs) +println("━ Polynomial Degree: ", degree) +println("━ Stamp Size: ", new_img_dim) +println("━ Signal to Noise Ratio Cutoff: ", snrCutoff) +println("━ Save YAML: ", YAMLSAVE) +println("━ Stopping Analytic Fit Gradient Star Vignets: ", minAnalyticGradientModel) +println("━ Stopping Analytic Fit Gradient Learned Vignets: ", minAnalyticGradientLearned) +println("━ Analytic Fit Stamp Size: ", AnalyticStampSize) +println("━ Stopping Pixel Fit Gradient: ", minPixelGradient) +println("━ Print Unicode Plots: ", UnicodePlotsPrint) +println("━ Python Plots: ", pythonPlots) +println("━ Parameters Histogram: ", parametersHistogram) +println("━ Parameters Scatterplot: ", parametersScatterplot) +println("━ Stream Plots: ", cairomakiePlots) +println("━ s Lower Bound: ", sLowerBound) +println("━ s Upper Bound: ", sUpperBound) +println("━ Training Ratio: ", training_ratio) +println("━ Sum Pixel Grid and Inputs to Unity: ", unity_sum) +println("━ Truncate Summary File: ", truncate_summary_file) +println("━ Outlier Sigma Clipping: ", outlier_sigma_clipping) +println("━ Comments: ", comments) + +#= +Utility Function used to count NaNs in an image, was used in testing and may be useful for future debugging in the data preprocessing process +=# +function countNaNs(arr) + count = 0 + + for element in arr + if isnan(element) + count += 1 + end + end + + return count +end + +function sinc(x) + if x == 0 + return 1 + else + return sin(x) / x + end +end + +function lanczos_kernel(x, a) + if abs(x) >= a + return 0 + elseif x == 0 + return 1 + else + return a * sinc(pi*x) * sinc(pi*x/a) + end +end + +function smooth(data, a) + kernel = [lanczos_kernel(i, a) * lanczos_kernel(j, a) for i in -a:a, j in -a:a] + kernel = kernel ./ sum(kernel) # Normalize the kernel + + nx, ny = size(data) + smoothed_data = copy(data) # Initialize with original data to keep NaN positions + + for i in 1:nx + for j in 1:ny + if isnan(data[i, j]) + continue # Skip smoothing for NaN values + end + norm_factor = 0.0 # This is to normalize considering non-NaN values + smoothed_value = 0.0 + for di in max(1, i-a):min(nx, i+a) + for dj in max(1, j-a):min(ny, j+a) + if !isnan(data[di, dj]) + smoothed_value += kernel[di-i+a+1, dj-j+a+1] * data[di, dj] + norm_factor += kernel[di-i+a+1, dj-j+a+1] + end + end + end + smoothed_data[i, j] = smoothed_value / norm_factor + end + end + + return smoothed_data +end + +function oversample_image(image, new_dim) + oversampled_image = zeros(new_dim, new_dim) + scale_factor = 1/(new_dim / size(image)[1]) + for y in 1:new_dim + for x in 1:new_dim + src_x = (x - 0.5) * scale_factor + 0.5 + src_y = (y - 0.5) * scale_factor + 0.5 + x0 = max(1, floor(Int, src_x)) + x1 = min(size(image)[1], x0 + 1) + y0 = max(1, floor(Int, src_y)) + y1 = min(size(image)[2], y0 + 1) + dx = src_x - x0 + dy = src_y - y0 + oversampled_image[y, x] = + (1 - dx) * (1 - dy) * image[y0, x0] + + dx * (1 - dy) * image[y0, x1] + + (1 - dx) * dy * image[y1, x0] + + dx * dy * image[y1, x1] + end + end + return oversampled_image +end + +function undersample_image(image, new_dim) + undersampled_image = imresize(image, new_dim, new_dim, algorithm =:bicubic) + return undersampled_image +end + +function sample_image(image, new_dim) + if new_dim > size(image)[1] + return oversample_image(image, new_dim) + else + return undersample_image(image, new_dim) + end +end +#= +function undersample_image(image, new_dim) + # Get the original dimensions of the image + height, width = size(image) + + # Create Lanczos interpolation object + interp = interpolate(image, Lanczos(4)) + + # Calculate the scaling factors + scale_height = (height - 1) / (new_dim - 1) + scale_width = (width - 1) / (new_dim - 1) + + # Generate coordinates for the new image + x_coords = 1:new_dim + y_coords = 1:new_dim + + # Perform downsampling with Lanczos interpolation + undersampled_image = zeros(new_dim, new_dim) + + for y in 1:new_dim + for x in 1:new_dim + undersampled_image[y, x] = interp(scale_height * (y - 1) + 1, scale_width * (x - 1) + 1) + end + end + + return undersampled_image +end +=# +#= +Supply an array and the new dimension you want and get the middle nxn of that array +=# +function get_middle_nxn(array, n) + rows, cols = size(array) + + # Calculate the starting row and column indices + row_start = div(rows - n, 2) + 1 + col_start = div(cols - n, 2) + 1 + + # Return the sub-matrix + return array[row_start:(row_start + n - 1), col_start:(col_start + n - 1)] +end + +#= +function get_middle_nxn(array, n) + rows, cols = size(array) + if isodd(rows) + if isodd(n) + row_start = div(rows,2) - (n ÷ 2) + col_start = div(cols,2) - (n ÷ 2) + return array[row_start:(row_start + (2*(n ÷ 2))), col_start:(col_start + (2*(n ÷ 2)))] + else + array = array[1:(rows - 1), 1:(cols - 1)] + rows, cols = size(array) + row_start = div(rows,2) - (n ÷ 2) + col_start = div(cols,2) - (n ÷ 2) + return array[(1 + row_start):(row_start + (2*(n ÷ 2))), ( 1 + col_start):(col_start + (2*(n ÷ 2)))] + end + else + if isodd(n) + row_start = div(rows,2) - (n ÷ 2) + col_start = div(cols,2) - (n ÷ 2) + return array[row_start:(row_start + (2*(n ÷ 2))), col_start:(col_start + (2*(n ÷ 2)))] + else + row_start = div(rows,2) - (n ÷ 2) + col_start = div(cols,2) - (n ÷ 2) + return array[(1 + row_start):(row_start + (2*(n ÷ 2))), (1 + col_start):(col_start + (2*(n ÷ 2)))] + end + end +end +=# + +#= +SnR function used to calculate the signal to noise ratio of an image +=# + +function signal_to_noise(I, e; nm=nanMask, nz=nanToZero) + signal_power = sum(nz(nm(I)).^2) + noise_power = sum(e.^2) + snr = 10*log10(signal_power/noise_power) + return snr +end + +#= +Consider +In [10]: len(f[2].data['DELTAWIN_J2000']) +Out[10]: 290 +=# + +function cataloging(args; nm=nanMask, nz=nanToZero, snr=signal_to_noise, dout=outliers_filter, alpha=alpha) + catalog = args[3] + + py""" + import numpy as np + from astropy.io import fits + + python_datadir = $catalog + mode = $mode + print(python_datadir) + alpha_default = $alpha + print(alpha_default) + if mode == "chisq": + print("Making error cutouts...") + f = fits.open(python_datadir) + def find_extension_with_colname(fits_file, target_colname): + with fits.open(fits_file) as hdulist: + matching_extensions = [] + + for idx, hdu in enumerate(hdulist): + if isinstance(hdu, fits.BinTableHDU) or isinstance(hdu, fits.TableHDU): + if target_colname in hdu.columns.names: + matching_extensions.append(idx) + + return matching_extensions + idx = find_extension_with_colname(python_datadir, 'VIGNET')[0] + + vignets = f[idx].data['VIGNET'] + if mode == "chisq": + backgrounds = f[idx].data['BACKGROUND'] + sigma_bs = [np.sqrt(np.var(background)) for background in backgrounds] + + + def sigma_squared_i(pi, g, sigma_b, alpha): + return sigma_b + pi / g + (alpha * pi)**2 + + def compute_sigma_i_for_vignet(vignet, sigma_b, g=4, alpha=alpha_default, eta=1.0): + sigma_i_cutout = np.zeros_like(vignet, dtype=float) + for i in range(vignet.shape[0]): + for j in range(vignet.shape[1]): + pi = vignet[i, j] + sigma_i_cutout[i, j] = sigma_squared_i(pi, g, sigma_b, alpha) + return sigma_i_cutout + + # Compute sigma_i for each vignette using its respective sigma_b + if mode == "chisq": + err_vignets = [compute_sigma_i_for_vignet(vignet, sigma_b) for vignet, sigma_b in zip(vignets, sigma_bs)] + #err_vignets = f[idx].data['ERR_VIGNET'] + l = len(vignets) + + u = f[idx].data['ALPHAWIN_J2000'] + v = f[idx].data['DELTAWIN_J2000'] + + snr = f[idx].data['SNR_WIN'] + median_img = np.nanmedian(f[idx].data['VIGNET'], axis=0) + """ + + datadir = py"python_datadir" + v = py"vignets" + median_img = py"median_img" + + if mode == "chisq" + err = py"err_vignets" + end +# print("━ Number of vignets: ", length(v)) + catalog = py"list(map(np.array, $v))" + + if mode == "chisq" + errVignets = py"list(map(np.array, $err))" + end + + u_coords = convert(Array{Float64,1}, py"u") + v_coords = convert(Array{Float64,1}, py"v") + snr = convert(Array{Float32,1}, py"snr") + + r = size(catalog[1], 1) + c = size(catalog[1], 2) + + catalogNew = [] + signal2noiseRatios = [] + for i in 1:length(catalog) + if unity_sum + push!(catalogNew, nm(catalog[i])./sum(nz(nm(catalog[i])))) + else + push!(catalogNew, nm(catalog[i])) + end + #push!(catalogNew, nm(catalog[i])./sum(nz(nm(catalog[i])))) + #push!(signal2noiseRatios, snr(catalog[i], errVignets[i])) + end + + new_snr_arr = Array{Float64}(undef, length(snr)) + for (i, element) in enumerate(snr) + new_snr_arr[i] = element + end + new_snr_arr = new_snr_arr[.!isnan.(new_snr_arr)] + + println(UnicodePlots.boxplot(["snr"], [new_snr_arr], title="signal to noise ratio")) + + + catalogNew, outlier_indices = dout(snr, catalogNew, snrCutoff) + + if mode == "chisq" + errVignets = [errVignets[i] for i in 1:length(errVignets) if i ∉ outlier_indices] + end + + println("━ Number of vignets: ", length(catalog)) + println("━ Removed $(length(catalog) - length(catalogNew)) outliers on the basis of Signal to Noise Ratio") + #= + for i in 1:length(catalogNew) + catalogNew[i] = get_middle_nxn(catalogNew[i], new_img_dim) + errVignets[i] = get_middle_nxn(errVignets[i], new_img_dim) + end + println("━ Sampled all vignets to $(new_img_dim) x $(new_img_dim) from $(r) x $(c) via cropping") + r = size(catalogNew[1], 1) + c = size(catalogNew[1], 2) + =# + + if new_img_dim/size(catalogNew[1], 1) !=1 + if new_img_dim/size(catalogNew[1], 1) < 1 + for i in 1:length(catalogNew) + catalogNew[i] = smooth(nm(get_middle_nxn(catalogNew[i], new_img_dim))/sum(nz(nm(get_middle_nxn(catalogNew[i], new_img_dim)))), lanczos) + if mode == "chisq" + errVignets[i] = get_middle_nxn(errVignets[i], new_img_dim) + end + end + else + for i in 1:length(catalogNew) + catalogNew[i] = smooth(nm(oversample_image(catalogNew[i], new_img_dim))/sum(nz(nm(oversample_image(catalogNew[i], new_img_dim)))), lanczos) + if mode == "chisq" + errVignets[i] = oversample_image(errVignets[i], new_img_dim) + end + end + end + end + + println("━ Sampled all vignets to $(size(catalogNew[1], 1)) x $(size(catalogNew[1], 2)) from $(r) x $(c) via cropping") + r = size(catalogNew[1], 1) + c = size(catalogNew[1], 2) + k = rand(1:length(catalogNew)) + if UnicodePlotsPrint + println(UnicodePlots.heatmap(nz(nm(get_middle_nxn(catalogNew[k],15))), colormap=:inferno, title="Sampled Vignet $k")) + end + + if mode == "chisq" + return catalogNew, r, c, length(catalogNew), u_coords, v_coords, outlier_indices, median_img, errVignets + else + return catalogNew, r, c, length(catalogNew), u_coords, v_coords, outlier_indices, median_img + end +end + + diff --git a/shopt/dependencies.jl b/shopt/dependencies.jl new file mode 100644 index 0000000..fe187c0 --- /dev/null +++ b/shopt/dependencies.jl @@ -0,0 +1,42 @@ +using Pkg + +function is_installed(package::AbstractString) + try + Pkg.installed(package) + return true + catch + return false + end +end + +function download_packages(file_path::AbstractString) + # Read the file and retrieve the import statements + file = open(file_path, "r") + imports = readlines(file) + close(file) + + # Download and import the packages + for package in imports + package = strip(package) + if !is_installed(package) + try + Pkg.add(package) + println("Successfully downloaded and added $package") + catch err + println("Failed to download and add $package: $err") + continue + end + end + + try + eval(Meta.parse("import $package")) + println("Successfully imported $package") + catch err + println("Failed to import $package: $err") + end + end +end + +# Provide the file path of the imports.txt file +download_packages("imports.txt") + diff --git a/shopt/fancyPrint.jl b/shopt/fancyPrint.jl new file mode 100644 index 0000000..fc663e9 --- /dev/null +++ b/shopt/fancyPrint.jl @@ -0,0 +1,27 @@ +#= +Making Log Statements Look Nicer, fancyPrints signify stepping stones in the program + +=# + +function fancyPrint(printStatement) + for i in 1:(length(printStatement) + 4) + print("-") + end + println("\n| "*printStatement*" |") + for i in 1:(length(printStatement) + 3) + print("-") + end + println("-") +end + +function fancyPrintTwo(printStatement) + for i in 1:(length(printStatement) + 4) + print("#") + end + println("\n# "*printStatement*" #") + for i in 1:(length(printStatement) + 3) + print("#") + end + println("#") +end + diff --git a/shopt/interpolate.jl b/shopt/interpolate.jl new file mode 100644 index 0000000..bdb47db --- /dev/null +++ b/shopt/interpolate.jl @@ -0,0 +1,149 @@ +#= +Helper functions to interpolate s(u,v), g1(u,v), g2(u,v) across the field of view +Reminder to make s, g1, g2 (u,v) of arbitrary degree! +=# + +#f(u, v) = a1u^3 + a2v^3 + a3u^2v + a4v^2u + a5u^2 + a6v^2 + a7uv + a8u + a9v + a10 +function interpCostS(p; truth=s_tuples) + I(u, v) = p[1]*u^3 + p[2]*v^3 + p[3]*u^2*v + p[4]*v^2*u + p[5]*u^2 + p[6]*v^2 + p[7]*u*v + p[8]*u + p[9]*v + p[10] + t = truth + function sumLoss(f, t) + totalLoss = 0 + for i in 1:length(t) #t = [(u,v, I), ... ] + totalLoss += (f(t[i][1], t[i][2]) - t[i][3])^2 + end + return totalLoss + end + return sumLoss(I, t) +end + +function polyG_s!(storage, p) + grad_cost = ForwardDiff.gradient(interpCostS, p) #∇ + storage[1] = grad_cost[1] + storage[2] = grad_cost[2] + storage[3] = grad_cost[3] + storage[4] = grad_cost[4] + storage[5] = grad_cost[5] + storage[6] = grad_cost[6] + storage[7] = grad_cost[7] + storage[8] = grad_cost[8] + storage[9] = grad_cost[9] + storage[10] = grad_cost[10] +end + +function interpCostg1(p; truth=g1_tuples) + I(u, v) = p[1]*u^3 + p[2]*v^3 + p[3]*u^2*v + p[4]*v^2*u + p[5]*u^2 + p[6]*v^2 + p[7]*u*v + p[8]*u + p[9]*v + p[10] + t = truth + function sumLoss(f, t) + totalLoss = 0 + for i in 1:length(t) #t = [(u,v, I), ... ] + totalLoss += (f(t[i][1], t[i][2]) - t[i][3])^2 + end + return totalLoss + end + return sumLoss(I, t) +end + +function polyG_g1!(storage, p) + grad_cost = ForwardDiff.gradient(interpCostg1, p) + storage[1] = grad_cost[1] + storage[2] = grad_cost[2] + storage[3] = grad_cost[3] + storage[4] = grad_cost[4] + storage[5] = grad_cost[5] + storage[6] = grad_cost[6] + storage[7] = grad_cost[7] + storage[8] = grad_cost[8] + storage[9] = grad_cost[9] + storage[10] = grad_cost[10] +end + +function interpCostg2(p; truth=g2_tuples) + + I(u, v) = p[1]*u^3 + p[2]*v^3 + p[3]*u^2*v + p[4]*v^2*u + p[5]*u^2 + p[6]*v^2 + p[7]*u*v + p[8]*u + p[9]*v + p[10] + t = truth + function sumLoss(f, t) + totalLoss = 0 + for i in 1:length(t) #t = [(u,v, I), ... ] + totalLoss += (f(t[i][1], t[i][2]) - t[i][3])^2 + end + return totalLoss + end + return sumLoss(I, t) +end + +function polyG_g2!(storage, p) + grad_cost = ForwardDiff.gradient(interpCostg2, p) + storage[1] = grad_cost[1] + storage[2] = grad_cost[2] + storage[3] = grad_cost[3] + storage[4] = grad_cost[4] + storage[5] = grad_cost[5] + storage[6] = grad_cost[6] + storage[7] = grad_cost[7] + storage[8] = grad_cost[8] + storage[9] = grad_cost[9] + storage[10] = grad_cost[10] +end + +function objective_function(p, x, y, degree) + num_coefficients = (degree + 1) * (degree + 2) ÷ 2 + value = 0 + counter = 0 + for i in 1:(degree + 1) + for j in 1:(degree + 1) + if (i - 1) + (j - 1) <= degree + counter += 1 + value += p[counter] * x^(i - 1) * y^(j - 1) + end + end + end + return value +end + + +function polynomial_optimizer(degree, x_data, y_data, z_data; median_constraint=median_constraint) + function generate_polynomial_terms(n) + terms = [] + for i in 0:n + for j in 0:(n-i) + push!(terms, (i, j)) + end + end + terms + end + polynomial_terms = generate_polynomial_terms(degree) + A = zeros(length(z_data), length(polynomial_terms)) + #replace nans with median + if median_constraint == true + A[:, 1] .= median(filter(!isnan, z_data)) + for i in 1:length(z_data) + for (index, (p, q)) in enumerate(polynomial_terms[2:end]) # Skip the constant term + A[i, index+1] = x_data[i]^p * y_data[i]^q + end + end + else + for i in 1:length(z_data) + if isnan(z_data[i]) + z_data[i] = median(filter(!isnan, z_data)) + end + end + for i in 1:length(z_data) + for (index, (p, q)) in enumerate(polynomial_terms) + A[i, index] = x_data[i]^p * y_data[i]^q + end + end + end + + coefficients = A \ z_data + if median_constraint == true + coefficients[1] = median(filter(!isnan, z_data)) + end + if length(z_data) < 5 + coefficients = zeros(length(coefficients)) + coefficients[1] = median(filter(!isnan, z_data)) + end + return coefficients +end + +# if statement that makes a new polynomial_optimizer function that has a median constraint. Objective function will only be a function of p-1 terms, will add median(z_data) for constant term diff --git a/shopt/kaisserSquires.jl b/shopt/kaisserSquires.jl new file mode 100644 index 0000000..c766251 --- /dev/null +++ b/shopt/kaisserSquires.jl @@ -0,0 +1,49 @@ +function truncate_to_square(mat1::Matrix{T}, mat2::Matrix{T}) where {T} + # Step 1: Get the size of the matrices + size1 = size(mat1) + size2 = size(mat2) + + # Step 2: Determine the minimum dimension + min_dim = min(size1[1], size1[2], size2[1], size2[2]) + + # Step 3: Create new square matrices + square_mat1 = Matrix{T}(undef, min_dim, min_dim) + square_mat2 = Matrix{T}(undef, min_dim, min_dim) + + # Step 4: Assign corresponding elements + for i in 1:min_dim, j in 1:min_dim + square_mat1[i, j] = mat1[i, j] + square_mat2[i, j] = mat2[i, j] + end + + return square_mat1, square_mat2 +end + + +function meshgrid(x, y) + X = [i for i in x, j in 1:length(y)] + Y = [j for i in 1:length(x), j in y] + return X, Y +end + +function ks(g1Map, g2Map; meshgrid = meshgrid, ts = truncate_to_square) + g1Map, g2Map = ts(g1Map, g2Map) + x, y = size(g1Map) + k1, k2 = meshgrid(fftfreq(y), fftfreq(x)) + + g1_hat = fft(g1Map) + g2_hat = fft(g2Map) + + p1 = k1 * k1 - k2 * k2 + p2 = 2 * k1 * k2 + k2 = k1 * k1 + k2 * k2 + k2[1,1] = 1 + + kEhat = (p1 * g1_hat + p2 * g2_hat) / k2 + kBhat = -(p2 * g1_hat - p1 * g2_hat) / k2 + + kE = real.(ifft(kEhat)) + kB = real.(ifft(kBhat)) + + return kE, kB +end diff --git a/shopt/masks.jl b/shopt/masks.jl new file mode 100644 index 0000000..7ecad62 --- /dev/null +++ b/shopt/masks.jl @@ -0,0 +1,79 @@ +#= +Masks for NaN Handeling, used in dataPreprocessing and for Pixel Grid Fitting +=# + +function nanMask(arr) + dummyArr = zeros(size(arr,1),size(arr,2)) + for i in 1:size(arr,1) + for j in 1:size(arr,2) + if arr[i,j] < -1000 + dummyArr[i,j] = NaN + else + dummyArr[i,j] = arr[i,j] + end + end + end + return dummyArr +end + +#For Purposes of Summing to unity +function nanToGaussian(arr, s, g1, g2, uc, vc, func=fGaussian) + dummyArr = zeros(size(arr,1),size(arr,2)) + for i in 1:size(arr,1) + for j in 1:size(arr,2) + if isnan(arr[i,j]) + mu = 0.0 + sigma = 1 + normal_dist = Normal(mu, sigma) + truncated_dist = Truncated(normal_dist, -0.1*func(i, j, s, g1, g2, uc, vc), 0.1*func(i, j, s, g1, g2, uc, vc)) + dummyArr[i,j] = func(i, j, s, g1, g2, uc, vc) + rand(truncated_dist) + else + dummyArr[i,j] = arr[i,j] + end + end + end + return dummyArr +end + +function nanToZero(arr) + dummyArr = zeros(size(arr,1),size(arr,2)) + for i in 1:size(arr,1) + for j in 1:size(arr,2) + if isnan(arr[i,j]) + dummyArr[i,j] = 0.001 + else + dummyArr[i,j] = arr[i,j] + end + end + end + return dummyArr +end +function nanToInf(arr) + dummyArr = zeros(size(arr,1),size(arr,2)) + for i in 1:size(arr,1) + for j in 1:size(arr,2) + if isnan(arr[i,j]) + dummyArr[i,j] = -1*Inf + else + dummyArr[i,j] = arr[i,j] + end + end + end + return dummyArr +end + +function nanMask2(arr) + dummyArr = zeros(size(arr,1),size(arr,2)) + for i in 1:size(arr,1) + for j in 1:size(arr,2) + if arr[i,j] < 0 + dummyArr[i,j] = NaN + else + dummyArr[i,j] = arr[i,j] + end + end + end + return dummyArr +end + +#note to self add nan to sky background diff --git a/shopt/outliers.jl b/shopt/outliers.jl new file mode 100644 index 0000000..42ab666 --- /dev/null +++ b/shopt/outliers.jl @@ -0,0 +1,46 @@ +#= +Functions for detecting and removing outliers from data, used in dataPreprocessing.jl +=# + +function detect_outliers(data::AbstractVector{T}; k::Float64=1.5) where T<:Real + q1 = quantile(data, 0.25) + q3 = quantile(data, 0.75) + iqr = q3 - q1 + lower_fence = q1 - k * iqr + upper_fence = q3 + k * iqr + filter = (data .< lower_fence) .| (data .> upper_fence) + outliers = data[filter] + return outliers +end + +function outliers_filter(snr::Vector{Float32}, img::Vector{Any}, k::Float64) + #snr_cleaned = filter(!isnan, snr) + q1 = quantile(snr, k) + #snr = [isnan(x) ? -1000 : x for x in snr] + println("━ Cutting off Stars below the $k Percentile of Signal to Noise Ratio: $q1 ") + q3 = quantile(snr, 0.75) + iqr = q3 - q1 + lower_fence = q1 - k * iqr + upper_fence = q3 + k * iqr + outlier_indices = findall(x -> x < q1, snr) #outlier_indices = findall(x -> x < lower_fence || x > upper_fence, snr) + println("━ outlier indices: $outlier_indices") + img_snr_cleaned = img + #wht_snr_cleaned = wht + for i in sort(outlier_indices, rev=true) + splice!(img_snr_cleaned, i) + #splice!(wht_snr_cleaned, i) + end + return img_snr_cleaned, outlier_indices +end + +function remove_outliers(data::AbstractVector{T}; k::Float64=1.5) where T<:Real + q1 = quantile(data, 0.25) + q3 = quantile(data, 0.75) + iqr = q3 - q1 + lower_fence = q1 - k * iqr + upper_fence = q3 + k * iqr + filter = (data .> lower_fence) .& (data .< upper_fence) + nonOutliers = data[filter] + return nonOutliers +end + diff --git a/shopt/oversample.jl b/shopt/oversample.jl new file mode 100644 index 0000000..d62dde0 --- /dev/null +++ b/shopt/oversample.jl @@ -0,0 +1,35 @@ +using Images +function oversample_image(image, new_dim) + # Create an empty array to store the oversampled image + oversampled_image = zeros(new_dim, new_dim) + + scale_factor = 1/(new_dim / size(image)[1]) + + for y in 1:new_dim + for x in 1:new_dim + src_x = (x - 0.5) * scale_factor + 0.5 + src_y = (y - 0.5) * scale_factor + 0.5 + + x0 = max(1, floor(Int, src_x)) + x1 = min(size(image)[1], x0 + 1) + y0 = max(1, floor(Int, src_y)) + y1 = min(size(image)[2], y0 + 1) + + dx = src_x - x0 + dy = src_y - y0 + + oversampled_image[y, x] = + (1 - dx) * (1 - dy) * image[y0, x0] + + dx * (1 - dy) * image[y0, x1] + + (1 - dx) * dy * image[y1, x0] + + dx * dy * image[y1, x1] + end + end + + return oversampled_image +end + +function undersample_image(image, new_dim) + undersampled_image = imresize(image, new_dim, new_dim, algorithm = :bilinear) + return undersampled_image +end diff --git a/shopt/pca.jl b/shopt/pca.jl new file mode 100644 index 0000000..2ba42bc --- /dev/null +++ b/shopt/pca.jl @@ -0,0 +1,20 @@ +function pca_image(image, ncomponents) + #Load img Matrix + img_matrix = image + + # Perform PCA + M = fit(PCA, img_matrix; maxoutdim=ncomponents) + + # Transform the image into the PCA space + transformed = MultivariateStats.transform(M, img_matrix) + + # Reconstruct the image + reconstructed = reconstruct(M, transformed) + + # Reshape the image back to its original shape + reconstructed_image = reshape(reconstructed, size(img_matrix)...) + + return reconstructed_image +end + + diff --git a/shopt/pixelGridAutoencoder.jl b/shopt/pixelGridAutoencoder.jl new file mode 100644 index 0000000..8378a88 --- /dev/null +++ b/shopt/pixelGridAutoencoder.jl @@ -0,0 +1,22 @@ +# Encoder +encoder = Chain( + Dense(r*c, 128, leakyrelu), + Dense(128, 64, leakyrelu), + Dense(64, 32, leakyrelu), + ) +#Decoder +decoder = Chain( + Dense(32, 64, leakyrelu), + Dense(64, 128, leakyrelu), + Dense(128, r*c, tanh), + ) +#Full autoencoder +autoencoder = Chain(encoder, decoder) + +#x̂ = autoencoder(x) +loss(x) = mse(autoencoder(x), x) +# + msle(shift_positive(x), shift_positive(autoencoder(x))) + +# Define the optimizer +optimizer = ADAM() + diff --git a/shopt/plot.jl b/shopt/plot.jl new file mode 100644 index 0000000..1000aaf --- /dev/null +++ b/shopt/plot.jl @@ -0,0 +1,455 @@ +#= +Julia Plots +=# + +function error_plot(model, learned, errorModel, errorLearned, t) + # Generate example data + x = [1, 5, 9] # Three points on the x-axis + labels = ["s", "g1", "g2"] # Labels for the points + model_data = model + learned_data = learned + errM = errorModel + errD = errorLearned + + # Plot scatter plot with error bars + Plots.scatter(x, model_data, yerr=errM, label="Model Data", color="blue", alpha=0.6) + Plots.xticks!(x, labels) + Plots.scatter!(x, learned_data, yerr=errD, label="Learned Data", color="red", markersize=4) + Plots.xticks!(x, labels) + + Plots.title!(t) + + # Show the plot + Plots.plot!(legend=:topright, titlefontsize=30, xguidefontsize=30, yguidefontsize=30, margin = 25mm) +end + +function plot_err(s_model=s_model, g1_model=g1_model, g2_model=g2_model, s_data=s_data, g1_data=g1_data, g2_data=g2_data, parametersScatterplot=parametersScatterplot) + + n = size(s_model, 1) + if parametersScatterplot + ps1 = error_plot([mean(s_model), mean(g1_model), mean(g2_model)], + [mean(s_data), mean(g1_data), mean(g2_data)], + [std(s_model)/sqrt(n), std(g1_model)/sqrt(n), std(g2_model)/sqrt(n)], + [std(s_data)/sqrt(n), std(g1_data)/sqrt(n), std(g2_data)/sqrt(n)], + "Learned vs True Parameters") + Plots.savefig(Plots.plot(ps1, size = (1920, 1080)), "parametersScatterplot.pdf") + Plots.savefig(Plots.plot(ps1, size = (1920, 1080)), "parametersScatterplot.png") + end +end + +function hist_1(x::Array{T, 1}, y::Array{T, 1}, t1, t2) where T<:Real + bin_edges = range(-1, stop=1, step=0.1) + Plots.histogram(x, + alpha=0.5, + bins=bin_edges, + label=t1, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30) + Plots.histogram!(y, + alpha=0.5, + bins=bin_edges, + label=t2, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30) + Plots.plot!(margin=15mm) +end +function hist_2(x::Array{T, 1}, y::Array{T, 1}, t1, t2) where T<:Real + bin_edges = range(-1, stop=1, step=0.1) + Plots.histogram(x, + alpha=0.5, + bins=bin_edges, + label=t1, + xticks = -1:0.1:1, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30) + Plots.histogram!(y, + alpha=0.5, + bins=bin_edges, + label=t2, + xticks = -1:0.1:1, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30) + Plots.plot!(margin=15mm) +end + +function plot_hist(s_model=s_model, g1_model=g1_model, g2_model=g2_model, s_data=s_data, g1_data=g1_data, g2_data=g2_data, hist_1=hist_1, hist_2=hist_2, parametersHistogram=parametersHistogram) + if parametersHistogram + hist1 = hist_1(s_model, s_data, "s Model", "s Data") + hist2 = hist_2(g1_model, g1_data, "g1 Model", "g1 Data") + hist3 = hist_2(g2_model, g2_data, "g2 Model", "g2 Data") + bin_edges = range(-1, stop=1, step=0.1) + hist4 = Plots.histogram(s_model - s_data, + alpha=0.5, + label="S Model and Data Residuals") + hist5 = Plots.histogram(g1_model - g1_data, + alpha=0.5, + closed =:both, + bins=range(-2, stop=2, step=0.1), + label="g1 Model and Data Residuals", + xticks = -1:0.1:1) + hist6 = Plots.histogram(g2_model - g2_data, + alpha=0.5, + closed =:both, + bins=range(-2, stop=2, step=0.1), + label="g2 Model and Data Residuals", + xticks = -1:0.1:1) + + Plots.savefig(Plots.plot(hist1, + hist2, + hist3, + hist4, + hist5, + hist6, + layout = (2,3), + size = (1920,1080)), + "parametersHistogram.pdf") + + Plots.savefig(Plots.plot(hist1, + hist2, + hist3, + hist4, + hist5, + hist6, + layout = (2,3), + size = (1920,1080)), + "parametersHistogram.png") + end +end + +function generate_heatmap(data::Array{T, 2}, cbmin, cbmax) where T<:Real + Plots.heatmap(data, + aspect_ratio=:equal, # ensure square cells + color=:winter, # use the Winter color map + cbar=:true, # add the color bar + xlabel="u", + ylabel="v", + clims=(cbmin, cbmax), # set the color limits + xlims=(0.5, size(data, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(data, 1) + 0.5), # set the y-axis limits to include the full cells + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30, + margin = 25mm + ) +end + +function generate_heatmap_sp_t(data::Array{T, 2}, t, cbmin, cbmax) where T<:Real + Plots.heatmap(data, + aspect_ratio=:equal, # ensure square cells + color=:winter, # use the Winter color map + cbar=:true, # add the color bar + title = t, + clims=(cbmin, cbmax), # set the color limits + xlims=(0.5, size(data, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(data, 1) + 0.5), # set the y-axis limits to include the full cells + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + margin = 15mm + ) +end + +function generate_heatmap_sp_t2(data::Array{T, 2}, t, cbmin, cbmax) where T<:Real + Plots.heatmap(data, + aspect_ratio=:equal, # ensure square cells + color=:winter, # use the Winter color map + cbar=:true, # add the color bar + title = t, + clims=(cbmin, cbmax), # set the color limits + xlims=(0.5, size(data, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(data, 1) + 0.5), # set the y-axis limits to include the full cells + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + margin = 15mm, + size=(1920, 1080) + ) +end + +function generate_heatmap_sp(data::Array{T, 2}, cbmin, cbmax) where T<:Real + Plots.heatmap(data, + aspect_ratio=:equal, # ensure square cells + color=:winter, # use the Winter color map + cbar=:true, # add the color bar + clims=(cbmin, cbmax), # set the color limits + xlims=(0.5, size(data, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(data, 1) + 0.5), # set the y-axis limits to include the full cells + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + margin = 15mm + ) +end + +function generate_heatmap_titled(data::Array{T, 2}, t, cbmin, cbmax) where T<:Real + Plots.heatmap(data, + aspect_ratio=:equal, # ensure square cells + color=:winter, # use the Winter color map + cbar=:true, # add the color bar + xlabel="u", + ylabel="v", + title = t, + clims=(cbmin, cbmax), # set the color limits + xlims=(0.5, size(data, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(data, 1) + 0.5), # set the y-axis limits to include the full cells + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30, + margin = 25mm + ) +end + +#= +function plot_hm(r=r, c=c, outdir=outdir, sampled_indices=sampled_indices, starCatalog=starCatalog, pixelGridFits=pixelGridFits, get_middle_15x15=get_middle_15x15, generate_heatmap_sp_t=generate_heatmap_sp_t, generate_heatmap_sp=generate_heatmap_sp,generate_heatmap_titled=generate_heatmap_titled)#p + star1 = sampled_indices[1] + star2 = sampled_indices[2] + star3 = sampled_indices[3] + #p_value = p + chiSquareTemplate = zeros(r, c) + chiSquare = [] + + for i in 1:3 + starCatalog[sampled_indices[i]] = Float64.(starCatalog[sampled_indices[i]]) + end + + for i in 1:3 + chiSquareTemplate = Float64.(starCatalog[sampled_indices[i]]) - pixelGridFits[sampled_indices[i]].^2 + chiSquareTemplate = chiSquareTemplate./(errVignets[sampled_indices[i]].^2) + #err cutout , divide matrices pointwise + push!(chiSquare, chiSquareTemplate) + end + + + hm11 = generate_heatmap_sp_t(get_middle_15x15(starCatalog[star1]), + "Model", + minimum([minimum(starCatalog[star1]), minimum(pixelGridFits[star1])]), + maximum([maximum(starCatalog[star1]), maximum(pixelGridFits[star1])])) + hm12 = generate_heatmap_sp_t(get_middle_15x15(pixelGridFits[star1]), + "Data", + minimum([minimum(starCatalog[star1]), minimum(pixelGridFits[star1])]), + maximum([maximum(starCatalog[star1]), maximum(pixelGridFits[star1])])) + hm13 = generate_heatmap_sp_t(get_middle_15x15(chiSquare[1]), + "ChiSquare", + minimum(chiSquare[1]), + maximum(chiSquare[1])) + + hm21 = generate_heatmap_sp(get_middle_15x15(starCatalog[star2]), + minimum([minimum(starCatalog[star2]), minimum(pixelGridFits[star2])]), + maximum([maximum(starCatalog[star2]), maximum(pixelGridFits[star2])])) + + hm22 = generate_heatmap_sp(get_middle_15x15(pixelGridFits[star2]), + minimum([minimum(starCatalog[star2]), minimum(pixelGridFits[star2])]), + maximum([maximum(starCatalog[star2]), maximum(pixelGridFits[star2])])) + + hm23 = generate_heatmap_sp(get_middle_15x15(chiSquare[2]), + minimum(chiSquare[2]), + maximum(chiSquare[2])) + + hm31 = generate_heatmap_sp(get_middle_15x15(starCatalog[star3]), + minimum([minimum(starCatalog[star3]), minimum(pixelGridFits[star3])]), + maximum([maximum(starCatalog[star3]), maximum(pixelGridFits[star3])])) + + hm32 = generate_heatmap_sp(get_middle_15x15(pixelGridFits[star3]), + minimum([minimum(starCatalog[star3]), minimum(pixelGridFits[star3])]), + maximum([maximum(starCatalog[star3]), maximum(pixelGridFits[star3])])) + + hm33 = generate_heatmap_sp(get_middle_15x15(chiSquare[3]), + minimum(chiSquare[3]), + maximum(chiSquare[3])) +#= + hm41 = generate_heatmap_sp(starCatalog[star4], minimum([minimum(starCatalog[star1]), minimum(pixelGridFits[star4])]), maximum([maximum(starCatalog[star4]), maximum(pixelGridFits[star4])])) + hm42 = generate_heatmap_sp(pixelGridFits[star4], minimum([minimum(starCatalog[star4]), minimum(pixelGridFits[star4])]), maximum([maximum(starCatalog[star4]), maximum(pixelGridFits[star4])])) + hm43 = generate_heatmap_sp(chiSquare[4], minimum(chiSquare[4]), maximum(chiSquare[4])) + + hm51 = generate_heatmap_sp(starCatalog[star5], minimum([minimum(starCatalog[star5]), minimum(pixelGridFits[star5])]), maximum([maximum(starCatalog[star5]), maximum(pixelGridFits[star5])])) + hm52 = generate_heatmap_sp(pixelGridFits[star5], minimum([minimum(starCatalog[star5]), minimum(pixelGridFits[star5])]), maximum([maximum(starCatalog[star5]), maximum(pixelGridFits[star5])])) + hm53 = generate_heatmap_sp(chiSquare[5], minimum(chiSquare[5]), maximum(chiSquare[5])) + + hm61 = generate_heatmap_sp(starCatalog[star6], minimum([minimum(starCatalog[star6]), minimum(pixelGridFits[star6])]), maximum([maximum(starCatalog[star6]), maximum(pixelGridFits[star6])])) + hm62 = generate_heatmap_sp(pixelGridFits[star6], minimum([minimum(starCatalog[star6]), minimum(pixelGridFits[star6])]), maximum([maximum(starCatalog[star6]), maximum(pixelGridFits[star6])])) + hm63 = generate_heatmap_sp(chiSquare[6], minimum(chiSquare[6]), maximum(chiSquare[6])) +=# + #= + fft_image1 = fft(complex.(starCatalog[star1] .- pixelGridFits[star1])) + fft_image1 = abs2.(fft_image1) + + pk1 = [] + for i in 1:10 + radius = range(1, max(r/2, c/2) - 1, length=10) + push!(pk1, powerSpectrum(fft_image1, radius[i])) + end + + fftmin1 = minimum(fft_image1) + fftmax1 = maximum(fft_image1) + + fft_image2 = fft(complex.(starCatalog[star2] .- pixelGridFits[star2])) + fft_image2 = abs2.(fft_image2) + + pk2 = [] + for i in 1:10 + radius = range(1, max(r/2, c/2) - 1, length=10) + push!(pk2, powerSpectrum(fft_image2, radius[i])) + end + + fftmin2 = minimum(fft_image2) + fftmax2 = maximum(fft_image2) + + fft_image3 = fft(complex.(starCatalog[star3] .- pixelGridFits[star3])) + fft_image3 = abs2.(fft_image3) + + pk3 = [] + for i in 1:10 + radius = range(1, max(r/2, c/2) - 1, length=10) + push!(pk3, powerSpectrum(fft_image3, radius[i])) + end + + fftmin3 = minimum(fft_image3) + fftmax3 = maximum(fft_image3) + + hm8 = generate_heatmap_sp_t(fft_image1, "FFT Residuals", fftmin1, fftmax1) + hm9 = generate_heatmap_sp(fft_image2, fftmin2, fftmax2) + hm10 = generate_heatmap_sp(fft_image3, fftmin3, fftmax3) + =# + zmn = minimum([minimum(starCatalog[star3]), minimum(pixelGridFits[star3])]) + zmx = maximum([maximum(starCatalog[star3]), maximum(pixelGridFits[star3])]) + + s1 = Plots.surface(title="Model Star", titlefontsize=30, starCatalog[star3], zlim=(zmn, zmx), colorbar = false) + s2 = Plots.surface(title="pixelGridFit", titlefontsize=30, pixelGridFits[star3], zlim=(zmn, zmx), colorbar = false) + s3 = Plots.surface(title="Residuals", titlefontsize=30, starCatalog[star3] .- pixelGridFits[star3], zlim=(2*zmn, 2*zmx), colorbar = false) + + meanRelativeError = zeros(r,c) + for j in 1:size(starCatalog[star1], 1) + for k in 1:size(starCatalog[star1], 2) + RelativeError = [] + for i in 1:length(starCatalog) + push!(RelativeError, abs.(starCatalog[i][j,k] .- pixelGridFits[i][j,k]) ./ abs.(starCatalog[i][j,k] .+ 1e-10)) + end + meanRelativeError[j,k] = mean(RelativeError) + end + end + + RelativeErrorHM = generate_heatmap_sp_t2(meanRelativeError, "Mean Relative Error at Each Pixel", minimum(meanRelativeError), maximum(meanRelativeError)) + + + + pk_1 = Plots.plot(pk1, + xlabel = "k", + linewidth=5, + tickfontsize=16, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30, + ylabel = "P(k)", + title = "Power Spectrum", + margin = 25mm) + pk_2 = Plots.plot(pk2, + xlabel = "k", + linewidth=5, + tickfontsize=16, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30, + ylabel = "P(k)", + margin = 25mm) + pk_3 = Plots.plot(pk3, + xlabel = "k", + linewidth=5, + tickfontsize=16, + titlefontsize=30, + xguidefontsize=30, + yguidefontsize=30, + ylabel = "P(k)", + margin = 25mm) + + #title1 = Plots.plot(title = "I(u,v) Model", titlefontsize=30, grid = false, showaxis = false, bottom_margin = -25Plots.px) + #title2 = Plots.plot(title = "I(u,v) Data", titlefontsize=30, grid = false, showaxis = false, bottom_margin = -25Plots.px) + #title3 = Plots.plot(title = "Chi-Square Residuals", titlefontsize=30, grid = false, showaxis = false, bottom_margin = -25Plots.px) + + Plots.savefig(Plots.plot(hm11, hm12, hm13, + hm21, hm22, hm23, + hm31, hm32, hm33, + layout = (3,3), + widths=ones(9), + size = (1920,1080)), + joinpath("outdir", "pixelGridFit.pdf")) + + Plots.savefig(Plots.plot(hm11, hm12, hm13, + hm21, hm22, hm23, + hm31, hm32, hm33, + layout = (3,3), + widths=ones(9), + size = (1920,1080)), + joinpath("outdir", "pixelGridFit.png")) + + Plots.savefig(RelativeErrorHM, joinpath("outdir", "relativeErrorHM.pdf")) + Plots.savefig(RelativeErrorHM, joinpath("outdir", "relativeErrorHM.png")) + Plots.savefig(Plots.plot(s1, + s2, + s3, + layout = (1,3), + size = (1920,1080)), + joinpath(outdir, "3dPixelGridFit.png")) + + Plots.savefig(Plots.plot(pk_1, + pk_2, + pk_3, + layout = (3,1), + size = (1920,1080)), + joinpath(outdir, "powerspectra.png")) + + Plots.savefig(Plots.plot(hm8, + hm9, + hm10, + layout = (3,1), + size = (1920,1080)), + joinpath(outdir, "fftresiduals.png")) + Plots.savefig(Plots.plot(pk_1, + pk_2, + pk_3, + layout = (3,1), + size = (1920,1080)), + joinpath(outdir, "powerspectra.pdf")) + + Plots.savefig(Plots.plot(hm8, + hm9, + hm10, + layout = (3,1), + size = (1920,1080)), + joinpath(outdir, "fftresiduals.pdf")) + +end +=# +#= +Revisit Later? + +=# +#= +#= +661 scale = 1/0.29 +662 ks93, k0 = ks(g1_map, g2_map) +663 ksCosmos = get_middle_15x15(imfilter(ks93, Kernel.gaussian(scale))) +664 kshm = Plots.heatmap(ksCosmos, +665 title="Kaisser-Squires", +666 xlabel="u", +667 ylabel="v", +668 xlims=(0.5, size(ksCosmos, 2) + 0.5), # set the x-axis limits to include the full cells +669 ylims=(0.5, size(ksCosmos, 1) + 0.5), # set the y-axis limits to include the full cells +670 aspect_ratio=:equal, +671 ticks=:none, # remove the ticks +672 frame=:box, # draw a box around the plot +673 grid=:none, # remove the grid lines +674 size=(1920,1080)) +675 +676 Plots.savefig(kshm, joinpath("outdir","kaisserSquires.png")) +677 =# +=# diff --git a/shopt/powerSpectrum.jl b/shopt/powerSpectrum.jl new file mode 100644 index 0000000..1db751e --- /dev/null +++ b/shopt/powerSpectrum.jl @@ -0,0 +1,16 @@ +#= +Computers the Power Spectra at a given annulus, called iteratively for a Power Spectra map +=# + +function powerSpectrum(data::Array{T, 2}, radius) where T<:Real + radiusPixels = [] + for u in 1:size(data,1) + for v in 1:size(data,2) + if round(sqrt((u - size(data,1)/2)^2 + (v - size(data,2)/2)^2) - radius) == 0 1 + push!(radiusPixels, data[u,v]) + end + end + end + return mean(radiusPixels) + +end diff --git a/shopt/radialProfiles.jl b/shopt/radialProfiles.jl new file mode 100644 index 0000000..09da5e7 --- /dev/null +++ b/shopt/radialProfiles.jl @@ -0,0 +1,27 @@ +#= +Gaussian and kolmogorov radial profiles +=# + +function fGaussian(u, v, g1, g2, s, uc, vc) + squareRadiusFactor = ([(1 + g1) g2; g2 (1 - g1)]*[(u - uc); (v- vc)])'*([(1 + g1) g2; g2 (1 - g1)]*[(u-uc);(v-vc)]) + matrixScalar = s^2/(1 - g1^2 - g2^2) + return exp(-(matrixScalar*squareRadiusFactor)) +end + +function Fkolmogorov(k, λ) + r0 = 1 + return exp(-(24*gamma(6.5)*0.2)^(5/6) * (λ*k/(2*π*r0)^(5/3))) +end + +function fkolmogorov(u, v, g1, g2, s, uc, vc, Fk=Fkolmogorov) + λ = 1 + squareRadiusFactor = ([(1 + g1) g2; g2 (1 - g1)]*[(u - uc); (v- vc)])'*([(1 + g1) g2; g2 (1 - g1)]*[(u-uc);(v-vc)]) + matrixScalar = s^2/(1 - g1^2 - g2^2) + squareRadius = squareRadiusFactor*matrixScalar + integrand(k) = Fk(k, λ) * besselj(0, k*(squareRadius^(0.5)) )* k + result, err = quadgk(integrand, 0, Inf) + return (1/(2*π)) * result +end + + + diff --git a/shopt/reader.jl b/shopt/reader.jl new file mode 100644 index 0000000..83f8b79 --- /dev/null +++ b/shopt/reader.jl @@ -0,0 +1,87 @@ +using PyCall + +#= +Functions for reading in summary.shopt file and returning the PSF at an arbitrary (u,v) +=# + +#= +Example: +polyMatrix, deg = read_shopt("summary.shopt") +p(0.5,0.5,polyMatrix, deg) +=# + +function read_shopt(shoptFile) + shoptFile = string(shoptFile) + py""" + shoptFile = $shoptFile + from astropy.io import fits + f = fits.open(shoptFile) + polyMatrix = f[0].data + degree = f[1].data['POLYNOMIAL_DEGREE'][0] + try: + s_matrix = f[6].data['s_MATRIX'] + g1_matrix = f[6].data['g1_MATRIX'] + g2_matrix = f[6].data['g2_MATRIX'] + except: + g1_matrix = f[2].data['g1_MATRIX'] + g2_matrix = f[2].data['g2_MATRIX'] + s_matrix = f[2].data['s_MATRIX'] + f.close() + """ + polynomialMatrix = convert(Array{Float64,3}, py"polyMatrix") + degree = convert(Int64, py"degree") + s_matrix = convert(Array{Float64,1}, py"s_matrix") + g1_matrix = convert(Array{Float64,1}, py"g1_matrix") + g2_matrix = convert(Array{Float64,1}, py"g2_matrix") + + return polynomialMatrix, degree, s_matrix, g1_matrix, g2_matrix +end + +function p(u,v, polMatrix, degree) + #Read in degree of Polynomial + #Read in Size of Star Vignet + r,c = size(polMatrix,1), size(polMatrix,2) + star = zeros(r,c) + + for a in 1:r + for b in 1:c + + function objective_function(p, x, y, degree) + value = 0 + counter = 0 + + for i in 1:(degree + 1) + for j in 1:(degree + 1) + if (i - 1) + (j - 1) <= degree + counter += 1 + value += p[counter] * x^(i - 1) * y^(j - 1) + end + end + end + return value + end + star[a,b] = objective_function(polMatrix[a,b,:] ,u,v,degree) # -> Filename Polynomial Matrix + end + end + star = star/sum(star) + return star +end + +function analytic_profile(u,v, s_matrix, g1_matrix, g2_matrix, radial_function) + s = s_matrix[1]*u^3 + s_matrix[2]*v^3 + s_matrix[3]*u^2*v + s_matrix[4]*v^2*u + s_matrix[5]*u^2 + s_matrix[6]*v^2 + s_matrix[7]*u*v + s_matrix[8]*u + s_matrix[9]*v + s_matrix[10] + g1 = g1_matrix[1]*u^3 + g1_matrix[2]*v^3 + g1_matrix[3]*u^2*v + g1_matrix[4]*v^2*u + g1_matrix[5]*u^2 + g1_matrix[6]*v^2 + g1_matrix[7]*u*v + g1_matrix[8]*u + g1_matrix[9]*v + g1_matrix[10] + g2 = g2_matrix[1]*u^3 + g2_matrix[2]*v^3 + g2_matrix[3]*u^2*v + g2_matrix[4]*v^2*u + g2_matrix[5]*u^2 + g2_matrix[6]*v^2 + g2_matrix[7]*u*v + g2_matrix[8]*u + g2_matrix[9]*v + g2_matrix[10] + return radial_function(s,g1,g2) +end + +function export_bibtex() + println("@misc{berman2023shoptjl,\n + title={ShOpt.jl: A Julia Package for Empirical Point Spread Function Characterization of JWST NIRCam Data},\n + author={Edward Berman and Jacqueline McCleary},\n + year={2023},\n + eprint={2310.00071},\n + archivePrefix={arXiv},\n + primaryClass={astro-ph.IM}\n +} +") +end diff --git a/shopt/shopt.jl b/shopt/shopt.jl new file mode 100644 index 0000000..599de8f --- /dev/null +++ b/shopt/shopt.jl @@ -0,0 +1,820 @@ +# ---------------------------------------------------------# +@time begin + include("argparser.jl") + include("fancyPrint.jl") + + try + process_arguments(ARGS) + catch err + println("Error: ", err) + println("Usage: julia shopt.jl ") + exit(1) + end + + configdir = ARGS[1] + outdir = ARGS[2] + catalog = ARGS[3] + + if isdir(outdir) + println("━ Outdir found") + else + println("━ Outdir not found, creating...") + mkdir(outdir) + end +end + +# ---------------------------------------------------------# +fancyPrint("Handling Imports") +@time begin + using Base: initarray! + using YAML + using BenchmarkTools + using Plots + using ForwardDiff + using LinearAlgebra + using PyCall + using Random + using Distributions + using SpecialFunctions + using Optim + using IterativeSolvers + using QuadGK + using DataFrames + using FFTW + using CSV + using Images, ImageFiltering + using ImageTransformations + using Measures + using ProgressBars + using UnicodePlots + using Flux + using Flux.Optimise + using Flux.Losses + using Flux: onehotbatch, throttle, mse, msle, mae + using CairoMakie + using Dates + using MultivariateStats + using Base.Threads + using MultivariatePolynomials + #using Interpolations +end +println("━ Start Time ", Dates.now()) +start = Dates.now() +# ---------------------------------------------------------# +fancyPrint("Reading .jl Files") +@time begin + include("plot.jl") + include("analyticLBFGS.jl") + include("radialProfiles.jl") + include("masks.jl") + include("outliers.jl") + include("dataOutprocessing.jl") + include("powerSpectrum.jl") + include("kaisserSquires.jl") + include("dataPreprocessing.jl") + include("interpolate.jl") + include("pixelGridAutoencoder.jl") + include("pca.jl") + include("chisq.jl") + include("reader.jl") + #include("lk.jl") +end +# ---------------------------------------------------------# +#fancyPrint("Running Source Extractor") +# ---------------------------------------------------------# + +fancyPrint("Processing Data for Fit") +@time begin + + if mode == "chisq" + starCatalog, r, c, itr, u_coordinates, v_coordinates, outlier_indices, median_img, errVignets = cataloging(ARGS) + else + starCatalog, r, c, itr, u_coordinates, v_coordinates, outlier_indices, median_img = cataloging(ARGS) + end + #starCatalog = starCatalog + #errVignets = errVignets + #u_coordinates = u_coordinates + #v_coordinates = v_coordinates + itr = length(starCatalog) + + + #= + starCatalog, r,c, itr = catalogingWEBBPSF() + u_coordinates = rand(2) + v_coordinates = rand(2) + + itr = length(starCatalog) + =# + + #= + starCatalog, r, c, u_coordinates, v_coordinates = gridPSFS() #return catalogNew, rows, cols, u_coords, v_coords + errVignets = starCatalog + itr = length(starCatalog) + =# +end + +starData = zeros(r, c) + +#A_model = zeros(itr) +s_model = zeros(itr) +g1_model = zeros(itr) +g2_model = zeros(itr) + +#A_data = zeros(itr) +s_data = zeros(itr) +g1_data = zeros(itr) +g2_data = zeros(itr) + +failedStars = [] +# ---------------------------------------------------------# +fancyPrint("Analytic Profile Fit for Model Star") +@time begin + pb = tqdm(1:itr) + for i in pb + initial_guess = rand(3) #println("\t initial guess [σ, e1, e2]: ", initial_guess) + set_description(pb, "Star $i/$itr Complete") + + global iteration = i + try + global x_cg = optimize(cost, + g!, + initial_guess, + LBFGS(),#ConjugateGradient() + Optim.Options(g_tol = minAnalyticGradientModel))#Optim.Options(callback = cb) #Optim.Options(g_tol = 1e-6)) + + s_model[i] = x_cg.minimizer[1]^2 + e1_guess = x_cg.minimizer[2] + e2_guess = x_cg.minimizer[3] + + ellipticityData = sqrt((e1_guess)^2 + (e2_guess)^2) + normGdata = sqrt(1 + 0.5*( (1/ellipticityData^2) - sqrt( (4/ellipticityData^2) + (1/ellipticityData^4) ) )) + ratioData = ellipticityData/normGdata + g1_model[i] = e1_guess/ratioData + g2_model[i] = e2_guess/ratioData + + catch + println("Star $i failed") + #push!(failedStars, i) + s_model[i] = 0 + g1_model[i] = 0 + g2_model[i] = 0 + continue + end + end +end + +s_blacklist = [] +for i in 1:length(s_model) + if (s_model[i] < sLowerBound || s_model[i] > sUpperBound) #i in failedStars is optional Since Failed Stars are assigned s=0 + push!(s_blacklist, i) + end +end + +println("\n━ Blacklisted Stars: ", s_blacklist) +println("\n━ Blacklisted $(length(s_blacklist)) stars on the basis of s < $sLowerBound or s > $sUpperBound (Failed Stars Assigned 0)." ) +println("\n━ NB: These blacklisted stars are being indexed after the initial removal on the basis of signal to noise, not based off of their original location in the star catalog.") +for i in sort(s_blacklist, rev=true) + splice!(starCatalog, i) + #splice!(errVignets, i) + splice!(s_model, i) + splice!(g1_model, i) + splice!(g2_model, i) + splice!(u_coordinates, i) + splice!(v_coordinates, i) +end + +#prePixelGridFits = starCatalog +#println(size(starCatalog[1])) +failedStars = [] + +# ---------------------------------------------------------# +fancyPrint("Pixel Grid Fit") +pixelGridFits = [] + +if mode == "chisq" + println("━ χ2 Mode...\n") + @time begin + pb = tqdm(1:itr) + for i in pb + set_description(pb, "Star $i/$itr Complete") + global iteration = i + initial_guess = rand(r*c) + try + global chisq_cg = optimize(chisq_cost, + chisq_g!, + initial_guess, + LBFGS())#,#ConjugateGradient() + + #Optim.Options(g_tol = chisq_stopping_gradient))#Optim.Options(callback = cb) #Optim.Options(g_tol = 1e-6)) + catch ex + println(ex) + println("Star $i failed") + push!(failedStars, i) + push!(pixelGridFits, zeros(r,c)) + continue + end + if unity_sum + pgf_current = reshape(chisq_cg.minimizer, (r, c))./sum(reshape(chisq_cg.minimizer, (r, c))) + else + pgf_current = reshape(chisq_cg.minimizer, (r, c)) + end + push!(pixelGridFits, pgf_current) + end + end +end + + +if mode == "autoencoder" + println("━ Autoencoder Mode...\n") + println(autoencoder,"\n") + @time begin + pb = tqdm(1:length(starCatalog)) + for i in pb + set_description(pb, "Star $i/$(length(starCatalog)) Complete") + global iteration = i + + # Format some random image data + #data = nanToGaussian(starCatalog[i], s_model[i], g1_model[i], g2_model[i], r/2, c/2) + #data = reshape(data, length(data)) + data = nanToZero(reshape(starCatalog[i], length(starCatalog[i]))) + #data = Float32.(reshape(nanToZero(starCatalog[i]), r, c, 1, 1)) + + # Train the autoencoder + #data_x̂ = sample_image(autoencoder(data_x),r) + try + min_gradient = minPixelGradient + for epoch in 1:epochs + Flux.train!(loss, Flux.params(autoencoder), [(data, )], optimizer) #loss#Flux.params(autoencoder)) + + grad = Flux.gradient(Flux.params(autoencoder)) do + loss(data) + end + grad_norm = norm(grad) + if (grad_norm < min_gradient) #min_gradient + #println(epoch) + break + end + end + # Take a sample input image + #input_image = reshape(starCatalog[i], length(starCatalog[i])) + + # Pass the input image through the autoencoder to get the reconstructed image + reconstructed_image = autoencoder(data) #autoencoder(input_image) + + if unity_sum + pgf_current = reshape(reconstructed_image, (r, c))./sum(reshape(reconstructed_image, (r, c))) + else + pgf_current = reshape(reconstructed_image, (r, c)) + end + push!(pixelGridFits, pgf_current) + catch ex + println(ex) + println("Star $i failed") + push!(failedStars, i) + push!(pixelGridFits, zeros(r,c)) + continue + end + end + end +end + +if mode == "PCA" + println("━ PCA Mode...") + @time begin + pb = tqdm(1:length(starCatalog)) + for i in pb + set_description(pb, "Star $i/$(length(starCatalog)) Complete") + global iteration = i + data = nanToZero(starCatalog[i]) + try + if unity_sum + try + push!(pixelGridFits, smooth(pca_image(data,PCAterms), lanczos)./sum(smooth(pca_image(data,PCAterms)), lanczos)) + catch + #println("Smoothing failed") + push!(pixelGridFits, pca_image(data,PCAterms)./sum(pca_image(data,PCAterms))) + end + else + try + push!(pixelGridFits, smooth(pca_image(data,PCAterms), lanczos)) + catch + #println("Smoothing failed") + push!(pixelGridFits, pca_image(data,PCAterms)) + end + end + catch + println("Star $i failed") + push!(failedStars, i) + push!(pixelGridFits, zeros(r,c)) + continue + end + + end + end +end + +if mode == "smooth" + println("━ Smooth Mode...") + @time begin + pb = tqdm(1:length(starCatalog)) + for i in pb + set_description(pb, "Star $i/$(length(starCatalog)) Complete") + global iteration = i + data = nanToZero(starCatalog[i]) + try + if unity_sum + try + push!(pixelGridFits, smooth(data, lanczos)./sum(smooth(data), lanczos)) + catch + #println("Smoothing failed") + push!(pixelGridFits, data./sum(data)) + end + else + try + push!(pixelGridFits, smooth(data, lanczos)) + catch + #println("Smoothing failed") + push!(pixelGridFits, data) + end + end + catch + println("Star $i failed") + push!(failedStars, i) + push!(pixelGridFits, zeros(r,c)) + continue + end + end + end +end + + + +GC.gc() + + +println("━ failed stars:", failedStars) +# ---------------------------------------------------------# +fancyPrint("Analytic Profile Fit for Learned Star") +#Copy Star Catalog then replace it with the learned pixel grid stars +@time begin + pb = tqdm(1:length(starCatalog)) + for i in pb + initial_guess = rand(3) #println("\t initial guess [σ, e1, e2]: ", initial_guess) + set_description(pb, "Star $i/$(length(starCatalog)) Complete") + + global iteration = i + try + global y_cg = optimize(costD, + gD!, + initial_guess, + LBFGS(),#ConjugateGradient()ConjugateGradient(), + Optim.Options(g_tol = minAnalyticGradientLearned)) #Optim.Options(callback = cb) + + s_data[i] = y_cg.minimizer[1]^2 + e1_guess = y_cg.minimizer[2] + e2_guess = y_cg.minimizer[3] + + ellipticityData = sqrt((e1_guess)^2 + (e2_guess)^2) + normGdata = sqrt(1 + 0.5*( (1/ellipticityData^2) - sqrt( (4/ellipticityData^2) + (1/ellipticityData^4) ) )) + ratioData = ellipticityData/normGdata + g1_data[i] = e1_guess/ratioData + g2_data[i] = e2_guess/ratioData + + + if s_data[i] < sLowerBound || s_data[i] > sUpperBound + push!(failedStars, i) + end + catch + println("Star $i failed") + s_data[i] = 0 + g1_data[i] = 0 + g2_data[i] = 0 + continue + end + #println("\t Found A: ", A_data[i], "\t s: ", s_data[i]^2, "\t g1: ", g1_data[i], "\t g2: ", g2_data[i]) + + end +end + + +println("━ failed stars: ", unique(failedStars)) +println("\n━ Rejected $(length(unique(failedStars))) more stars for failing or having either s < $sLowerBound or s > $sUpperBound when fitting an analytic profile to an autoencoded image.") +println("\n━ NB: These failed stars are being indexed after both the screening of signal to noise and the blacklisting of s values.") +failedStars = unique(failedStars) + +for i in sort(failedStars, rev=true) + splice!(pixelGridFits, i) + splice!(s_data, i) + splice!(g1_data, i) + splice!(g2_data, i) + splice!(s_model, i) + splice!(g1_model, i) + splice!(g2_model, i) + splice!(u_coordinates, i) + splice!(v_coordinates, i) + splice!(starCatalog, i) + #splice!(errVignets, i) +end + +GC.gc() + +# ---------------------------------------------------------# +fancyPrint("Transforming (x,y) -> (u,v) | Interpolation [s, g1, g2] Across the Field of View") + +s_data = s_data[1:length(pixelGridFits)] +g1_data = g1_data[1:length(pixelGridFits)] +g2_data = g2_data[1:length(pixelGridFits)] + +s_tuples = [] +for i in 1:length(starCatalog) + push!(s_tuples, (u_coordinates[i], v_coordinates[i], s_data[i])) +end + +s_fov = optimize(interpCostS, polyG_s!, rand(10), LBFGS()) +sC = s_fov.minimizer +println("\ns(u,v) = \n$(sC[1]) u³ \n+ $(sC[2]) v³ \n+ $(sC[3]) u²v \n+ $(sC[4]) v²u \n+ $(sC[5]) u² \n+ $(sC[6]) v² \n+ $(sC[7]) uv \n+ $(sC[8]) u \n+ $(sC[9]) v \n+ $(sC[10])\n") + +s(u,v) = sC[1]*u^3 + sC[2]*v^3 + sC[3]*u^2*v + sC[4]*v^2*u + sC[5]*u^2 + sC[6]*v^2 + sC[7]*u*v + sC[8]*u + sC[9]*v + sC[10] +ds_du(u,v) = sC[1]*3*u^2 + sC[3]*2*u*v + sC[4]*v^2 + sC[5]*2*u + sC[7]*v + sC[8] +ds_dv(u,v) = sC[2]*3*v^2 + sC[3]*u^2 + sC[4]*2*u*v + sC[6]*2*v + sC[7]*u + sC[9] + +g1_tuples = [] +for i in 1:length(starCatalog) + push!(g1_tuples, (u_coordinates[i], v_coordinates[i], g1_data[i])) +end + +g1_fov = optimize(interpCostg1, polyG_g1!, rand(10), LBFGS()) +g1C = g1_fov.minimizer +println("\ng1(u,v) = \n$(g1C[1]) u³ \n+ $(g1C[2]) v³ \n+ $(g1C[3]) u²v \n+ $(g1C[4]) v²u \n+ $(g1C[5]) u² \n+ $(g1C[6]) v² \n+ $(g1C[7]) uv \n+ $(g1C[8]) u \n+ $(g1C[9]) v \n+ $(g1C[10])\n") + +g1(u,v) = g1C[1]*u^3 + g1C[2]*v^3 + g1C[3]*u^2*v + g1C[4]*v^2*u + g1C[5]*u^2 + g1C[6]*v^2 + g1C[7]*u*v + g1C[8]*u + g1C[9]*v + g1C[10] +dg1_du(u,v) = g1C[1]*3*u^2 + g1C[3]*2*u*v + g1C[4]*v^2 + g1C[5]*2*u + g1C[7]*v + g1C[8] +dg1_dv(u,v) = g1C[2]*3*v^2 + g1C[3]*u^2 + g1C[4]*2*u*v + g1C[6]*2*v + g1C[7]*u + g1C[9] + +g2_tuples = [] +for i in 1:length(starCatalog) + push!(g2_tuples, (u_coordinates[i], v_coordinates[i], g2_data[i])) +end +h_uv_data = g2_tuples + +g2_fov = optimize(interpCostg2, polyG_g2!, rand(10), LBFGS()) +g2C = g2_fov.minimizer +println("\ng2(u,v) = \n$(g2C[1]) u³ \n+ $(g2C[2]) v³ \n+ $(g2C[3]) u²v \n+ $(g2C[4]) v²u \n+ $(g2C[5]) u² \n+ $(g2C[6]) v² \n+ $(g2C[7]) uv \n+ $(g2C[8]) u \n+ $(g2C[9]) v \n+ $(g2C[10])\n") + +g2(u,v) = g2C[1]*u^3 + g2C[2]*v^3 + g2C[3]*u^2*v + g2C[4]*v^2*u + g2C[5]*u^2 + g2C[6]*v^2 + g2C[7]*u*v + g2C[8]*u + g2C[9]*v + g2C[10] +dg2_du(u,v) = g2C[1]*3*u^2 + g2C[3]*2*u*v + g2C[4]*v^2 + g2C[5]*2*u + g2C[7]*v + g2C[8] +dg2_dv(u,v) = g2C[2]*3*v^2 + g2C[3]*u^2 + g2C[4]*2*u*v + g2C[6]*2*v + g2C[7]*u + g2C[9] + +#println("\n** Adding a Progress Bar Dramatically Increases the Run Time, but note that Interpolation across the FOV is taking place! **\n") + +#PolynomialMatrix = ones(r,c, 10) + +function sample_indices(array, k) + indices = collect(1:length(array)) # Create an array of indices + return randperm(length(indices))[1:k] #sample(indices, k, replace = false) +end + +total_samples = length(pixelGridFits) +#training_ratio = 0.8 +training_samples = round(Int, training_ratio * total_samples) + +training_indices = sample_indices(pixelGridFits, training_samples) +training_stars = pixelGridFits[training_indices] +training_u_coords = u_coordinates[training_indices] +training_v_coords = v_coordinates[training_indices] + +validation_indices = setdiff(1:total_samples, training_indices) +validation_stars = pixelGridFits[validation_indices] +validation_u_coords = u_coordinates[validation_indices] +validation_v_coords = v_coordinates[validation_indices] +validation_star_catalog = starCatalog[validation_indices] + +PolynomialMatrix = ones(r,c, (degree + 1) * (degree + 2) ÷ 2 ) + +fancyPrint("Transforming (x,y) -> (u,v) | Interpolation Pixel by Pixel Across the Field of View") + +function compute_single_star_reconstructed_value(PolynomialMatrix, x, y, degree) + r, c, _ = size(PolynomialMatrix) + reconstructed_star = zeros(r, c) + for i in 1:r + for j in 1:c + p = PolynomialMatrix[i, j, :] + reconstructed_star[i, j] = objective_function(p, x, y, degree) + end + end + return reconstructed_star +end + + +function compute_mse(reconstructed_matrix, true_matrix) #err_map + return mean((reconstructed_matrix .- true_matrix) .^ 2 ) #./(err_map.^2) +end + +function worst_10_percent(errors) + n = length(errors) + star_errors = [(i, errors[i]) for i in 1:n] + sort!(star_errors, by = x->x[2], rev=true) + threshold_idx = Int(ceil(0.10 * n)) + worst_stars = star_errors[1:threshold_idx] + return [star[1] for star in worst_stars] +end + +if outlier_sigma_clipping + function worst_10_percent(errors) + n = length(errors) + mean_error = mean(errors) + std_error = std(errors) + threshold = mean_error + 3 * std_error + star_errors = [(i, errors[i]) for i in 1:n] + outliers = filter(x -> x[2] > threshold, star_errors) + return [outlier[1] for outlier in outliers] + end +end + +@time begin + global training_stars, training_u_coords, training_v_coords + for loop in 1:iterationsPolyfit + #print(length(iterationsPolyfit)) + if loop == 1 + println("$(length(training_stars)) training stars") + end + println("━ Iteration: $loop") + + @threads for i in 1:r + for j in 1:c + z_data = [star[i, j] for star in training_stars] + pC = polynomial_optimizer(degree, training_u_coords, training_v_coords, z_data) + PolynomialMatrix[i,j,:] .= pC + end + end + #println("Here") + training_errors = [] + for idx in 1:length(training_stars) + reconstructed_star = compute_single_star_reconstructed_value(PolynomialMatrix, training_u_coords[idx], training_v_coords[idx], degree) + push!(training_errors, compute_mse(reconstructed_star, training_stars[idx])) #errVignets[idx] + end + + bad_indices = worst_10_percent(training_errors) + + if loop != iterationsPolyfit + new_training_stars = [] + new_training_u_coords = [] + new_training_v_coords = [] + for i in 1:length(training_stars) + if i ∉ bad_indices + push!(new_training_stars, training_stars[i]) + push!(new_training_u_coords, training_u_coords[i]) + push!(new_training_v_coords, training_v_coords[i]) + end + end + + global training_stars = new_training_stars + println("$(length(new_training_stars)) training stars") + global training_u_coords = new_training_u_coords + global training_v_coords = new_training_v_coords + end + end +end + +#= Future Work for iterative refinement +global itr_count = 1 + +@time begin + global training_stars, training_u_coords, training_v_coords + + # Initialize a flag for the while loop + global outliers_exist = true + while outliers_exist # Keep iterating as long as there are outliers + #println("here") + current_itr = itr_count + println("━ Iteration $current_itr") + global itr_count += 1 + @threads for i in 1:r + for j in 1:c + z_data = [star[i, j] for star in training_stars] + pC = polynomial_optimizer(degree, training_u_coords, training_v_coords, z_data) + PolynomialMatrix[i,j,:] .= pC + end + end + + training_errors = [] + for idx in 1:length(training_stars) + reconstructed_star = compute_single_star_reconstructed_value(PolynomialMatrix, training_u_coords[idx], training_v_coords[idx], degree) + push!(training_errors, compute_mse(reconstructed_star, training_stars[idx])) + end + + # Calculate standard deviation of training_errors + error_std_dev = std(training_errors) + + # Find indices of outliers which have an error greater than 2 times the std deviation + bad_indices = findall(x -> x > 3 * error_std_dev, training_errors) + + if isempty(bad_indices) # If no outliers found, set the flag to false + #println("here") + global outliers_exist = false + else + new_training_stars = [] + new_training_u_coords = [] + new_training_v_coords = [] + for i in 1:length(training_stars) + if i ∉ bad_indices + push!(new_training_stars, training_stars[i]) + push!(new_training_u_coords, training_u_coords[i]) + push!(new_training_v_coords, training_v_coords[i]) + end + end + + global training_stars = new_training_stars + println("Number of training stars after iteration $current_itr: ", length(new_training_stars)) + if length(new_training_stars) < 30 + global outliers_exist = false + println("Training Stars < 30, terminating...") + end + + global training_u_coords = new_training_u_coords + global training_v_coords = new_training_v_coords + end + end +end +=# + +GC.gc() + +try + global sampled_indices = sort(sample_indices(validation_indices, 3)) +catch + global sampled_indices = rand(3) +end +#= +println("Sampled indices: ", sampled_indices) +meanRelativeError = [] +for i in 1:length(starCatalog) + RelativeError = [] + for j in 1:size(starCatalog[i], 1) + for k in 1:size(starCatalog[i], 2) + push!(RelativeError, abs.(starCatalog[i][j,k] .- pixelGridFits[i][j,k]) ./ abs.(starCatalog[i][j,k] .+ 1e-10)) + end + end + push!(meanRelativeError, mean(RelativeError)) +end +=# + +# ---------------------------------------------------------# +fancyPrint("Plotting") +@time begin + + plot_hist() + plot_err() + + starSample = rand(1:length(starCatalog)) + a = starCatalog[starSample] + b = pixelGridFits[starSample] + + a = nanToZero(a) + b = nanToZero(b) + + cmx = maximum([maximum(a), maximum(b)]) + cmn = minimum([minimum(a), minimum(b)]) + + function symlog(x, linthresh) + sign_x = sign(x) + abs_x = abs(x) + scaled = linthresh * log10(abs_x / linthresh + 1) + return sign_x * scaled + end + + if UnicodePlotsPrint + println(UnicodePlots.heatmap(symlog.(get_middle_nxn(a,75),0.0001), cmax = cmx, cmin = cmn, colormap=:inferno, title="Heatmap of star $starSample")) + println(UnicodePlots.heatmap(symlog.(get_middle_nxn(b,75),0.0001), cmax = cmx, cmin = cmn, colormap=:inferno, title="Heatmap of Pixel Grid Fit $starSample")) + println(UnicodePlots.heatmap(get_middle_nxn(a - b, 75), colormap=:inferno, title="Heatmap of Residuals")) + end + + if cairomakiePlots + testField(u, v) = Point2f(ds_du(u,v), ds_dv(u,v)) # x'(t) = -x, y'(t) = 2y + u = range(minimum(u_coordinates), stop=maximum(u_coordinates), step=0.0001) + v = range(minimum(v_coordinates), stop=maximum(v_coordinates), step=0.0001) + + s_map = [s(u,v) for u in u, v in v] + + fig1 = Figure(resolution = (1920, 1080), fontsize = 60, fonts = (;regular="CMU Serif")) + ax1 = fig1[1, 1] = CairoMakie.Axis(fig1, xlabel = L"u", ylabel = L"v",xticklabelsize = 80, yticklabelsize = 80, xlabelsize = 80, ylabelsize = 80) + fs1 = CairoMakie.heatmap!(ax1, u, v, s_map, colormap = Reverse(:plasma)) + CairoMakie.streamplot!(ax1, + testField, + u, + v, + colormap = Reverse(:plasma), + gridsize = (32, 32), + density = 0.25, + arrow_size = 10) + + CairoMakie.Colorbar(fig1[1, 2], + fs1, + label = L"s(u,v)", + width = 20, + labelsize = 100, + ticklabelsize = 40) + + CairoMakie.colgap!(fig1.layout, 5) + + save("s_uv.png", fig1) + + testField(u, v) = Point2f(dg1_du(u,v), dg1_dv(u,v)) # x'(t) = -x, y'(t) = 2y + u = range(minimum(u_coordinates), stop=maximum(u_coordinates), step=0.0001) + v = range(minimum(v_coordinates), stop=maximum(v_coordinates), step=0.0001) + + g1_map = [g1(u,v) for u in u, v in v] + + fig2 = Figure(resolution = (1920, 1080), fontsize = 60, fonts = (;regular="CMU Serif")) + ax2 = fig2[1, 1] = CairoMakie.Axis(fig2, xlabel = L"u", ylabel = L"v",xticklabelsize = 80, yticklabelsize = 80, xlabelsize = 80, ylabelsize = 80) + fs2 = CairoMakie.heatmap!(ax2, u, v, g1_map, colormap = Reverse(:plasma)) + CairoMakie.streamplot!(ax2, + testField, + u, + v, + colormap = Reverse(:plasma), + gridsize = (32, 32), + density = 0.25, + arrow_size = 10) + + CairoMakie.Colorbar(fig2[1, 2], + fs2, + label = L"g1(u,v)", + width = 20, + labelsize = 100, + ticklabelsize = 40) + + CairoMakie.colgap!(fig2.layout, 5) + + save("g1_uv.png", fig2) + + testField(u, v) = Point2f(dg2_du(u,v), dg2_dv(u,v)) # x'(t) = -x, y'(t) = 2y + u = range(minimum(u_coordinates), stop=maximum(u_coordinates), step=0.0001) + v = range(minimum(v_coordinates), stop=maximum(v_coordinates), step=0.0001) + + g2_map = [g2(u,v) for u in u, v in v] + + fig3 = Figure(resolution = (1920, 1080), fontsize = 60, fonts = (;regular="CMU Serif")) + ax3 = fig3[1, 1] = CairoMakie.Axis(fig3, xlabel = L"u", ylabel = L"v", xticklabelsize = 80, yticklabelsize = 80, xlabelsize = 80, ylabelsize = 80) + fs3 = CairoMakie.heatmap!(ax3, u, v, g2_map, colormap = Reverse(:plasma)) + CairoMakie.streamplot!(ax3, + testField, + u, + v, + colormap = Reverse(:plasma), + gridsize = (32, 32), + density = 0.25, + arrow_size = 10) + + CairoMakie.Colorbar(fig3[1, 2], + fs3, + label = L"g2(u,v)", + width = 20, + labelsize = 100, + ticklabelsize = 40) + + CairoMakie.colgap!(fig3.layout, 5) + + save("g2_uv.png", fig3) + end + + #= + scale = 1/0.29 + ks93, k0 = ks(g1_map, g2_map) + ksCosmos = get_middle_15x15(imfilter(ks93, Kernel.gaussian(scale))) + kshm = Plots.heatmap(ksCosmos, + title="Kaisser-Squires", + xlabel="u", + ylabel="v", + xlims=(0.5, size(ksCosmos, 2) + 0.5), # set the x-axis limits to include the full cells + ylims=(0.5, size(ksCosmos, 1) + 0.5), # set the y-axis limits to include the full cells + aspect_ratio=:equal, + ticks=:none, # remove the ticks + frame=:box, # draw a box around the plot + grid=:none, # remove the grid lines + size=(1920,1080)) + + Plots.savefig(kshm, joinpath("outdir","kaisserSquires.png")) + =# +end + +if UnicodePlotsPrint + println(UnicodePlots.boxplot(["s model", "s data", "g1 model", "g1 data", "g2 model", "g2 data"], + [s_model, s_data, g1_model, g1_data, g2_model, g2_data], + title="Boxplot of df.shopt")) +end + +GC.gc() + +# ---------------------------------------------------------# +fancyPrint("Saving Data to summary.shopt") +writeFitsData() + +GC.gc() + +# ---------------------------------------------------------# +fancyPrint("Done! =]") +end_time = Dates.now() +println("━ Total Time: ", (end_time - start) / Dates.Millisecond(1000), " seconds") +#println("━ Total Time: ", Dates.format(now() - start, "HH:MM:SS")) +println("For more on ShOpt.jl, see https://github.com/EdwardBerman/shopt")