Crate Development - System Module

Structure of the system module

I am in the process of defining the first System module. The initial structure I tried to create was using Generic types that reference each other, as shown below:

trait Point<N : Neighbor<P>>{}
trait Neighbor<P : Point<N>>{}

However, this approach was not compatible with Rust syntax, so I tried to code using method steps with generic types, as shown below:

trait Point{
  fn neighbor<N : Neighbor>(&self) -> N;
}

trait Neighbor{
  fn inclusion<P : Point>(&self, pos : &P) -> bool;
}

However, even with this approach, it was impossible to achieve the structure I wanted. My desired structure was to create pairs between specific structures P with Point trait and N with Neighbor trait. With the above implementation, P needs to write functions for generic neighbor N, and N needs to write functions for generic point P, making it impossible to define relationships between two specific structures. In other words, when there is another structure P2 with Point trait, Neighbor’s inclusion must be written to encompass not only P but also P2. This means we cannot use any of the unique functionalities that the structures inherently possess.

The solution was surprisingly simple. We just need to follow the actual structure of Topology. Topology is originally defined as a pair of a set of points and a basis set. Point and Neighbor must form a pair to be defined as a Topology. Point and Neighbor, which don’t yet require any specific functionality, can be defined as empty Traits, and Topology can be defined as follows:

/// Point in the system
///
/// A trait representing points in a topology.
/// Each point must be able to identify its corresponding neighbors.
pub trait Point{
}

/// Neighbor of a point
///
/// A trait representing neighbors of a point in a topology.
/// Each neighbor must be able to check whether a point is in its neighborhood or not.
pub trait Neighbor<'a>{
}


/// Topology on the system
///
/// A trait representing the topology of a system,
/// which is defined using Point and Neighbor traits.
/// The relationship between points and neighbors is defined only within the topology.
pub trait Topology<'a, P, N>
    where P : Point,
          N : Neighbor<'a>{
    /// Return a neighbor of the point
    ///
    /// A function that returns the neighbor of a point.
    /// In the simulation, the practical meaning of a neighbor is a set of points that can be reached in one step.
    fn neighbor<'b : 'a>(&self, pos : &'b P) -> N;

    /// Check whether a point is in the neighbor or not
    ///
    /// A function that checks whether a point is in the neighbor or not.
    /// It is used to test whether a particle can move to a specific position or not.
    fn inclusion(&self, pos : &P, neigh : &N) -> bool;
}

Here, the lifetime issue becomes important. In the actual space, a neighbor is usually defined as an open ball centered on a specific point, which is not an efficient way in terms of memory and speed if the coordinates of the point are cloned directly. So we chose to store the reference of coordinates in the neighbor, and because Rust does not allow implicit lifetimes, a lifetime identifier is also attached to the Neighbor.

Continuous Topology

The first thing I organized was open-ball topology in continuous space. This is a topology in which each point is defined, and an open ball centered on that point with a specific distance is defined as a neighbor. Here, specific distance represents the maximum step size that a particle can move in the system. Since we will approximate the motion of a particle moving continuously as a superposition of infinitesimal changes, we need to check whether a step is too large to reduce the error before it occurs within the code range.

The code structure written from this perspective is as follows.

/// Cartessian coordinate of n-dimensional continuous system
#[derive(Clone, Debug, PartialEq, PartialOrd, Hash)]
pub struct CartessianND<T>{
    /// Coordinate of Point
    pub coord : Vec<T>,
}

impl Point for CartessianND<f64>{}

impl CartessianND<f64>{
    /// return a dimension of system
    pub fn dim(&self) -> usize{
        self.coord.len()
    }

    /// return a distance to other point
    pub fn distance(&self, other : &Self) -> f64{
        let dim = self.dim();
        if dim != other.dim(){
            panic!("{}", ErrorCode::InvalidDimension);
        }

        let mut r = 0f64;
        for i in 0..dim{
            let dx = self.coord[i] - other.coord[i];
            r += dx * dx;
        }

        return r.sqrt();
    }
}

/// Open-Ball neighborhood corresponds to [`Cartessian1D`](struct@Cartessian1D)
#[derive(Copy, Clone, Debug, PartialEq, PartialOrd)]
pub struct OpenBallND<'a>{
    /// Point of ball center
    pub center : &'a CartessianND<f64>,
}

impl Neighbor<'_> for OpenBallND<'_>{}

#[derive(Copy, Clone, Debug, PartialEq, PartialOrd)]
pub struct ContinuousTopology{
    /// Maximum step size available in the system
    max_step : f64,
}

impl<'a> Topology<'a, CartessianND<f64>, OpenBallND<'a>> for ContinuousTopology{
    /// Return a openball of center pos.
    fn neighbor<'b : 'a>(&self, pos : &'b CartessianND<f64>) -> OpenBallND<'b>{
        OpenBallND{
            center : &pos,
        }
    }

    /// Check whether the point is in a neighbor or not.
    fn inclusion(&self, pos : &CartessianND<f64>, neigh : &OpenBallND) -> bool{
        let r = pos.distance(neigh.center);
        r < self.max_step
    }
}

The next change is the dimension part. If the position is represented as a vector as above, every time two vectors are operated, the dimension must be checked, which is inconvenient. Since an if statement is added, the speed is also slower. Since we often use systems up to 4D, we thought it would be better to define coordinates as fixed-size arrays rather than vectors. If they are defined like this, vectors of different dimensions will be filtered out at the compile stage, so there is no need to use if statements. However, writing repetitive code multiple times is also a hassle, so we defined multiple coordinate systems through the following macro.

