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    for bank in &lane.modal {
50        for a in sva_formula::modal::atoms(bank, sva_formula::Origin::UNKNOWN) {
51            sum = sum + eval_atom(&a, t)?;
52        }
53    }
54    Ok(sum)
55}
56
57pub fn eval_spectral_sum(n: &SpectralSum, c: usize, t: f64) -> Result<C64, CollapseError> {
58    eval_lane(&n.lanes[c.min(n.lanes.len() - 1)], t)
59}
60
61/// The fallback for a closed form with no spectral sum at all: `tanh`, `sat`, `abs`, `log`, `sqrt`,
62/// a non-integer power, a non-affine `sin`. Nothing here is claimed exact.
63pub fn eval_body(
64    fm: &Body,
65    component: usize,
66    t: f64,
67    refs: &dyn Refs,
68) -> Result<C64, CollapseError> {
69    let of = |p: &sva_formula::Part| eval_body(&p.body, component, t, refs);
70    let value = match fm {
71        Body::Const(c) => *c,
72        Body::Line => C64::real(t),
73        Body::Add(parts) => parts.iter().try_fold(C64::ZERO, |a, p| Ok(a + of(p)?))?,
74        Body::Mul(parts) => parts.iter().try_fold(C64::ONE, |a, p| Ok(a * of(p)?))?,
75        Body::Div(a, b) => of(a)? / of(b)?,
76        Body::Pow(a, n) => power(of(a)?, *n),
77        Body::Apply(op, a) => unary(*op, of(a)?),
78        Body::Fold(op, parts) => fold(*op, parts, component, t, refs)?,
79        Body::Shift { by, of: inner } => eval_body(&inner.body, component, t - by, refs)?,
80        Body::Warp { at, of: inner } => {
81            let when = eval_body(&at.body, component, t, refs)?.re;
82            eval_body(&inner.body, component, when, refs)?
83        }
84        Body::Crop {
85            of: inner,
86            l,
87            r,
88            rise,
89            fall,
90        } => match raised_cosine(t, l.value(), r.value(), *rise, *fall) {
91            0.0 => C64::ZERO,
92            gain => of(inner)?.scale(gain),
93        },
94        Body::Channel(inner, k) => eval_body(&inner.body, usize::from(*k), t, refs)?,
95        Body::Delta { order, .. } => {
96            return Err(CollapseError::SingularInCt {
97                at: t,
98                order: i32::from(*order),
99            });
100        }
101        Body::Pv(_) => return Err(CollapseError::SingularInCt { at: t, order: -1 }),
102        Body::Keyed { seed, of: key } => C64::real(sva_formula::draw(*seed, of(key)?.re)),
103        Body::Join(parts) => {
104            let widths: Vec<usize> = parts.iter().map(|p| width_of(&p.body, refs)).collect();
105            let (at, inner) = lane_of(&widths, component)
106                .ok_or(CollapseError::NotEvaluable("a component past the width"))?;
107            eval_body(&parts[at].body, inner, t, refs)?
108        }
109        Body::Modal(bank) => {
110            let mut sum = C64::ZERO;
111            for a in sva_formula::modal::atoms(bank, sva_formula::Origin::UNKNOWN) {
112                sum = sum + eval_atom(&a, t)?;
113            }
114            sum
115        }
116        Body::Node(id) => refs.value(*id, component, t)?,
117        other => return Err(CollapseError::NotEvaluable(sketch(other))),
118    };
119    match value.is_finite() {
120        true => Ok(value),
121        false => Err(CollapseError::NotEvaluable(
122            "a division or a remainder by zero",
123        )),
124    }
125}
126
127/// FORMAT 15.6's window: one over the plateau, a raised cosine over each shoulder, zero out.
128fn raised_cosine(t: f64, l: f64, r: f64, rise: f64, fall: f64) -> f64 {
129    if t < l || t >= r {
130        return 0.0;
131    }
132    let opening = shoulder(t - l, rise);
133    let closing = shoulder(r - t, fall);
134    opening.min(closing)
135}
136
137fn shoulder(into: f64, span: f64) -> f64 {
138    if span <= 0.0 || into >= span {
139        return 1.0;
140    }
141    0.5 - 0.5 * (std::f64::consts::PI * into / span).cos()
142}
143
144/// Which operand of a `join` holds one component of the joined value, and which of its own
145/// components that is.
146pub fn lane_of(widths: &[usize], component: usize) -> Option<(usize, usize)> {
147    let mut left = component;
148    for (at, width) in widths.iter().enumerate() {
149        if left < *width {
150            return Some((at, left));
151        }
152        left -= width;
153    }
154    None
155}
156
157/// The component count a written subterm answers for. Only `join` widens, and only `ch` and
158/// a node narrow back.
159pub(crate) fn width_of(f: &Body, refs: &dyn Refs) -> usize {
160    match f {
161        Body::Join(parts) => parts.iter().map(|p| width_of(&p.body, refs)).sum(),
162        Body::Channel(..) => 1,
163        Body::Node(id) => refs.width(*id).max(1),
164        other => children(other)
165            .iter()
166            .map(|p| width_of(&p.body, refs))
167            .max()
168            .unwrap_or(1),
169    }
170}
171
172fn power(x: C64, n: i32) -> C64 {
173    match n {
174        0.. => x.powi(n as u32),
175        _ => x.powi(n.unsigned_abs()).inv(),
176    }
177}
178
179pub fn unary(op: Unary, x: C64) -> C64 {
180    match op {
181        Unary::Exp => x.exp(),
182        Unary::Sin => C64::new(x.re.sin() * x.im.cosh(), x.re.cos() * x.im.sinh()),
183        Unary::Cos => C64::new(x.re.cos() * x.im.cosh(), -x.re.sin() * x.im.sinh()),
184        Unary::Tanh => C64::real(x.re.tanh()),
185        Unary::Sat => C64::real(x.re.clamp(-1.0, 1.0)),
186        Unary::Abs => C64::real(x.abs()),
187        Unary::Log => C64::real(x.re.ln()),
188        Unary::Sqrt => C64::real(x.re.sqrt()),
189    }
190}
191
192fn fold(
193    op: Fold,
194    parts: &[sva_formula::Part],
195    component: usize,
196    t: f64,
197    refs: &dyn Refs,
198) -> Result<C64, CollapseError> {
199    let mut it = parts.iter();
200    let head = it.next().expect("a fold holds one part");
201    let first = eval_body(&head.body, component, t, refs)?;
202    it.try_fold(first, |acc, p| {
203        let v = eval_body(&p.body, component, t, refs)?;
204        Ok(C64::real(match op {
205            Fold::Max => acc.re.max(v.re),
206            Fold::Min => acc.re.min(v.re),
207            Fold::Mod => acc.re.rem_euclid(v.re),
208        }))
209    })
210}
211
212fn sketch(f: &Body) -> &'static str {
213    match f {
214        Body::Param(_) => "an unsubstituted parameter",
215        Body::Index(_) => "a free series index",
216        Body::Deriv { .. } => "a derivative",
217        Body::Rational(_) => "a rational",
218        Body::Series(_) => "a series",
219        _ => "this subterm",
220    }
221}
sva_samples/collapse/point.rs

sva_samples/collapse/
point.rs
