//! Defines lattices and lattice component.
//!
//! [`LatticeCyclic`] is the structure that encode the lattice information like
//! the lattice spacing, the number of point and the dimension.
//! It is used to do operation on [`LatticePoint`], [`LatticeLink`] and
//! [`LatticeLinkCanonical`].
//! Or to get an iterator over these elements.
//!
//! [`LatticePoint`], [`LatticeLink`] and [`LatticeLinkCanonical`] are elements on the lattice.
//! They encode where a field element is situated.
use std::cmp::Ordering;
use std::convert::TryInto;
use std::iter::FusedIterator;
use std::ops::{Index, IndexMut, Neg};
use lattice_qcd_rs_procedural_macro::{implement_direction_from, implement_direction_list};
use na::{SVector, Vector4};
#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Serialize};
use super::field::{EField, LinkMatrix};
use super::{error::LatticeInitializationError, Real};
/// A Cyclic lattice in space. Does not store point and links but is used to generate them.
///
/// The generic parameter `D` is the dimension.
///
/// This lattice is Cyclic more precisely if the lattice has N points in each direction.
/// Then we can move alongside a direction going though point 0, 1, ... N-1. The next step in
/// the same direction goes back to the point at 0.
///
/// This structure is used for operation on [`LatticePoint`], [`LatticeLink`] and
/// [`LatticeLinkCanonical`].
// For example, theses three structures are abstract and are in general use to
/// access data on the lattice. These data are stored [`LinkMatrix`] and [`EField`] which are just
/// a wrapper around a [`Vec`]. `LatticeCyclic` is used to convert the lattice element to
/// an index to access these data.
///
/// This contain very few data and can be cloned at almost no cost even though
/// it does not implement [`Copy`].
#[derive(Clone, Debug, PartialEq)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct LatticeCyclic<const D: usize> {
/// The lattice spacing.
size: Real,
/// The number of point *per* dimension.
dim: usize,
}
impl<const D: usize> LatticeCyclic<D> {
/// Number space + time dimension, this is the `D` parameter.
///
/// Not to confuse with [`LatticeCyclic::dim`] which is the number of point per dimension.
pub const fn dim_st() -> usize {
D
}
/// see [`LatticeLinkCanonical`], a conical link is a link whose direction is always positive.
/// That means that a link form `[x, y, z, t]` with direction `-x`
/// the link return is `[x - 1, y, z, t]` (modulo the `lattice::dim()``) with direction `+x`
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{LatticeCyclic, DirectionEnum, LatticePoint, LatticeLinkCanonical};
/// # use lattice_qcd_rs::error::ImplementationError;
/// fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
/// let point = LatticePoint::<4>::new([1, 0, 2, 0].into());
/// assert_eq!(
/// lattice.link_canonical(point, DirectionEnum::XNeg.into()),
/// LatticeLinkCanonical::new(LatticePoint::new([0, 0, 2, 0].into()), DirectionEnum::XPos.into()).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// lattice.link_canonical(point, DirectionEnum::XPos.into()),
/// LatticeLinkCanonical::new(LatticePoint::new([1, 0, 2, 0].into()), DirectionEnum::XPos.into()).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// lattice.link_canonical(point, DirectionEnum::YNeg.into()),
/// LatticeLinkCanonical::new(LatticePoint::new([1, 3, 2, 0].into()), DirectionEnum::YPos.into()).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// # Ok(())
/// # }
/// ```
pub fn link_canonical(
&self,
pos: LatticePoint<D>,
dir: Direction<D>,
) -> LatticeLinkCanonical<D> {
let mut pos_link = pos;
if !dir.is_positive() {
pos_link = self.add_point_direction(pos_link, &dir);
}
for i in 0..pos.len() {
pos_link[i] %= self.dim();
}
LatticeLinkCanonical::new(pos_link, dir.to_positive()).unwrap()
}
/// Return a link build from `pos` and `dir`.
///
/// It is similar to [`LatticeLink::new`]. It however enforce that the point is inside the bounds.
/// If it is not, it will use the modulus of the bound.
/// # Example
/// ```
/// # use lattice_qcd_rs::{lattice::{LatticeCyclic, Direction, LatticePoint}, error::ImplementationError};
/// # use nalgebra::SVector;
/// fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let l = LatticeCyclic::<3>::new(1_f64, 4)?;
/// let dir = Direction::new(0, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let pt = LatticePoint::new(SVector::<_, 3>::new(1, 2, 5));
/// let link = l.link(pt, dir);
/// assert_eq!(
/// *link.pos(),
/// LatticePoint::new(SVector::<_, 3>::new(1, 2, 1))
/// );
/// # Ok(())
/// # }
/// ```
pub fn link(&self, pos: LatticePoint<D>, dir: Direction<D>) -> LatticeLink<D> {
let mut pos_link = LatticePoint::new_zero();
for i in 0..pos.len() {
pos_link[i] = pos[i] % self.dim();
}
LatticeLink::new(pos_link, dir)
}
/// Transform a [`LatticeLink`] into a [`LatticeLinkCanonical`].
///
/// See [`LatticeCyclic::link_canonical`] and [`LatticeLinkCanonical`].
pub fn into_canonical(&self, l: LatticeLink<D>) -> LatticeLinkCanonical<D> {
self.link_canonical(*l.pos(), *l.dir())
}
/// Get the number of points in a single direction.
///
/// use [`LatticeCyclic::number_of_points`] for the total number of points.
/// Not to confuse with [`LatticeCyclic::dim_st`] which is the dimension of space-time.
pub const fn dim(&self) -> usize {
self.dim
}
/// Get an Iterator over all points of the lattice.
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticeCyclic;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// for i in [4, 8, 16, 32, 64].into_iter() {
/// let l = LatticeCyclic::<4>::new(1_f64, i)?;
/// assert_eq!(l.get_points().size_hint().0, l.number_of_points());
/// }
/// # Ok(())
/// # }
/// ```
pub fn get_points(&self) -> IteratorLatticePoint<'_, D> {
IteratorLatticePoint::new(self, LatticePoint::new_zero())
}
/// Get an Iterator over all canonical link of the lattice.
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticeCyclic;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// for i in [4, 8, 16, 32, 64].into_iter() {
/// let l = LatticeCyclic::<4>::new(1_f64, i)?;
/// assert_eq!(
/// l.get_links().size_hint().0,
/// l.number_of_canonical_links_space()
/// );
/// }
/// # Ok(())
/// # }
/// ```
pub fn get_links(&self) -> IteratorLatticeLinkCanonical<'_, D> {
return IteratorLatticeLinkCanonical::new(
self,
self.link_canonical(
LatticePoint::new_zero(),
*Direction::positive_directions().first().unwrap(),
),
);
}
/// create a new lattice with `size` the lattice size parameter, and `dim` the number of
/// points in each spatial dimension.
///
/// # Errors
/// Size should be greater than 0 and dim greater or equal to 2, otherwise return an error.
pub fn new(size: Real, dim: usize) -> Result<Self, LatticeInitializationError> {
if D == 0 {
return Err(LatticeInitializationError::ZeroDimension);
}
if size <= 0_f64 || size.is_nan() || size.is_infinite() {
return Err(LatticeInitializationError::NonPositiveSize);
}
if dim < 2 {
return Err(LatticeInitializationError::DimTooSmall);
}
Ok(Self { size, dim })
}
/// Total number of canonical links oriented in space for a set time.
///
/// Basically the number of element return by [`LatticeCyclic::get_links`].
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticeCyclic;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let l = LatticeCyclic::<4>::new(1_f64, 8)?;
/// assert_eq!(l.number_of_canonical_links_space(), 8_usize.pow(4) * 4);
///
/// let l = LatticeCyclic::<4>::new(1_f64, 16)?;
/// assert_eq!(l.number_of_canonical_links_space(), 16_usize.pow(4) * 4);
///
/// let l = LatticeCyclic::<3>::new(1_f64, 8)?;
/// assert_eq!(l.number_of_canonical_links_space(), 8_usize.pow(3) * 3);
/// # Ok(())
/// # }
/// ```
pub fn number_of_canonical_links_space(&self) -> usize {
self.number_of_points() * D
}
/// Total number of point in the lattice for a set time.
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticeCyclic;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let l = LatticeCyclic::<4>::new(1_f64, 8)?;
/// assert_eq!(l.number_of_points(), 8_usize.pow(4));
///
/// let l = LatticeCyclic::<4>::new(1_f64, 16)?;
/// assert_eq!(l.number_of_points(), 16_usize.pow(4));
///
/// let l = LatticeCyclic::<3>::new(1_f64, 8)?;
/// assert_eq!(l.number_of_points(), 8_usize.pow(3));
/// # Ok(())
/// # }
/// ```
pub fn number_of_points(&self) -> usize {
self.dim().pow(D.try_into().unwrap())
}
/// Return the lattice size factor.
pub const fn size(&self) -> Real {
self.size
}
/// Get the next point in the lattice following the direction `dir`.
/// It follows the Cyclic property of the lattice.
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{LatticeCyclic, DirectionEnum, LatticePoint};
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
/// let point = LatticePoint::<4>::from([1, 0, 2, 0]);
/// assert_eq!(
/// lattice.add_point_direction(point, &DirectionEnum::XPos.into()),
/// LatticePoint::from([2, 0, 2, 0])
/// );
/// // In the following case we get [_, 3, _, _] because `dim = 4`, and this lattice is Cyclic.
/// assert_eq!(
/// lattice.add_point_direction(point, &DirectionEnum::YNeg.into()),
/// LatticePoint::from([1, 3, 2, 0])
/// );
/// # Ok(())
/// # }
/// ```
pub fn add_point_direction(
&self,
point: LatticePoint<D>,
dir: &Direction<D>,
) -> LatticePoint<D> {
self.add_point_direction_n(point, dir, 1)
}
/// Returns the point given y moving `shift_number` times in direction `dir` from position `point`.
/// It follows the Cyclic property of the lattice.
///
/// It is equivalent of doing [`Self::add_point_direction`] n times.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{LatticeCyclic, DirectionEnum, LatticePoint};
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
/// let point = LatticePoint::<4>::from([1, 0, 2, 0]);
/// assert_eq!(
/// lattice.add_point_direction_n(point, &DirectionEnum::XPos.into(), 2),
/// LatticePoint::from([3, 0, 2, 0])
/// );
/// // In the following case we get [_, 1, _, _] because `dim = 4`, and this lattice is Cyclic.
/// assert_eq!(
/// lattice.add_point_direction_n(point, &DirectionEnum::YNeg.into(), 3),
/// LatticePoint::from([1, 1, 2, 0])
/// );
/// # Ok(())
/// # }
/// ```
pub fn add_point_direction_n(
&self,
mut point: LatticePoint<D>,
dir: &Direction<D>,
shift_number: usize,
) -> LatticePoint<D> {
let shift_number = shift_number % self.dim(); // we ensure that shift_number < % self.dim()
if dir.is_positive() {
point[dir.index()] = (point[dir.index()] + shift_number) % self.dim();
}
else {
let dir_pos = dir.to_positive();
if point[dir_pos.index()] < shift_number {
point[dir_pos.index()] = self.dim() - (shift_number - point[dir_pos.index()]);
}
else {
point[dir_pos.index()] = (point[dir_pos.index()] - shift_number) % self.dim();
}
}
point
}
/// Returns whether the number of canonical link is the same as the length of `links`.
pub fn has_compatible_length_links(&self, links: &LinkMatrix) -> bool {
self.number_of_canonical_links_space() == links.len()
}
/// Returns wether the number of point is the same as the length of `e_field`.
pub fn has_compatible_length_e_field(&self, e_field: &EField<D>) -> bool {
self.number_of_points() == e_field.len()
}
/// Returns the length is compatible for both `links` and `e_field`.
/// See [`Self::has_compatible_length_links`] and [`Self::has_compatible_length_e_field`].
pub fn has_compatible_length(&self, links: &LinkMatrix, e_field: &EField<D>) -> bool {
self.has_compatible_length_links(links) && self.has_compatible_length_e_field(e_field)
}
}
impl<const D: usize> std::fmt::Display for LatticeCyclic<D> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(
f,
"Cyclic lattice with {}^{} points and spacing {}",
self.dim, D, self.size
)
}
}
/// Iterator over [`LatticeLinkCanonical`] associated to a particular [`LatticeCyclic`].
#[derive(Clone, Debug, PartialEq)]
pub struct IteratorLatticeLinkCanonical<'a, const D: usize> {
lattice: &'a LatticeCyclic<D>,
element: Option<LatticeLinkCanonical<D>>,
}
impl<'a, const D: usize> IteratorLatticeLinkCanonical<'a, D> {
/// create a new iterator. The first [`IteratorLatticeLinkCanonical::next()`] will return `first_el`.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorLatticeLinkCanonical, LatticeCyclic, LatticeLinkCanonical, LatticePoint, DirectionEnum};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
/// let first_el = LatticeLinkCanonical::<4>::new(LatticePoint::from([1, 0, 2, 0]), DirectionEnum::YPos.into()).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let mut iter = IteratorLatticeLinkCanonical::new(&lattice, first_el);
/// assert_eq!(iter.next().ok_or(ImplementationError::OptionWithUnexpectedNone)?, first_el);
/// # Ok(())
/// # }
/// ```
pub const fn new(lattice: &'a LatticeCyclic<D>, first_el: LatticeLinkCanonical<D>) -> Self {
Self {
lattice,
element: Some(first_el),
}
}
}
impl<'a, const D: usize> Iterator for IteratorLatticeLinkCanonical<'a, D> {
type Item = LatticeLinkCanonical<D>;
// TODO improve
fn next(&mut self) -> Option<Self::Item> {
let previous_el = self.element;
if let Some(ref mut element) = self.element {
let mut iter_dir = IteratorDirection::<D, true>::new(Some(element.dir)).unwrap();
let new_dir = iter_dir.next();
match new_dir {
Some(dir) => element.set_dir(dir),
None => {
element.set_dir(Direction::new(0, true).unwrap());
let mut iter = IteratorLatticePoint::new(self.lattice, *element.pos());
match iter.nth(1) {
// get the second element
Some(array) => *element.pos_mut() = array,
None => {
self.element = None;
return previous_el;
}
}
}
}
}
previous_el
}
fn size_hint(&self) -> (usize, Option<usize>) {
match self.element {
None => (0, Some(0)),
Some(element) => {
let val =
self.lattice.number_of_canonical_links_space() - element.to_index(self.lattice);
(val, Some(val))
}
}
}
}
impl<'a, const D: usize> FusedIterator for IteratorLatticeLinkCanonical<'a, D> {}
impl<'a, const D: usize> ExactSizeIterator for IteratorLatticeLinkCanonical<'a, D> {}
/// Enum for internal use of iterator. It store the previous element returned by `next`
#[derive(Clone, Debug, Copy, Hash, PartialOrd, Ord, PartialEq, Eq)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub enum IteratorElement<T> {
/// First element of the iterator
FirstElement,
/// An element of the iterator
Element(T),
/// The Iterator is exhausted
LastElement,
}
impl<T: std::fmt::Display> std::fmt::Display for IteratorElement<T> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Self::FirstElement => write!(f, "first element"),
Self::Element(t) => write!(f, "element {}", t),
Self::LastElement => write!(f, "last element"),
}
}
}
impl<T> Default for IteratorElement<T> {
fn default() -> Self {
Self::FirstElement
}
}
impl<T> From<IteratorElement<T>> for Option<T> {
fn from(element: IteratorElement<T>) -> Self {
match element {
IteratorElement::Element(el) => Some(el),
IteratorElement::LastElement | IteratorElement::FirstElement => None,
}
}
}
/// Iterator over [`Direction`] with the same sign.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorDirection, Direction, IteratorElement};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let mut iter = IteratorDirection::<4, true>::new(None)
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(0, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(1, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(2, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(3, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(iter.next(), None);
/// assert_eq!(iter.next(), None);
/// # Ok(())
/// # }
/// ```
#[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Hash)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct IteratorDirection<const D: usize, const IS_POSITIVE_DIRECTION: bool> {
element: IteratorElement<Direction<D>>,
}
impl<const D: usize, const IS_POSITIVE_DIRECTION: bool>
IteratorDirection<D, IS_POSITIVE_DIRECTION>
{
/// Create an iterator where the first element upon calling [`Self::next`] is the direction
/// after `element`. Giving `None` results in the first element being the direction with index 0
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorDirection, Direction, IteratorElement};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let mut iter = IteratorDirection::<4, true>::new(None)
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(0, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(1, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// let element =
/// Direction::<4>::new(0, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let mut iter = IteratorDirection::<4, false>::new(Some(element))
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(2, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// # Ok(())
/// # }
/// ```
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorDirection, Direction, IteratorElement};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let iter = IteratorDirection::<0, true>::new(None);
/// // 0 is invalid
/// assert!(iter.is_none());
/// let iter = IteratorDirection::<4, true>::new(Some(
/// Direction::new(2, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// ));
/// // the sign of the direction is invalid
/// assert!(iter.is_none());
/// # Ok(())
/// # }
/// ```
pub const fn new(element: Option<Direction<D>>) -> Option<Self> {
match element {
None => Self::new_from_element(IteratorElement::FirstElement),
Some(dir) => Self::new_from_element(IteratorElement::Element(dir)),
}
}
/// create a new iterator. The first [`IteratorLatticeLinkCanonical::next()`] the element just after the one given
/// or the first element if [`IteratorElement::FirstElement`] is given.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorDirection, Direction, IteratorElement};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let mut iter = IteratorDirection::<4, true>::new_from_element(IteratorElement::FirstElement)
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(0, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(1, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// let element =
/// Direction::<4>::new(0, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let mut iter =
/// IteratorDirection::<4, false>::new_from_element(IteratorElement::Element(element))
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// iter.next()
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?,
/// Direction::new(2, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// # Ok(())
/// # }
/// ```
pub const fn new_from_element(element: IteratorElement<Direction<D>>) -> Option<Self> {
if D == 0 {
return None;
}
if let IteratorElement::Element(ref el) = element {
if el.is_positive() != IS_POSITIVE_DIRECTION {
return None;
}
}
Some(Self { element })
}
}
impl<const D: usize, const IS_POSITIVE_DIRECTION: bool> Iterator
for IteratorDirection<D, IS_POSITIVE_DIRECTION>
{
type Item = Direction<D>;
fn next(&mut self) -> Option<Self::Item> {
let next_element = match self.element {
IteratorElement::FirstElement => {
IteratorElement::Element(Direction::new(0, IS_POSITIVE_DIRECTION).unwrap())
}
IteratorElement::Element(ref dir) => {
if let Some(dir) = Direction::new(dir.index() + 1, IS_POSITIVE_DIRECTION) {
IteratorElement::Element(dir)
}
else {
IteratorElement::LastElement
}
}
IteratorElement::LastElement => IteratorElement::LastElement,
};
self.element = next_element;
next_element.into()
}
fn size_hint(&self) -> (usize, Option<usize>) {
let size = match self.element {
IteratorElement::FirstElement => D,
IteratorElement::Element(ref dir) => D - (dir.index() + 1),
IteratorElement::LastElement => 0,
};
(size, Some(size))
}
}
impl<const D: usize, const IS_POSITIVE_DIRECTION: bool> FusedIterator
for IteratorDirection<D, IS_POSITIVE_DIRECTION>
{
}
impl<const D: usize, const IS_POSITIVE_DIRECTION: bool> ExactSizeIterator
for IteratorDirection<D, IS_POSITIVE_DIRECTION>
{
}
/// Iterator over [`LatticePoint`]
#[derive(Clone, Debug, PartialEq)]
pub struct IteratorLatticePoint<'a, const D: usize> {
lattice: &'a LatticeCyclic<D>,
element: Option<LatticePoint<D>>,
}
impl<'a, const D: usize> IteratorLatticePoint<'a, D> {
/// create a new iterator. The first [`IteratorLatticePoint::next()`] will return `first_el`.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{IteratorLatticePoint, LatticeCyclic, LatticePoint, Direction};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let lattice = LatticeCyclic::new(1_f64, 4)?;
/// let first_el = LatticePoint::from([1, 0, 2, 0]);
/// let mut iter = IteratorLatticePoint::new(&lattice, first_el);
/// assert_eq!(iter.next().ok_or(ImplementationError::OptionWithUnexpectedNone)?, first_el);
/// # Ok(())
/// # }
/// ```
pub const fn new(lattice: &'a LatticeCyclic<D>, first_el: LatticePoint<D>) -> Self {
Self {
lattice,
element: Some(first_el),
}
}
}
impl<'a, const D: usize> Iterator for IteratorLatticePoint<'a, D> {
type Item = LatticePoint<D>;
// TODO improve
fn next(&mut self) -> Option<Self::Item> {
let previous_el = self.element;
if let Some(ref mut el) = &mut self.element {
el[0] += 1;
for i in 0..el.len() {
while el[i] >= self.lattice.dim() {
if i < el.len() - 1 {
// every element except the last one
el[i + 1] += 1;
}
else {
self.element = None;
return previous_el;
}
el[i] -= self.lattice.dim();
}
}
}
previous_el
}
fn size_hint(&self) -> (usize, Option<usize>) {
match self.element {
None => (0, Some(0)),
Some(element) => {
let val = self.lattice.number_of_points() - element.to_index(self.lattice);
(val, Some(val))
}
}
}
}
impl<'a, const D: usize> FusedIterator for IteratorLatticePoint<'a, D> {}
impl<'a, const D: usize> ExactSizeIterator for IteratorLatticePoint<'a, D> {}
/// Represents point on a (any) lattice.
#[derive(Clone, Debug, Copy, PartialEq, Eq, PartialOrd, Hash)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct LatticePoint<const D: usize> {
data: na::SVector<usize, D>,
}
impl<const D: usize> LatticePoint<D> {
/// Create a new lattice point.
///
/// It can be outside a lattice.
pub const fn new(data: SVector<usize, D>) -> Self {
Self { data }
}
/// Create a point at the origin
pub fn new_zero() -> Self {
Self {
data: SVector::zeros(),
}
}
/// Create a point using the closure generate elements with the index as input.
///
/// See [`nalgebra::base::Matrix::from_fn`].
pub fn from_fn<F>(mut f: F) -> Self
where
F: FnMut(usize) -> usize,
{
Self::new(SVector::from_fn(|index, _| f(index)))
}
/// Number of elements in [`LatticePoint`]. This is `D`.
#[allow(clippy::unused_self)]
pub const fn len(&self) -> usize {
// this is in order to have a const function.
// we could have called self.data.len()
D
}
/// Return if LatticePoint contain no data. True when the dimension is 0, false otherwise.
#[allow(clippy::unused_self)]
pub const fn is_empty(&self) -> bool {
D == 0
}
/// Get an iterator on the data.
pub fn iter(&self) -> impl Iterator<Item = &usize> + ExactSizeIterator + FusedIterator {
self.data.iter()
}
/// Get an iterator on the data as mutable.
pub fn iter_mut(
&mut self,
) -> impl Iterator<Item = &mut usize> + ExactSizeIterator + FusedIterator {
self.data.iter_mut()
}
/// Get the point as as [`nalgebra::SVector<usize, D>`]
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticePoint;
/// #
/// # let point = LatticePoint::<4>::default();
/// let max = point.as_svector().max();
/// let min = point.as_ref().min();
/// ```
pub const fn as_svector(&self) -> &SVector<usize, D> {
&self.data
}
/// Get the point as a mut ref to [`nalgebra::SVector<usize, D>`]
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::LatticePoint;
/// #
/// # let mut point = LatticePoint::<4>::new_zero();
/// point.as_svector_mut()[2] = 2;
/// point.as_mut()[0] = 1;
/// ```
pub fn as_svector_mut(&mut self) -> &mut SVector<usize, D> {
&mut self.data
}
}
impl<const D: usize> Default for LatticePoint<D> {
fn default() -> Self {
Self::new_zero()
}
}
impl<const D: usize> std::fmt::Display for LatticePoint<D> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "{}", self.data)
}
}
impl<'a, const D: usize> IntoIterator for &'a LatticePoint<D> {
type IntoIter = <&'a SVector<usize, D> as IntoIterator>::IntoIter;
type Item = &'a usize;
fn into_iter(self) -> Self::IntoIter {
self.data.iter()
}
}
impl<'a, const D: usize> IntoIterator for &'a mut LatticePoint<D> {
type IntoIter = <&'a mut SVector<usize, D> as IntoIterator>::IntoIter;
type Item = &'a mut usize;
fn into_iter(self) -> Self::IntoIter {
self.data.iter_mut()
}
}
impl<const D: usize> Index<usize> for LatticePoint<D> {
type Output = usize;
/// Get the element at position `pos`
/// # Panic
/// Panics if the position is out of bound
/// ```should_panic
/// # use lattice_qcd_rs::lattice::LatticePoint;
/// let point = LatticePoint::new([0; 4].into());
/// let _ = point[4];
/// ```
fn index(&self, pos: usize) -> &Self::Output {
&self.data[pos]
}
}
impl<const D: usize> IndexMut<usize> for LatticePoint<D> {
/// Get the element at position `pos`
/// # Panic
/// Panics if the position is out of bound
/// ```should_panic
/// # use lattice_qcd_rs::lattice::LatticePoint;
/// let mut point = LatticePoint::new([0; 4].into());
/// point[4] += 1;
/// ```
fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
&mut self.data[pos]
}
}
impl<T, const D: usize> From<T> for LatticePoint<D>
where
SVector<usize, D>: From<T>,
{
fn from(data: T) -> Self {
LatticePoint::new(SVector::from(data))
}
}
impl<const D: usize> From<LatticePoint<D>> for [usize; D]
where
SVector<usize, D>: Into<[usize; D]>,
{
fn from(data: LatticePoint<D>) -> [usize; D] {
data.data.into()
}
}
impl<const D: usize> AsRef<SVector<usize, D>> for LatticePoint<D> {
fn as_ref(&self) -> &SVector<usize, D> {
self.as_svector()
}
}
impl<const D: usize> AsMut<SVector<usize, D>> for LatticePoint<D> {
fn as_mut(&mut self) -> &mut SVector<usize, D> {
self.as_svector_mut()
}
}
/// Trait to convert an element on a lattice to an [`usize`].
///
/// Used mainly to index field on the lattice using [`std::vec::Vec`]
// TODO change name ? make it less confusing
pub trait LatticeElementToIndex<const D: usize> {
/// Given a lattice return an index from the element
fn to_index(&self, l: &LatticeCyclic<D>) -> usize;
}
impl<const D: usize> LatticeElementToIndex<D> for LatticePoint<D> {
fn to_index(&self, l: &LatticeCyclic<D>) -> usize {
self.iter()
.enumerate()
.map(|(index, pos)| (pos % l.dim()) * l.dim().pow(index.try_into().unwrap()))
.sum()
}
}
impl<const D: usize> LatticeElementToIndex<D> for Direction<D> {
/// equivalent to [`Direction::to_index()`]
fn to_index(&self, _: &LatticeCyclic<D>) -> usize {
self.index()
}
}
impl<const D: usize> LatticeElementToIndex<D> for LatticeLinkCanonical<D> {
fn to_index(&self, l: &LatticeCyclic<D>) -> usize {
self.pos().to_index(l) * D + self.dir().index()
}
}
/// This is just the identity.
///
/// It is implemented for compatibility reason
/// such that function that require a [`LatticeElementToIndex`] can also accept [`usize`].
impl<const D: usize> LatticeElementToIndex<D> for usize {
/// return self
fn to_index(&self, _l: &LatticeCyclic<D>) -> usize {
*self
}
}
/// A canonical link of a lattice. It contain a position and a direction.
///
/// The direction should always be positive.
/// By itself the link does not store data about the lattice. Hence most function require a [`LatticeCyclic`].
/// It also means that there is no guarantee that the object is inside a lattice.
/// You can use modulus over the elements to use inside a lattice.
///
/// This object can be used to safely index in a [`std::collections::HashMap`]
#[derive(Clone, Debug, PartialEq, Eq, Hash, Copy)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct LatticeLinkCanonical<const D: usize> {
from: LatticePoint<D>,
dir: Direction<D>,
}
impl<const D: usize> LatticeLinkCanonical<D> {
/// Try create a LatticeLinkCanonical. If the dir is negative it fails.
///
/// To guaranty creating an element see [LatticeCyclic::link_canonical].
/// The creation of an element this ways does not guaranties that the element is inside a lattice.
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{LatticeLinkCanonical, LatticePoint, DirectionEnum};
/// let l = LatticeLinkCanonical::new(LatticePoint::new([0; 4].into()), DirectionEnum::XNeg.into());
/// assert_eq!(l, None);
///
/// let l = LatticeLinkCanonical::new(LatticePoint::new([0; 4].into()), DirectionEnum::XPos.into());
/// assert!(l.is_some());
/// ```
pub const fn new(from: LatticePoint<D>, dir: Direction<D>) -> Option<Self> {
if dir.is_negative() {
return None;
}
Some(Self { from, dir })
}
/// Position of the link.
pub const fn pos(&self) -> &LatticePoint<D> {
&self.from
}
/// Get a mutable reference on the position of the link.
pub fn pos_mut(&mut self) -> &mut LatticePoint<D> {
&mut self.from
}
/// Direction of the link.
pub const fn dir(&self) -> &Direction<D> {
&self.dir
}
/// Set the direction to dir
///
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::{LatticeLinkCanonical, LatticePoint, DirectionEnum};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let mut lattice_link_canonical =
/// LatticeLinkCanonical::new(LatticePoint::new([0; 4].into()), DirectionEnum::YPos.into())
/// .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// lattice_link_canonical.set_dir(DirectionEnum::XPos.into());
/// assert_eq!(*lattice_link_canonical.dir(), DirectionEnum::XPos.into());
/// # Ok(())
/// # }
/// ```
///
/// # Panic
/// panic if a negative direction is given.
/// ```should_panic
/// # use lattice_qcd_rs::lattice::{LatticeLinkCanonical, LatticePoint, DirectionEnum};
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// # let mut lattice_link_canonical = LatticeLinkCanonical::new(LatticePoint::new([0; 4].into()), DirectionEnum::XPos.into()).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// lattice_link_canonical.set_dir(DirectionEnum::XNeg.into()); // Panics !
/// # Ok(())
/// # }
/// ```
pub fn set_dir(&mut self, dir: Direction<D>) {
if dir.is_negative() {
panic!("Cannot set a negative direction to a canonical link.");
}
self.dir = dir;
}
/// Set the direction using positive direction. i.e. if a direction `-x` is passed
/// the direction assigned will be `+x`.
///
/// This is equivalent to `link.set_dir(dir.to_positive())`.
pub fn set_dir_positive(&mut self, dir: Direction<D>) {
self.dir = dir.to_positive();
}
}
impl<const D: usize> std::fmt::Display for LatticeLinkCanonical<D> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(
f,
"canonical link [position {}, direction {}]",
self.from, self.dir
)
}
}
impl<const D: usize> From<LatticeLinkCanonical<D>> for LatticeLink<D> {
fn from(l: LatticeLinkCanonical<D>) -> Self {
LatticeLink::new(l.from, l.dir)
}
}
impl<const D: usize> From<&LatticeLinkCanonical<D>> for LatticeLink<D> {
fn from(l: &LatticeLinkCanonical<D>) -> Self {
LatticeLink::new(l.from, l.dir)
}
}
/// A lattice link, contrary to [`LatticeLinkCanonical`] the direction can be negative.
///
/// This means that multiple link can be equivalent but does not have the same data
/// and therefore hash (hopefully).
///
/// By itself the link does not store data about the lattice. Hence most function require a [`LatticeCyclic`].
/// It also means that there is no guarantee that the object is inside a lattice.
/// You can use modulus over the elements to use inside a lattice.
#[derive(Clone, Debug, PartialEq, Eq, Hash, Copy)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct LatticeLink<const D: usize> {
from: LatticePoint<D>,
dir: Direction<D>,
}
impl<const D: usize> LatticeLink<D> {
/// Create a link from position `from` and direction `dir`.
pub const fn new(from: LatticePoint<D>, dir: Direction<D>) -> Self {
Self { from, dir }
}
/// Get the position of the link.
pub const fn pos(&self) -> &LatticePoint<D> {
&self.from
}
/// Get a mutable reference to the position of the link.
pub fn pos_mut(&mut self) -> &mut LatticePoint<D> {
&mut self.from
}
/// Get the direction of the link.
pub const fn dir(&self) -> &Direction<D> {
&self.dir
}
/// Get a mutable reference to the direction of the link.
pub fn dir_mut(&mut self) -> &mut Direction<D> {
&mut self.dir
}
/// Get if the direction of the link is positive.
pub const fn is_dir_positive(&self) -> bool {
self.dir.is_positive()
}
/// Get if the direction of the link is negative.
pub const fn is_dir_negative(&self) -> bool {
self.dir.is_negative()
}
}
impl<const D: usize> std::fmt::Display for LatticeLink<D> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "link [position {}, direction {}]", self.from, self.dir)
}
}
/// Represent a cardinal direction
#[derive(Clone, Debug, PartialEq, Eq, Hash, Copy)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub struct Direction<const D: usize> {
index_dir: usize,
is_positive: bool,
}
impl<const D: usize> Direction<D> {
/// New direction from a direction as an idex and wether it is in the positive direction.
pub const fn new(index_dir: usize, is_positive: bool) -> Option<Self> {
if index_dir >= D {
return None;
}
Some(Self {
index_dir,
is_positive,
})
}
/// List of all positives directions.
/// This is very slow use [`Self::positive_directions`] instead.
pub fn positives_vec() -> Vec<Self> {
let mut x = Vec::with_capacity(D);
for i in 0..D {
x.push(Self::new(i, true).unwrap());
}
x
}
/// List all directions.
/// This is very slow use [`DirectionList::directions`] instead.
pub fn directions_vec() -> Vec<Self> {
let mut x = Vec::with_capacity(2 * D);
for i in 0..D {
x.push(Self::new(i, true).unwrap());
x.push(Self::new(i, false).unwrap());
}
x
}
// TODO add const function for all direction once operation on const generic are added
/// Get all direction with the sign `IS_POSITIVE`.
pub const fn directions_array<const IS_POSITIVE: bool>() -> [Self; D] {
// TODO use unsafe code to avoid the allocation
let mut i = 0_usize;
let mut array = [Direction {
index_dir: 0,
is_positive: IS_POSITIVE,
}; D];
while i < D {
array[i] = Direction {
index_dir: i,
is_positive: IS_POSITIVE,
};
i += 1;
}
array
}
/// Get all negative direction
pub const fn negative_directions() -> [Self; D] {
Self::directions_array::<false>()
}
/// Get all positive direction
pub const fn positive_directions() -> [Self; D] {
Self::directions_array::<true>()
}
/// Get if the position is positive.
pub const fn is_positive(&self) -> bool {
self.is_positive
}
/// Get if the position is Negative. see [`Direction::is_positive`]
pub const fn is_negative(&self) -> bool {
!self.is_positive()
}
/// Return the positive direction associated, for example `-x` gives `+x`
/// and `+x` gives `+x`.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::Direction;
/// assert_eq!(
/// Direction::<4>::new(1, false).unwrap().to_positive(),
/// Direction::<4>::new(1, true).unwrap()
/// );
/// assert_eq!(
/// Direction::<4>::new(1, true).unwrap().to_positive(),
/// Direction::<4>::new(1, true).unwrap()
/// );
/// ```
pub const fn to_positive(mut self) -> Self {
self.is_positive = true;
self
}
/// Get a index associated to the direction.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::Direction;
/// assert_eq!(Direction::<4>::new(1, false).unwrap().index(), 1);
/// ```
pub const fn index(&self) -> usize {
self.index_dir
}
/// Convert the direction into a vector of norm `a`;
pub fn to_vector(self, a: f64) -> SVector<Real, D> {
self.to_unit_vector() * a
}
/// Returns the dimension
pub const fn dim() -> usize {
D
}
/// Convert the direction into a vector of norm `1`;
pub fn to_unit_vector(self) -> SVector<Real, D> {
let mut v = SVector::zeros();
v[self.index_dir] = 1_f64;
v
}
/// Find the direction the vector point the most.
/// For a zero vector return [`DirectionEnum::XPos`].
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::Direction;
/// # use lattice_qcd_rs::error::ImplementationError;
/// # fn main() -> Result<(), Box< dyn std::error::Error>> {
/// assert_eq!(
/// Direction::from_vector(&nalgebra::Vector4::new(1_f64, 0_f64, 0_f64, 0_f64)),
/// Direction::<4>::new(0, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// Direction::from_vector(&nalgebra::Vector4::new(0_f64, -1_f64, 0_f64, 0_f64)),
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// assert_eq!(
/// Direction::from_vector(&nalgebra::Vector4::new(0.5_f64, 1_f64, 0_f64, 2_f64)),
/// Direction::<4>::new(3, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?
/// );
/// # Ok(())
/// # }
/// ```
pub fn from_vector(v: &SVector<Real, D>) -> Self {
let mut max = 0_f64;
let mut index_max: usize = 0;
let mut is_positive = true;
for (i, dir) in Self::positive_directions().iter().enumerate() {
let scalar_prod = v.dot(&dir.to_vector(1_f64));
if scalar_prod.abs() > max {
max = scalar_prod.abs();
index_max = i;
is_positive = scalar_prod > 0_f64;
}
}
Self::new(index_max, is_positive).expect("Unreachable")
}
}
// TODO default when condition on const generic is available
/*
impl<const D: usize> Default for LatticeLinkCanonical<D>
where
Direction<D>: Default,
{
fn default() -> Self {
Self {
from: LatticePoint::default(),
dir: Direction::default(),
}
}
}
impl<const D: usize> Default for LatticeLink<D>
where
Direction<D>: Default,
{
fn default() -> Self {
Self {
from: LatticePoint::default(),
dir: Direction::default(),
}
}
}
*/
impl<const D: usize> std::fmt::Display for Direction<D> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(
f,
"[index {}, is positive {}]",
self.index(),
self.is_positive()
)
}
}
/// List all possible direction
pub trait DirectionList: Sized {
/// List all directions.
fn directions() -> &'static [Self];
/// List all positive directions.
fn positive_directions() -> &'static [Self];
}
implement_direction_list!();
implement_direction_from!();
/// Partial ordering is set as follows: two directions can be compared if they have the same index
/// or the same direction sign. In the first case a positive direction is greater than a negative direction
/// In the latter case the ordering is done on the index.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::Direction;
/// # use lattice_qcd_rs::error::ImplementationError;
/// use std::cmp::Ordering;
/// # fn main() -> Result<(), Box< dyn std::error::Error>> {
///
/// let dir_1 =
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let dir_2 =
/// Direction::<4>::new(1, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert!(dir_1 < dir_2);
/// assert_eq!(dir_1.partial_cmp(&dir_2), Some(Ordering::Less));
/// //--------
/// let dir_3 =
/// Direction::<4>::new(3, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let dir_4 =
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert!(dir_3 > dir_4);
/// assert_eq!(dir_3.partial_cmp(&dir_4), Some(Ordering::Greater));
/// //--------
/// let dir_5 =
/// Direction::<4>::new(3, true).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let dir_6 =
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(dir_5.partial_cmp(&dir_6), None);
/// //--------
/// let dir_5 =
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// let dir_6 =
/// Direction::<4>::new(1, false).ok_or(ImplementationError::OptionWithUnexpectedNone)?;
/// assert_eq!(dir_5.partial_cmp(&dir_6), Some(Ordering::Equal));
/// # Ok(())
/// # }
/// ```
impl<const D: usize> PartialOrd for Direction<D> {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
if self == other {
Some(Ordering::Equal)
}
else if self.is_positive() == other.is_positive() {
self.index().partial_cmp(&other.index())
}
else if self.index() == other.index() {
self.is_positive().partial_cmp(&other.is_positive())
}
else {
None
}
}
}
impl<const D: usize> Neg for Direction<D> {
type Output = Self;
fn neg(mut self) -> Self::Output {
self.is_positive = !self.is_positive;
self
}
}
impl<const D: usize> Neg for &Direction<D> {
type Output = Direction<D>;
fn neg(self) -> Self::Output {
-*self
}
}
/// Return [`Direction::to_index`].
impl<const D: usize> From<Direction<D>> for usize {
fn from(d: Direction<D>) -> Self {
d.index()
}
}
/// Return [`Direction::to_index`].
impl<const D: usize> From<&Direction<D>> for usize {
fn from(d: &Direction<D>) -> Self {
d.index()
}
}
/// Return [`DirectionEnum::from_vector`].
impl<const D: usize> From<SVector<Real, D>> for Direction<D> {
fn from(v: SVector<Real, D>) -> Self {
Direction::from_vector(&v)
}
}
/// Return [`DirectionEnum::from_vector`].
impl<const D: usize> From<&SVector<Real, D>> for Direction<D> {
fn from(v: &SVector<Real, D>) -> Self {
Direction::<D>::from_vector(v)
}
}
/// Return [`Direction::to_unit_vector`].
impl<const D: usize> From<Direction<D>> for SVector<Real, D> {
fn from(d: Direction<D>) -> Self {
d.to_unit_vector()
}
}
/// Return [`Direction::to_unit_vector`].
impl<const D: usize> From<&Direction<D>> for SVector<Real, D> {
fn from(d: &Direction<D>) -> Self {
d.to_unit_vector()
}
}
impl From<DirectionEnum> for Direction<4> {
fn from(d: DirectionEnum) -> Self {
Self::new(d.to_index(), d.is_positive()).expect("unreachable")
}
}
impl From<&DirectionEnum> for Direction<4> {
fn from(d: &DirectionEnum) -> Self {
Self::new(d.to_index(), d.is_positive()).expect("unreachable")
}
}
/// Represent a cardinal direction
#[derive(Clone, Debug, PartialEq, Eq, Hash, Copy)]
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
pub enum DirectionEnum {
/// Positive x direction.
XPos,
/// Negative x direction.
XNeg,
/// Positive y direction.
YPos,
/// Negative y direction.
YNeg,
/// Positive z direction.
ZPos,
/// Negative z direction.
ZNeg,
/// Positive t direction.
TPos,
/// Negative t direction.
TNeg,
}
impl DirectionEnum {
/// List all directions.
pub const DIRECTIONS: [Self; 8] = [
DirectionEnum::XPos,
DirectionEnum::YPos,
DirectionEnum::ZPos,
DirectionEnum::TPos,
DirectionEnum::XNeg,
DirectionEnum::YNeg,
DirectionEnum::ZNeg,
DirectionEnum::TNeg,
];
/// List all spatial directions.
pub const DIRECTIONS_SPACE: [Self; 6] = [
DirectionEnum::XPos,
DirectionEnum::YPos,
DirectionEnum::ZPos,
DirectionEnum::XNeg,
DirectionEnum::YNeg,
DirectionEnum::ZNeg,
];
/// List of all positives directions.
pub const POSITIVES: [Self; 4] = [
DirectionEnum::XPos,
DirectionEnum::YPos,
DirectionEnum::ZPos,
DirectionEnum::TPos,
];
/// List spatial positive direction.
pub const POSITIVES_SPACE: [Self; 3] = [
DirectionEnum::XPos,
DirectionEnum::YPos,
DirectionEnum::ZPos,
];
/// Convert the direction into a vector of norm `a`;
pub fn to_vector(self, a: f64) -> Vector4<Real> {
self.to_unit_vector() * a
}
/// Convert the direction into a vector of norm `1`;
pub const fn to_unit_vector(self) -> Vector4<Real> {
match self {
DirectionEnum::XPos => Vector4::<Real>::new(1_f64, 0_f64, 0_f64, 0_f64),
DirectionEnum::XNeg => Vector4::<Real>::new(-1_f64, 0_f64, 0_f64, 0_f64),
DirectionEnum::YPos => Vector4::<Real>::new(0_f64, 1_f64, 0_f64, 0_f64),
DirectionEnum::YNeg => Vector4::<Real>::new(0_f64, -1_f64, 0_f64, 0_f64),
DirectionEnum::ZPos => Vector4::<Real>::new(0_f64, 0_f64, 1_f64, 0_f64),
DirectionEnum::ZNeg => Vector4::<Real>::new(0_f64, 0_f64, -1_f64, 0_f64),
DirectionEnum::TPos => Vector4::<Real>::new(0_f64, 0_f64, 0_f64, 1_f64),
DirectionEnum::TNeg => Vector4::<Real>::new(0_f64, 0_f64, 0_f64, -1_f64),
}
}
/// Get if the position is positive.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::DirectionEnum;
/// assert_eq!(DirectionEnum::XPos.is_positive(), true);
/// assert_eq!(DirectionEnum::TPos.is_positive(), true);
/// assert_eq!(DirectionEnum::YNeg.is_positive(), false);
/// ```
pub const fn is_positive(self) -> bool {
match self {
DirectionEnum::XPos
| DirectionEnum::YPos
| DirectionEnum::ZPos
| DirectionEnum::TPos => true,
DirectionEnum::XNeg
| DirectionEnum::YNeg
| DirectionEnum::ZNeg
| DirectionEnum::TNeg => false,
}
}
/// Get if the position is Negative. see [`DirectionEnum::is_positive`]
pub const fn is_negative(self) -> bool {
!self.is_positive()
}
/// Find the direction the vector point the most.
/// For a zero vector return [`DirectionEnum::XPos`].
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::DirectionEnum;
/// # extern crate nalgebra;
/// assert_eq!(
/// DirectionEnum::from_vector(&nalgebra::Vector4::new(1_f64, 0_f64, 0_f64, 0_f64)),
/// DirectionEnum::XPos
/// );
/// assert_eq!(
/// DirectionEnum::from_vector(&nalgebra::Vector4::new(0_f64, -1_f64, 0_f64, 0_f64)),
/// DirectionEnum::YNeg
/// );
/// assert_eq!(
/// DirectionEnum::from_vector(&nalgebra::Vector4::new(0.5_f64, 1_f64, 0_f64, 2_f64)),
/// DirectionEnum::TPos
/// );
/// ```
#[allow(clippy::needless_return)] // for readability
pub fn from_vector(v: &Vector4<Real>) -> Self {
let mut max = 0_f64;
let mut index_max: usize = 0;
let mut is_positive = true;
for i in 0..DirectionEnum::POSITIVES.len() {
let scalar_prod = v.dot(&DirectionEnum::POSITIVES[i].to_vector(1_f64));
if scalar_prod.abs() > max {
max = scalar_prod.abs();
index_max = i;
is_positive = scalar_prod > 0_f64;
}
}
match index_max {
0 => {
if is_positive {
return DirectionEnum::XPos;
}
else {
return DirectionEnum::XNeg;
}
}
1 => {
if is_positive {
return DirectionEnum::YPos;
}
else {
return DirectionEnum::YNeg;
}
}
2 => {
if is_positive {
return DirectionEnum::ZPos;
}
else {
return DirectionEnum::ZNeg;
}
}
3 => {
if is_positive {
return DirectionEnum::TPos;
}
else {
return DirectionEnum::TNeg;
}
}
_ => {
// the code should attain this code. and therefore panicking is not expected.
unreachable!("Implementation error : invalid index")
}
}
}
/// Return the positive direction associated, for example `-x` gives `+x`
/// and `+x` gives `+x`.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::DirectionEnum;
/// assert_eq!(DirectionEnum::XNeg.to_positive(), DirectionEnum::XPos);
/// assert_eq!(DirectionEnum::YPos.to_positive(), DirectionEnum::YPos);
/// ```
pub const fn to_positive(self) -> Self {
match self {
DirectionEnum::XNeg => DirectionEnum::XPos,
DirectionEnum::YNeg => DirectionEnum::YPos,
DirectionEnum::ZNeg => DirectionEnum::ZPos,
DirectionEnum::TNeg => DirectionEnum::TPos,
_ => self,
}
}
/// Get a index associated to the direction.
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::DirectionEnum;
/// assert_eq!(DirectionEnum::XPos.to_index(), 0);
/// assert_eq!(DirectionEnum::XNeg.to_index(), 0);
/// assert_eq!(DirectionEnum::YPos.to_index(), 1);
/// assert_eq!(DirectionEnum::YNeg.to_index(), 1);
/// assert_eq!(DirectionEnum::ZPos.to_index(), 2);
/// assert_eq!(DirectionEnum::ZNeg.to_index(), 2);
/// assert_eq!(DirectionEnum::TPos.to_index(), 3);
/// assert_eq!(DirectionEnum::TNeg.to_index(), 3);
/// ```
pub const fn to_index(self) -> usize {
match self {
DirectionEnum::XPos | DirectionEnum::XNeg => 0,
DirectionEnum::YPos | DirectionEnum::YNeg => 1,
DirectionEnum::ZPos | DirectionEnum::ZNeg => 2,
DirectionEnum::TPos | DirectionEnum::TNeg => 3,
}
}
}
/// Return [`DirectionEnum::XPos`]
impl Default for DirectionEnum {
///Return [`DirectionEnum::XPos`]
fn default() -> Self {
Self::XPos
}
}
impl std::fmt::Display for DirectionEnum {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
DirectionEnum::XPos => write!(f, "positive X direction"),
DirectionEnum::XNeg => write!(f, "negative X direction"),
DirectionEnum::YPos => write!(f, "positive Y direction"),
DirectionEnum::YNeg => write!(f, "negative Y direction"),
DirectionEnum::ZPos => write!(f, "positive Z direction"),
DirectionEnum::ZNeg => write!(f, "negative Z direction"),
DirectionEnum::TPos => write!(f, "positive T direction"),
DirectionEnum::TNeg => write!(f, "negative T direction"),
}
}
}
impl DirectionList for DirectionEnum {
fn directions() -> &'static [Self] {
&Self::DIRECTIONS
}
fn positive_directions() -> &'static [Self] {
&Self::POSITIVES
}
}
impl PartialOrd for DirectionEnum {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Direction::<4>::from(self).partial_cmp(&other.into())
}
}
/// Return the negative of a direction
/// # Example
/// ```
/// # use lattice_qcd_rs::lattice::DirectionEnum;
/// assert_eq!(-DirectionEnum::XNeg, DirectionEnum::XPos);
/// assert_eq!(-DirectionEnum::YPos, DirectionEnum::YNeg);
/// ```
impl Neg for DirectionEnum {
type Output = Self;
fn neg(self) -> Self::Output {
match self {
DirectionEnum::XPos => DirectionEnum::XNeg,
DirectionEnum::XNeg => DirectionEnum::XPos,
DirectionEnum::YPos => DirectionEnum::YNeg,
DirectionEnum::YNeg => DirectionEnum::YPos,
DirectionEnum::ZPos => DirectionEnum::ZNeg,
DirectionEnum::ZNeg => DirectionEnum::ZPos,
DirectionEnum::TPos => DirectionEnum::TNeg,
DirectionEnum::TNeg => DirectionEnum::TPos,
}
}
}
impl Neg for &DirectionEnum {
type Output = Self;
fn neg(self) -> Self::Output {
match self {
DirectionEnum::XPos => &DirectionEnum::XNeg,
DirectionEnum::XNeg => &DirectionEnum::XPos,
DirectionEnum::YPos => &DirectionEnum::YNeg,
DirectionEnum::YNeg => &DirectionEnum::YPos,
DirectionEnum::ZPos => &DirectionEnum::ZNeg,
DirectionEnum::ZNeg => &DirectionEnum::ZPos,
DirectionEnum::TPos => &DirectionEnum::TNeg,
DirectionEnum::TNeg => &DirectionEnum::TPos,
}
}
}
/// Return [`DirectionEnum::to_index`].
impl From<DirectionEnum> for usize {
fn from(d: DirectionEnum) -> Self {
d.to_index()
}
}
/// Return [`DirectionEnum::to_index`].
impl From<&DirectionEnum> for usize {
fn from(d: &DirectionEnum) -> Self {
d.to_index()
}
}
/// Return [`DirectionEnum::from_vector`].
impl From<Vector4<Real>> for DirectionEnum {
fn from(v: Vector4<Real>) -> Self {
DirectionEnum::from_vector(&v)
}
}
/// Return [`DirectionEnum::from_vector`].
impl From<&Vector4<Real>> for DirectionEnum {
fn from(v: &Vector4<Real>) -> Self {
DirectionEnum::from_vector(v)
}
}
/// Return [`DirectionEnum::to_unit_vector`].
impl From<DirectionEnum> for Vector4<Real> {
fn from(d: DirectionEnum) -> Self {
d.to_unit_vector()
}
}
/// Return [`DirectionEnum::to_unit_vector`].
impl From<&DirectionEnum> for Vector4<Real> {
fn from(d: &DirectionEnum) -> Self {
d.to_unit_vector()
}
}
#[cfg(test)]
mod test {
use super::*;
#[test]
fn directions_list_cst() {
let dirs = Direction::<3>::positive_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, true).unwrap());
}
let dirs = Direction::<3>::negative_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, false).unwrap());
}
let dirs = Direction::<8>::positive_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, true).unwrap());
}
let dirs = Direction::<8>::negative_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, false).unwrap());
}
let dirs = Direction::<0>::positive_directions();
assert_eq!(dirs.len(), 0);
let dirs = Direction::<0>::negative_directions();
assert_eq!(dirs.len(), 0);
let dirs = Direction::<128>::positive_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, true).unwrap());
}
let dirs = Direction::<128>::negative_directions();
for (index, dir) in dirs.iter().enumerate() {
assert_eq!(*dir, Direction::new(index, false).unwrap());
}
let dirs = Direction::<128>::positive_directions();
let dir_vec = Direction::<128>::positives_vec();
assert_eq!(dirs.to_vec(), dir_vec);
let dir_vec = Direction::<128>::directions_vec();
assert_eq!(dir_vec.len(), 256);
}
#[test]
#[allow(clippy::explicit_counter_loop)]
fn lattice_pt() {
let array = [1, 2, 3, 4];
let mut pt = LatticePoint::from(array);
assert_eq!(LatticePoint::from(array), LatticePoint::new(array.into()));
assert_eq!(<[usize; 4]>::from(pt), array);
assert!(!pt.is_empty());
assert_eq!(pt.iter().count(), 4);
assert_eq!(pt.iter_mut().count(), 4);
assert_eq!(pt.iter().copied().collect::<Vec<_>>(), vec![1, 2, 3, 4]);
assert_eq!(pt.iter_mut().collect::<Vec<_>>(), vec![&1, &2, &3, &4]);
let mut j = 0;
for i in &pt {
assert_eq!(*i, array[j]);
j += 1;
}
let mut j = 0;
for i in &mut pt {
assert_eq!(*i, array[j]);
j += 1;
}
pt.iter_mut().for_each(|el| *el = 0);
assert_eq!(pt, LatticePoint::new_zero());
assert_eq!(LatticePoint::<4>::default(), LatticePoint::new_zero());
let mut pt_2 = LatticePoint::<0>::new_zero();
assert!(pt_2.is_empty());
assert_eq!(pt_2.iter().count(), 0);
assert_eq!(pt_2.iter_mut().count(), 0);
assert_eq!(pt.to_string(), pt.as_ref().to_string());
}
#[test]
#[should_panic]
fn set_dir_neg() {
let mut lattice_link_canonical =
LatticeLinkCanonical::new(LatticePoint::new([0; 4].into()), DirectionEnum::XPos.into())
.unwrap();
lattice_link_canonical.set_dir(DirectionEnum::XNeg.into());
}
#[test]
#[allow(clippy::cognitive_complexity)]
fn dir() {
assert!(Direction::<4>::new(4, true).is_none());
assert!(Direction::<4>::new(32, true).is_none());
assert!(Direction::<4>::new(3, true).is_some());
assert!(Direction::<127>::new(128, true).is_none());
assert_eq!(Direction::<4>::dim(), 4);
assert_eq!(Direction::<124>::dim(), 124);
assert!(DirectionEnum::TNeg.is_negative());
assert!(!DirectionEnum::ZPos.is_negative());
assert_eq!(-&DirectionEnum::TNeg, &DirectionEnum::TPos);
assert_eq!(-&DirectionEnum::TPos, &DirectionEnum::TNeg);
assert_eq!(-&DirectionEnum::ZNeg, &DirectionEnum::ZPos);
assert_eq!(-&DirectionEnum::ZPos, &DirectionEnum::ZNeg);
assert_eq!(-&DirectionEnum::YNeg, &DirectionEnum::YPos);
assert_eq!(-&DirectionEnum::YPos, &DirectionEnum::YNeg);
assert_eq!(-&DirectionEnum::XNeg, &DirectionEnum::XPos);
assert_eq!(-&DirectionEnum::XPos, &DirectionEnum::XNeg);
assert_eq!(-DirectionEnum::TNeg, DirectionEnum::TPos);
assert_eq!(-DirectionEnum::TPos, DirectionEnum::TNeg);
assert_eq!(-DirectionEnum::ZNeg, DirectionEnum::ZPos);
assert_eq!(-DirectionEnum::ZPos, DirectionEnum::ZNeg);
assert_eq!(-DirectionEnum::YNeg, DirectionEnum::YPos);
assert_eq!(-DirectionEnum::YPos, DirectionEnum::YNeg);
assert_eq!(-DirectionEnum::XNeg, DirectionEnum::XPos);
assert_eq!(-DirectionEnum::XPos, DirectionEnum::XNeg);
assert_eq!(DirectionEnum::directions().len(), 8);
assert_eq!(DirectionEnum::positive_directions().len(), 4);
assert_eq!(Direction::<4>::directions().len(), 8);
assert_eq!(Direction::<4>::positive_directions().len(), 4);
let l = LatticeCyclic::new(1_f64, 4).unwrap();
for dir in Direction::<3>::directions() {
assert_eq!(
<Direction<3> as LatticeElementToIndex<3>>::to_index(dir, &l),
dir.index()
);
}
let array_dir_name = ["X", "Y", "Z", "T"];
let array_pos = ["positive", "negative"];
for (i, dir) in DirectionEnum::directions().iter().enumerate() {
assert_eq!(
dir.to_string(),
format!("{} {} direction", array_pos[i / 4], array_dir_name[i % 4])
);
}
}
#[test]
fn lattice_link() {
let pt = [0, 0, 0, 0].into();
let dir = DirectionEnum::XNeg.into();
let mut link = LatticeLinkCanonical::new(pt, DirectionEnum::YPos.into()).unwrap();
link.set_dir_positive(dir);
assert_eq!(
link,
LatticeLinkCanonical::new(pt, dir.to_positive()).unwrap()
);
let mut link = LatticeLink::new(pt, DirectionEnum::YPos.into());
let pos = *link.pos();
let dir = *link.dir();
assert_eq!(pos, *link.pos_mut());
assert_eq!(dir, *link.dir_mut());
assert!(link.is_dir_positive());
*link.dir_mut() = DirectionEnum::YNeg.into();
assert!(!link.is_dir_positive());
assert_eq!(
link.to_string(),
format!(
"link [position {}, direction [index 1, is positive false]]",
SVector::<usize, 4>::zeros()
)
);
let vector = SVector::<usize, 5>::from([1, 0, 0, 0, 0]);
let canonical_link = LatticeLinkCanonical::<5>::new(
LatticePoint::new(vector),
Direction::new(0, true).unwrap(),
)
.unwrap();
println!("{}", canonical_link);
assert_eq!(
canonical_link.to_string(),
format!(
"canonical link [position {}, direction [index 0, is positive true]]",
vector
)
);
}
#[allow(clippy::cognitive_complexity)]
#[test]
fn iterator() {
let l = LatticeCyclic::<2>::new(1_f64, 4).unwrap();
let mut iterator = l.get_points();
assert_eq!(
iterator.size_hint(),
(l.number_of_points(), Some(l.number_of_points()))
);
assert_eq!(iterator.size_hint(), (16, Some(16)));
iterator.nth(9);
assert_eq!(
iterator.size_hint(),
(
l.number_of_points() - 10, // 6
Some(l.number_of_points() - 10) // 6
)
);
assert!(iterator.nth(4).is_some());
assert_eq!(iterator.size_hint(), (1, Some(1)));
assert!(iterator.next().is_some());
assert_eq!(iterator.size_hint(), (0, Some(0)));
assert!(iterator.next().is_none());
assert_eq!(iterator.size_hint(), (0, Some(0)));
let mut iterator = l.get_links();
assert_eq!(
iterator.size_hint(),
(
l.number_of_canonical_links_space(),
Some(l.number_of_canonical_links_space())
)
);
assert_eq!(iterator.size_hint(), (16 * 2, Some(16 * 2)));
iterator.nth(9);
assert_eq!(
iterator.size_hint(),
(
l.number_of_canonical_links_space() - 10, // 6
Some(l.number_of_canonical_links_space() - 10) // 6
)
);
assert!(iterator.nth(20).is_some());
assert_eq!(iterator.size_hint(), (1, Some(1)));
assert!(iterator.next().is_some());
assert_eq!(iterator.size_hint(), (0, Some(0)));
assert!(iterator.next().is_none());
assert_eq!(iterator.size_hint(), (0, Some(0)));
let mut iterator = IteratorDirection::<2, true>::new(None).unwrap();
assert_eq!(iterator.size_hint(), (2, Some(2)));
assert!(iterator.next().is_some());
assert_eq!(iterator.size_hint(), (1, Some(1)));
assert!(iterator.next().is_some());
assert_eq!(iterator.size_hint(), (0, Some(0)));
assert!(iterator.next().is_none());
assert!(iterator.next().is_none());
//----
assert_eq!(
IteratorElement::<i32>::FirstElement.to_string(),
"first element"
);
assert_eq!(IteratorElement::Element(0_i32).to_string(), "element 0");
assert_eq!(
IteratorElement::<i32>::LastElement.to_string(),
"last element"
);
assert_eq!(
IteratorElement::<i32>::default(),
IteratorElement::<i32>::FirstElement,
);
}
#[test]
fn lattice() {
let lattice = LatticeCyclic::<3>::new(1_f64, 8).unwrap();
assert_eq!(
lattice.to_string(),
"Cyclic lattice with 8^3 points and spacing 1"
);
}
}