rbf_interpolation/
builder.rs

1use std::iter::{repeat, Product};
2
3use na::{
4    ArrayStorage, Const, DefaultAllocator, DimMin, DimName, Matrix, RealField, Scalar,
5    SquareMatrix, Storage, StorageMut, U1, Vector, allocator::Allocator,
6};
7
8use crate::{powers::monomial_exponents, rbf::RBFInterpolator};
9
10/// You must provide the number of monomials terms
11/// for your choice of polynomial degree and points
12/// dimension.
13///
14/// i.e. A degree 1 polynomial in two variables will
15/// incur the following monomials: 1, x, y
16/// Degree 2 will incur: 1, x, y, xy, x^2, y^2
17///
18/// In general, monomials = (degree + dimension) choose (degree)
19pub enum RBFInterpolatorBuilder<
20    T,
21    const DEGREE: usize,
22    const MONOMIALS: usize,
23    const POINTS: usize,
24    const DIM: usize,
25> {
26    /// r
27    Linear,
28    /// r^2 * log(r)
29    ThinPlateSpline,
30    /// r^3
31    Cubic,
32    /// r^5
33    Quintic,
34    /// -sqrt(1 + r^2)
35    Multiquadratic { epsilon: T },
36    /// 1/sqrt(1 + r^2)
37    InverseMultiquadratic { epsilon: T },
38    /// 1/(1 + r^2)
39    InverseQuadratic { epsilon: T },
40    /// exp(-r^2)
41    Gaussian { epsilon: T },
42}
43
44impl<T, const DEGREE: usize, const MONOMIALS: usize, const POINTS: usize, const DIM: usize>
45    RBFInterpolatorBuilder<T, DEGREE, MONOMIALS, POINTS, DIM>
46where
47    T: Scalar + RealField + Copy + Product,
48{
49    /// Constructs NxN matrix of phi(||a - b||) for each combination of points a, b
50    fn construct_phi<SP>(
51        &self,
52        points: &Matrix<T, Const<DIM>, Const<POINTS>, SP>,
53    ) -> SquareMatrix<
54        T,
55        Const<POINTS>,
56        <DefaultAllocator as Allocator<Const<POINTS>, Const<POINTS>>>::Buffer<T>,
57    >
58    where
59        SP: Storage<T, Const<DIM>, Const<POINTS>>,
60    {
61        let iter = points
62            .column_iter()
63            .flat_map(|point_a| points.column_iter().zip(repeat(point_a)))
64            .map(|(point_a, point_b)| {
65                self.kernel((point_a - point_b).map(|e| e.powi(2)).sum().sqrt())
66            });
67
68        SquareMatrix::<
69            T,
70            Const<POINTS>,
71            <DefaultAllocator as Allocator<Const<POINTS>, Const<POINTS>>>::Buffer<T>,
72        >::from_iterator(iter)
73    }
74
75    /// Add polynomial constraints to phi matrix (POINTS + MONOMIALS)x(POINTS + MONOMIALS)
76    fn embelish_phi<SP, SH>(
77        &self,
78        points: &Matrix<T, Const<DIM>, Const<POINTS>, SP>,
79        phi: SquareMatrix<T, Const<POINTS>, SH>,
80    ) -> SquareMatrix<
81        T,
82        Const<{ POINTS + MONOMIALS }>,
83        <DefaultAllocator as Allocator<
84            Const<{ POINTS + MONOMIALS }>,
85            Const<{ POINTS + MONOMIALS }>,
86        >>::Buffer<T>,
87    >
88    where
89        SP: Storage<T, Const<DIM>, Const<POINTS>>,
90        SH: StorageMut<T, Const<POINTS>, Const<POINTS>>,
91        Const<{ POINTS + MONOMIALS }>: DimName,
92        DefaultAllocator: Allocator<Const<POINTS>, Const<POINTS>>
93            + Allocator<Const<{ POINTS + MONOMIALS }>, Const<{ POINTS + MONOMIALS }>>,
94    {
95        let mut phi = phi.resize_generic(
96            Const::<{ POINTS + MONOMIALS }>,
97            Const::<{ POINTS + MONOMIALS }>,
98            na::zero(),
99        );
100        let exponents = monomial_exponents::<DIM, DEGREE>();
101
102        assert_eq!(
103            exponents.len(),
104            MONOMIALS,
105            "The choices of generics polynomial DEGREE and number of variables DIM did not match the number of MONONOMIALS."
106        );
107
108        for (i, point) in points.column_iter().enumerate() {
109            for (j, exponent) in exponents.iter().enumerate() {
110                debug_assert_eq!(exponent.len(), point.len());
111                let value: T = exponent
112                    .iter()
113                    .zip(point.row_iter())
114                    .map(|(&exponent, ordinate)| ordinate[(0, 0)].powi(exponent))
115                    .product();
116
117                phi[(i, POINTS + j)] = value;
118                phi[(POINTS + j, i)] = value;
119            }
120        }
121
122        phi
123    }
124
125    /// Add polynomial constraints to the values vector (POINTS + MONOMIALS)
126    fn embelish_values<SV>(
127        &self,
128        values: Vector<T, Const<POINTS>, SV>,
129    ) -> Vector<
130        T,
131        Const<{ POINTS + MONOMIALS }>,
132        <DefaultAllocator as Allocator<Const<{ POINTS + MONOMIALS }>, U1>>::Buffer<T>,
133    >
134    where
135        SV: StorageMut<T, Const<POINTS>, U1>,
136        Const<{ POINTS + MONOMIALS }>: DimName,
137    {
138        values.resize_generic(Const::<{ POINTS + MONOMIALS }>, U1, na::zero())
139    }
140
141    /// Solves (phi) * (weights) = (values), storing the solved weights into values
142    /// Introduces new generic M to simplify the traits
143    fn solve_phi<const M: usize, SH, SV>(
144        &self,
145        phi: SquareMatrix<T, Const<M>, SH>,
146        values: &mut Vector<T, Const<M>, SV>,
147    ) -> bool
148    where
149        SH: Storage<T, Const<M>, Const<M>>,
150        SV: StorageMut<T, Const<M>, U1>,
151        Const<M>: DimMin<Const<M>, Output = Const<M>>,
152        DefaultAllocator: Allocator<<Const<M> as DimMin<Const<M>>>::Output>,
153    {
154        phi.lu().solve_mut(values)
155    }
156
157    /// Number of points with dimension provided as a matrix of
158    /// collum vectors, with their values in a seperate vector.
159    ///
160    /// Returns None when the linear system was not solveable.
161    /// Will panic when the choice of added polynomial DEGREE
162    /// and number of corresponding MONOMIAL terms are incompatible.
163    /// Should satisfy: MONOMIAL = (DIM+DEGREE) choose DEGREE.
164    pub fn build<SP, SV>(
165        self,
166        points: Matrix<T, Const<DIM>, Const<POINTS>, SP>,
167        values: Vector<T, Const<POINTS>, SV>,
168    ) -> Option<
169        RBFInterpolator<
170            T,
171            DEGREE,
172            MONOMIALS,
173            POINTS,
174            DIM,
175            SP,
176            ArrayStorage<T, { POINTS + MONOMIALS }, 1>,
177        >,
178    >
179    where
180        SP: Storage<T, Const<DIM>, Const<POINTS>>,
181        SV: StorageMut<T, Const<POINTS>, U1>,
182        //Convince compiler that phi is infact square
183        Const<{ POINTS + MONOMIALS }>:
184            DimMin<Const<{ POINTS + MONOMIALS }>, Output = Const<{ POINTS + MONOMIALS }>>,
185    {
186        let phi = self.construct_phi::<SP>(&points);
187        let phi = self.embelish_phi::<SP, ArrayStorage<T, POINTS, POINTS>>(&points, phi);
188        let mut values = self.embelish_values(values);
189        if self.solve_phi(phi, &mut values) {
190            Some(
191                RBFInterpolator::<T, _, _, _, _, SP, ArrayStorage<T, { POINTS + MONOMIALS }, 1>> {
192                    kernel: self,
193                    points,
194                    weights: values,
195                },
196            )
197        } else {
198            None
199        }
200    }
201}
rbf_interpolation/builder.rs

rbf_interpolation/
builder.rs