macro_rules! impl_cartessian_nD{
    ($type_name:ident, $dim:expr) =>{
        #[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Eq, Ord, Hash)]
        pub struct $type_name<T>{
            /// Coordinate of Point
            pub coord : [T; $dim],
        }

        impl Point for $type_name<f64>{}

        impl $type_name<f64>{
            pub fn dim(&self) -> usize{
                $dim
            }

            pub fn distance(&self, other : &Self) -> f64{
                let mut r = 0f64;
                for i in 0..$dim{
                    let dx = self.coord[i] - other.coord[i];
                    r += dx * dx;
                }

                return r.sqrt();
            }
        }
    }
}

impl_cartessian_nD!(Cartessian2D, 2);

macro_rules! impl_neighbor_cartessian_nD{
    ($neighbor_name:ident, $cartessian_name:ident, $dim : expr) =>{
        #[derive(Copy, Clone, Debug, PartialEq, PartialOrd)]
        pub struct $neighbor_name<'a>{
            /// Point of ball center
            pub center : &'a $cartessian_name<f64>,
        }

        impl Neighbor<'_> for $neighbor_name<'_>{}

        impl<'a> $neighbor_name<'a>{
            pub fn new(center : &'a $cartessian_name<f64>) -> Self{
                $neighbor_name{
                    center,
                }
            }
        }
    }
}

impl_neighbor_cartessian_nD!(OpenBall2D, Cartessian2D, 2);

macro_rules! impl_continuous_topology_nD{
    ($cartessian_name:ident, $neighbor_name:ident, $dim:expr) =>{
        impl<'a> Topology<'a, $cartessian_name<f64>, $neighbor_name<'a>> for ContinuousTopology{
            fn neighbor<'b : 'a>(&self, pos : &'b $cartessian_name<f64>) -> $neighbor_name<'b>{
                $neighbor_name{
                    center : &pos,
                }
            }

            fn inclusion(&self, pos : &$cartessian_name<f64>, neigh : &$neighbor_name) -> bool{
                let r = pos.distance(neigh.center);
                r < self.max_step
            }
        }
    }
}

impl_continuous_topology_nD!(Cartessian2D, OpenBall2D, 2);

During this process, it was necessary to add comments for doc test, but the macro does not parse the value inside the comment. Therefore, we added a macro provided by the doc_comment crate to create documents.

Remove neighbor trait

However, I thought it was meaningless to define neighbor like this. It seemed like a variable that just stored a reference was not particularly meaningful. The most important role of neighbor is to determine whether a specific movement is possible or not. However, in continuous systems or lattices, it is enough to test the distance between two points, and in networks, we need to define an adjoint matrix in the system, so we can test it based on that. Then, there is no need for neighbor in this case. So I removed the neighbor trait and rewrote the code.

/// Point in the system
pub trait Point{}

/// Cartessian coordinate of 1-dimensional continuous system
#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Eq, Ord, Hash)]
pub struct Cartessian1D<T>{
    /// Coordinate of Point
    pub coord : T,
}

impl Point for Cartessian1D<f64>{}

impl Cartessian1D<f64>{
    /// return a dimension of system
    pub fn dim(&self) -> usize{
        1
    }

    /// return a distance to other point
    pub fn distance(&self, other : &Self) -> f64{
        (self.coord - other.coord).abs()
    }

    /// return a norm of vector
    pub fn norm(&self) -> f64{
        return self.coord.abs();
    }
}

/// Topology on the system
pub trait Topology<P>
    where P : Point{
    /// Check whether a movement is valid or not.
    fn check_move(&self, pos : &P, movement : &P) -> bool;
}

/// Normal openball topology for continuous system
#[derive(Copy, Clone, Debug, PartialEq, PartialOrd)]
pub struct ContinuousTopology{
    /// Maximum step size available in the system
    max_step : f64,
}

impl Topology<Cartessian1D<f64>> for ContinuousTopology{
    /// Check whether the movement is valid or not
    fn check_move(&self, _pos : &Cartessian1D<f64>, movement : &Cartessian1D<f64>) -> bool{
        let r = movement.norm();
        r < self.max_step
    }
}

Lattice Topology

Lattice is almost fully developed. I developed Cartessian<i32> together with Cartessian<f64>. Here, I defined taxi_distance and taxi_norm, and based on that, I was able to define the topology.

impl Topology<CartessianND<i32>> for LatticeTopology{
    /// Check whether the movement is valid or not
    fn check_move(&self, pos : &CartessianND<i32>, movement : &CartessianND<i32>) -> bool{
        if pos.dim() != movement.dim(){
            panic!("{}", ErrorCode::InvalidDimension);
        }
        let r = movement.taxi_norm();
        r <= self.max_step
    }
}

Network Topology

Network topology needs to be developed from point. Here, Point only needs to be given an index, and the connection information is placed in the topology structure. At this time, the network may have a weight on the edge, so I defined the network as a generic type with multiple variables, and used default for 0 or false, which is disconnected.

/// node of a network
#[derive(Copy, Clone, Debug, PartialEq, PartialOrd, Eq, Ord, Hash)]
pub struct NodeIndex{
    /// Index of node
    pub index : usize,
}

impl Point for NodeIndex{}

/// Network topology
#[derive(Clone, Debug, PartialEq, PartialOrd)]
pub struct NetworkTopology<T>{
    /// The number of nodes
    num_node : usize,
    /// Adjoint matrix,
    adjoint : Vec<Vec<T>>,
}

impl<T : Copy + Default + PartialEq> Topology<NodeIndex> for NetworkTopology<T>{
    /// Check whether the movement is valid or not
    fn check_move(&self, depart : &NodeIndex, arrival : &NodeIndex) -> bool{
        let connect = self.adjoint[depart.index][arrival.index];
        !(connect == T::default())
    }
}