developments

dataPreprocessing.jl

@@ -34,7 +34,7 @@ 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"] 
 #=
 Log these config choices
 =#
@@ -48,6 +48,7 @@ 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)
@@ -312,10 +313,12 @@ function cataloging(args; nm=nanMask, nz=nanToZero, snr=signal_to_noise, dout=ou
   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"
@@ -400,9 +403,9 @@ function cataloging(args; nm=nanMask, nz=nanToZero, snr=signal_to_noise, dout=ou
   end
 
   if mode == "chisq"
-    return catalogNew, r, c, length(catalogNew), u_coords, v_coords, outlier_indices, errVignets
+    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
+    return catalogNew, r, c, length(catalogNew), u_coords, v_coords, outlier_indices, median_img
   end
 end
 

shopt.jl

@@ -55,6 +55,7 @@ fancyPrint("Handling Imports")
   using Dates
   using MultivariateStats
   using Base.Threads
+  using MultivariatePolynomials
   #using Interpolations
 end
 println("━ Start Time ", Dates.now())
@@ -86,9 +87,9 @@ fancyPrint("Processing Data for Fit")
 @time begin
 
   if mode == "chisq"
-    starCatalog, r, c, itr, u_coordinates, v_coordinates, outlier_indices, errVignets = cataloging(ARGS)
+    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 = cataloging(ARGS)
+    starCatalog, r, c, itr, u_coordinates, v_coordinates, outlier_indices, median_img = cataloging(ARGS)
   end
   #starCatalog = starCatalog
   #errVignets = errVignets
@@ -182,6 +183,8 @@ for i in sort(s_blacklist, rev=true)
   splice!(v_coordinates, i)
 end
 
+#prePixelGridFits = starCatalog
+#println(size(starCatalog[1]))
 failedStars = []
 
 # ---------------------------------------------------------#
@@ -310,6 +313,40 @@ if mode == "PCA"
   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()
@@ -468,6 +505,23 @@ function compute_single_star_reconstructed_value(PolynomialMatrix, x, y, degree)
     return reconstructed_star
 end
 
+if median_constraint == true
+    function compute_single_star_reconstructed_value(PolynomialMatrix, x, y, degree, ppgf=median_img)
+        #println("Found pre Pixel Grid Fits")
+        r, c, _ = size(PolynomialMatrix)
+        reconstructed_star = zeros(r, c)
+        for i in 1:r
+            for j in 1:c
+                p = PolynomialMatrix[i, j, :]
+                z_data = median_img[i, j]
+                reconstructed_star[i, j] = objective_function(p, x, y, degree, median_z=z_data)
+            end
+        end
+        return reconstructed_star
+    end
+end
+
+
 function compute_mse(reconstructed_matrix, true_matrix) #err_map
   return mean((reconstructed_matrix .- true_matrix) .^ 2 ) #./(err_map.^2)
 end
@@ -493,26 +547,20 @@ end
     @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
-    #=
-    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
+        if median_constraint == true
+            pC = polynomial_optimizer(degree, training_u_coords, training_v_coords, z_data, median_img[i,j])
+            PolynomialMatrix[i,j,:] = vcat(median_img[i,j], pC)
+        else
+            pC = polynomial_optimizer(degree, training_u_coords, training_v_coords, z_data)
+            PolynomialMatrix[i,j,:] .= pC
+        end
       end
     end
-    =#
-
-    #training_errors = Threads.@spawn [compute_mse(compute_single_star_reconstructed_value(PolynomialMatrix, training_u_coords[idx], training_v_coords[idx], degree), training_stars[idx]) for idx in 1:length(training_stars)]
-
+    #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]
+        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)