1use num::{Float, Zero};
40use std::iter::Sum;
41use std::ops::Mul;
42
43#[inline(always)]
45pub fn inner_product<T>(a: &[T], b: &[T]) -> T
46where
47 T: Float + Sum<T> + Mul<T, Output = T>,
48{
49 assert!(a.len() == b.len());
50
51 a.iter().zip(b.iter()).map(|(x, y)| (*x) * (*y)).sum()
52}
53
54#[inline(always)]
56pub fn norm1<T>(a: &[T]) -> T
57where
58 T: Float + Sum<T>,
59{
60 a.iter().map(|x| x.abs()).sum()
61}
62
63#[inline(always)]
65pub fn norm2<T>(a: &[T]) -> T
66where
67 T: Float + Sum<T> + Mul<T, Output = T>,
68{
69 let norm: T = norm2_squared(a);
70 norm.sqrt()
71}
72
73#[inline(always)]
75pub fn norm2_squared_diff<T>(a: &[T], b: &[T]) -> T
76where
77 T: Float + Sum<T> + Mul<T, Output = T> + std::ops::AddAssign,
78{
79 a.iter().zip(b.iter()).fold(T::zero(), |mut sum, (&x, &y)| {
80 sum += (x - y).powi(2);
81 sum
82 })
83}
84
85#[inline(always)]
87pub fn norm2_squared<T>(a: &[T]) -> T
88where
89 T: Float + Sum<T> + Mul<T, Output = T>,
90{
91 let norm: T = a.iter().map(|x| (*x) * (*x)).sum();
92 norm
93}
94
95#[inline(always)]
97pub fn sum<T>(a: &[T]) -> T
98where
99 T: Float + Sum<T> + Mul<T, Output = T>,
100{
101 let norm: T = a.iter().copied().sum();
102 norm
103}
104
105#[inline(always)]
107pub fn norm_inf<T>(a: &[T]) -> T
108where
109 T: Float + Zero,
110{
111 a.iter()
112 .fold(T::zero(), |current_max, x| x.abs().max(current_max))
113}
114
115#[inline(always)]
117pub fn norm_inf_diff<T>(a: &[T], b: &[T]) -> T
118where
119 T: Float + Zero,
120{
121 assert_eq!(a.len(), b.len());
122 a.iter()
123 .zip(b.iter())
124 .fold(T::zero(), |current_max, (x, y)| {
125 (*x - *y).abs().max(current_max)
126 })
127}
128
129#[inline(always)]
136pub fn is_finite<T>(a: &[T]) -> bool
137where
138 T: Float,
139{
140 !a.iter().any(|&xi| !xi.is_finite())
141}
142
143#[cfg(test)]
147mod tests {
148 use crate::*;
149
150 #[test]
151 fn t_inner_product_test() {
152 unit_test_utils::assert_nearly_equal(
153 14.0,
154 matrix_operations::inner_product(&[1.0, 2.0, 3.0], &[1.0, 2.0, 3.0]),
155 1e-10,
156 1e-16,
157 "inner product",
158 );
159 }
160
161 #[test]
162 #[should_panic]
163 fn t_inner_product_test_panic() {
164 matrix_operations::inner_product(&[2.0, 3.0], &[1.0, 2.0, 3.0]);
165 }
166
167 #[test]
168 fn t_norm1_test() {
169 unit_test_utils::assert_nearly_equal(
170 6.0,
171 matrix_operations::norm1(&[1.0, -2.0, -3.0]),
172 1e-10,
173 1e-16,
174 "norm1",
175 );
176 }
177
178 #[test]
179 fn t_norm2_test() {
180 unit_test_utils::assert_nearly_equal(
181 5.0,
182 matrix_operations::norm2(&[3.0, 4.0]),
183 1e-10,
184 1e-16,
185 "norm2",
186 );
187 }
188
189 #[test]
190 fn t_norm_inf_test() {
191 unit_test_utils::assert_nearly_equal(
192 3.0,
193 matrix_operations::norm_inf(&[1.0, -2.0, -3.0]),
194 1e-10,
195 1e-16,
196 "norm infinity of vector",
197 );
198 unit_test_utils::assert_nearly_equal(
199 8.0,
200 matrix_operations::norm_inf(&[1.0, -8.0, -3.0, 0.0]),
201 1e-10,
202 1e-16,
203 "infinity norm",
204 );
205 }
206
207 #[test]
208 fn t_norm_inf_diff() {
209 let x = [1.0, 2.0, 1.0];
210 let y = [-4.0, 0.0, 3.0];
211 let norm_diff = matrix_operations::norm_inf_diff(&x, &y);
212 unit_test_utils::assert_nearly_equal(5.0f64, norm_diff, 1e-10, 1e-9, "norm of difference");
213 unit_test_utils::assert_nearly_equal(
214 0.0f64,
215 matrix_operations::norm_inf_diff(&x, &x),
216 1e-10,
217 1e-16,
218 "difference of same vector",
219 );
220 unit_test_utils::assert_nearly_equal(
221 0.0f64,
222 matrix_operations::norm_inf_diff(&[], &[]),
223 1e-10,
224 1e-16,
225 "difference of empty vectors",
226 );
227 }
228
229 #[test]
230 #[should_panic]
231 fn t_norm_inf_diff_panic() {
232 let x = [1.0, 2.0, 3.0];
233 let y = [0.0, 3.0];
234 let _ = matrix_operations::norm_inf_diff(&x, &y);
235 }
236
237 #[test]
238 fn t_norm2_squared_diff_test() {
239 let x = [2.0, 5.0, 7.0, -1.0];
240 let y = [4.0, 1.0, 0.0, 10.0];
241 let norm2sq = matrix_operations::norm2_squared_diff(&x, &y);
242 unit_test_utils::assert_nearly_equal(190., norm2sq, 1e-10, 1e-12, "norm sq diff");
243 }
244}