Constraint Examples¶
Example Code¶
Representative script:
"""IMSL CONFT example 1: Monotone constrained B-spline fit.
Reproduces Example 1 from the IMSL CONFT documentation:
- Function: f(x) = 0.5*x + sin(0.5*x)
- Data: 15 points with random noise in [0, 10]
- Unconstrained fit (BSLSQ) vs Monotone constrained fit (CONFT)
- Constraint: f'(x) >= 0 at 15 points
Outputs:
- Average error comparison (BSLSQ vs CONFT)
- SVG plot: tests/test_output/demo_imsl_conft_monotone.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Any, Dict
import sys
import matplotlib.pyplot as plt
import numpy as np
# Allow running this file directly without requiring external PYTHONPATH setup.
REPO_ROOT = Path(__file__).resolve().parent.parent
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
from interpolation import (
BSplineRepresentation,
bspline_evaluate,
least_squares_bspline_constrained,
least_squares_bspline_fixed_knots,
)
from interpolation.core import _bspline_basis # private internal helper
def run_demo_imsl_conft_monotone() -> Dict[str, Any]:
"""Run CONFT monotone fit demo and save comparison plot."""
# Seed for reproducibility (aligns roughly with IMSL RNSET if possible,
# but we use standard numpy seed)
np.random.seed(234579)
korder = 4
ncoef = 8
ndata = 15
grdsiz = 10.0
def f_true(x):
return 0.5 * x + np.sin(0.5 * x)
# Generate data with noise
x_data = np.linspace(0.0, grdsiz, ndata)
# IMSL: FDATA(I) = F1(XDATA(I)) + (FDATA(I)-.5) where FDATA is RNUNF
noise = np.random.uniform(0.0, 1.0, ndata) - 0.5
f_data = f_true(x_data) + noise
# Knots: IMSL computes interior knots based on NCOEF - KORDER + 2 breakpoints
# For NCOEF=8, KORDER=4, this is 8-4+2=6 breakpoints (0 to 10)
# result: [0, 2, 4, 6, 8, 10]
n_breakpoints = ncoef - korder + 2
inner_breaks = np.linspace(0.0, grdsiz, n_breakpoints)
knots = np.concatenate(
([inner_breaks[0]] * (korder - 1), inner_breaks, [inner_breaks[-1]] * (korder - 1))
)
# 1. Unconstrained B-spline LSQ fit
# Note: least_squares_bspline_fixed_knots should return BSplineRepresentation
bs_lsq = least_squares_bspline_fixed_knots(x_data, f_data, knots, korder)
# 2. Constrained fit (Monotone: f' >= 0)
nxval = 15
x_val = np.linspace(0.0, grdsiz, nxval)
n_coeff = len(knots) - korder
# In this implementation, we must provide the constraint matrix and RHS directly
# for least_squares_bspline_constrained.
# Derivative order 1 constraints at x_val
def bspline_deriv_basis(i, k, t, x):
eps = 1e-6
return (_bspline_basis(i, k, t, x + eps) - _bspline_basis(i, k, t, x - eps)) / (2 * eps)
c_mat = np.zeros((nxval, n_coeff))
for i, xi in enumerate(x_val):
for j in range(n_coeff):
c_mat[i, j] = bspline_deriv_basis(j, korder, knots, xi)
# This CONFT analogue enforces equalities only, so we approximate monotonicity
# by steering derivatives toward the nonnegative part of the unconstrained fit.
# A moderate penalty avoids over-constraining the LSQ problem.
slope_eps = 1e-4
c_rhs = np.asarray(
[
max(
0.0,
(
bspline_evaluate(bs_lsq, float(xi + slope_eps))
- bspline_evaluate(bs_lsq, float(xi - slope_eps))
)
/ (2.0 * slope_eps),
)
for xi in x_val
],
dtype=float,
)
bs_conft = least_squares_bspline_constrained(
x_data, f_data, knots, korder, constraints_matrix=c_mat, constraints_rhs=c_rhs, penalty=10.0
)
# Evaluation for error comparison
x_eval = np.linspace(0.0, grdsiz, 100)
y_true = f_true(x_eval)
# Handle bspline_evaluate returning array or scalar
y_lsq = np.asarray([bspline_evaluate(bs_lsq, xi) for xi in x_eval])
y_conft = np.asarray([bspline_evaluate(bs_conft, xi) for xi in x_eval])
avg_err_lsq = np.mean(np.abs(y_true - y_lsq))
avg_err_conft = np.mean(np.abs(y_true - y_conft))
print(f"\nIMSL CONFT Example 1 Parity")
print(f"Average error with unconstrained fit: {avg_err_lsq:.5f}")
print(f"Average error with monotone fit: {avg_err_conft:.5f}")
# Plotting
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
output_path = output_dir / "demo_imsl_conft_monotone.svg"
fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(x_eval, y_true, "k--", alpha=0.6, label="True f(x)")
ax.plot(x_eval, y_lsq, color="red", label="Unconstrained (BSLSQ)")
ax.plot(x_eval, y_conft, color="blue", label="Constrained (CONFT)")
ax.scatter(x_data, f_data, color="black", s=20, label="Noisy Data")
ax.set_title("IMSL CONFT Example 1: Monotone B-Spline Fit")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig(output_path)
plt.close(fig)
return {
"avg_err_lsq": avg_err_lsq,
"avg_err_conft": avg_err_conft,
"plot_path": str(output_path),
}
if __name__ == "__main__":
run_demo_imsl_conft_monotone()
Input (Console)¶
Run the constraint scripts from the package root:
python examples/example_imsl_conft_1.py
python examples/example_imsl_conft_2.py
Plot Output¶
Generated SVG plots:
Output Console¶
Summary console output:
script input setup representative numeric outputs
example_imsl_conft_1.py monotone constraints on [0,10] average unconstrained error = 0.1100303928219149; average constrained error = 0.1147826851898706
example_imsl_conft_2.py periodic and flat constraints on [-7,7] average unconstrained error = 0.01893213473141974; average constrained error = 0.01390170332275927
Detailed integration result tables:
CONFT integration output
example avg_error_unconstrained avg_error_constrained
CONFT-1 monotone 0.1100303928219149 0.1147826851898706
CONFT-2 periodic-flat 0.01893213473141974 0.01390170332275927