Fitting Curves from Partial Data: Historical Morphing & Boundary Smoothing
TL;DR: Estimating full-season trajectory profiles (such as crop growth or energy demand) when only the initial 25% of data is observed causes runaway polynomial extrapolation or discontinuous historical splicing. In this post, we develop a 3-tier mathematical framework: 4-parameter affine profile morphing to adapt historical archetypes, a \(C^1\)-continuous Hermite blending bridge to eliminate boundary jump discontinuities, and a constrained quadratic program (QP) to strictly enforce cumulative macro targets with minimal curvature roughness. Try the interactive visual simulator and collapsed production SciPy code below.
"When predicting the future from the past, the true art lies not in following the historical average blindly, but in stretching, shifting, and smoothly welding the template to match today's unfolding reality."
1. The Core Problem: The 25% Incomplete Curve Dilemma
In agricultural analytics, remote sensing, and industrial forecasting, practitioners frequently encounter the partial curve completion problem:
- Known Reality (\(0 \le t \le t_{\text{cut}}\)): We possess high-precision measurements for only the first quarter of the season (e.g., Day of Year \(1\) to \(35\) for a \(140\)-day corn growing cycle).
- Historical Prior (\(0 \le t \le T\)): We have a multi-year historical average growth curve \(f_{\text{hist}}(t)\) describing standard sigmoidal biomass accumulation.
- External Constraint (Optional): We may have an independent macro-level prediction (e.g., satellite yield forecast or seasonal nitrogen budget) dictating the target total biomass \(Y_{\text{total}} = \int_0^T y(t) \, dt\) or final maturity value \(y(T)\).
flowchart LR
subgraph Known["Known Region (0 to 25%)"]
Obs["Noisy Observations y_obs(t)<br/>Real-time In-Season Data"]
end
subgraph Interface["The Boundary (t = 25%)"]
Jump["Discontinuity Hazard!<br/>Step-jump & Derivative Mismatch"]
end
subgraph Unknown["Forecast Region (25% to 100%)"]
Hist["Historical Baseline f_hist(t)<br/>Provides Phenological Shape"]
Target["Macro Target Constraint Y_total<br/>Enforces Physical Mass Balance"]
end
Obs --> Interface
Hist --> Interface
Target --> Interface
Interface --> SmoothCurve["Seamless C^1 Continuous Full-Season Trajectory"]
style Known fill:#065f46,stroke:#047857,stroke-width:2px,color:#fff
style Interface fill:#b91c1c,stroke:#ef4444,stroke-width:2px,color:#fff
style Unknown fill:#1e3a8a,stroke:#1e40af,stroke-width:2px,color:#fff
style SmoothCurve fill:#0284c7,stroke:#0369a1,stroke-width:2px,color:#fffWhy Standard Methods Fail
- Direct Polynomial / Spline Extrapolation: Fitting a high-degree polynomial or unconstrained spline to the first \(25\%\) of points causes explosive Runge oscillations as \(t \to T\).
- Raw Historical Splicing: Sticking the unadjusted historical curve onto the last observed point creates an abrupt angle or vertical cliff (discontinuity in \(y(t)\) or \(y'(t)\)).
- Pure Proportional Rescaling: Multiplying the entire historical curve by \(\alpha = y_{\text{obs}}(t_{\text{cut}}) / f_{\text{hist}}(t_{\text{cut}})\) ignores timing delays (phenological shifts) and cannot satisfy cumulative target constraints.
2. Approach 1: 4-Parameter Affine Profile Morphing
The simplest, most robust baseline is to assume that current season dynamics represent an affine transformation (translation and linear stretching in time and amplitude) of the historical archetype \(f_{\text{hist}}(t)\)[^1]:
where the parameter vector \(\mathbf{\theta} = [s_x, s_y, \Delta t, \Delta y]^T\) represents: - \(s_x\) (Time Dilation / Phenology Speed): \(s_x > 1\) represents accelerated growing degree days (GDD); \(s_x < 1\) represents delayed crop development. - \(\Delta t\) (Time Shift): Accounts for late or early planting dates. - \(s_y\) (Amplitude Scaling): Reflects higher or lower overall biomass vigor. - \(\Delta y\) (Baseline Offset): Initial background offset (usually bounded near \(0\)).
Affine Parameter Morphing
y ↑ Historical f_hist(t)
│ ╭───────────
│ ╭─────╯
│ ╭─────╯ (Scaled by s_y, Shifted by Δt)
│ ╭─────╯ • • • • Fitted Forecast y(t)
│ ╭────────╯ •
│ ╭──╯ ▲ Known 25% •
│ ╭──╯ │ Obs Points
└──────────────────────────────────────────────────► t
Optimization Formulation
We estimate \(\mathbf{\theta}^*\) by solving a bounded nonlinear least-squares problem over the observed window \(t \in [0, t_{\text{cut}}]\):
where \(w_i = \exp\left( \frac{t_i - t_{\text{cut}}}{\sigma} \right)\) optionally gives greater weight to recent points, and \(\mathbf{\theta}_{\text{prior}} = [1, 1, 0, 0]^T\) prevents unphysical distortions.
View Python Implementation: Affine Profile Morphing
import numpy as np
from scipy.interpolate import CubicSpline
from scipy.optimize import minimize
def fit_affine_curve(t_obs, y_obs, t_hist, y_hist, t_full):
"""
Fits 4-parameter affine transformation: y(t) = s_y * f_hist(s_x * (t - dt)) + dy
"""
# Create continuous historical template function
f_hist = CubicSpline(t_hist, y_hist, extrapolate=True)
def objective(params):
sx, sy, dt, dy = params
# Transformed time coordinates
t_warped = sx * (t_obs - dt)
y_pred = sy * f_hist(t_warped) + dy
# Weighted least squares with prior regularization
residuals = y_obs - y_pred
prior_penalty = 10.0 * ((sx - 1.0)**2 + (dt / 10.0)**2 + (dy)**2)
return np.sum(residuals**2) + prior_penalty
# Initial guess [sx, sy, dt, dy]
p0 = [1.0, np.mean(y_obs[-3:]) / (f_hist(t_obs[-1]) + 1e-6), 0.0, 0.0]
bounds = [(0.7, 1.4), (0.2, 3.0), (-20.0, 20.0), (-5.0, 5.0)]
res = minimize(objective, p0, bounds=bounds, method='L-BFGS-B')
sx_opt, sy_opt, dt_opt, dy_opt = res.x
# Project across full season
y_forecast = sy_opt * f_hist(sx_opt * (t_full - dt_opt)) + dy_opt
return y_forecast, res.x
3. Approach 2: Boundary-Continuous Hermite Blending Bridge
Even with an optimal affine fit, the fitted curve \(\hat{y}(t)\) evaluated at \(t_{\text{cut}}\) rarely matches the exact final observed point \(y_{\text{obs}}(t_{\text{cut}})\) and its instantaneous slope \(y'_{\text{obs}}(t_{\text{cut}})\).
To achieve \(C^1\) continuity (smooth value and slope transition with zero jump), we introduce a Hermite Blending Bridge[^2]:
where the residual correction \(\delta(t)\) decays smoothly from the boundary discrepancy \(\delta_0 = y_{\text{obs}}(t_{\text{cut}}) - \hat{y}(t_{\text{cut}})\) to zero over a relaxation timescale \(\tau\):
where \(\delta'_0 = y'_{\text{obs}}(t_{\text{cut}}) - \hat{y}'(t_{\text{cut}})\) matches the instantaneous first derivative.
Discontinuity Without Bridge Smooth Hermite Blending Bridge
y ↑ y ↑
│ / (Fitted Forecast) │ / (Smooth Transition)
│ / │ ╭╯
│ / │ ╭╯
│ ● ◄── Discontinuous Jump Δy │ ● ◄── Zero Jump (C¹ Continuous)
│ / │ /
│ / (Known 25%) │ / (Known 25%)
└──┴───────────────────────► t └──┴───────────────────────► t
t_cut t_cut
View Python Implementation: Hermite Blending Bridge
def apply_hermite_bridge(t_full, y_forecast, t_cut, y_cut, dy_cut, tau=15.0):
"""
Applies C1-continuous residual decay to eliminate jump at t_cut.
"""
idx_cut = np.searchsorted(t_full, t_cut)
y_model_cut = y_forecast[idx_cut]
# Numerical derivative of forecast at cut
dy_model_cut = (y_forecast[idx_cut + 1] - y_forecast[idx_cut - 1]) / (t_full[idx_cut + 1] - t_full[idx_cut - 1])
delta_0 = y_cut - y_model_cut
delta_prime_0 = dy_cut - dy_model_cut
y_smooth = np.copy(y_forecast)
t_after = t_full[idx_cut:] - t_cut
# Hermite gaussian-damped bridge kernel
bridge = (delta_0 + delta_prime_0 * t_after) * np.exp(-(t_after**2) / (2 * tau**2))
y_smooth[idx_cut:] += bridge
return y_smooth
4. Approach 3: Constrained Target Updating (Quadratic Programming)
In advanced applications, we are given an external cumulative total target \(Y_{\text{total}}\) (e.g., predicted end-of-season total dry biomass of \(18.5\text{ t ha}^{-1}\) or total cumulative water volume):
The remaining integral to be allocated across the unknown future \(t > t_{\text{cut}}\) is:
The Optimal Distribution Problem
How do we distribute \(Y_{\text{remaining}}\) across \(t \in (t_{\text{cut}}, T]\) such that: 1. The curve mirrors the historical distribution shape \(p(t) \propto f_{\text{hist}}(t)\). 2. The total sum strictly equals \(Y_{\text{remaining}}\). 3. The curve joins seamlessly at \(t_{\text{cut}}\) with minimum curvature energy (\(\int [y''(t)]^2 dt\))[^3].
We formulate this as a Constrained Quadratic Program (QP):
flowchart TD
TotalTarget["External Macro Target: Y_total"] --> RemCalc["Calculate Remaining Mass: Y_rem = Y_total - ∑y_known"]
HistShape["Historical Template Shape: f_hist(t)"] --> QP["Quadratic Program Formulation"]
BoundaryC1["Boundary C^1 Conditions: y(t_cut), y'(t_cut)"] --> QP
QP --> MatrixSolve["Solve KKT Linear System: [H A^T; A 0] [y; λ] = [d; b]"]
MatrixSolve --> OptimalCurve["Globally Optimal Future Trajectory:<br/>- Perfect Total Target Match<br/>- Zero Boundary Discontinuity<br/>- Maximum Geometric Smoothness"]
style TotalTarget fill:#f59e0b,stroke:#d97706,stroke-width:2px,color:#fff
style QP fill:#1e293b,stroke:#3b82f6,stroke-width:2px,color:#fff
style MatrixSolve fill:#0284c7,stroke:#0369a1,stroke-width:2px,color:#fff
style OptimalCurve fill:#059669,stroke:#10b981,stroke-width:2px,color:#fffView Python Implementation: Target-Constrained QP Solver
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
def fit_target_constrained_curve(t_future, f_hist_future, y_cut, dy_cut, target_remaining, dt, alpha=0.1):
"""
Solves QP for smooth trajectory matching historical shape and exact target integral.
"""
M = len(t_future)
# Second-order difference matrix (curvature)
D2 = sp.diags([1, -2, 1], [0, 1, 2], shape=(M - 2, M)).tocsc()
H = 2 * (D2.T @ D2) + 2 * alpha * sp.eye(M).tocsc()
# Scale historical template to approximate remaining target
gamma = target_remaining / (np.sum(f_hist_future) * dt)
d = 2 * alpha * gamma * f_hist_future
# Equality constraints:
# 1. Sum constraint: dt * sum(y) = target_remaining
# 2. Boundary value: y[0] = y_cut
# 3. Boundary slope: y[1] - y[0] = dy_cut * dt
A_sum = np.ones((1, M)) * dt
A_b0 = np.zeros((1, M)); A_b0[0, 0] = 1.0
A_b1 = np.zeros((1, M)); A_b1[0, 0] = -1.0; A_b1[0, 1] = 1.0
A = np.vstack([A_sum, A_b0, A_b1])
b = np.array([target_remaining, y_cut, dy_cut * dt])
# Build KKT saddle-point system
KKT_top = sp.hstack([H, A.T])
KKT_bot = sp.hstack([A, sp.csc_matrix((3, 3))])
KKT = sp.vstack([KKT_top, KKT_bot]).tocsc()
rhs = np.concatenate([d, b])
sol = spsolve(KKT, rhs)
y_optimal = sol[:M]
return y_optimal
6. Production-Ready Python Pipeline: PartialCurveFitter
Here is a unified, object-oriented Python implementation using numpy and scipy that encapsulates all three methods into a single class:
View Production Python Module: PartialGrowthCurveFitter
import numpy as np
import scipy.sparse as sp
from scipy.interpolate import CubicSpline
from scipy.optimize import minimize
from scipy.sparse.linalg import spsolve
class PartialGrowthCurveFitter:
"""
Hierarchical framework for completing partial growth time series.
Methods:
1. 'affine': 4-parameter least-squares shape morphing.
2. 'bridge': Affine morphing + C1 Hermite boundary smoothing.
3. 'constrained_target': Quadratic program for exact target total allocation.
"""
def __init__(self, t_hist: np.ndarray, y_hist: np.ndarray):
self.t_hist = np.asarray(t_hist, dtype=float)
self.y_hist = np.asarray(y_hist, dtype=float)
self.spline_hist = CubicSpline(self.t_hist, self.y_hist, extrapolate=True)
def fit_affine(self, t_obs: np.ndarray, y_obs: np.ndarray, t_full: np.ndarray):
"""Tier 1: Affine Least-Squares Morphing"""
t_cut = t_obs[-1]
def loss(params):
sx, sy, dt, dy = params
pred = sy * self.spline_hist(sx * (t_obs - dt)) + dy
# Recent points receive exponential weighting
weights = np.exp((t_obs - t_cut) / 10.0)
res = weights * (y_obs - pred)
reg = 5.0 * ((sx - 1.0)**2 + (dt / 10.0)**2 + dy**2)
return np.sum(res**2) + reg
p0 = [1.0, np.mean(y_obs[-3:]) / (self.spline_hist(t_cut) + 1e-5), 0.0, 0.0]
bounds = [(0.6, 1.6), (0.2, 3.0), (-25.0, 25.0), (-5.0, 5.0)]
res = minimize(loss, p0, bounds=bounds, method='L-BFGS-B')
sx, sy, dt, dy = res.x
y_full = sy * self.spline_hist(sx * (t_full - dt)) + dy
return y_full, res.x
def fit_with_bridge(self, t_obs: np.ndarray, y_obs: np.ndarray, t_full: np.ndarray, tau: float = 15.0):
"""Tier 2: Affine Morphing + C1 Hermite Boundary Bridge"""
y_affine, params = self.fit_affine(t_obs, y_obs, t_full)
t_cut = t_obs[-1]
y_cut = y_obs[-1]
# Estimate observed boundary slope using last 3 points
dy_obs_cut = np.polyfit(t_obs[-3:], y_obs[-3:], deg=1)[0]
idx_cut = np.searchsorted(t_full, t_cut)
y_model_cut = y_affine[idx_cut]
dy_model_cut = (y_affine[idx_cut + 1] - y_affine[idx_cut - 1]) / (t_full[idx_cut + 1] - t_full[idx_cut - 1])
delta_0 = y_cut - y_model_cut
delta_p0 = dy_obs_cut - dy_model_cut
y_bridged = np.copy(y_affine)
dt = t_full[idx_cut:] - t_cut
bridge = (delta_0 + delta_p0 * dt) * np.exp(-(dt**2) / (2.0 * tau**2))
y_bridged[idx_cut:] += bridge
return y_bridged, params
def fit_constrained_target(self, t_obs: np.ndarray, y_obs: np.ndarray, t_full: np.ndarray,
target_total: float, alpha: float = 0.05):
"""Tier 3: Constrained Quadratic Program for Exact Target Total Integral"""
t_cut = t_obs[-1]
y_cut = y_obs[-1]
dy_obs_cut = np.polyfit(t_obs[-3:], y_obs[-3:], deg=1)[0]
idx_cut = np.searchsorted(t_full, t_cut)
dt_step = t_full[1] - t_full[0]
# Known integral up to t_cut
known_integral = np.trapz(y_obs, t_obs)
target_remaining = max(1e-3, target_total - known_integral)
t_future = t_full[idx_cut:]
M = len(t_future)
f_future = self.spline_hist(t_future)
# 1. Curvature penalty matrix D2
D2 = sp.diags([1.0, -2.0, 1.0], [0, 1, 2], shape=(M - 2, M)).tocsc()
H = 2.0 * (D2.T @ D2) + 2.0 * alpha * sp.eye(M).tocsc()
# 2. Linear gradient term
gamma = target_remaining / (np.sum(f_future) * dt_step)
d = 2.0 * alpha * gamma * f_future
# 3. Constraints: [Sum, Value, Slope]
A_sum = np.ones((1, M)) * dt_step
A_b0 = np.zeros((1, M)); A_b0[0, 0] = 1.0
A_b1 = np.zeros((1, M)); A_b1[0, 0] = -1.0; A_b1[0, 1] = 1.0
A = np.vstack([A_sum, A_b0, A_b1])
b = np.array([target_remaining, y_cut, dy_obs_cut * dt_step])
# 4. KKT Linear System Solve
KKT = sp.vstack([
sp.hstack([H, A.T]),
sp.hstack([A, sp.csc_matrix((3, 3))])
]).tocsc()
rhs = np.concatenate([d, b])
sol = spsolve(KKT, rhs)
y_future = sol[:M]
y_full = np.zeros_like(t_full)
y_full[:idx_cut] = np.interp(t_full[:idx_cut], t_obs, y_obs)
y_full[idx_cut:] = y_future
return y_full
7. Interactive Visual Playground
Experiment with all three algorithms below! Drag the controls to adjust in-season crop vigor, simulate late planting shifts, set an end-of-season target total, and observe the live curve fitting algorithms in action:
🌾 Partial Curve Morphing & Target Fitting Sandbox
7. Practical Engineering Guidelines & Best Practices
- Avoid Over-Fitting the 25% Window: Do not use high-degree polynomials (e.g. degree \(\ge 3\)) on early points. Crop growth is inherently sigmoidal (logistic / Gompertz / Richards curves); fitting a polynomial to the initial exponential phase will inevitably explode to infinity.
- Prior Regularization is Mandatory: Always regularize parameter deviations (\(\|\mathbf{\theta} - \mathbf{\theta}_{\text{prior}}\|^2\)) in nonlinear optimization. In early spring, 5 consecutive warm days can artificially double \(s_y\) if unconstrained.
- Decay Window (\(\tau\)) Selection: Set the Hermite bridge decay timescale \(\tau\) proportional to approximately \(10–20\%\) of the remaining season length (\(15–25\text{ days}\)). A timescale that is too short produces an abrupt kink, while a timescale that is too long overrides the historical biological inflection point.
Data and Code Availability
- Python Reference Implementations: Complete Jupyter Notebooks and reusable SciPy scripts are hosted on GitHub:
kamingfung/partial-curve-fitting. - Interactive Sandbox: The live browser visualizer above is implemented in vanilla HTML5 Canvas and runs entirely client-side without external dependencies.
References
[^1]: Yin, X., Goudriaan, J., Lantinga, E. A., Vos, J. & Spiertz, H. J. A flexible sigmoid function of determinate growth. Ann. Bot. 91, 361–371 (2003). https://doi.org/10.1093/aob/mcg029 [^2]: Farin, G. Curves and Surfaces for CAGD: A Practical Guide. (Morgan Kaufmann, San Francisco, 2002). https://doi.org/10.1016/B978-1-55860-737-8.X5000-5 [^3]: Boyd, S. & Vandenberghe, L. Convex Optimization. (Cambridge University Press, Cambridge, 2004). https://doi.org/10.1017/CBO9780511804441