edits

dataPreprocessing.jl

@@ -394,7 +394,7 @@ function cataloging(args; nm=nanMask, nz=nanToZero, snr=signal_to_noise, dout=ou
     end
   end
 
-  println("━ Sampled all vignets to $(size(catalogNew[1], 1)) x $(size(catalogNew[1], 2)) from $(r) x $(c) via over/under sampling")
+  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))

interpolate.jl

@@ -102,7 +102,7 @@ function objective_function(p, x, y, degree)
 end
 
 
-function polynomial_optimizer(degree, x_data, y_data, z_data)
+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
@@ -115,55 +115,31 @@ function polynomial_optimizer(degree, x_data, y_data, z_data)
     polynomial_terms = generate_polynomial_terms(degree)
     A = zeros(length(z_data), length(polynomial_terms))
     #replace nans with median 
-    for i in 1:length(z_data)
-        if isnan(z_data[i])
-            z_data[i] = median(filter(!isnan, z_data))
+    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
-    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
+    else 
+        for i in 1:length(z_data)
+            if isnan(z_data[i])
+                z_data[i] = median(filter(!isnan, z_data))
+            end
         end
-    end
-    coefficients = A \ z_data
-    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
-if median_constraint == true
-    function objective_function(p, x, y, degree; median_z=median_z)
-        num_coefficients = ((degree + 1) * (degree + 2) ÷ 2 ) - 1
-        value = 0
-        counter = 0
-        for i in 1:(degree + 1)  
-            for j in 1:(degree + 1)
-                if (i - 1) + (j - 1) <= degree
-                    if i == 1 && j == 1
-                        #print(size(z_data))
-                        value += median_z 
-                    else
-                        counter += 1
-                        value += p[counter] * x^(i - 1) * y^(j - 1) 
-                    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
-        return value
     end
-    function polynomial_optimizer(degree, x_data, y_data, z_data, median_z)
-        #println("Median constraint")
-        num_coefficients = ((degree + 1) * (degree + 2) ÷ 2) - 1
-        initial_guess = zeros(num_coefficients)  
-        initial_guess[1] = mean(z_data) 
-        median_z = median_z
-        function objective(p)
-            loss = sum((objective_function(p, x_val, y_val, degree, median_z=median_z ) - z_actual)^2  for ((x_val, y_val), z_actual) in zip(zip(x_data, y_data), z_data) if !isnan(z_actual))
-            loss += 1e-4*sum(p.^2) #L2 regularization
-            return loss
-        end
-        result = optimize(objective, initial_guess, autodiff=:forward, LBFGS())
-        #prepend median(z_data) to the coefficients
-        return Optim.minimizer(result)
+
+    coefficients = A \ z_data
+    if median_constraint == true
+        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

shopt.jl

@@ -505,22 +505,6 @@ 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)
@@ -547,13 +531,8 @@ end
     @threads for i in 1:r
       for j in 1:c
         z_data = [star[i, j] for star in training_stars]
-        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
+        pC = polynomial_optimizer(degree, training_u_coords, training_v_coords, z_data)
+        PolynomialMatrix[i,j,:] .= pC
       end
     end
     #println("Here")