fastsim_core/utils/
interp.rs

1use crate::imports::*;
2
3pub(crate) trait InterpolatorScanValues {
4    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, f: F) -> Result<(), E>;
5}
6
7impl InterpolatorScanValues for Interp0D<f64> {
8    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, mut f: F) -> Result<(), E> {
9        f(self.0)
10    }
11}
12
13impl InterpolatorScanValues for Interp1D<f64, strategy::enums::Strategy1DEnum<f64>> {
14    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, mut f: F) -> Result<(), E> {
15        self.data.values.iter().copied().try_for_each(&mut f)
16    }
17}
18
19impl InterpolatorScanValues for Interp2D<f64, strategy::enums::Strategy2DEnum<f64>> {
20    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, mut f: F) -> Result<(), E> {
21        self.data.values.iter().copied().try_for_each(&mut f)
22    }
23}
24
25impl InterpolatorScanValues for Interp3D<f64, strategy::enums::Strategy3DEnum<f64>> {
26    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, mut f: F) -> Result<(), E> {
27        self.data.values.iter().copied().try_for_each(&mut f)
28    }
29}
30
31impl InterpolatorScanValues for InterpND<f64, strategy::enums::StrategyNDEnum<f64>> {
32    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, mut f: F) -> Result<(), E> {
33        self.data.values.iter().copied().try_for_each(&mut f)
34    }
35}
36
37impl InterpolatorScanValues for InterpolatorEnum<f64> {
38    fn try_for_each_value<E, F: FnMut(f64) -> Result<(), E>>(&self, f: F) -> Result<(), E> {
39        match self {
40            Self::Interp0D(interp) => interp.try_for_each_value(f),
41            Self::Interp1D(interp) => interp.try_for_each_value(f),
42            Self::Interp2D(interp) => interp.try_for_each_value(f),
43            Self::Interp3D(interp) => interp.try_for_each_value(f),
44            Self::InterpND(interp) => interp.try_for_each_value(f),
45        }
46    }
47}
48
49/// Methods for mutating interpolator data, e.g. proportionally scaling
50/// interpolator function data
51pub trait InterpolatorMutMethods {
52    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()>;
53    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()>;
54    fn set_range(&mut self, range: f64) -> anyhow::Result<()>;
55}
56
57impl InterpolatorMutMethods for Interp0D<f64> {
58    fn set_min(&mut self, min: f64, _scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
59        self.0 = min;
60        Ok(())
61    }
62
63    fn set_max(&mut self, max: f64, _scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
64        self.0 = max;
65        Ok(())
66    }
67
68    fn set_range(&mut self, _range: f64) -> anyhow::Result<()> {
69        bail!("Cannot set range for 0D interpolator")
70    }
71}
72
73impl<S> InterpolatorMutMethods for Interp1D<f64, S>
74where
75    S: ninterp::strategy::traits::Strategy1D<ndarray::OwnedRepr<f64>> + Clone,
76{
77    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
78        let old_min = *self.min()?;
79        match scaling.unwrap_or_default() {
80            utils::interp::ScalingMethods::Proportional => {
81                ensure!(
82                    old_min != 0.,
83                    "Cannot modify min proportionally when old_min == 0."
84                );
85                self.data.values.map_inplace(|v| *v *= min / old_min);
86            }
87            utils::interp::ScalingMethods::AnchoredProportional => {
88                todo!()
89            }
90            utils::interp::ScalingMethods::Offset => {
91                todo!()
92            }
93        }
94        self.validate()?;
95        Ok(())
96    }
97
98    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
99        let old_max = *self.max()?;
100        match scaling.unwrap_or_default() {
101            utils::interp::ScalingMethods::Proportional => {
102                ensure!(
103                    old_max != 0.,
104                    "Cannot modify max proportionally when old_max == 0."
105                );
106                self.data.values.map_inplace(|v| *v *= max / old_max);
107            }
108            utils::interp::ScalingMethods::AnchoredProportional => {
109                todo!()
110            }
111            utils::interp::ScalingMethods::Offset => {
112                todo!()
113            }
114        }
115        self.validate()?;
116        Ok(())
117    }
118
119    fn set_range(&mut self, range: f64) -> anyhow::Result<()> {
120        let old_max = *self.max()?;
121        let old_range = old_max - self.min()?;
122        ensure!(old_range != 0., "Cannot modify range when min == max");
123        // if the new range is 0., chooses the max as the value for all elements of the array
124        if range == 0. {
125            self.data.values = self.data.values.map(|_| old_max);
126        } else {
127            self.data.values = self
128                .data
129                .values
130                .map(|x| old_max + (x - old_max) * range / old_range);
131        }
132        self.validate()?;
133        Ok(())
134    }
135}
136
137impl<S> InterpolatorMutMethods for Interp2D<f64, S>
138where
139    S: ninterp::strategy::traits::Strategy2D<ndarray::OwnedRepr<f64>> + Clone,
140{
141    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
142        let old_min = *self.min()?;
143        match scaling.unwrap_or_default() {
144            utils::interp::ScalingMethods::Proportional => {
145                ensure!(
146                    old_min != 0.,
147                    "Cannot modify min proportionally when old_min == 0."
148                );
149                self.data.values.map_inplace(|v| *v *= min / old_min);
150            }
151            utils::interp::ScalingMethods::AnchoredProportional => {
152                todo!()
153            }
154            utils::interp::ScalingMethods::Offset => {
155                todo!()
156            }
157        }
158        self.validate()?;
159        Ok(())
160    }
161
162    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
163        let old_max = *self.max()?;
164        match scaling.unwrap_or_default() {
165            utils::interp::ScalingMethods::Proportional => {
166                ensure!(
167                    old_max != 0.,
168                    "Cannot modify max proportionally when old_max == 0."
169                );
170                self.data.values.map_inplace(|v| *v *= max / old_max);
171            }
172            utils::interp::ScalingMethods::AnchoredProportional => {
173                todo!()
174            }
175            utils::interp::ScalingMethods::Offset => {
176                todo!()
177            }
178        }
179        self.validate()?;
180        Ok(())
181    }
182
183    fn set_range(&mut self, range: f64) -> anyhow::Result<()> {
184        let old_max = *self.max()?;
185        let old_range = old_max - self.min()?;
186        ensure!(old_range != 0., "Cannot modify range when min == max");
187        // if the new range is 0., chooses the max as the value for all elements of the array
188        if range == 0. {
189            self.data.values = self.data.values.map(|_| old_max);
190        } else {
191            self.data.values = self
192                .data
193                .values
194                .map(|x| old_max + (x - old_max) * range / old_range);
195        }
196        self.validate()?;
197        Ok(())
198    }
199}
200
201impl<S> InterpolatorMutMethods for Interp3D<f64, S>
202where
203    S: ninterp::strategy::traits::Strategy3D<ndarray::OwnedRepr<f64>> + Clone,
204{
205    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
206        let old_min = *self.min()?;
207        match scaling.unwrap_or_default() {
208            utils::interp::ScalingMethods::Proportional => {
209                ensure!(
210                    old_min != 0.,
211                    "Cannot modify min proportionally when old_min == 0."
212                );
213                self.data.values.map_inplace(|v| *v *= min / old_min);
214            }
215            utils::interp::ScalingMethods::AnchoredProportional => {
216                todo!()
217            }
218            utils::interp::ScalingMethods::Offset => {
219                todo!()
220            }
221        }
222        self.validate()?;
223        Ok(())
224    }
225
226    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
227        let old_max = *self.max()?;
228        match scaling.unwrap_or_default() {
229            utils::interp::ScalingMethods::Proportional => {
230                ensure!(
231                    old_max != 0.,
232                    "Cannot modify max proportionally when old_max == 0."
233                );
234                self.data.values.map_inplace(|v| *v *= max / old_max);
235            }
236            utils::interp::ScalingMethods::AnchoredProportional => {
237                todo!()
238            }
239            utils::interp::ScalingMethods::Offset => {
240                todo!()
241            }
242        }
243        self.validate()?;
244        Ok(())
245    }
246
247    fn set_range(&mut self, range: f64) -> anyhow::Result<()> {
248        let old_max = *self.max()?;
249        let old_range = old_max - self.min()?;
250        ensure!(old_range != 0., "Cannot modify range when min == max");
251        // if the new range is 0., chooses the max as the value for all elements of the array
252        if range == 0. {
253            self.data.values = self.data.values.map(|_| old_max);
254        } else {
255            self.data.values = self
256                .data
257                .values
258                .map(|x| old_max + (x - old_max) * range / old_range);
259        }
260        self.validate()?;
261        Ok(())
262    }
263}
264
265impl<S> InterpolatorMutMethods for InterpND<f64, S>
266where
267    S: ninterp::strategy::traits::StrategyND<ndarray::OwnedRepr<f64>> + Clone,
268{
269    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
270        let old_min = *self.min()?;
271        match scaling.unwrap_or_default() {
272            utils::interp::ScalingMethods::Proportional => {
273                ensure!(
274                    old_min != 0.,
275                    "Cannot modify min proportionally when old_min == 0."
276                );
277                self.data.values.map_inplace(|v| *v *= min / old_min);
278            }
279            utils::interp::ScalingMethods::AnchoredProportional => {
280                todo!()
281            }
282            utils::interp::ScalingMethods::Offset => {
283                todo!()
284            }
285        }
286        self.validate()?;
287        Ok(())
288    }
289
290    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
291        let old_max = *self.max()?;
292        match scaling.unwrap_or_default() {
293            utils::interp::ScalingMethods::Proportional => {
294                ensure!(
295                    old_max != 0.,
296                    "Cannot modify max proportionally when old_max == 0."
297                );
298                self.data.values.map_inplace(|v| *v *= max / old_max);
299            }
300            utils::interp::ScalingMethods::AnchoredProportional => {
301                todo!()
302            }
303            utils::interp::ScalingMethods::Offset => {
304                todo!()
305            }
306        }
307        self.validate()?;
308        Ok(())
309    }
310
311    fn set_range(&mut self, range: f64) -> anyhow::Result<()> {
312        let old_max = *self.max()?;
313        let old_range = old_max - self.min()?;
314        ensure!(old_range != 0., "Cannot modify range when min == max");
315        // if the new range is 0., chooses the max as the value for all elements of the array
316        if range == 0. {
317            self.data.values = self.data.values.map(|_| old_max);
318        } else {
319            self.data.values = self
320                .data
321                .values
322                .map(|x| old_max + (x - old_max) * range / old_range);
323        }
324        self.validate()?;
325        Ok(())
326    }
327}
328
329// This can be made more generic by using a `ninterp::num_traits` bound instead of f64
330// If there are future methods that *do not* mutate the interpolator,
331// we should define a new trait and impl it for `InterpolatorEnum<D> where D: ndarray::Data`
332impl InterpolatorMutMethods for InterpolatorEnum<f64> {
333    // scale all values so that the min is the new min
334    // (Note: this may change the max, depending on what scaling method is chosen)
335    fn set_min(&mut self, min: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
336        match self {
337            Self::Interp0D(interp) => interp.set_min(min, scaling),
338            Self::Interp1D(interp) => interp.set_min(min, scaling),
339            Self::Interp2D(interp) => interp.set_min(min, scaling),
340            Self::Interp3D(interp) => interp.set_min(min, scaling),
341            Self::InterpND(interp) => interp.set_min(min, scaling),
342        }
343    }
344
345    // scale all values so that the max is the new max
346    // (Note: may change the min, depending on what scaling method is chosen)
347    fn set_max(&mut self, max: f64, scaling: Option<ScalingMethods>) -> anyhow::Result<()> {
348        match self {
349            Self::Interp0D(interp) => interp.set_max(max, scaling),
350            Self::Interp1D(interp) => interp.set_max(max, scaling),
351            Self::Interp2D(interp) => interp.set_max(max, scaling),
352            Self::Interp3D(interp) => interp.set_max(max, scaling),
353            Self::InterpND(interp) => interp.set_max(max, scaling),
354        }
355    }
356
357    fn set_range(&mut self, range: f64) -> anyhow::Result<()> {
358        match self {
359            Self::Interp0D(interp) => interp.set_range(range),
360            Self::Interp1D(interp) => interp.set_range(range),
361            Self::Interp2D(interp) => interp.set_range(range),
362            Self::Interp3D(interp) => interp.set_range(range),
363            Self::InterpND(interp) => interp.set_range(range),
364        }
365    }
366}
367
368impl<D> Init for InterpolatorEnumBase<D>
369where
370    D: ndarray::Data + ndarray::RawDataClone + Clone,
371    D::Elem: ninterp::num_traits::Float + ninterp::num_traits::Euclid + std::fmt::Debug,
372{
373    fn init(&mut self) -> Result<(), Error> {
374        self.validate()
375            .map_err(|e| Error::NinterpError(e.to_string()))
376    }
377}
378impl<D> SerdeAPI for InterpolatorEnumBase<D>
379where
380    D: ndarray::Data + ndarray::RawDataClone + Clone + ndarray::DataOwned,
381    D::Elem: ninterp::num_traits::Float
382        + ninterp::num_traits::Euclid
383        + std::fmt::Debug
384        + Serialize
385        + serde::de::DeserializeOwned,
386{
387    #[cfg(feature = "resources")]
388    const RESOURCES_SUBDIR: &'static str = "interpolators";
389}
390
391#[derive(Default, Debug, Serialize, Deserialize, Clone)]
392pub enum ScalingMethods {
393    #[default]
394    /// Scales everything by the same factor -- e.g. setting min of [1, 2, 3] to 0.5 yields [0.5, 1, 1.5]
395    Proportional,
396    /// Scales proportionally to distance from min/max -- e.g. setting min of [1, 2, 3] to 0.5 yields [0.5, 1.5, 3]
397    AnchoredProportional,
398    /// Scaling by sliding all values up or down by the same offset -- e.g. setting min of [1, 2, 3] to 0.5 yields [0.5, 1.5, 2.5]
399    Offset,
400}
401
402#[cfg(test)]
403mod tests {
404    use super::*;
405
406    #[test]
407    fn test_min() {
408        let x = array![0.05, 0.10, 0.15];
409        let y = array![0.10, 0.20, 0.30];
410        let z = array![0.20, 0.40, 0.60];
411        let f_xy = array![[0.1, 1., 2.], [3., 4., 5.], [6., 7., 8.]];
412        let f_xyz = array![
413            [[0.1, 1., 2.], [3., 4., 5.], [6., 7., 8.]],
414            [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]],
415            [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.],],
416        ];
417        let mut interp_1d = InterpolatorEnum::new_1d(
418            array![85.0, 90.0],
419            array![0.2, 1.0],
420            strategy::Linear,
421            Extrapolate::Clamp,
422        )
423        .unwrap();
424        let mut interp_2d = InterpolatorEnum::new_2d(
425            x.clone(),
426            y.clone(),
427            f_xy.clone(),
428            strategy::Linear,
429            Extrapolate::Clamp,
430        )
431        .unwrap();
432        let mut interp_3d = InterpolatorEnum::new_3d(
433            x.clone(),
434            y.clone(),
435            z.clone(),
436            f_xyz.clone(),
437            strategy::Linear,
438            Extrapolate::Clamp,
439        )
440        .unwrap();
441        assert_eq!(interp_1d.min().unwrap(), &0.2);
442        assert_eq!(interp_2d.min().unwrap(), &0.1);
443        assert_eq!(interp_3d.min().unwrap(), &0.1);
444        interp_1d.set_min(0.1, None).unwrap();
445        interp_2d.set_min(0.3, None).unwrap();
446        interp_3d.set_min(0.3, None).unwrap();
447        println!("{:?}", interp_1d.min().unwrap());
448        println!("{:?}", interp_2d.min().unwrap());
449        println!("{:?}", interp_3d.min().unwrap());
450        assert!(almost_eq(*interp_1d.min().unwrap(), 0.1, Some(1e-3)));
451        assert!(almost_eq(*interp_2d.min().unwrap(), 0.3, Some(1e-3)));
452        assert!(almost_eq(*interp_3d.min().unwrap(), 0.3, Some(1e-3)));
453    }
454
455    #[test]
456    fn test_max() {
457        let x = array![0.05, 0.10, 0.15];
458        let y = array![0.10, 0.20, 0.30];
459        let z = array![0.20, 0.40, 0.60];
460        let f_xy = array![[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]];
461        let f_xyz = array![
462            [[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]],
463            [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]],
464            [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.]],
465        ];
466        let mut interp_1d = InterpolatorEnum::new_1d(
467            array![85.0, 90.0],
468            array![0.2, 1.0],
469            strategy::Linear,
470            Extrapolate::Clamp,
471        )
472        .unwrap();
473        let mut interp_2d = InterpolatorEnum::new_2d(
474            x.clone(),
475            y.clone(),
476            f_xy.clone(),
477            strategy::Linear,
478            Extrapolate::Clamp,
479        )
480        .unwrap();
481        let mut interp_3d = InterpolatorEnum::new_3d(
482            x.clone(),
483            y.clone(),
484            z.clone(),
485            f_xyz.clone(),
486            strategy::Linear,
487            Extrapolate::Clamp,
488        )
489        .unwrap();
490        assert_eq!(interp_1d.max().unwrap(), &1.0);
491        assert_eq!(interp_2d.max().unwrap(), &8.);
492        assert_eq!(interp_3d.max().unwrap(), &26.);
493        interp_1d.set_max(2., None).unwrap();
494        interp_2d.set_max(7., None).unwrap();
495        interp_3d.set_max(5., None).unwrap();
496        println!("{:?}", interp_1d.max().unwrap());
497        println!("{:?}", interp_2d.max().unwrap());
498        println!("{:?}", interp_3d.max().unwrap());
499        assert!(almost_eq(*interp_1d.max().unwrap(), 2., Some(1e-3)));
500        assert!(almost_eq(*interp_2d.max().unwrap(), 7., Some(1e-3)));
501        assert!(almost_eq(*interp_3d.max().unwrap(), 5., Some(1e-3)));
502    }
503
504    #[test]
505    fn test_range() {
506        let x = array![0.05, 0.10, 0.15];
507        let y = array![0.10, 0.20, 0.30];
508        let z = array![0.20, 0.40, 0.60];
509        let f_xy = array![[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]];
510        let f_xyz = array![
511            [[0., 1., 2.], [3., 4., 5.], [6., 7., 8.]],
512            [[9., 10., 11.], [12., 13., 14.], [15., 16., 17.]],
513            [[18., 19., 20.], [21., 22., 23.], [24., 25., 26.]],
514        ];
515        let mut interp_1d = InterpolatorEnum::new_1d(
516            array![85.0, 90.0],
517            array![0.2, 1.0],
518            strategy::Linear,
519            Extrapolate::Clamp,
520        )
521        .unwrap();
522        let mut interp_2d = InterpolatorEnum::new_2d(
523            x.clone(),
524            y.clone(),
525            f_xy.clone(),
526            strategy::Linear,
527            Extrapolate::Clamp,
528        )
529        .unwrap();
530        let mut interp_3d = InterpolatorEnum::new_3d(
531            x.clone(),
532            y.clone(),
533            z.clone(),
534            f_xyz.clone(),
535            strategy::Linear,
536            Extrapolate::Clamp,
537        )
538        .unwrap();
539        assert_eq!(interp_1d.range().unwrap(), 0.8);
540        assert_eq!(interp_2d.range().unwrap(), 8.);
541        assert_eq!(interp_3d.range().unwrap(), 26.);
542        interp_1d.set_range(2.).unwrap();
543        interp_2d.set_range(7.).unwrap();
544        interp_3d.set_range(5.).unwrap();
545        println!("{:?}", interp_1d.range().unwrap());
546        println!("{:?}", interp_2d.range().unwrap());
547        println!("{:?}", interp_3d.range().unwrap());
548        assert!(almost_eq(interp_1d.range().unwrap(), 2., Some(1e-3)));
549        assert!(almost_eq(interp_2d.range().unwrap(), 7., Some(1e-3)));
550        assert!(almost_eq(interp_3d.range().unwrap(), 5., Some(1e-3)));
551    }
552
553    type StructWithResources = InterpolatorEnum<f64>;
554
555    #[test]
556    fn test_resources() {
557        let resource_list = StructWithResources::list_resources().unwrap();
558        assert!(!resource_list.is_empty());
559
560        // verify that resources can all load
561        for resource in resource_list {
562            StructWithResources::from_resource(resource.clone(), false)
563                .with_context(|| format_dbg!(resource))
564                .unwrap();
565        }
566    }
567}
fastsim_core/utils/interp.rs

fastsim_core/utils/
interp.rs
