diff --git a/CHANGELOG.md b/CHANGELOG.md index 71846f87..e5b81265 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,37 @@ documented here. This project adheres to [Semantic Versioning](https://semver.org/). +## [0.32.4] (19 Feb 2023) +- Add the `glam-0.25` feature to enable conversion from/to types from `glam` v0.25. + +## [0.32.3] (09 July 2023) + +### Modified +- Statically sized matrices are now serialized as tuples to match how serde + serialized plain arrays. +- Don’t require `Scalar` for matrix `PartialEq` and `Eq`. + +### Added +- Allow trailing punctuation in macros `vector!`, `matrix!`, `point!`, etc. +- Add the methods `Matrix1::as_scalar`, `::as_scalar_mut`, `::to_scalar`, `::into_scalar`. +- Add `Rotation3::euler_angles_ordered`, a generalized euler angles calculation. +- Add the `glam-0.24` feature to enable conversion from/to types from `glam` v0.24. +- Add the `lerp` method to points. +- Implement `Clone` for `MatrixIter`. + +### Fixed +- Fixed severe catastrophic cancellation issue in variance calculation. + +## [0.32.2] (07 March 2023) + +### Added +- Add the `glam-0.23` to enable conversion from/to type from `glam` v0.23. + +## [0.32.1] (14 Jan. 2023) + +### Modified +- Updated `nalgebra-macros` to use the new `Dyn`, avoiding macro-generated deprecation warnings. + ## [0.32.0] (14 Jan. 2023) ### Modified diff --git a/Cargo.toml b/Cargo.toml index 03217548..b62c4579 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra" -version = "0.32.0" +version = "0.32.4" authors = [ "Sébastien Crozet " ] description = "General-purpose linear algebra library with transformations and statically-sized or dynamically-sized matrices." @@ -47,6 +47,9 @@ convert-glam019 = [ "glam019" ] convert-glam020 = [ "glam020" ] convert-glam021 = [ "glam021" ] convert-glam022 = [ "glam022" ] +convert-glam023 = [ "glam023" ] +convert-glam024 = [ "glam024" ] +convert-glam025 = [ "glam025" ] # Serialization ## To use serde in a #[no-std] environment, enable the @@ -56,7 +59,7 @@ convert-glam022 = [ "glam022" ] serde-serialize-no-std = [ "serde", "num-complex/serde" ] serde-serialize = [ "serde-serialize-no-std", "serde/std" ] rkyv-serialize-no-std = [ "rkyv/size_32" ] -rkyv-serialize = [ "rkyv-serialize-no-std", "rkyv/std", "rkyv/validation", "bytecheck" ] +rkyv-serialize = [ "rkyv-serialize-no-std", "rkyv/std", "rkyv/validation" ] # Randomness ## To use rand in a #[no-std] environment, enable the @@ -71,7 +74,7 @@ slow-tests = [] rkyv-safe-deser = [ "rkyv-serialize", "rkyv/validation" ] [dependencies] -nalgebra-macros = { version = "0.1", path = "nalgebra-macros", optional = true } +nalgebra-macros = { version = "0.2.1", path = "nalgebra-macros", optional = true } typenum = "1.12" rand-package = { package = "rand", version = "0.8", optional = true, default-features = false } num-traits = { version = "0.2", default-features = false } @@ -83,8 +86,7 @@ alga = { version = "0.9", default-features = false, optional = true } rand_distr = { version = "0.4", default-features = false, optional = true } matrixmultiply = { version = "0.3", optional = true } serde = { version = "1.0", default-features = false, features = [ "derive" ], optional = true } -rkyv = { version = "0.7", default-features = false, optional = true } -bytecheck = { version = "~0.6.1", optional = true } +rkyv = { version = "0.7.41", default-features = false, optional = true } mint = { version = "0.5", optional = true } quickcheck = { version = "1", optional = true } pest = { version = "2", optional = true } @@ -101,6 +103,9 @@ glam019 = { package = "glam", version = "0.19", optional = true } glam020 = { package = "glam", version = "0.20", optional = true } glam021 = { package = "glam", version = "0.21", optional = true } glam022 = { package = "glam", version = "0.22", optional = true } +glam023 = { package = "glam", version = "0.23", optional = true } +glam024 = { package = "glam", version = "0.24", optional = true } +glam025 = { package = "glam", version = "0.25", optional = true } cust_core = { version = "0.1", optional = true } rayon = { version = "1.6", optional = true } @@ -109,6 +114,7 @@ serde_json = "1.0" rand_xorshift = "0.3" rand_isaac = "0.3" criterion = { version = "0.4", features = ["html_reports"] } +nalgebra = { path = ".", features = ["debug", "compare", "rand", "macros"]} # For matrix comparison macro matrixcompare = "0.3.0" diff --git a/README.md b/README.md index 62ab4759..857c84f2 100644 --- a/README.md +++ b/README.md @@ -29,22 +29,3 @@

----- - -## Acknowledgements -nalgebra is supported by our **platinum** sponsors: -

- - - -

- -And our gold sponsors: - -

- - - - - - -

\ No newline at end of file diff --git a/nalgebra-glm/src/gtc/epsilon.rs b/nalgebra-glm/src/gtc/epsilon.rs index fae29981..efe6ddd6 100644 --- a/nalgebra-glm/src/gtc/epsilon.rs +++ b/nalgebra-glm/src/gtc/epsilon.rs @@ -7,24 +7,24 @@ use na::DefaultAllocator; use crate::traits::{Alloc, Number, Dimension}; use crate::aliases::TVec; -/// Component-wise approximate equality beween two vectors. +/// Component-wise approximate equality between two vectors. pub fn epsilon_equal(x: &TVec, y: &TVec, epsilon: T) -> TVec where DefaultAllocator: Alloc { x.zip_map(y, |x, y| abs_diff_eq!(x, y, epsilon = epsilon)) } -/// Component-wise approximate equality beween two scalars. +/// Component-wise approximate equality between two scalars. pub fn epsilon_equal2>(x: T, y: T, epsilon: T) -> bool { abs_diff_eq!(x, y, epsilon = epsilon) } -/// Component-wise approximate non-equality beween two vectors. +/// Component-wise approximate non-equality between two vectors. pub fn epsilon_not_equal(x: &TVec, y: &TVec, epsilon: T) -> TVec where DefaultAllocator: Alloc { x.zip_map(y, |x, y| abs_diff_ne!(x, y, epsilon = epsilon)) } -/// Component-wise approximate non-equality beween two scalars. +/// Component-wise approximate non-equality between two scalars. pub fn epsilon_not_equal2>(x: T, y: T, epsilon: T) -> bool { abs_diff_ne!(x, y, epsilon = epsilon) } diff --git a/nalgebra-glm/src/gtx/quaternion.rs b/nalgebra-glm/src/gtx/quaternion.rs index d4f82af2..736d3bbb 100644 --- a/nalgebra-glm/src/gtx/quaternion.rs +++ b/nalgebra-glm/src/gtx/quaternion.rs @@ -80,7 +80,7 @@ pub fn quat_to_mat3(x: &Qua) -> TMat3 { .into_inner() } -/// Converts a quaternion to a rotation matrix in homogenous coordinates. +/// Converts a quaternion to a rotation matrix in homogeneous coordinates. pub fn quat_to_mat4(x: &Qua) -> TMat4 { UnitQuaternion::new_unchecked(*x).to_homogeneous() } diff --git a/nalgebra-macros/Cargo.toml b/nalgebra-macros/Cargo.toml index 041e36a7..8124ec78 100644 --- a/nalgebra-macros/Cargo.toml +++ b/nalgebra-macros/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-macros" -version = "0.1.0" +version = "0.2.1" authors = [ "Andreas Longva", "Sébastien Crozet " ] edition = "2018" description = "Procedural macros for nalgebra" diff --git a/nalgebra-macros/src/lib.rs b/nalgebra-macros/src/lib.rs index e12139d2..827d6080 100644 --- a/nalgebra-macros/src/lib.rs +++ b/nalgebra-macros/src/lib.rs @@ -223,7 +223,7 @@ impl Parse for Vector { elements: Vec::new(), }) } else { - let elements = MatrixRowSyntax::parse_separated_nonempty(input)? + let elements = MatrixRowSyntax::parse_terminated(input)? .into_iter() .collect(); Ok(Self { elements }) diff --git a/nalgebra-macros/tests/tests.rs b/nalgebra-macros/tests/tests.rs index 0e52da1f..ed6353d0 100644 --- a/nalgebra-macros/tests/tests.rs +++ b/nalgebra-macros/tests/tests.rs @@ -94,6 +94,12 @@ fn dmatrix_small_dims_exhaustive() { DMatrix::from_row_slice(4, 4, &[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16])); } +#[test] +fn matrix_trailing_semi() { + matrix![1, 2;]; + dmatrix![1, 2;]; +} + // Skip rustfmt because it just makes the test bloated without making it more readable #[rustfmt::skip] #[test] @@ -151,6 +157,13 @@ fn dvector_small_dims_exhaustive() { assert_eq_and_type!(dvector![1, 2, 3, 4, 5, 6], DVector::from_column_slice(&[1, 2, 3, 4, 5, 6])); } +#[test] +fn vector_trailing_comma() { + vector![1, 2,]; + point![1, 2,]; + dvector![1, 2,]; +} + #[test] fn matrix_trybuild_tests() { let t = trybuild::TestCases::new(); diff --git a/nalgebra-sparse/src/cs.rs b/nalgebra-sparse/src/cs.rs index 474eb2c0..e000e2de 100644 --- a/nalgebra-sparse/src/cs.rs +++ b/nalgebra-sparse/src/cs.rs @@ -494,7 +494,7 @@ where assert_eq!(source_minor_indices.len(), values.len()); let nnz = values.len(); - // Count the number of occurences of each minor index + // Count the number of occurrences of each minor index let mut minor_counts = vec![0; minor_dim]; for minor_idx in source_minor_indices { minor_counts[*minor_idx] += 1; diff --git a/nalgebra-sparse/src/ops/serial/csc.rs b/nalgebra-sparse/src/ops/serial/csc.rs index 5cf8ab23..a18cca3c 100644 --- a/nalgebra-sparse/src/ops/serial/csc.rs +++ b/nalgebra-sparse/src/ops/serial/csc.rs @@ -98,7 +98,7 @@ where /// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`. /// This will not return an error even if the patterns don't match. -/// Should be used for situations where pattern creation immediately preceeds multiplication. +/// Should be used for situations where pattern creation immediately precedes multiplication. /// /// Panics if the dimensions of the matrices involved are not compatible with the expression. pub fn spmm_csc_prealloc_unchecked( diff --git a/nalgebra-sparse/src/ops/serial/csr.rs b/nalgebra-sparse/src/ops/serial/csr.rs index d69bc54c..6384f26d 100644 --- a/nalgebra-sparse/src/ops/serial/csr.rs +++ b/nalgebra-sparse/src/ops/serial/csr.rs @@ -89,7 +89,7 @@ where /// Faster sparse-sparse matrix multiplication, `C <- beta * C + alpha * op(A) * op(B)`. /// This will not return an error even if the patterns don't match. -/// Should be used for situations where pattern creation immediately preceeds multiplication. +/// Should be used for situations where pattern creation immediately precedes multiplication. /// /// Panics if the dimensions of the matrices involved are not compatible with the expression. pub fn spmm_csr_prealloc_unchecked( diff --git a/src/base/alias.rs b/src/base/alias.rs index e3ac40b0..a07ea9e7 100644 --- a/src/base/alias.rs +++ b/src/base/alias.rs @@ -81,32 +81,32 @@ pub type MatrixXx5 = Matrix>; #[cfg(any(feature = "std", feature = "alloc"))] pub type MatrixXx6 = Matrix>; -/// A heap-allocated, row-major, matrix with 1 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 1 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] pub type Matrix1xX = Matrix>; -/// A heap-allocated, row-major, matrix with 2 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 2 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] pub type Matrix2xX = Matrix>; -/// A heap-allocated, row-major, matrix with 3 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 3 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] pub type Matrix3xX = Matrix>; -/// A heap-allocated, row-major, matrix with 4 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 4 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] pub type Matrix4xX = Matrix>; -/// A heap-allocated, row-major, matrix with 5 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 5 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] pub type Matrix5xX = Matrix>; -/// A heap-allocated, row-major, matrix with 6 rows and a dynamic number of columns. +/// A heap-allocated, column-major, matrix with 6 rows and a dynamic number of columns. /// /// **Because this is an alias, not all its methods are listed here. See the [`Matrix`](crate::base::Matrix) type too.** #[cfg(any(feature = "std", feature = "alloc"))] diff --git a/src/base/allocator.rs b/src/base/allocator.rs index 6458b8cb..10c4bd31 100644 --- a/src/base/allocator.rs +++ b/src/base/allocator.rs @@ -20,9 +20,9 @@ use std::mem::MaybeUninit; /// Every allocator must be both static and dynamic. Though not all implementations may share the /// same `Buffer` type. pub trait Allocator: Any + Sized { - /// The type of buffer this allocator can instanciate. + /// The type of buffer this allocator can instantiate. type Buffer: StorageMut + IsContiguous + Clone + Debug; - /// The type of buffer with uninitialized components this allocator can instanciate. + /// The type of buffer with uninitialized components this allocator can instantiate. type BufferUninit: RawStorageMut, R, C> + IsContiguous; /// Allocates a buffer with the given number of rows and columns without initializing its content. diff --git a/src/base/array_storage.rs b/src/base/array_storage.rs index 5c165399..56e88d47 100644 --- a/src/base/array_storage.rs +++ b/src/base/array_storage.rs @@ -5,12 +5,15 @@ use std::ops::Mul; #[cfg(feature = "serde-serialize-no-std")] use serde::de::{Error, SeqAccess, Visitor}; #[cfg(feature = "serde-serialize-no-std")] -use serde::ser::SerializeSeq; +use serde::ser::SerializeTuple; #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Deserializer, Serialize, Serializer}; #[cfg(feature = "serde-serialize-no-std")] use std::marker::PhantomData; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + use crate::base::allocator::Allocator; use crate::base::default_allocator::DefaultAllocator; use crate::base::dimension::{Const, ToTypenum}; @@ -189,7 +192,7 @@ where where S: Serializer, { - let mut serializer = serializer.serialize_seq(Some(R * C))?; + let mut serializer = serializer.serialize_tuple(R * C)?; for e in self.as_slice().iter() { serializer.serialize_element(e)?; @@ -208,7 +211,7 @@ where where D: Deserializer<'a>, { - deserializer.deserialize_seq(ArrayStorageVisitor::new()) + deserializer.deserialize_tuple(R * C, ArrayStorageVisitor::new()) } } diff --git a/src/base/dimension.rs b/src/base/dimension.rs index 97616129..11743dd8 100644 --- a/src/base/dimension.rs +++ b/src/base/dimension.rs @@ -8,6 +8,8 @@ use std::fmt::Debug; use std::ops::{Add, Div, Mul, Sub}; use typenum::{self, Diff, Max, Maximum, Min, Minimum, Prod, Quot, Sum, Unsigned}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Deserializer, Serialize, Serializer}; diff --git a/src/base/edition.rs b/src/base/edition.rs index e482fa24..8994eed7 100644 --- a/src/base/edition.rs +++ b/src/base/edition.rs @@ -1077,7 +1077,7 @@ where } // Move the elements of `data` in such a way that the matrix with -// the rows `[i, i + nremove[` deleted is represented in a contigous +// the rows `[i, i + nremove[` deleted is represented in a contiguous // way in `data` after this method completes. // Every deleted element are manually dropped by this method. unsafe fn compress_rows( diff --git a/src/base/iter.rs b/src/base/iter.rs index 0e4aa8d4..ebffdb07 100644 --- a/src/base/iter.rs +++ b/src/base/iter.rs @@ -15,9 +15,9 @@ use crate::base::storage::{RawStorage, RawStorageMut}; use crate::base::{Matrix, MatrixView, MatrixViewMut, Scalar}; macro_rules! iterator { - (struct $Name:ident for $Storage:ident.$ptr: ident -> $Ptr:ty, $Ref:ty, $SRef: ty) => { + (struct $Name:ident for $Storage:ident.$ptr: ident -> $Ptr:ty, $Ref:ty, $SRef: ty, $($derives:ident),* $(,)?) => { /// An iterator through a dense matrix with arbitrary strides matrix. - #[derive(Debug)] + #[derive($($derives),*)] pub struct $Name<'a, T, R: Dim, C: Dim, S: 'a + $Storage> { ptr: $Ptr, inner_ptr: $Ptr, @@ -39,7 +39,7 @@ macro_rules! iterator { let ptr = storage.$ptr(); // If we have a size of 0, 'ptr' must be - // dangling. Howver, 'inner_offset' might + // dangling. However, 'inner_offset' might // not be zero if only one dimension is zero, so // we don't want to call 'offset'. // This pointer will never actually get used @@ -177,8 +177,8 @@ macro_rules! iterator { }; } -iterator!(struct MatrixIter for RawStorage.ptr -> *const T, &'a T, &'a S); -iterator!(struct MatrixIterMut for RawStorageMut.ptr_mut -> *mut T, &'a mut T, &'a mut S); +iterator!(struct MatrixIter for RawStorage.ptr -> *const T, &'a T, &'a S, Clone, Debug); +iterator!(struct MatrixIterMut for RawStorageMut.ptr_mut -> *mut T, &'a mut T, &'a mut S, Debug); /* * diff --git a/src/base/matrix.rs b/src/base/matrix.rs index d4875944..af5609df 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -13,6 +13,8 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer}; #[cfg(feature = "rkyv-serialize-no-std")] use super::rkyv_wrappers::CustomPhantom; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; #[cfg(feature = "rkyv-serialize-no-std")] use rkyv::{with::With, Archive, Archived}; @@ -1859,14 +1861,14 @@ where impl Eq for Matrix where - T: Scalar + Eq, + T: Eq, S: RawStorage, { } impl PartialEq> for Matrix where - T: Scalar + PartialEq, + T: PartialEq, C: Dim, C2: Dim, R: Dim, @@ -2244,3 +2246,102 @@ where Unit::new_unchecked(crate::convert_ref(self.as_ref())) } } + +impl Matrix +where + S: RawStorage, +{ + /// Returns a reference to the single element in this matrix. + /// + /// As opposed to indexing, using this provides type-safety + /// when flattening dimensions. + /// + /// # Example + /// ``` + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let inner_product: f32 = *(v.transpose() * v).as_scalar(); + /// ``` + /// + ///```compile_fail + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let inner_product = (v * v.transpose()).item(); // Typo, does not compile. + ///``` + pub fn as_scalar(&self) -> &T { + &self[(0, 0)] + } + /// Get a mutable reference to the single element in this matrix + /// + /// As opposed to indexing, using this provides type-safety + /// when flattening dimensions. + /// + /// # Example + /// ``` + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let mut inner_product = (v.transpose() * v); + /// *inner_product.as_scalar_mut() = 3.; + /// ``` + /// + ///```compile_fail + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let mut inner_product = (v * v.transpose()); + /// *inner_product.as_scalar_mut() = 3.; + ///``` + pub fn as_scalar_mut(&mut self) -> &mut T + where + S: RawStorageMut, + { + &mut self[(0, 0)] + } + /// Convert this 1x1 matrix by reference into a scalar. + /// + /// As opposed to indexing, using this provides type-safety + /// when flattening dimensions. + /// + /// # Example + /// ``` + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let mut inner_product: f32 = (v.transpose() * v).to_scalar(); + /// ``` + /// + ///```compile_fail + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let mut inner_product: f32 = (v * v.transpose()).to_scalar(); + ///``` + pub fn to_scalar(&self) -> T + where + T: Clone, + { + self.as_scalar().clone() + } +} + +impl super::alias::Matrix1 { + /// Convert this 1x1 matrix into a scalar. + /// + /// As opposed to indexing, using this provides type-safety + /// when flattening dimensions. + /// + /// # Example + /// ``` + /// # use nalgebra::{Vector3, Matrix2, U1}; + /// let v = Vector3::new(0., 0., 1.); + /// let inner_product: f32 = (v.transpose() * v).into_scalar(); + /// assert_eq!(inner_product, 1.); + /// ``` + /// + ///```compile_fail + /// # use nalgebra::Vector3; + /// let v = Vector3::new(0., 0., 1.); + /// let mut inner_product: f32 = (v * v.transpose()).into_scalar(); + ///``` + pub fn into_scalar(self) -> T { + let [[scalar]] = self.data.0; + scalar + } +} diff --git a/src/base/par_iter.rs b/src/base/par_iter.rs index af5e1cb7..c4af719a 100644 --- a/src/base/par_iter.rs +++ b/src/base/par_iter.rs @@ -11,7 +11,7 @@ use crate::{ use rayon::iter::plumbing::Producer; use rayon::{iter::plumbing::bridge, prelude::*}; -/// A rayon parallel iterator over the colums of a matrix. It is created +/// A rayon parallel iterator over the columns of a matrix. It is created /// using the [`par_column_iter`] method of [`Matrix`]. /// /// *Only available if compiled with the feature `rayon`.* @@ -89,7 +89,7 @@ pub struct ParColumnIterMut< } #[cfg_attr(doc_cfg, doc(cfg(feature = "rayon")))] -/// *only availabe if compiled with the feature `rayon`* +/// *only available if compiled with the feature `rayon`* impl<'a, T, R, Cols, S> ParColumnIterMut<'a, T, R, Cols, S> where R: Dim, @@ -161,7 +161,7 @@ where S: Sync, { /// Iterate through the columns of the matrix in parallel using rayon. - /// This iterates over *immutable* references ot the columns of the matrix, + /// This iterates over *immutable* references to the columns of the matrix, /// if *mutable* access to the columns is required, use [`par_column_iter_mut`] /// instead. /// diff --git a/src/base/statistics.rs b/src/base/statistics.rs index 9f0e0ee6..6007f8c7 100644 --- a/src/base/statistics.rs +++ b/src/base/statistics.rs @@ -335,12 +335,12 @@ impl> Matrix { if self.is_empty() { T::zero() } else { - let val = self.iter().cloned().fold((T::zero(), T::zero()), |a, b| { - (a.0 + b.clone() * b.clone(), a.1 + b) - }); - let denom = T::one() / crate::convert::<_, T>(self.len() as f64); - let vd = val.1 * denom.clone(); - val.0 * denom - vd.clone() * vd + let n_elements: T = crate::convert(self.len() as f64); + let mean = self.mean(); + + self.iter().cloned().fold(T::zero(), |acc, x| { + acc + (x.clone() - mean.clone()) * (x.clone() - mean.clone()) + }) / n_elements } } diff --git a/src/base/uninit.rs b/src/base/uninit.rs index ad2759eb..401e3336 100644 --- a/src/base/uninit.rs +++ b/src/base/uninit.rs @@ -34,7 +34,7 @@ pub unsafe trait InitStatus: Copy { /// A type implementing `InitStatus` indicating that the value is completely initialized. pub struct Init; #[derive(Copy, Clone, Debug, PartialEq, Eq)] -/// A type implementing `InitStatus` indicating that the value is completely unitialized. +/// A type implementing `InitStatus` indicating that the value is completely uninitialized. pub struct Uninit; unsafe impl InitStatus for Init { diff --git a/src/base/unit.rs b/src/base/unit.rs index 2fc51107..12c70963 100644 --- a/src/base/unit.rs +++ b/src/base/unit.rs @@ -9,6 +9,9 @@ use crate::base::DefaultAllocator; use crate::storage::RawStorage; use crate::{Dim, Matrix, OMatrix, RealField, Scalar, SimdComplexField, SimdRealField}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A wrapper that ensures the underlying algebraic entity has a unit norm. /// /// **It is likely that the only piece of documentation that you need in this page are:** diff --git a/src/base/vec_storage.rs b/src/base/vec_storage.rs index 4614598b..42c4511b 100644 --- a/src/base/vec_storage.rs +++ b/src/base/vec_storage.rs @@ -148,7 +148,7 @@ impl VecStorage { }; // Avoid double-free by forgetting `self` because its data buffer has - // been transfered to `new_data`. + // been transferred to `new_data`. std::mem::forget(self); new_data } diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs index bae04f46..2272b622 100644 --- a/src/geometry/dual_quaternion.rs +++ b/src/geometry/dual_quaternion.rs @@ -1,6 +1,9 @@ // The macros break if the references are taken out, for some reason. #![allow(clippy::op_ref)] +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + use crate::{ Isometry3, Matrix4, Normed, OVector, Point3, Quaternion, Scalar, SimdRealField, Translation3, Unit, UnitQuaternion, Vector3, Zero, U8, diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs old mode 100755 new mode 100644 index e3f44075..376af06f --- a/src/geometry/isometry.rs +++ b/src/geometry/isometry.rs @@ -14,6 +14,9 @@ use crate::base::storage::Owned; use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar, Unit}; use crate::geometry::{AbstractRotation, Point, Translation}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A direct isometry, i.e., a rotation followed by a translation (aka. a rigid-body motion). /// /// This is also known as an element of a Special Euclidean (SE) group. diff --git a/src/geometry/isometry_interpolation.rs b/src/geometry/isometry_interpolation.rs index 90f2c7ae..d6c20503 100644 --- a/src/geometry/isometry_interpolation.rs +++ b/src/geometry/isometry_interpolation.rs @@ -42,7 +42,7 @@ impl Isometry3 { /// Attempts to interpolate between two isometries using a linear interpolation for the translation part, /// and a spherical interpolation for the rotation part. /// - /// Retuns `None` if the angle between both rotations is 180 degrees (in which case the interpolation + /// Returns `None` if the angle between both rotations is 180 degrees (in which case the interpolation /// is not well-defined). /// /// # Examples: @@ -118,7 +118,7 @@ impl IsometryMatrix3 { /// Attempts to interpolate between two isometries using a linear interpolation for the translation part, /// and a spherical interpolation for the rotation part. /// - /// Retuns `None` if the angle between both rotations is 180 degrees (in which case the interpolation + /// Returns `None` if the angle between both rotations is 180 degrees (in which case the interpolation /// is not well-defined). /// /// # Examples: diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs index 06d0b471..2d809fe4 100644 --- a/src/geometry/orthographic.rs +++ b/src/geometry/orthographic.rs @@ -17,6 +17,9 @@ use crate::base::{Matrix4, Vector, Vector3}; use crate::geometry::{Point3, Projective3}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A 3D orthographic projection stored as a homogeneous 4x4 matrix. #[repr(C)] #[cfg_attr( diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs index a5fc19a8..383c99b9 100644 --- a/src/geometry/perspective.rs +++ b/src/geometry/perspective.rs @@ -18,6 +18,9 @@ use crate::base::{Matrix4, Vector, Vector3}; use crate::geometry::{Point3, Projective3}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A 3D perspective projection stored as a homogeneous 4x4 matrix. #[repr(C)] #[cfg_attr( diff --git a/src/geometry/point.rs b/src/geometry/point.rs index 71ace3a9..e936c3da 100644 --- a/src/geometry/point.rs +++ b/src/geometry/point.rs @@ -1,9 +1,11 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq}; -use num::One; +use num::{One, Zero}; use std::cmp::Ordering; use std::fmt; use std::hash; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; #[cfg(feature = "serde-serialize-no-std")] use serde::{Deserialize, Deserializer, Serialize, Serializer}; @@ -13,6 +15,7 @@ use crate::base::allocator::Allocator; use crate::base::dimension::{DimName, DimNameAdd, DimNameSum, U1}; use crate::base::iter::{MatrixIter, MatrixIterMut}; use crate::base::{Const, DefaultAllocator, OVector, Scalar}; +use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub}; use std::mem::MaybeUninit; /// A point in an euclidean space. @@ -221,6 +224,31 @@ where unsafe { res.assume_init() } } + /// Linear interpolation between two points. + /// + /// Returns `self * (1.0 - t) + rhs.coords * t`, i.e., the linear blend of the points + /// `self` and `rhs` using the scalar value `t`. + /// + /// The value for a is not restricted to the range `[0, 1]`. + /// + /// # Examples: + /// + /// ``` + /// # use nalgebra::Point3; + /// let a = Point3::new(1.0, 2.0, 3.0); + /// let b = Point3::new(10.0, 20.0, 30.0); + /// assert_eq!(a.lerp(&b, 0.1), Point3::new(1.9, 3.8, 5.7)); + /// ``` + #[must_use] + pub fn lerp(&self, rhs: &OPoint, t: T) -> OPoint + where + T: Scalar + Zero + One + ClosedAdd + ClosedSub + ClosedMul, + { + OPoint { + coords: self.coords.lerp(&rhs.coords, t), + } + } + /// Creates a new point with the given coordinates. #[deprecated(note = "Use Point::from(vector) instead.")] #[inline] diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs old mode 100755 new mode 100644 index 1b251b29..7dcbce51 --- a/src/geometry/quaternion.rs +++ b/src/geometry/quaternion.rs @@ -19,6 +19,9 @@ use crate::base::{ use crate::geometry::{Point3, Rotation}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A quaternion. See the type alias `UnitQuaternion = Unit` for a quaternion /// that may be used as a rotation. #[repr(C)] @@ -1577,7 +1580,7 @@ where #[inline] #[must_use] pub fn inverse_transform_point(&self, pt: &Point3) -> Point3 { - // TODO: would it be useful performancewise not to call inverse explicitly (i-e. implement + // TODO: would it be useful performance-wise not to call inverse explicitly (i-e. implement // the inverse transformation explicitly here) ? self.inverse() * pt } diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs old mode 100755 new mode 100644 index 5eceec21..67a242cb --- a/src/geometry/rotation.rs +++ b/src/geometry/rotation.rs @@ -17,6 +17,9 @@ use crate::base::dimension::{DimNameAdd, DimNameSum, U1}; use crate::base::{Const, DefaultAllocator, OMatrix, SMatrix, SVector, Scalar, Unit}; use crate::geometry::Point; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A rotation matrix. /// /// This is also known as an element of a Special Orthogonal (SO) group. diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs index 1ea5cd92..eecc9b21 100644 --- a/src/geometry/rotation_specialization.rs +++ b/src/geometry/rotation_specialization.rs @@ -478,9 +478,10 @@ where SB: Storage, SC: Storage, { + // Gram–Schmidt process let zaxis = dir.normalize(); let xaxis = up.cross(&zaxis).normalize(); - let yaxis = zaxis.cross(&xaxis).normalize(); + let yaxis = zaxis.cross(&xaxis); Self::from_matrix_unchecked(SMatrix::::new( xaxis.x.clone(), @@ -979,6 +980,176 @@ impl Rotation3 { ) } } + + /// Represent this rotation as Euler angles. + /// + /// Returns the angles produced in the order provided by seq parameter, along with the + /// observability flag. The Euler axes passed to seq must form an orthonormal basis. If the + /// rotation is gimbal locked, then the observability flag is false. + /// + /// # Panics + /// + /// Panics if the Euler axes in `seq` are not orthonormal. + /// + /// # Example 1: + /// ``` + /// use std::f64::consts::PI; + /// use approx::assert_relative_eq; + /// use nalgebra::{Matrix3, Rotation3, Unit, Vector3}; + /// + /// // 3-1-2 + /// let n = [ + /// Unit::new_unchecked(Vector3::new(0.0, 0.0, 1.0)), + /// Unit::new_unchecked(Vector3::new(1.0, 0.0, 0.0)), + /// Unit::new_unchecked(Vector3::new(0.0, 1.0, 0.0)), + /// ]; + /// + /// let r1 = Rotation3::from_axis_angle(&n[2], 20.0 * PI / 180.0); + /// let r2 = Rotation3::from_axis_angle(&n[1], 30.0 * PI / 180.0); + /// let r3 = Rotation3::from_axis_angle(&n[0], 45.0 * PI / 180.0); + /// + /// let d = r3 * r2 * r1; + /// + /// let (angles, observable) = d.euler_angles_ordered(n, false); + /// assert!(observable); + /// assert_relative_eq!(angles[0] * 180.0 / PI, 45.0, epsilon = 1e-12); + /// assert_relative_eq!(angles[1] * 180.0 / PI, 30.0, epsilon = 1e-12); + /// assert_relative_eq!(angles[2] * 180.0 / PI, 20.0, epsilon = 1e-12); + /// ``` + /// + /// # Example 2: + /// ``` + /// use std::f64::consts::PI; + /// use approx::assert_relative_eq; + /// use nalgebra::{Matrix3, Rotation3, Unit, Vector3}; + /// + /// let sqrt_2 = 2.0_f64.sqrt(); + /// let n = [ + /// Unit::new_unchecked(Vector3::new(1.0 / sqrt_2, 1.0 / sqrt_2, 0.0)), + /// Unit::new_unchecked(Vector3::new(1.0 / sqrt_2, -1.0 / sqrt_2, 0.0)), + /// Unit::new_unchecked(Vector3::new(0.0, 0.0, 1.0)), + /// ]; + /// + /// let r1 = Rotation3::from_axis_angle(&n[2], 20.0 * PI / 180.0); + /// let r2 = Rotation3::from_axis_angle(&n[1], 30.0 * PI / 180.0); + /// let r3 = Rotation3::from_axis_angle(&n[0], 45.0 * PI / 180.0); + /// + /// let d = r3 * r2 * r1; + /// + /// let (angles, observable) = d.euler_angles_ordered(n, false); + /// assert!(observable); + /// assert_relative_eq!(angles[0] * 180.0 / PI, 45.0, epsilon = 1e-12); + /// assert_relative_eq!(angles[1] * 180.0 / PI, 30.0, epsilon = 1e-12); + /// assert_relative_eq!(angles[2] * 180.0 / PI, 20.0, epsilon = 1e-12); + /// ``` + /// + /// Algorithm based on: + /// Malcolm D. Shuster, F. Landis Markley, “General formula for extraction the Euler + /// angles”, Journal of guidance, control, and dynamics, vol. 29.1, pp. 215-221. 2006, + /// and modified to be able to produce extrinsic rotations. + #[must_use] + pub fn euler_angles_ordered( + &self, + mut seq: [Unit>; 3], + extrinsic: bool, + ) -> ([T; 3], bool) + where + T: RealField + Copy, + { + let mut angles = [T::zero(); 3]; + let eps = T::from_subset(&1e-7); + let _2 = T::from_subset(&2.0); + + if extrinsic { + seq.reverse(); + } + + let [n1, n2, n3] = &seq; + assert_relative_eq!(n1.dot(n2), T::zero(), epsilon = eps); + assert_relative_eq!(n3.dot(n1), T::zero(), epsilon = eps); + + let n1_c_n2 = n1.cross(n2); + let s1 = n1_c_n2.dot(n3); + let c1 = n1.dot(n3); + let lambda = s1.atan2(c1); + + let mut c = Matrix3::zeros(); + c.column_mut(0).copy_from(n2); + c.column_mut(1).copy_from(&n1_c_n2); + c.column_mut(2).copy_from(n1); + c.transpose_mut(); + + let r1l = Matrix3::new( + T::one(), + T::zero(), + T::zero(), + T::zero(), + c1, + s1, + T::zero(), + -s1, + c1, + ); + let o_t = &c * self.matrix() * (c.transpose() * r1l); + angles[1] = o_t.m33.acos(); + + let safe1 = angles[1].abs() >= eps; + let safe2 = (angles[1] - T::pi()).abs() >= eps; + let observable = safe1 && safe2; + angles[1] += lambda; + + if observable { + angles[0] = o_t.m13.atan2(-o_t.m23); + angles[2] = o_t.m31.atan2(o_t.m32); + } else { + // gimbal lock detected + if extrinsic { + // angle1 is initialized to zero + if !safe1 { + angles[2] = (o_t.m12 - o_t.m21).atan2(o_t.m11 + o_t.m22); + } else { + angles[2] = -(o_t.m12 + o_t.m21).atan2(o_t.m11 - o_t.m22); + }; + } else { + // angle3 is initialized to zero + if !safe1 { + angles[0] = (o_t.m12 - o_t.m21).atan2(o_t.m11 + o_t.m22); + } else { + angles[0] = (o_t.m12 + o_t.m21).atan2(o_t.m11 - o_t.m22); + }; + }; + }; + + let adjust = if seq[0] == seq[2] { + // lambda = 0, so ensure angle2 -> [0, pi] + angles[1] < T::zero() || angles[1] > T::pi() + } else { + // lamda = + or - pi/2, so ensure angle2 -> [-pi/2, pi/2] + angles[1] < -T::frac_pi_2() || angles[1] > T::frac_pi_2() + }; + + // dont adjust gimbal locked rotation + if adjust && observable { + angles[0] += T::pi(); + angles[1] = _2 * lambda - angles[1]; + angles[2] -= T::pi(); + } + + // ensure all angles are within [-pi, pi] + for angle in angles.as_mut_slice().iter_mut() { + if *angle < -T::pi() { + *angle += T::two_pi(); + } else if *angle > T::pi() { + *angle -= T::two_pi(); + } + } + + if extrinsic { + angles.reverse(); + } + + (angles, observable) + } } #[cfg(feature = "rand-no-std")] diff --git a/src/geometry/scale.rs b/src/geometry/scale.rs old mode 100755 new mode 100644 index 36a68066..c95c0967 --- a/src/geometry/scale.rs +++ b/src/geometry/scale.rs @@ -15,6 +15,9 @@ use crate::ClosedMul; use crate::geometry::Point; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A scale which supports non-uniform scaling. #[repr(C)] #[cfg_attr( diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs old mode 100755 new mode 100644 index 989b721b..4d931947 --- a/src/geometry/similarity.rs +++ b/src/geometry/similarity.rs @@ -15,6 +15,9 @@ use crate::base::storage::Owned; use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar}; use crate::geometry::{AbstractRotation, Isometry, Point, Translation}; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation. #[repr(C)] #[derive(Debug, Copy, Clone)] diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs index 2a7ca112..73dc8d8a 100755 --- a/src/geometry/transform.rs +++ b/src/geometry/transform.rs @@ -122,7 +122,7 @@ macro_rules! category_mul_impl( )*} ); -// We require stability uppon multiplication. +// We require stability upon multiplication. impl TCategoryMul for T { type Representative = T; } diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs old mode 100755 new mode 100644 index 39fae3b6..4fc50777 --- a/src/geometry/translation.rs +++ b/src/geometry/translation.rs @@ -15,6 +15,9 @@ use crate::base::{Const, DefaultAllocator, OMatrix, SVector, Scalar}; use crate::geometry::Point; +#[cfg(feature = "rkyv-serialize")] +use rkyv::bytecheck; + /// A translation. #[repr(C)] #[cfg_attr( diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs index 8e44f71a..d6c0ade5 100755 --- a/src/geometry/unit_complex.rs +++ b/src/geometry/unit_complex.rs @@ -347,7 +347,7 @@ where #[inline] #[must_use] pub fn inverse_transform_point(&self, pt: &Point2) -> Point2 { - // TODO: would it be useful performancewise not to call inverse explicitly (i-e. implement + // TODO: would it be useful performance-wise not to call inverse explicitly (i-e. implement // the inverse transformation explicitly here) ? self.inverse() * pt } diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs index 2402bb3d..b0abe994 100644 --- a/src/linalg/convolution.rs +++ b/src/linalg/convolution.rs @@ -47,11 +47,11 @@ impl> Vector { let u_f = cmp::min(i, vec - 1); if u_i == u_f { - conv[i] += self[u_i].clone() * kernel[(i - u_i)].clone(); + conv[i] += self[u_i].clone() * kernel[i - u_i].clone(); } else { for u in u_i..(u_f + 1) { if i - u < ker { - conv[i] += self[u].clone() * kernel[(i - u)].clone(); + conv[i] += self[u].clone() * kernel[i - u].clone(); } } } diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 06bae4a3..39283e24 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -724,7 +724,7 @@ where /// Sort the estimated components of the SVD by its singular values in descending order. /// Such an ordering is often implicitly required when the decompositions are used for estimation or fitting purposes. - /// Using this function is only required if `new_unordered` or `try_new_unorderd` were used and the specific sorting is required afterward. + /// Using this function is only required if `new_unordered` or `try_new_unordered` were used and the specific sorting is required afterward. pub fn sort_by_singular_values(&mut self) { const VALUE_PROCESSED: usize = usize::MAX; diff --git a/src/sparse/cs_matrix.rs b/src/sparse/cs_matrix.rs index 5b63e537..9a240ff6 100644 --- a/src/sparse/cs_matrix.rs +++ b/src/sparse/cs_matrix.rs @@ -498,7 +498,7 @@ where } } - // Remove dupliate entries on a sorted CsMatrix. + // Remove duplicate entries on a sorted CsMatrix. pub(crate) fn dedup(&mut self) where T: Zero + ClosedAdd, diff --git a/src/third_party/glam/mod.rs b/src/third_party/glam/mod.rs index 811e88b2..680e65e5 100644 --- a/src/third_party/glam/mod.rs +++ b/src/third_party/glam/mod.rs @@ -16,3 +16,9 @@ mod v020; mod v021; #[cfg(feature = "glam022")] mod v022; +#[cfg(feature = "glam023")] +mod v023; +#[cfg(feature = "glam024")] +mod v024; +#[cfg(feature = "glam025")] +mod v025; diff --git a/src/third_party/glam/v023/mod.rs b/src/third_party/glam/v023/mod.rs new file mode 100644 index 00000000..45b1ad47 --- /dev/null +++ b/src/third_party/glam/v023/mod.rs @@ -0,0 +1,18 @@ +#[path = "../common/glam_isometry.rs"] +mod glam_isometry; +#[path = "../common/glam_matrix.rs"] +mod glam_matrix; +#[path = "../common/glam_point.rs"] +mod glam_point; +#[path = "../common/glam_quaternion.rs"] +mod glam_quaternion; +#[path = "../common/glam_rotation.rs"] +mod glam_rotation; +#[path = "../common/glam_similarity.rs"] +mod glam_similarity; +#[path = "../common/glam_translation.rs"] +mod glam_translation; +#[path = "../common/glam_unit_complex.rs"] +mod glam_unit_complex; + +pub(self) use glam023 as glam; diff --git a/src/third_party/glam/v024/mod.rs b/src/third_party/glam/v024/mod.rs new file mode 100644 index 00000000..3dae276f --- /dev/null +++ b/src/third_party/glam/v024/mod.rs @@ -0,0 +1,18 @@ +#[path = "../common/glam_isometry.rs"] +mod glam_isometry; +#[path = "../common/glam_matrix.rs"] +mod glam_matrix; +#[path = "../common/glam_point.rs"] +mod glam_point; +#[path = "../common/glam_quaternion.rs"] +mod glam_quaternion; +#[path = "../common/glam_rotation.rs"] +mod glam_rotation; +#[path = "../common/glam_similarity.rs"] +mod glam_similarity; +#[path = "../common/glam_translation.rs"] +mod glam_translation; +#[path = "../common/glam_unit_complex.rs"] +mod glam_unit_complex; + +pub(self) use glam024 as glam; diff --git a/src/third_party/glam/v025/mod.rs b/src/third_party/glam/v025/mod.rs new file mode 100644 index 00000000..b9f41e34 --- /dev/null +++ b/src/third_party/glam/v025/mod.rs @@ -0,0 +1,18 @@ +#[path = "../common/glam_isometry.rs"] +mod glam_isometry; +#[path = "../common/glam_matrix.rs"] +mod glam_matrix; +#[path = "../common/glam_point.rs"] +mod glam_point; +#[path = "../common/glam_quaternion.rs"] +mod glam_quaternion; +#[path = "../common/glam_rotation.rs"] +mod glam_rotation; +#[path = "../common/glam_similarity.rs"] +mod glam_similarity; +#[path = "../common/glam_translation.rs"] +mod glam_translation; +#[path = "../common/glam_unit_complex.rs"] +mod glam_unit_complex; + +pub(self) use glam025 as glam; diff --git a/tests/core/mod.rs b/tests/core/mod.rs index 0f7ee85b..f0484e4d 100644 --- a/tests/core/mod.rs +++ b/tests/core/mod.rs @@ -11,6 +11,7 @@ mod reshape; #[cfg(feature = "rkyv-serialize-no-std")] mod rkyv; mod serde; +mod variance; #[cfg(feature = "compare")] mod matrixcompare; diff --git a/tests/core/variance.rs b/tests/core/variance.rs new file mode 100644 index 00000000..eb08ea0f --- /dev/null +++ b/tests/core/variance.rs @@ -0,0 +1,18 @@ +use nalgebra::DVector; + +#[test] +fn test_variance_catastrophic_cancellation() { + let long_repeating_vector = DVector::repeat(10_000, 100000000.0); + assert_eq!(long_repeating_vector.variance(), 0.0); + + let short_vec = DVector::from_vec(vec![1., 2., 3.]); + assert_eq!(short_vec.variance(), 2.0 / 3.0); + + let short_vec = + DVector::::from_vec(vec![1.0e8 + 4.0, 1.0e8 + 7.0, 1.0e8 + 13.0, 1.0e8 + 16.0]); + assert_eq!(short_vec.variance(), 22.5); + + let short_vec = + DVector::::from_vec(vec![1.0e9 + 4.0, 1.0e9 + 7.0, 1.0e9 + 13.0, 1.0e9 + 16.0]); + assert_eq!(short_vec.variance(), 22.5); +} diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 162aad6a..a5dcf835 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -123,7 +123,7 @@ fn symmetric_eigen_singular_24x24() { // // /* // * NOTE: for the following tests, we use only upper-triangular matrices. -// * Thes ensures the schur decomposition will work, and allows use to test the eigenvector +// * This ensures the schur decomposition will work, and allows use to test the eigenvector // * computation. // */ // fn eigen(n: usize) -> bool { @@ -134,11 +134,11 @@ fn symmetric_eigen_singular_24x24() { // verify_eigenvectors(m, eig) // } // -// fn eigen_with_adjascent_duplicate_diagonals(n: usize) -> bool { +// fn eigen_with_adjacent_duplicate_diagonals(n: usize) -> bool { // let n = cmp::max(1, cmp::min(n, 10)); // let mut m = DMatrix::::new_random(n, n).upper_triangle(); // -// // Suplicate some adjascent diagonal elements. +// // Suplicate some adjacent diagonal elements. // for i in 0 .. n / 2 { // m[(i * 2 + 1, i * 2 + 1)] = m[(i * 2, i * 2)]; // } @@ -147,7 +147,7 @@ fn symmetric_eigen_singular_24x24() { // verify_eigenvectors(m, eig) // } // -// fn eigen_with_nonadjascent_duplicate_diagonals(n: usize) -> bool { +// fn eigen_with_nonadjacent_duplicate_diagonals(n: usize) -> bool { // let n = cmp::max(3, cmp::min(n, 10)); // let mut m = DMatrix::::new_random(n, n).upper_triangle(); //