rsl_interpolation/spline2d.rs
1use crate::Accelerator;
2use crate::DynInterp2dType;
3use crate::Interp2dType;
4use crate::Interpolation2d;
5use crate::{DomainError, InterpolationError};
6
7/// 2D Higher level interface.
8///
9/// A 2D Spline owns the data it is constructed with, and provides the same evaluation methods as the
10/// lower-level Interpolator object, without needing to provide the data arrays in every call.
11///
12/// # Example
13/// ```
14/// # use rsl_interpolation::*;
15/// #
16/// # fn main() -> Result<(), InterpolationError>{
17/// let mut xacc = Accelerator::new();
18/// let mut yacc = Accelerator::new();
19///
20/// let xa = [0.0, 1.0, 2.0, 3.0];
21/// let ya = [0.0, 2.0, 4.0, 6.0];
22/// // z = x + y, in column-major order
23/// let za = [
24/// 0.0, 1.0, 2.0, 3.0,
25/// 2.0, 3.0, 4.0, 5.0,
26/// 4.0, 5.0, 6.0, 7.0,
27/// 6.0, 7.0, 8.0, 9.0,
28/// ];
29///
30/// let interp = Bicubic.build(&xa, &ya, &za)?;
31///
32/// let typ = Bicubic;
33/// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
34///
35/// let (x, y) = (2.5, 4.1);
36/// let y_interp = interp.eval(&xa, &ya, &za, x, y, &mut xacc, &mut yacc)?;
37/// let y_spline = spline.eval(x, y, &mut xacc, &mut yacc)?;
38///
39/// assert_eq!(y_interp, y_spline);
40/// #
41/// # Ok(())
42/// # }
43/// ```
44pub struct Spline2d<I, T>
45where
46 I: Interp2dType<T> + Send + Sync + 'static,
47{
48 /// The lower-level [`2D Interpolator`].
49 ///
50 /// [`2D Interpolator`]: Interpolation2d#implementors
51 pub interp: I::Interpolation2d,
52 /// The owned x data.
53 pub xa: Box<[T]>,
54 /// The owned y data.
55 pub ya: Box<[T]>,
56 /// The owned z data.
57 pub za: Box<[T]>,
58 name: Box<str>,
59 min_size: usize,
60}
61
62impl<I, T> Spline2d<I, T>
63where
64 I: Interp2dType<T> + Send + Sync + 'static,
65{
66 /// Constructs a 2D Spline of a 2D Interpolation type `typ` from the data arrays `xa`, `ya` and
67 /// `za`.
68 ///
69 /// # Example
70 /// ```
71 /// # use rsl_interpolation::*;
72 /// #
73 /// # fn main() -> Result<(), InterpolationError>{
74 /// let mut xacc = Accelerator::new();
75 /// let mut yacc = Accelerator::new();
76 ///
77 /// let xa = [0.0, 1.0, 2.0, 3.0];
78 /// let ya = [0.0, 2.0, 4.0, 6.0];
79 /// // z = x + y, in column-major order
80 /// let za = [
81 /// 0.0, 1.0, 2.0, 3.0,
82 /// 2.0, 3.0, 4.0, 5.0,
83 /// 4.0, 5.0, 6.0, 7.0,
84 /// 6.0, 7.0, 8.0, 9.0,
85 /// ];
86 ///
87 /// let typ = Bicubic;
88 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
89 /// #
90 /// # Ok(())
91 /// # }
92 /// ```
93 #[doc(alias = "gsl_spline2d_init")]
94 pub fn new(typ: I, xa: &[T], ya: &[T], za: &[T]) -> Result<Self, InterpolationError>
95 where
96 T: Clone,
97 {
98 let spline = Self {
99 interp: typ.build(xa, ya, za)?,
100 xa: xa.into(),
101 ya: ya.into(),
102 za: za.into(),
103 name: typ.name().into(),
104 min_size: typ.min_size(),
105 };
106
107 Ok(spline)
108 }
109
110 /// Returns the interpolated value of `z` for a given point (`x`, `y`), using the
111 /// [`Accelerators`] `xacc` and `yacc`.
112 ///
113 /// # Example
114 ///
115 /// ```
116 /// # use rsl_interpolation::*;
117 /// #
118 /// # fn main() -> Result<(), InterpolationError>{
119 /// let mut xacc = Accelerator::new();
120 /// let mut yacc = Accelerator::new();
121 ///
122 /// let xa = [0.0, 1.0, 2.0];
123 /// let ya = [0.0, 2.0, 4.0];
124 /// // z = x + y
125 /// let za = [
126 /// 0.0, 1.0, 2.0,
127 /// 2.0, 3.0, 4.0,
128 /// 4.0, 5.0, 6.0,
129 /// ];
130 ///
131 /// let typ = Bilinear;
132 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
133 ///
134 /// let z = spline.eval(1.5, 3.0, &mut xacc, &mut yacc)?;
135 ///
136 /// assert_eq!(z, 4.5);
137 /// # Ok(())
138 /// # }
139 /// ```
140 ///
141 /// # Errors
142 ///
143 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
144 /// of `ya`.
145 ///
146 /// [`Accelerators`]: Accelerator
147 #[doc(alias = "gsl_spline2d_eval")]
148 #[doc(alias = "gsl_spline2d_eval_e")]
149 pub fn eval(
150 &self,
151 x: T,
152 y: T,
153 xacc: &mut Accelerator,
154 yacc: &mut Accelerator,
155 ) -> Result<T, DomainError>
156 where
157 T: crate::Num,
158 {
159 self.interp
160 .eval(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
161 }
162
163 /// Returns the interpolated value of `z` for a given point (`x`, `y`), using the
164 /// [`Accelerators`]` `xacc` and `yacc`.
165 ///
166 /// # Note
167 ///
168 /// This function performs *no bound checking*, so when `x` is outside the range of `xa` or y
169 /// is outside the range of `ya`, extrapolation is performed.
170 ///
171 /// # Example
172 ///
173 /// ```
174 /// # use rsl_interpolation::*;
175 /// #
176 /// # fn main() -> Result<(), InterpolationError>{
177 /// let mut xacc = Accelerator::new();
178 /// let mut yacc = Accelerator::new();
179 ///
180 /// let xa = [0.0, 1.0, 2.0];
181 /// let ya = [0.0, 2.0, 4.0];
182 /// // z = x + y
183 /// let za = [
184 /// 0.0, 1.0, 2.0,
185 /// 2.0, 3.0, 4.0,
186 /// 4.0, 5.0, 6.0,
187 /// ];
188 ///
189 /// let typ = Bilinear;
190 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
191 ///
192 /// let z = spline.eval_extrap(3.0, 6.0, &mut xacc, &mut yacc)?;
193 ///
194 /// assert_eq!(z, 9.0);
195 /// # Ok(())
196 /// # }
197 /// ```
198 ///
199 /// # Errors
200 ///
201 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
202 /// of `ya`.
203 ///
204 /// [`Accelerators`]: Accelerator
205 #[doc(alias = "gsl_spline2d_eval_extrap")]
206 #[doc(alias = "gsl_spline2d_eval_extrap_e")]
207 pub fn eval_extrap(
208 &self,
209 x: T,
210 y: T,
211 xacc: &mut Accelerator,
212 yacc: &mut Accelerator,
213 ) -> Result<T, DomainError> {
214 self.interp
215 .eval_extrap(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
216 }
217
218 /// Returns the interpolated value `d = ∂z/∂x` for a given point (`x`, `y`), using the
219 /// [`Accelerators`] `xacc` and `yacc`.
220 ///
221 /// # Example
222 ///
223 /// ```
224 /// # use rsl_interpolation::*;
225 /// #
226 /// # fn main() -> Result<(), InterpolationError>{
227 /// let mut xacc = Accelerator::new();
228 /// let mut yacc = Accelerator::new();
229 ///
230 /// let xa = [0.0, 1.0, 2.0];
231 /// let ya = [0.0, 2.0, 4.0];
232 /// // z = x² + y²
233 /// let za = [
234 /// 0.0, 1.0, 4.0,
235 /// 4.0, 5.0, 8.0,
236 /// 16.0, 17.0, 20.0,
237 /// ];
238 ///
239 /// let typ = Bilinear;
240 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
241 ///
242 /// let dzdx = spline.eval_deriv_x(1.5, 3.0, &mut xacc, &mut yacc)?;
243 ///
244 /// assert_eq!(dzdx, 3.0);
245 /// # Ok(())
246 /// # }
247 /// ```
248 ///
249 /// # Errors
250 ///
251 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
252 /// of `ya`.
253 ///
254 /// [`Accelerators`]: Accelerator
255 #[doc(alias = "gsl_spline2d_eval_deriv_x")]
256 #[doc(alias = "gsl_spline2d_eval_deriv_x_e")]
257 pub fn eval_deriv_x(
258 &self,
259 x: T,
260 y: T,
261 xacc: &mut Accelerator,
262 yacc: &mut Accelerator,
263 ) -> Result<T, DomainError> {
264 self.interp
265 .eval_deriv_x(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
266 }
267
268 /// Returns the interpolated value `d = ∂z/∂y` for a given point (`x`, `y`), using the
269 /// [`Accelerators`] `xacc` and `yacc`.
270 ///
271 /// # Example
272 ///
273 /// ```
274 /// # use rsl_interpolation::*;
275 /// #
276 /// # fn main() -> Result<(), InterpolationError>{
277 /// let mut xacc = Accelerator::new();
278 /// let mut yacc = Accelerator::new();
279 ///
280 /// let xa = [0.0, 1.0, 2.0];
281 /// let ya = [0.0, 2.0, 4.0];
282 /// // z = x² + y²
283 /// let za = [
284 /// 0.0, 1.0, 4.0,
285 /// 4.0, 5.0, 8.0,
286 /// 16.0, 17.0, 20.0,
287 /// ];
288 ///
289 /// let typ = Bilinear;
290 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
291 ///
292 /// let dzdx = spline.eval_deriv_y(1.5, 3.0, &mut xacc, &mut yacc)?;
293 ///
294 /// assert_eq!(dzdx, 6.0);
295 /// # Ok(())
296 /// # }
297 /// ```
298 ///
299 /// # Errors
300 ///
301 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
302 /// of `ya`.
303 ///
304 /// [`Accelerators`]: Accelerator
305 #[doc(alias = "gsl_spline2d_eval_deriv_y")]
306 #[doc(alias = "gsl_spline2d_eval_deriv_y_e")]
307 pub fn eval_deriv_y(
308 &self,
309 x: T,
310 y: T,
311 xacc: &mut Accelerator,
312 yacc: &mut Accelerator,
313 ) -> Result<T, DomainError> {
314 self.interp
315 .eval_deriv_y(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
316 }
317
318 /// Returns the interpolated value `d = 𝜕²z/𝜕x²` for a given point (`x`, `y`), using the
319 /// [`Accelerators`] `xacc` and `yacc`.
320 ///
321 /// # Example
322 ///
323 /// ```
324 /// # use rsl_interpolation::*;
325 /// #
326 /// # fn main() -> Result<(), InterpolationError>{
327 /// let mut xacc = Accelerator::new();
328 /// let mut yacc = Accelerator::new();
329 ///
330 /// let xa = [0.0, 1.0, 2.0];
331 /// let ya = [0.0, 2.0, 4.0];
332 /// // z = x² + y²
333 /// let za = [
334 /// 0.0, 1.0, 4.0,
335 /// 4.0, 5.0, 8.0,
336 /// 16.0, 17.0, 20.0,
337 /// ];
338 ///
339 /// let typ = Bilinear;
340 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
341 ///
342 /// let dzdx2 = spline.eval_deriv_xx(1.5, 3.0, &mut xacc, &mut yacc)?;
343 ///
344 /// assert_eq!(dzdx2, 0.0); // Linear Interpolation!
345 /// # Ok(())
346 /// # }
347 /// ```
348 ///
349 /// # Errors
350 ///
351 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
352 /// of `ya`.
353 ///
354 /// [`Accelerators`]: Accelerator
355 #[doc(alias = "gsl_interp2d_eval_deriv_xx")]
356 #[doc(alias = "gsl_interp2d_eval_deriv_xx_e")]
357 pub fn eval_deriv_xx(
358 &self,
359 x: T,
360 y: T,
361 xacc: &mut Accelerator,
362 yacc: &mut Accelerator,
363 ) -> Result<T, DomainError> {
364 self.interp
365 .eval_deriv_xx(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
366 }
367
368 /// Returns the interpolated value `d = 𝜕²z/𝜕x²` for a given point (`x`, `y`), using the
369 /// [`Accelerators`] `xacc` and `yacc`.
370 ///
371 /// # Example
372 ///
373 /// ```
374 /// # use rsl_interpolation::*;
375 /// #
376 /// # fn main() -> Result<(), InterpolationError>{
377 /// let mut xacc = Accelerator::new();
378 /// let mut yacc = Accelerator::new();
379 ///
380 /// let xa = [0.0, 1.0, 2.0];
381 /// let ya = [0.0, 2.0, 4.0];
382 /// // z = x² + y²
383 /// let za = [
384 ///
385 /// 0.0, 1.0, 4.0,
386 /// 4.0, 5.0, 8.0,
387 /// 16.0, 17.0, 20.0,
388 /// ];
389 /// let typ = Bilinear;
390 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
391 ///
392 /// let dzdx2 = spline.eval_deriv_yy(1.5, 3.0, &mut xacc, &mut yacc)?;
393 ///
394 /// assert_eq!(dzdx2, 0.0); // Linear Interpolation!
395 /// # Ok(())
396 /// # }
397 /// ```
398 ///
399 /// # Errors
400 ///
401 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
402 /// of `ya`.
403 ///
404 /// [`Accelerators`]: Accelerator
405 #[doc(alias = "gsl_interp2d_eval_deriv_yy")]
406 #[doc(alias = "gsl_interp2d_eval_deriv_yy_e")]
407 pub fn eval_deriv_yy(
408 &self,
409 x: T,
410 y: T,
411 xacc: &mut Accelerator,
412 yacc: &mut Accelerator,
413 ) -> Result<T, DomainError> {
414 self.interp
415 .eval_deriv_yy(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
416 }
417
418 /// Returns the interpolated value `d = 𝜕²z/𝜕x𝜕y` for a given point (`x`, `y`), using the
419 /// [`Accelerators`] `xacc` and `yacc`.
420 ///
421 /// # Example
422 ///
423 /// ```
424 /// # use rsl_interpolation::*;
425 /// #
426 /// # fn main() -> Result<(), InterpolationError>{
427 /// let mut xacc = Accelerator::new();
428 /// let mut yacc = Accelerator::new();
429 ///
430 /// let xa = [0.0, 1.0, 2.0];
431 /// let ya = [0.0, 2.0, 4.0];
432 /// // z = x² + y²
433 /// let za = [
434 /// 0.0, 1.0, 4.0,
435 /// 4.0, 5.0, 8.0,
436 /// 16.0, 17.0, 20.0,
437 /// ];
438 ///
439 /// let typ = Bilinear;
440 /// let spline = Spline2d::new(typ, &xa, &ya, &za)?;
441 ///
442 /// let dzdxy = spline.eval_deriv_xy(1.5, 3.0, &mut xacc, &mut yacc)?;
443 ///
444 /// assert_eq!(dzdxy, 0.0);
445 /// # Ok(())
446 /// # }
447 /// ```
448 ///
449 /// # Errors
450 ///
451 /// Returns a [`DomainError`] if `x` is outside the range of `xa` or `y` is outside the range
452 /// of `ya`.
453 ///
454 /// [`Accelerators`]: Accelerator
455 #[doc(alias = "gsl_interp2d_eval_deriv_xy")]
456 #[doc(alias = "gsl_interp2d_eval_deriv_xy_e")]
457 pub fn eval_deriv_xy(
458 &self,
459 x: T,
460 y: T,
461 xacc: &mut Accelerator,
462 yacc: &mut Accelerator,
463 ) -> Result<T, DomainError> {
464 self.interp
465 .eval_deriv_xy(&self.xa, &self.ya, &self.za, x, y, xacc, yacc)
466 }
467
468 /// Returns the name of the Interpolator.
469 #[doc(alias = "gsl_interp2d_name")]
470 pub fn name(&self) -> &str {
471 &self.name
472 }
473
474 /// Returns the minimum number of points required by the Interpolator.
475 #[doc(alias = "gsl_interp2d_min_size")]
476 pub fn min_size(&self) -> usize {
477 self.min_size
478 }
479}
480
481/// 2D Spline with runtime-determined Interpolation Type.
482pub type DynSpline2d<T> = Spline2d<DynInterp2dType<T>, T>;
483
484impl<T> DynSpline2d<T> {
485 /// Constructs a 2d Spline of a dynamic 2d Interpolation type `typ` from the data arrays `xa`,
486 /// `ya` and `za`.
487 ///
488 /// # Example
489 /// ```
490 /// # use rsl_interpolation::*;
491 /// #
492 /// # fn main() -> Result<(), InterpolationError> {
493 /// let xa = [0.0, 1.0, 2.0, 3.0];
494 /// let ya = [0.0, 2.0, 4.0, 6.0];
495 /// // z = x + y, in column-major order
496 /// let za = [
497 /// 0.0, 1.0, 2.0, 3.0,
498 /// 2.0, 3.0, 4.0, 5.0,
499 /// 4.0, 5.0, 6.0, 7.0,
500 /// 6.0, 7.0, 8.0, 9.0,
501 /// ];
502 /// let typ = "bicubic";
503 ///
504 /// let spline = match typ {
505 /// "bilinear" => Spline2d::new_dyn(Bilinear, &xa, &ya, &za)?,
506 /// "bicubic" => Spline2d::new_dyn(Bicubic, &xa, &ya, &za)?,
507 /// _ => unreachable!()
508 /// };
509 /// # Ok(())
510 /// # }
511 /// ```
512 #[doc(alias = "gsl_spline2d_init")]
513 pub fn new_dyn<I>(typ: I, xa: &[T], ya: &[T], za: &[T]) -> Result<Self, InterpolationError>
514 where
515 T: Clone,
516 I: Interp2dType<T> + Send + Sync + 'static,
517 I::Interpolation2d: Send + Sync + 'static,
518 {
519 Self::new(DynInterp2dType::new(typ), xa, ya, za)
520 }
521}
522
523/// Creates a [`DynSpline2d`] of `typ` type.
524///
525/// Useful when `typ` is not known at compile time.
526///
527/// # Example
528/// ```
529/// # use rsl_interpolation::*;
530/// #
531/// # fn main() -> Result<(), InterpolationError> {
532/// let xa = [0.0, 1.0, 2.0, 3.0];
533/// let ya = [0.0, 2.0, 4.0, 6.0];
534/// // z = x + y, in column-major order
535/// let za = [
536/// 0.0, 1.0, 2.0, 3.0,
537/// 2.0, 3.0, 4.0, 5.0,
538/// 4.0, 5.0, 6.0, 7.0,
539/// 6.0, 7.0, 8.0, 9.0,
540/// ];
541/// let typ = "bicubic";
542///
543/// let spline = make_spline2d(typ, &xa, &ya, &za)?;
544/// # Ok(())
545/// # }
546/// ```
547pub fn make_spline2d<T>(
548 typ: &str,
549 xa: &[T],
550 ya: &[T],
551 za: &[T],
552) -> Result<DynSpline2d<T>, InterpolationError>
553where
554 T: crate::Num + ndarray_linalg::Lapack,
555{
556 use crate::*;
557
558 match typ.to_lowercase().as_str() {
559 "bilinear" => Ok(DynSpline2d::new_dyn(Bilinear, xa, ya, za)?),
560 "bicubic" => Ok(DynSpline2d::new_dyn(Bicubic, xa, ya, za)?),
561 _ => Err(InterpolationError::InvalidType(typ.into())),
562 }
563}
564
565#[cfg(test)]
566mod test {
567 use crate::tests::build_comparator;
568 use crate::*;
569
570 #[test]
571 fn test_spline2d_creation() {
572 let xa = [0.0, 1.0, 2.0];
573 let ya = [0.0, 2.0, 4.0];
574 let za = [0.0, 1.0, 2.0, 2.0, 3.0, 4.0, 4.0, 5.0, 6.0];
575
576 let spline2d = Spline2d::new(Bilinear, &xa, &ya, &za).unwrap();
577 let _: &str = spline2d.name();
578 let _: usize = spline2d.min_size();
579 }
580
581 #[test]
582 fn test_spline2d_eval() {
583 let comp = build_comparator::<f64>();
584
585 let mut xacc = Accelerator::new();
586 let mut yacc = Accelerator::new();
587
588 let xa = [0.0, 1.0, 2.0, 3.0];
589 let ya = [0.0, 1.0, 2.0, 3.0];
590 #[rustfmt::skip]
591 let za = [
592 1.0, 1.1, 1.2, 1.3,
593 1.1, 1.2, 1.3, 1.4,
594 1.2, 1.3, 1.4, 1.5,
595 1.3, 1.4, 1.5, 1.6,
596 ];
597
598 let spline2d = Spline2d::new(Bicubic, &xa, &ya, &za).unwrap();
599
600 let (x, y) = (1.5, 1.5);
601 let z = spline2d.eval(x, y, &mut xacc, &mut yacc).unwrap();
602 let dzdx = spline2d.eval_deriv_x(x, y, &mut xacc, &mut yacc).unwrap();
603 let dzdy = spline2d.eval_deriv_y(x, y, &mut xacc, &mut yacc).unwrap();
604 let dzdx2 = spline2d.eval_deriv_xx(x, y, &mut xacc, &mut yacc).unwrap();
605 let dzdy2 = spline2d.eval_deriv_yy(x, y, &mut xacc, &mut yacc).unwrap();
606 let dzdxy = spline2d.eval_deriv_xy(x, y, &mut xacc, &mut yacc).unwrap();
607
608 assert!(comp.is_close(z, 1.3));
609 assert!(comp.is_close(dzdx, 0.1));
610 assert!(comp.is_close(dzdy, 0.1));
611 assert!(comp.is_close(dzdx2, 0.0));
612 assert!(comp.is_close(dzdy2, 0.0));
613 assert!(comp.is_close(dzdxy, 0.0));
614
615 let ze = spline2d
616 .eval_extrap(4.0, 4.0, &mut xacc, &mut yacc)
617 .unwrap();
618 assert!(comp.is_close(ze, 1.8));
619 }
620
621 #[test]
622 fn test_dyn_spline2d() {
623 let xa = [0.0, 1.0, 2.0, 3.0];
624 let ya = [0.0, 1.0, 2.0, 3.0];
625 #[rustfmt::skip]
626 let za = [
627 1.0, 1.1, 1.2, 1.3,
628 1.1, 1.2, 1.3, 1.4,
629 1.2, 1.3, 1.4, 1.5,
630 1.3, 1.4, 1.5, 1.6,
631 ];
632
633 Spline2d::new(Bicubic, &xa, &ya, &za).unwrap();
634 Spline2d::new_dyn(Bicubic, &xa, &ya, &za).unwrap();
635 }
636}