From dd054e901609ee3a8b14a633f94426a37861e236 Mon Sep 17 00:00:00 2001 From: numzero Date: Tue, 25 Jun 2024 11:23:37 +0300 Subject: [PATCH] Move stuff around --- src/bin/flat/fns.rs | 110 +++++++++++++++++++++++++++++++++++ src/bin/flat/main.rs | 134 ++++--------------------------------------- 2 files changed, 121 insertions(+), 123 deletions(-) create mode 100644 src/bin/flat/fns.rs diff --git a/src/bin/flat/fns.rs b/src/bin/flat/fns.rs new file mode 100644 index 0000000..c6d94d8 --- /dev/null +++ b/src/bin/flat/fns.rs @@ -0,0 +1,110 @@ +use crate::FloatExt2; + +pub struct LinearLimiter { + pub min: f32, + pub max: f32, +} + +impl LinearLimiter { + pub fn value(&self, x: f32) -> f32 { (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0) } + pub fn derivative(&self, x: f32) -> f32 { if x.abs() > self.min && x.abs() < self.max { x.signum() / (self.max - self.min) } else { 0.0 } } +} + +pub struct QuadraticAccelerator { + pub internal: f32, + pub external: f32, +} + +/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1. +fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { + if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) } +} + +/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0. +fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { + if t.abs() <= lim { f(t) } else { val } +} + +impl QuadraticAccelerator { + fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) } + fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 } + fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() } + + pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) } + pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) } + pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) } + pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) } + pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) } +} + +#[cfg(test)] +mod test { + use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq}; + + #[test] + fn test_linear_limiter() { + let testee = super::LinearLimiter { min: 20.0, max: 30.0 }; + let ε = 1.0e-4f32; + let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t. + let margin = 1.0 / 16.0; + let mul = 1.0 + margin; + for x in itertools_num::linspace(0., testee.min, 10) { + assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε); + assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε); + } + for x in itertools_num::linspace(testee.max, mul * testee.max, 10) { + assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε); + assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε); + } + for x in itertools_num::linspace(-mul * testee.max, mul * testee.max, 100) { + // Currently, the derivative is discontinuous at ±min and ±max... let’s just skip these for now. + if x.abs().abs_diff_eq(&testee.min, δ) || x.abs().abs_diff_eq(&testee.max, δ) { + continue; + } + let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ); + let df_expl = testee.derivative(x); + assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n"); + } + } + + #[test] + fn test_quadratic_accelerator() { + let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 }; + let ε = 1.0e-4f32; + let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t. + let margin = 1.0 / 16.0; + let mul = 1.0 + margin; + assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε); + assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε); + assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε); + assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε); + for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) { + let ux = testee.u(x); + let xux = testee.x(ux); + assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n"); + + let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ); + let du_expl = testee.du(x); + assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n"); + + let dudx = du_expl * testee.dx(ux); + assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n"); + + let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ); + let d2u_expl = testee.d2u(x); + assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n"); + } + for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) { + let xu = testee.x(u); + let uxu = testee.u(xu); + assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n"); + + let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ); + let dx_expl = testee.dx(u); + assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n"); + + let dudx = testee.du(xu) * dx_expl; + assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n"); + } + } +} diff --git a/src/bin/flat/main.rs b/src/bin/flat/main.rs index 12e4181..244b2f7 100644 --- a/src/bin/flat/main.rs +++ b/src/bin/flat/main.rs @@ -4,6 +4,7 @@ use flo_canvas::*; use glam::*; mod riemann; +mod fns; use riemann::{Tens2, Decomp2, Metric, trace_iter}; use shape::Shape; @@ -624,129 +625,6 @@ impl Tube { pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) } } -mod fns { - use crate::FloatExt2; - - pub struct LinearLimiter { - pub min: f32, - pub max: f32, - } - - impl LinearLimiter { - pub fn value(&self, x: f32) -> f32 { (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0) } - pub fn derivative(&self, x: f32) -> f32 { if x.abs() > self.min && x.abs() < self.max { x.signum() / (self.max - self.min) } else { 0.0 } } - } - - pub struct QuadraticAccelerator { - pub internal: f32, - pub external: f32, - } - - /// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1. - fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { - if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) } - } - - /// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0. - fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 { - if t.abs() <= lim { f(t) } else { val } - } - - impl QuadraticAccelerator { - fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) } - fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 } - fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() } - - pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) } - pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) } - pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) } - pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) } - pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) } - } - - #[cfg(test)] - mod test { - use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq}; - - #[test] - fn test_linear_limiter() { - let testee = super::LinearLimiter { min: 20.0, max: 30.0 }; - let ε = 1.0e-4f32; - let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t. - let margin = 1.0 / 16.0; - let mul = 1.0 + margin; - for x in itertools_num::linspace(0., testee.min, 10) { - assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε); - assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε); - } - for x in itertools_num::linspace(testee.max, mul * testee.max, 10) { - assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε); - assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε); - } - for x in itertools_num::linspace(-mul * testee.max, mul * testee.max, 100) { - // Currently, the derivative is discontinuous at ±min and ±max... let’s just skip these for now. - if x.abs().abs_diff_eq(&testee.min, δ) || x.abs().abs_diff_eq(&testee.max, δ) { - continue; - } - let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ); - let df_expl = testee.derivative(x); - assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit: {df_expl}\n"); - } - } - - #[test] - fn test_quadratic_accelerator() { - let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 }; - let ε = 1.0e-4f32; - let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t. - let margin = 1.0 / 16.0; - let mul = 1.0 + margin; - assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε); - assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε); - assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε); - assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε); - for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) { - let ux = testee.u(x); - let xux = testee.x(ux); - assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x): {ux}\nx(u(x)): {xux}\n"); - - let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ); - let du_expl = testee.du(x); - assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit: {du_expl}\n"); - - let dudx = du_expl * testee.dx(ux); - assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n"); - - let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ); - let d2u_expl = testee.d2u(x); - assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit: {d2u_expl}\n"); - } - for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) { - let xu = testee.x(u); - let uxu = testee.u(xu); - assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u): {xu}\nu(x(u)): {uxu}\n"); - - let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ); - let dx_expl = testee.dx(u); - assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit: {dx_expl}\n"); - - let dudx = testee.du(xu) * dx_expl; - assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n"); - } - } - } -} - -trait FloatExt2 { - fn lerp(&self, t: f32) -> f32; - fn inverse_lerp(&self, y: f32) -> f32; -} - -impl FloatExt2 for (f32, f32) { - fn lerp(&self, t: f32) -> f32 { f32::lerp(self.0, self.1, t) } - fn inverse_lerp(&self, y: f32) -> f32 { f32::inverse_lerp(self.0, self.1, y) } -} - impl Metric for Tube { fn sqrt_at(&self, pos: Vec2) -> Decomp2 { let sx = self.fx().value(pos.x); @@ -802,3 +680,13 @@ fn test_tube_metric_derivs() { } } } + +trait FloatExt2 { + fn lerp(&self, t: f32) -> f32; + fn inverse_lerp(&self, y: f32) -> f32; +} + +impl FloatExt2 for (f32, f32) { + fn lerp(&self, t: f32) -> f32 { f32::lerp(self.0, self.1, t) } + fn inverse_lerp(&self, y: f32) -> f32 { f32::inverse_lerp(self.0, self.1, y) } +}