Skip to main content

sva_samples/collapse/
truncate.rs

1// Concern: truncates every series to the terms the profile leaves, once per collapse | Non-concern: evaluating what is left (point.rs) | IO: (&SpectralSum or &Body) -> the same, series-free
2
3use std::f64::consts::TAU;
4
5use sva_formula::affine::{Axis, axis, exact_constant};
6use sva_formula::closed_form::{Bound, Series, children, map_children};
7use sva_formula::series::{mentions_line, ratio, substitute};
8use sva_formula::spectral_sum::atom::{Exp, Factors, Singular, SpectralAtom};
9use sva_formula::table::series::{Shape, read};
10use sva_formula::{
11    Body, C64, Codomain, Env, Lane, NodeId, ParamId, Part, Run, SpectralSum, Ty, Unary, Var, lines,
12};
13
14use crate::error::CollapseError;
15use crate::profile::Profile;
16
17/// What a series is truncated against: the observation's own ceiling, the amplitude a
18/// dropped tail provably stays within, and the floor a tail no bound sums falls under.
19#[derive(Clone, Copy, Debug, PartialEq)]
20pub struct Audible {
21    ceiling: f64,
22    precision: f64,
23    floor_db: f64,
24}
25
26impl Audible {
27    pub fn of(profile: &Profile, rate: u32) -> Audible {
28        let ceiling = profile.ceiling(rate);
29        Audible {
30            ceiling,
31            precision: profile.half_lsb(),
32            floor_db: profile.floor(ceiling),
33        }
34    }
35}
36
37/// A line series is placed analytically under `sva_formula`'s own bound; an expanded one
38/// becomes a formula the sample loop walks, which is what this caps.
39const MAX_EXPANDED_TERMS: usize = 1 << 13;
40
41/// FORMAT 6.2 truncates a series once, at collapse: every term the ceiling and the precision
42/// leave becomes an ordinary atom before the first sample is read.
43pub fn spectral_sum(n: &SpectralSum, band: Audible) -> Result<SpectralSum, CollapseError> {
44    if n.lanes.iter().all(|l| l.series.is_empty()) {
45        return Ok(n.clone());
46    }
47    let mut lanes = Vec::with_capacity(n.lanes.len());
48    for lane in &n.lanes {
49        let mut held = Lane {
50            series: Vec::new(),
51            ..lane.clone()
52        };
53        for s in &lane.series {
54            held.atoms.extend(atoms(s, band)?);
55        }
56        lanes.push(held);
57    }
58    Ok(SpectralSum::of(n.var, lanes))
59}
60
61/// The same truncation over a written closed form, whose series a spectral sum never reached.
62pub fn written(f: &Body, band: Audible) -> Result<Body, CollapseError> {
63    if let Body::Series(s) = f {
64        return expanded(s, band);
65    }
66    let mut found = None;
67    let out = map_children(f, |p| match written(&p.body, band) {
68        Ok(body) => Part::new(p.origin, body),
69        Err(e) => {
70            found = Some(e);
71            p.clone()
72        }
73    });
74    match found {
75        Some(e) => Err(e),
76        None => Ok(out),
77    }
78}
79
80fn atoms(s: &Series, band: Audible) -> Result<Vec<SpectralAtom>, CollapseError> {
81    let (body, window) = sva_formula::crop_peeled(&s.term.body);
82    let bare = Series {
83        term: Part::new(s.term.origin, body),
84        ..s.clone()
85    };
86    let taken = enumerated(&bare, band);
87    if !taken.is_empty() {
88        return Ok(taken
89            .into_iter()
90            .map(|l| {
91                SpectralAtom::new(
92                    l.amp,
93                    Factors {
94                        exp: Some(Exp::at(0.0, TAU * l.hz)),
95                        ind: window,
96                        ..Factors::NONE
97                    },
98                    Singular::Regular,
99                    s.term.origin,
100                )
101            })
102            .collect());
103    }
104    let body = expanded(s, band)?;
105    let held = sva_formula::normalize(&body, Var::T)
106        .map_err(|_| CollapseError::NotEvaluable("a series term"))?;
107    Ok(held.lanes.into_iter().flat_map(|l| l.atoms).collect())
108}
109
110/// A line series keeps its terms as runs; anything else is summed term by term.
111fn expanded(s: &Series, band: Audible) -> Result<Body, CollapseError> {
112    let taken = enumerated(s, band);
113    if !taken.is_empty() {
114        let runs = Run::of(&taken).into_iter();
115        return Ok(sum(runs.map(|r| Body::Run(Box::new(r))).collect()));
116    }
117    let count = terms(s, band)?;
118    let mut parts = Vec::with_capacity(count);
119    for i in 0..count {
120        let form = substitute(&s.term.body, s.index, (s.lo + i as i64) as f64);
121        parts.push(written(&form, band)?);
122    }
123    Ok(sum(parts))
124}
125
126fn enumerated(s: &Series, band: Audible) -> Vec<sva_formula::Line> {
127    match read(&s.term.body) {
128        Some(Shape::Lines(_)) => lines(s, band.ceiling, band.floor_db, band.precision).taken,
129        _ => Vec::new(),
130    }
131}
132
133fn sum(parts: Vec<Body>) -> Body {
134    match parts.len() {
135        0 => Body::Const(C64::ZERO),
136        1 => parts.into_iter().next().expect("one part"),
137        _ => Body::Add(parts.into_iter().map(Part::bare).collect()),
138    }
139}
140
141/// A nesting is priced whole before it expands: a written bound counts like the floor's.
142fn terms(s: &Series, band: Audible) -> Result<usize, CollapseError> {
143    let cost = cost(s, band).ok_or(CollapseError::NotEvaluable(
144        "a series whose term count no coefficient bounds",
145    ))?;
146    if cost.depth == 1 && cost.own > MAX_EXPANDED_TERMS {
147        return Err(CollapseError::NotEvaluable(
148            "a series of more terms than one instant expands",
149        ));
150    }
151    match cost.terms > MAX_EXPANDED_TERMS {
152        true => Err(CollapseError::NestedSeries {
153            depth: cost.depth,
154            terms: cost.terms,
155            bound: MAX_EXPANDED_TERMS,
156        }),
157        false => Ok(cost.own),
158    }
159}
160
161/// Series deep, terms the whole nesting takes, terms this level alone takes.
162struct Cost {
163    depth: usize,
164    terms: usize,
165    own: usize,
166}
167
168fn cost(s: &Series, band: Audible) -> Option<Cost> {
169    let own = counted(s, band)?;
170    let inside = substitute(&s.term.body, s.index, s.lo as f64);
171    let (depth, inner) = price(&inside, band)?;
172    Some(Cost {
173        depth: depth + 1,
174        terms: own.saturating_mul(inner),
175        own,
176    })
177}
178
179/// A sum of series costs their counts added; every other spelling multiplies.
180fn price(f: &Body, band: Audible) -> Option<(usize, usize)> {
181    if let Body::Series(s) = f {
182        let cost = cost(s, band)?;
183        return Some((cost.depth, cost.terms));
184    }
185    let summed = matches!(f, Body::Add(_));
186    let mut depth = 0;
187    let mut terms: usize = if summed { 0 } else { 1 };
188    for p in children(f) {
189        let (d, n) = price(&p.body, band)?;
190        depth = depth.max(d);
191        terms = match summed {
192            true => terms.saturating_add(n),
193            false => terms.saturating_mul(n),
194        };
195    }
196    Some((depth, terms.max(1)))
197}
198
199/// A geometric magnitude bound stops where its whole tail rounds away; any other stand-in at
200/// the profile's floor.
201fn counted(s: &Series, band: Audible) -> Option<usize> {
202    let hi = match s.hi {
203        Bound::Finite(n) => return usize::try_from((n - s.lo + 1).max(0)).ok(),
204        Bound::Infinite => MAX_EXPANDED_TERMS,
205    };
206    let bound = bound(&s.term.body, band)?;
207    let ratio = match bound.exact {
208        true => ratio(&bound.body, s.index).filter(|r| *r < 1.0),
209        false => None,
210    };
211    let floor = 10f64.powf(band.floor_db / 20.0);
212    let mut peak = 0.0f64;
213    for i in 0..hi {
214        let at = substitute(&bound.body, s.index, (s.lo + i as i64) as f64);
215        let held = exact_constant(&at)?.abs();
216        peak = peak.max(held);
217        let gone = match ratio {
218            Some(r) => held / (1.0 - r) <= band.precision,
219            None => peak > 0.0 && held < peak * floor,
220        };
221        if gone {
222            return Some(i.max(1));
223        }
224    }
225    None
226}
227
228/// A term's coefficient in the index. `exact` where it bounds the term's magnitude; elsewhere
229/// it only stands in, a turning factor or a name read as one.
230struct Bounded {
231    body: Body,
232    exact: bool,
233}
234
235fn bound(f: &Body, band: Audible) -> Option<Bounded> {
236    let held = |body: Body, exact: bool| Some(Bounded { body, exact });
237    let under = |p: &Part| bound(&p.body, band).map(|b| (Part::new(p.origin, b.body), b.exact));
238    let all = |parts: &[Part], wrap: &dyn Fn(Part) -> Part| {
239        let mut exact = true;
240        let mut kept = Vec::with_capacity(parts.len());
241        for p in parts {
242            let (part, known) = under(p)?;
243            exact &= known;
244            kept.push(wrap(part));
245        }
246        Some((kept, exact))
247    };
248    match f {
249        Body::Series(s) => match total(s, band) {
250            Some(sum) => held(Body::Const(C64::real(sum)), true),
251            None => held(Body::Const(C64::ONE), false),
252        },
253        Body::Line | Body::Node(_) => held(Body::Const(C64::ONE), false),
254        Body::Mul(parts) => {
255            let (kept, exact) = all(parts, &|p| p)?;
256            held(Body::Mul(kept), exact)
257        }
258        Body::Add(parts) => {
259            let (kept, exact) = all(parts, &|p| Part::bare(Body::Apply(Unary::Abs, p)))?;
260            held(Body::Add(kept), exact)
261        }
262        Body::Div(a, b) if !mentions_line(&b.body) => {
263            let (num, exact) = under(a)?;
264            held(Body::Div(num, b.clone()), exact)
265        }
266        Body::Div(a, b) => {
267            let ((num, _), (den, _)) = (under(a)?, under(b)?);
268            held(Body::Div(num, den), false)
269        }
270        Body::Apply(Unary::Sin | Unary::Cos | Unary::Tanh | Unary::Sat, arg) => {
271            held(Body::Const(C64::ONE), real(&arg.body))
272        }
273        Body::Crop { of, .. } | Body::Shift { of, .. } | Body::Warp { of, .. } => {
274            bound(&of.body, band)
275        }
276        Body::Pow(base, n) => {
277            let (part, exact) = under(base)?;
278            held(Body::Pow(part, *n), exact && *n >= 0)
279        }
280        other if !mentions_line(other) => held(other.clone(), true),
281        _ => None,
282    }
283}
284
285/// A geometric series sums to its first bound over `1 - ratio`; any other to what it keeps.
286fn total(s: &Series, band: Audible) -> Option<f64> {
287    let bound = bound(&s.term.body, band).filter(|b| b.exact)?.body;
288    let at = |i: i64| Some(exact_constant(&substitute(&bound, s.index, i as f64))?.abs());
289    if let (Bound::Infinite, Some(r)) = (s.hi, ratio(&bound, s.index).filter(|r| *r < 1.0)) {
290        return Some(at(s.lo)? / (1.0 - r));
291    }
292    let kept = i64::try_from(counted(s, band).filter(|n| *n <= MAX_EXPANDED_TERMS)?).ok()?;
293    (s.lo..s.lo + kept).try_fold(0.0, |held, i| Some(held + at(i)?))
294}
295
296fn real(f: &Body) -> bool {
297    axis(f, &Unread) == Axis::Real
298}
299
300struct Unread;
301
302impl Env for Unread {
303    fn node(&self, _: NodeId) -> Ty {
304        Ty::form(Var::T, false, Codomain::Complex)
305    }
306
307    fn param(&self, _: ParamId) -> Ty {
308        Ty::form(Var::T, false, Codomain::Complex)
309    }
310}