sepehr/Entropyk
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

Complete Story 0.2 phantom-gradient regularization standard.

Browse Source 
Add analytic derivatives and C1 tests for smooth_max/min/clamp, document epsilon/width guidance, and declare flow_regularization as the HX reference pattern without changing component physics.

Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
sepehr
2026-07-18 16:33:48 +02:00
parent f88cd7f7d8
commit 88620790d6
7 changed files with 889 additions and 290 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
crates/components/src/heat_exchanger/flow_regularization.rs
View File

@@ -1,5 +1,9 @@
//! C¹/C∞ flow regularization for zero-flow-safe heat-exchanger residuals.
//!
//! **Reference specialized pattern** (Story 0.2): keep these HX helpers here;
//! the canonical general toolkit is [`entropyk_core::smoothing`]. See
//! `docs/components/phantom-gradient-regularization.md`.
//!
//! Zero mass flow is a **valid** operating state (circuit staging, isolation,
//! Newton trial steps). Hard branches `if |m| < ε { Q = 0 }` create Jacobian
//! discontinuities and can leave energy residuals under-ranked.

220
crates/core/src/smoothing.rs
View File

@@ -24,6 +24,28 @@
//! - [`smooth_abs`] — C^∞ approximation of `|x|`.
//! - [`smooth_clamp`] — C¹ clamp into `[lo, hi]`.
//!
//! Analytic first derivatives are exported for every Newton-facing primitive
//! ([`smooth_max_derivative`], [`smooth_min_derivative`], [`smooth_abs_derivative`],
//! [`smooth_clamp_derivative`]) so component Jacobians stay exact.
//!
//! # ε / width selection (physical scale of the blend)
//!
//! | Primitive | Parameter | Meaning | Rule of thumb |
//! |---|---|---|---|
//! | [`smooth_abs`] | `eps` | value at `x = 0`; blend half-width | ~1% of characteristic `\|x\|`, or smaller if the kink is far from the operating region |
//! | [`smooth_max`] / [`smooth_min`] | `k` | overshoot at equality = `k / 2` | `k` much smaller than `\|a - b\|` in the region you care about; never so large it shifts equilibrium |
//! | [`smooth_clamp`] | `width` | transition length at each bound | `width <= (hi - lo) / 2`; keep the band physically small (e.g. opening `1e-3`..`1e-2` on `[0, 1]`) |
//!
//! Heat-exchanger mass-flow activity (`flow_activity`, `m_eps`) lives in
//! `entropyk-components` as a specialized reference pattern — see
//! `docs/components/phantom-gradient-regularization.md`.
//!
//! # Degenerate parameters (never panic)
//!
//! - `k <= 0` → hard [`f64::max`] / [`f64::min`]
//! - `eps <= 0` → hard `|x|` / `sign(x)`
//! - `width <= 0` or inverted `[lo, hi]` → hard [`f64::clamp`] (bounds swapped if needed)
//!
//! # Example: smoothing a friction-factor transition
//!
//! ```rust
@@ -250,6 +272,8 @@ pub fn quintic_blend_derivative(a: f64, b: f64, x: f64, edge0: f64, edge1: f64)
/// arguments. A small `k` (relative to the scale of `a - b`) gives a sharp but
/// still infinitely differentiable maximum.
///
/// When `k <= 0`, returns the hard `a.max(b)` (no panic).
///
/// # Example
///
/// ```rust
@@ -259,17 +283,51 @@ pub fn quintic_blend_derivative(a: f64, b: f64, x: f64, edge0: f64, edge1: f64)
/// assert!((smooth_max(5.0, 1.0, 1e-3) - 5.0).abs() < 1e-2);
/// // Always an upper bound on the true max.
/// assert!(smooth_max(2.0, 2.0, 0.1) >= 2.0);
/// assert_eq!(smooth_max(5.0, 1.0, 0.0), 5.0);
/// ```
#[inline]
pub fn smooth_max(a: f64, b: f64, k: f64) -> f64 {
if k <= 0.0 {
return a.max(b);
}
0.5 * (a + b + ((a - b) * (a - b) + k * k).sqrt())
}
/// Analytic derivative of [`smooth_max`] with respect to `a` (`b`, `k` fixed).
///
/// Equals `0.5 * (1 + (a - b) / sqrt((a - b)² + k²))` for `k > 0`. When
/// `k <= 0`, returns the hard-max subgradient pick: `1` if `a > b`, `0` if
/// `a < b`, and `0.5` if `a == b`.
///
/// # Example
///
/// ```rust
/// use entropyk_core::smoothing::smooth_max_derivative;
///
/// assert!((smooth_max_derivative(5.0, 1.0, 1e-3) - 1.0).abs() < 1e-2);
/// assert!((smooth_max_derivative(1.0, 5.0, 1e-3)).abs() < 1e-2);
/// ```
#[inline]
pub fn smooth_max_derivative(a: f64, b: f64, k: f64) -> f64 {
if k <= 0.0 {
return if a > b {
1.0
} else if a < b {
0.0
} else {
0.5
};
}
0.5 * (1.0 + (a - b) / ((a - b) * (a - b) + k * k).sqrt())
}
/// C^∞ smooth minimum of `a` and `b` with sharpness controlled by `k > 0`.
///
/// Uses `0.5 * (a + b - sqrt((a - b)² + k²))`. Always `<= min(a, b)`, converging
/// to `min(a, b)` as `k -> 0`.
///
/// When `k <= 0`, returns the hard `a.min(b)` (no panic).
///
/// # Example
///
/// ```rust
@@ -277,12 +335,44 @@ pub fn smooth_max(a: f64, b: f64, k: f64) -> f64 {
///
/// assert!((smooth_min(5.0, 1.0, 1e-3) - 1.0).abs() < 1e-2);
/// assert!(smooth_min(2.0, 2.0, 0.1) <= 2.0);
/// assert_eq!(smooth_min(5.0, 1.0, -1.0), 1.0);
/// ```
#[inline]
pub fn smooth_min(a: f64, b: f64, k: f64) -> f64 {
if k <= 0.0 {
return a.min(b);
}
0.5 * (a + b - ((a - b) * (a - b) + k * k).sqrt())
}
/// Analytic derivative of [`smooth_min`] with respect to `a` (`b`, `k` fixed).
///
/// Equals `0.5 * (1 - (a - b) / sqrt((a - b)² + k²))` for `k > 0`. When
/// `k <= 0`, returns the hard-min subgradient pick: `1` if `a < b`, `0` if
/// `a > b`, and `0.5` if `a == b`.
///
/// # Example
///
/// ```rust
/// use entropyk_core::smoothing::smooth_min_derivative;
///
/// assert!((smooth_min_derivative(1.0, 5.0, 1e-3) - 1.0).abs() < 1e-2);
/// assert!((smooth_min_derivative(5.0, 1.0, 1e-3)).abs() < 1e-2);
/// ```
#[inline]
pub fn smooth_min_derivative(a: f64, b: f64, k: f64) -> f64 {
if k <= 0.0 {
return if a < b {
1.0
} else if a > b {
0.0
} else {
0.5
};
}
0.5 * (1.0 - (a - b) / ((a - b) * (a - b) + k * k).sqrt())
}
/// C^∞ approximation of the absolute value `|x|`.
///
/// Returns `sqrt(x² + eps²)`. Always `>= |x|`, converging to `|x|` as
@@ -290,6 +380,8 @@ pub fn smooth_min(a: f64, b: f64, k: f64) -> f64 {
/// `x / sqrt(x² + eps²)` is smooth through the origin, replacing the
/// non-differentiable corner of `|x|` that would otherwise break Newton.
///
/// When `eps <= 0`, returns the hard `|x|` (no panic).
///
/// # Example
///
/// ```rust
@@ -298,16 +390,21 @@ pub fn smooth_min(a: f64, b: f64, k: f64) -> f64 {
/// assert!((smooth_abs(0.0, 1e-3) - 1e-3).abs() < 1e-12);
/// assert!((smooth_abs(5.0, 1e-3) - 5.0).abs() < 1e-3);
/// assert!(smooth_abs(-3.0, 0.1) > 3.0);
/// assert_eq!(smooth_abs(-3.0, 0.0), 3.0);
/// ```
#[inline]
pub fn smooth_abs(x: f64, eps: f64) -> f64 {
if eps <= 0.0 {
return x.abs();
}
(x * x + eps * eps).sqrt()
}
/// Analytic derivative of [`smooth_abs`] with respect to `x`.
///
/// Equals `x / sqrt(x² + eps²)`, a smooth approximation of `sign(x)` that
/// passes continuously through `0` at the origin.
/// passes continuously through `0` at the origin. When `eps <= 0`, returns
/// hard `sign(x)` with `0.0` at the origin.
///
/// # Example
///
@@ -320,6 +417,15 @@ pub fn smooth_abs(x: f64, eps: f64) -> f64 {
/// ```
#[inline]
pub fn smooth_abs_derivative(x: f64, eps: f64) -> f64 {
if eps <= 0.0 {
return if x > 0.0 {
1.0
} else if x < 0.0 {
-1.0
} else {
0.0
};
}
x / (x * x + eps * eps).sqrt()
}
@@ -333,6 +439,10 @@ pub fn smooth_abs_derivative(x: f64, eps: f64) -> f64 {
/// always within `[lo, hi]`. For a flat identity region to exist, pass
/// `width <= (hi - lo) / 2`; larger widths still produce a bounded, C¹ result.
///
/// Degenerate cases (never panic):
/// - `hi < lo` → bounds are swapped
/// - `width <= 0` → hard `x.clamp(lo, hi)`
///
/// # Example
///
/// ```rust
@@ -343,9 +453,14 @@ pub fn smooth_abs_derivative(x: f64, eps: f64) -> f64 {
/// // Values well outside saturate to the bounds.
/// assert!((smooth_clamp(5.0, -1.0, 1.0, 0.1) - 1.0).abs() < 1e-9);
/// assert!((smooth_clamp(-5.0, -1.0, 1.0, 0.1) + 1.0).abs() < 1e-9);
/// assert_eq!(smooth_clamp(5.0, -1.0, 1.0, 0.0), 1.0);
/// ```
#[inline]
pub fn smooth_clamp(x: f64, lo: f64, hi: f64, width: f64) -> f64 {
let (lo, hi) = if hi < lo { (hi, lo) } else { (lo, hi) };
if width <= 0.0 {
return x.clamp(lo, hi);
}
// Lower ramp: equals `lo` for x <= lo, eases to identity `x` by x = lo + width.
// Stays within [lo, x] and is C¹ at both joins (zero slope at lo, unit slope
// at lo + width).
@@ -354,6 +469,36 @@ pub fn smooth_clamp(x: f64, lo: f64, hi: f64, width: f64) -> f64 {
hi + (lower - hi) * smoothstep(hi, hi - width, lower)
}
/// Analytic derivative of [`smooth_clamp`] with respect to `x`.
///
/// Chain-rules the lower then upper [`smoothstep`] ramps. Outside the
/// transition bands the slope is `0` (saturated) or `1` (identity interior).
/// Degenerate parameters match [`smooth_clamp`]: hard-clamp derivative is `1`
/// on the open interior and `0` elsewhere (including the kink points).
///
/// # Example
///
/// ```rust
/// use entropyk_core::smoothing::smooth_clamp_derivative;
///
/// // Flat exterior, unit slope deep interior.
/// assert!((smooth_clamp_derivative(5.0, -1.0, 1.0, 0.1)).abs() < 1e-9);
/// assert!((smooth_clamp_derivative(0.0, -1.0, 1.0, 0.1) - 1.0).abs() < 1e-9);
/// ```
#[inline]
pub fn smooth_clamp_derivative(x: f64, lo: f64, hi: f64, width: f64) -> f64 {
let (lo, hi) = if hi < lo { (hi, lo) } else { (lo, hi) };
if width <= 0.0 {
return if x > lo && x < hi { 1.0 } else { 0.0 };
}
let s_lo = smoothstep(lo, lo + width, x);
let d_lower = s_lo + (x - lo) * smoothstep_derivative(lo, lo + width, x);
let lower = lo + (x - lo) * s_lo;
let s_hi = smoothstep(hi, hi - width, lower);
let ds_hi = smoothstep_derivative(hi, hi - width, lower);
d_lower * (s_hi + (lower - hi) * ds_hi)
}
#[cfg(test)]
mod tests {
use super::*;
@@ -528,15 +673,84 @@ mod tests {
#[test]
fn test_smooth_max_derivative_matches_fd() {
// d/da smooth_max = 0.5 * (1 + (a-b)/sqrt((a-b)^2 + k^2))
let (b, k): (f64, f64) = (1.0, 0.3);
for a in [-1.0, 0.5, 1.0, 1.5, 3.0] {
let analytic = 0.5 * (1.0 + (a - b) / ((a - b) * (a - b) + k * k).sqrt());
let analytic = smooth_max_derivative(a, b, k);
let numeric = fd(|a| smooth_max(a, b, k), a, 1e-6);
assert!((analytic - numeric).abs() < 1e-5, "at a={}", a);
}
}
#[test]
fn test_smooth_min_derivative_matches_fd() {
let (b, k): (f64, f64) = (1.0, 0.3);
for a in [-1.0, 0.5, 1.0, 1.5, 3.0] {
let analytic = smooth_min_derivative(a, b, k);
let numeric = fd(|a| smooth_min(a, b, k), a, 1e-6);
assert!((analytic - numeric).abs() < 1e-5, "at a={}", a);
}
}
#[test]
fn test_smooth_clamp_derivative_matches_fd() {
let (lo, hi, width) = (-1.0, 1.0, 0.1);
// Exterior, both ramps, deep interior.
for x in [-5.0, -1.05, -0.95, -0.5, 0.0, 0.5, 0.95, 1.05, 5.0] {
let analytic = smooth_clamp_derivative(x, lo, hi, width);
let numeric = fd(|x| smooth_clamp(x, lo, hi, width), x, 1e-6);
assert!(
(analytic - numeric).abs() < 1e-5,
"smooth_clamp deriv mismatch at {}: {} vs {}",
x,
analytic,
numeric
);
}
}
#[test]
fn test_smooth_clamp_c1_at_junctions() {
// Slope continuity across the four joins: lo, lo+width, hi-width, hi.
let (lo, hi, width) = (-1.0, 1.0, 0.2);
let junctions = [lo, lo + width, hi - width, hi];
let h = 1e-5;
for j in junctions {
let slope_left = fd(|x| smooth_clamp(x, lo, hi, width), j - h, h);
let slope_right = fd(|x| smooth_clamp(x, lo, hi, width), j + h, h);
let analytic = smooth_clamp_derivative(j, lo, hi, width);
assert!(
(slope_left - slope_right).abs() < 5e-3,
"slope jump at {}: left={} right={}",
j,
slope_left,
slope_right
);
assert!(
(analytic - 0.5 * (slope_left + slope_right)).abs() < 5e-3,
"analytic slope mismatch at {}: analytic={} fd_avg={}",
j,
analytic,
0.5 * (slope_left + slope_right)
);
}
}
#[test]
fn test_degenerate_parameters_do_not_panic() {
assert_eq!(smooth_max(3.0, 1.0, 0.0), 3.0);
assert_eq!(smooth_max(3.0, 1.0, -0.5), 3.0);
assert_eq!(smooth_min(3.0, 1.0, -1.0), 1.0);
assert_eq!(smooth_abs(-2.0, 0.0), 2.0);
assert_eq!(smooth_abs_derivative(-2.0, 0.0), -1.0);
assert_eq!(smooth_clamp(5.0, -1.0, 1.0, 0.0), 1.0);
assert_eq!(smooth_clamp_derivative(0.0, -1.0, 1.0, 0.0), 1.0);
assert_eq!(smooth_clamp_derivative(5.0, -1.0, 1.0, 0.0), 0.0);
// Inverted bounds swap.
assert!((smooth_clamp(0.0, 1.0, -1.0, 0.1) - 0.0).abs() < 1e-9);
assert_eq!(smooth_max_derivative(3.0, 1.0, 0.0), 1.0);
assert_eq!(smooth_min_derivative(3.0, 1.0, 0.0), 0.0);
}
#[test]
fn test_smooth_abs_properties_and_derivative() {
assert!((smooth_abs(0.0, 1e-3) - 1e-3).abs() < 1e-12);