sva_samples/collapse/
point.rs

1// Concern: evaluates one closed form at one instant and routes the components it holds | Non-concern: what the samples are labelled (collapse.rs) | IO: (&SpectralSum or &Body, t, component) -> C64
2
3use sva_formula::closed_form::{Fold, Unary, children};
4use sva_formula::spectral_sum::atom::{Singular, SpectralAtom};
5use sva_formula::{Body, C64, Lane, NodeId, SpectralSum};
6
7use crate::error::CollapseError;
8
9/// How a `Body::Node` reaches a value at one instant, and how many components it holds: a
10/// closed form on its own holds no node, so the caller holding the graph answers both.
11pub trait Refs {
12    fn value(&self, id: NodeId, component: usize, t: f64) -> Result<C64, CollapseError>;
13    fn width(&self, id: NodeId) -> usize;
14}
15
16/// What a caller with no graph behind it answers a node with.
17pub(crate) struct NoRefs;
18
19impl Refs for NoRefs {
20    fn value(&self, _: NodeId, _: usize, _: f64) -> Result<C64, CollapseError> {
21        Err(CollapseError::NotEvaluable("a node"))
22    }
23
24    fn width(&self, _: NodeId) -> usize {
25        1
26    }
27}
28
29/// The six atom factors in closed form. A delta has no ordinary value, and a principal
30/// value has none at its own pole.
31pub fn eval_atom(a: &SpectralAtom, t: f64) -> Result<C64, CollapseError> {
32    if let Singular::Delta { at, order } = a.sing {
33        return Err(CollapseError::SingularInCt {
34            at,
35            order: i32::from(order),
36        });
37    }
38    a.smooth_at(t).ok_or(CollapseError::SingularInCt {
39        at: a.pole.map_or(t, |p| p.at.re),
40        order: -1,
41    })
42}
43
44pub fn eval_lane(lane: &Lane, t: f64) -> Result<C64, CollapseError> {
45    let mut sum = C64::ZERO;
46    for a in &lane.atoms {
47        sum = sum + eval_atom(a, t)?;
48    }
49    with_modal(lane, sum, t)
50}
51
52/// Every atom `live` leaves out is exactly zero here.
53pub(crate) fn eval_among(lane: &Lane, live: &[usize], t: f64) -> Result<C64, CollapseError> {
54    let mut sum = C64::ZERO;
55    for &i in live {
56        sum = sum + eval_atom(&lane.atoms[i], t)?;
57    }
58    with_modal(lane, sum, t)
59}
60
61fn with_modal(lane: &Lane, mut sum: C64, t: f64) -> Result<C64, CollapseError> {
62    for bank in &lane.modal {
63        for a in sva_formula::modal::atoms(bank, sva_formula::Origin::UNKNOWN) {
64            sum = sum + eval_atom(&a, t)?;
65        }
66    }
67    Ok(sum)
68}
69
70pub fn eval_spectral_sum(n: &SpectralSum, c: usize, t: f64) -> Result<C64, CollapseError> {
71    eval_lane(&n.lanes[c.min(n.lanes.len() - 1)], t)
72}
73
74/// A shut crop zeroes a factor beside it that passes a double.
75fn product(
76    factors: impl Iterator<Item = Result<C64, CollapseError>>,
77) -> Result<C64, CollapseError> {
78    let mut held = Ok(C64::ONE);
79    for factor in factors {
80        match (factor, &held) {
81            (Ok(v), _) if v.is_zero() => return Ok(C64::ZERO),
82            (Ok(v), Ok(acc)) => held = Ok(*acc * v),
83            (Err(e), Ok(_)) => held = Err(e),
84            (_, Err(_)) => {}
85        }
86    }
87    held
88}
89
90/// The fallback for a closed form with no spectral sum at all: `tanh`, `sat`, `abs`, `log`, `sqrt`,
91/// a non-integer power, a non-affine `sin`. Nothing here is claimed exact.
92pub fn eval_body(
93    fm: &Body,
94    component: usize,
95    t: f64,
96    refs: &dyn Refs,
97) -> Result<C64, CollapseError> {
98    let of = |p: &sva_formula::Part| eval_body(&p.body, component, t, refs);
99    let value = match fm {
100        Body::Const(c) => *c,
101        Body::Line => C64::real(t),
102        Body::Add(parts) => parts.iter().try_fold(C64::ZERO, |a, p| Ok(a + of(p)?))?,
103        Body::Mul(parts) => product(parts.iter().map(of))?,
104        Body::Div(a, b) => of(a)? / of(b)?,
105        Body::Pow(a, n) => power(of(a)?, *n),
106        Body::Apply(op, a) => unary(*op, of(a)?),
107        Body::Fold(op, parts) => fold(*op, parts, component, t, refs)?,
108        Body::Shift { by, of: inner } => eval_body(&inner.body, component, t - by, refs)?,
109        Body::Warp { at, of: inner } => {
110            let when = eval_body(&at.body, component, t, refs)?.re;
111            eval_body(&inner.body, component, when, refs)?
112        }
113        Body::Crop {
114            of: inner,
115            l,
116            r,
117            rise,
118            fall,
119        } => match crop_gain(t, l.value(), r.value(), *rise, *fall) {
120            0.0 => C64::ZERO,
121            gain => of(inner)?.scale(gain),
122        },
123        Body::Channel(inner, k) => eval_body(&inner.body, usize::from(*k), t, refs)?,
124        Body::Delta { order, .. } => {
125            return Err(CollapseError::SingularInCt {
126                at: t,
127                order: i32::from(*order),
128            });
129        }
130        Body::Pv(_) => return Err(CollapseError::SingularInCt { at: t, order: -1 }),
131        Body::Keyed { seed, of: key } => C64::real(
132            sva_formula::draw_nearest(*seed, of(key)?.re)
133                .ok_or(CollapseError::NotEvaluable("a key past any step"))?,
134        ),
135        Body::Join(parts) => {
136            let widths: Vec<usize> = parts.iter().map(|p| width_of(&p.body, refs)).collect();
137            let (at, inner) = lane_of(&widths, component)
138                .ok_or(CollapseError::NotEvaluable("a component past the width"))?;
139            eval_body(&parts[at].body, inner, t, refs)?
140        }
141        Body::Modal(bank) => {
142            let mut sum = C64::ZERO;
143            for a in sva_formula::modal::atoms(bank, sva_formula::Origin::UNKNOWN) {
144                sum = sum + eval_atom(&a, t)?;
145            }
146            sum
147        }
148        Body::Node(id) => refs.value(*id, component, t)?,
149        Body::Run(run) => super::run::at(run, t),
150        other => return Err(CollapseError::NotEvaluable(sketch(other))),
151    };
152    finite(value)
153}
154
155/// `eval_body`'s sum from +0 over only `live`: one left out is exact zero at `t`.
156pub(crate) fn eval_addends(
157    parts: &[&sva_formula::Part],
158    live: &[usize],
159    component: usize,
160    t: f64,
161) -> Result<C64, CollapseError> {
162    let sum = live.iter().try_fold(C64::ZERO, |held, &i| {
163        Ok(held + eval_body(&parts[i].body, component, t, &NoRefs)?)
164    })?;
165    finite(sum)
166}
167
168fn finite(value: C64) -> Result<C64, CollapseError> {
169    match value.is_finite() {
170        true => Ok(value),
171        false => Err(CollapseError::NotEvaluable(
172            "a division or a remainder by zero, or an infinite value",
173        )),
174    }
175}
176
177/// FORMAT 15.6's window: one over the plateau, a raised cosine over each shoulder, zero out.
178pub fn crop_gain(t: f64, l: f64, r: f64, rise: f64, fall: f64) -> f64 {
179    if t < l || t >= r {
180        return 0.0;
181    }
182    shoulders(t, l, r, rise, fall)
183}
184
185/// The window's gain at an instant already known to lie inside it.
186pub fn shoulders(t: f64, l: f64, r: f64, rise: f64, fall: f64) -> f64 {
187    let opening = shoulder(t - l, rise);
188    let closing = shoulder(r - t, fall);
189    opening.min(closing)
190}
191
192fn shoulder(into: f64, span: f64) -> f64 {
193    if span <= 0.0 || into >= span {
194        return 1.0;
195    }
196    0.5 - 0.5 * (std::f64::consts::PI * into / span).cos()
197}
198
199/// Which operand of a `join` holds one component of the joined value, and which of its own
200/// components that is.
201pub fn lane_of(widths: &[usize], component: usize) -> Option<(usize, usize)> {
202    let mut left = component;
203    for (at, width) in widths.iter().enumerate() {
204        if left < *width {
205            return Some((at, left));
206        }
207        left -= width;
208    }
209    None
210}
211
212/// The component count a written subterm answers for. Only `join` widens, and only `ch` and
213/// a node narrow back.
214pub(crate) fn width_of(f: &Body, refs: &dyn Refs) -> usize {
215    match f {
216        Body::Join(parts) => parts.iter().map(|p| width_of(&p.body, refs)).sum(),
217        Body::Channel(..) => 1,
218        Body::Node(id) => refs.width(*id).max(1),
219        other => children(other)
220            .iter()
221            .map(|p| width_of(&p.body, refs))
222            .max()
223            .unwrap_or(1),
224    }
225}
226
227fn power(x: C64, n: i32) -> C64 {
228    match n {
229        0.. => x.powi(n as u32),
230        _ => x.powi(n.unsigned_abs()).inv(),
231    }
232}
233
234pub fn unary(op: Unary, x: C64) -> C64 {
235    match op {
236        Unary::Exp => x.exp(),
237        Unary::Sin => C64::new(x.re.sin() * x.im.cosh(), x.re.cos() * x.im.sinh()),
238        Unary::Cos => C64::new(x.re.cos() * x.im.cosh(), -x.re.sin() * x.im.sinh()),
239        Unary::Tanh => C64::real(x.re.tanh()),
240        Unary::Sat => C64::real(x.re.clamp(-1.0, 1.0)),
241        Unary::Abs => C64::real(x.abs()),
242        Unary::Log => C64::real(x.re.ln()),
243        Unary::Sqrt => C64::real(x.re.sqrt()),
244        Unary::Step => C64::real(sva_formula::affine::step(x.re)),
245    }
246}
247
248fn fold(
249    op: Fold,
250    parts: &[sva_formula::Part],
251    component: usize,
252    t: f64,
253    refs: &dyn Refs,
254) -> Result<C64, CollapseError> {
255    let mut it = parts.iter();
256    let head = it.next().expect("a fold holds one part");
257    let first = eval_body(&head.body, component, t, refs)?;
258    it.try_fold(first, |acc, p| {
259        let v = eval_body(&p.body, component, t, refs)?;
260        Ok(C64::real(match op {
261            Fold::Max => acc.re.max(v.re),
262            Fold::Min => acc.re.min(v.re),
263            Fold::Mod => acc.re.rem_euclid(v.re),
264        }))
265    })
266}
267
268fn sketch(f: &Body) -> &'static str {
269    match f {
270        Body::Param(_) => "an unsubstituted parameter",
271        Body::Index(_) => "a free series index",
272        Body::Deriv { .. } => "a derivative",
273        Body::Rational(_) => "a rational",
274        Body::Series(_) => "a series",
275        _ => "this subterm",
276    }
277}
sva_samples/collapse/point.rs

sva_samples/collapse/
point.rs
