From fe7f6da7ccf0ebef2ce6c6adc302843a01c11102 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?=
 <luiseduardoramirezmontoya@ciencias.unam.mx>
Date: Fri, 5 Jan 2024 08:25:58 -0600
Subject: [PATCH] Polynomial fit

---
 Project.toml                          |  2 +-
 src/NEOs.jl                           |  2 +-
 src/observations/observation_night.jl | 55 ++++++++++++++++++++++-----
 3 files changed, 48 insertions(+), 11 deletions(-)

diff --git a/Project.toml b/Project.toml
index 5b800100..841c0c7c 100755
--- a/Project.toml
+++ b/Project.toml
@@ -12,7 +12,6 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
-GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 HORIZONS = "5a3ac768-beb4-554a-9c98-3342fe3377f5"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
@@ -21,6 +20,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
diff --git a/src/NEOs.jl b/src/NEOs.jl
index f81e23e0..a48fb4ae 100755
--- a/src/NEOs.jl
+++ b/src/NEOs.jl
@@ -23,7 +23,7 @@ using PlanetaryEphemeris: daysec, yr, TaylorInterpolant, auday2kmsec, su, ea, au
       meanmotion, meananomaly, selecteph
 using Healpix: ang2pixRing, Resolution
 using StatsBase: mean, std
-using GLM: @formula, lm, coef
+using LsqFit: curve_fit
 using Roots: find_zeros
 using Clustering: kmeans
 
diff --git a/src/observations/observation_night.jl b/src/observations/observation_night.jl
index c5a9569f..b8a2fa76 100644
--- a/src/observations/observation_night.jl
+++ b/src/observations/observation_night.jl
@@ -59,6 +59,45 @@ end
 # Order in ObservationNight is given by date
 isless(a::ObservationNight{T}, b::ObservationNight{T}) where {T <: AbstractFloat} = a.date < b.date
 
+# Evaluate a polynomial with coefficients p in every element of x
+polymodel(x, p) = map(y -> evalpoly(y, p), x)
+
+# Normalized mean square residual for polynomial fit
+polyerror(x) = sum(x .^ 2) / length(x)
+
+@doc raw"""
+    polyfit(x::Vector{T}, y::Vector{T}; tol::T = 1e-4) where {T <: AbstractFloat}
+
+Fit a polynomial to points `(x, y)`. The order of the polynomial is increased
+until `polyerror` is less than `tol`.
+"""
+function polyfit(x::Vector{T}, y::Vector{T}; tol::T = 1e-4) where {T <: AbstractFloat}
+    # Avoid odd and high orders (to have a defined concavity and avoid overfit)
+    for order in [1, 2, 4, 6]
+        # Initial guess for coefficients
+        coeffs = ones(T, order+1)
+        # Polynomial fit
+        fit = curve_fit(polymodel, x, y, coeffs)
+        # Convergence condition
+        if polyerror(fit.resid) < tol || order == 6
+            return fit.param
+        end
+    end
+end
+
+@doc raw"""
+    diffcoeffs(x::Vector{T}) where {T <: AbstractFloat}
+
+Return the coefficients of the derivative of a polynomial with coefficients `x`.
+"""
+function diffcoeffs(x::Vector{T}) where {T <: AbstractFloat}
+    y = Vector{T}(undef, length(x)-1)
+    for i in eachindex(y)
+        y[i] = i * x[i+1]
+    end
+    return y
+end
+
 # Outer constructor
 function ObservationNight(radec::Vector{RadecMPC{T}}, df::AbstractDataFrame) where {T <: AbstractFloat}
     # Defining quantities of an ObservationNight
@@ -98,17 +137,15 @@ function ObservationNight(radec::Vector{RadecMPC{T}}, df::AbstractDataFrame) whe
         df.α = map(x -> x < π ? x + 2π : x, df.α)
     end
 
-    # Linear regression
-    α_p = lm(@formula(α ~ t_rel), df)
-    α_coef = coef(α_p)
-    δ_p = lm(@formula(δ ~ t_rel), df)
-    δ_coef = coef(δ_p)
+    # Polynomial regression
+    α_coef = polyfit(df.t_rel, df.α)
+    δ_coef = polyfit(df.t_rel, df.δ)
 
     # Evaluate polynomials at mean date 
-    α = mod2pi(α_coef[1] + α_coef[2] * t_mean)
-    δ = δ_coef[1] + δ_coef[2] * t_mean
-    v_α = α_coef[2]
-    v_δ = δ_coef[2]
+    α = mod2pi(polymodel(t_mean, α_coef))
+    δ = polymodel(t_mean, δ_coef)
+    v_α = polymodel(t_mean, diffcoeffs(α_coef))
+    v_δ = polymodel(t_mean, diffcoeffs(δ_coef))
 
     # Mean apparent magnitude
     mag = mean(filter(!isnan, df.mag))