use ndarray::{self, ArrayBase};
use std::cmp;
use std::collections::HashSet;
use std::convert::AsRef;
use std::hash::Hash;
use std::iter::{Enumerate, FilterMap, IntoIterator, Peekable, Sum, Zip};
use std::marker::PhantomData;
use std::ops::{Add, Deref, DerefMut, Index, IndexMut, Mul, Neg, Sub};
use std::slice::{self, Iter, IterMut};
use Ix1;
use num_traits::{Num, Zero};
use array_backend::Array2;
use errors::SprsError;
use indexing::SpIndex;
use sparse::csmat::CompressedStorage::{CSC, CSR};
use sparse::permutation::PermViewI;
use sparse::prelude::*;
use sparse::utils;
use sparse::{binop, prod};
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub struct NnzIndex(pub usize);
pub trait VecDim<N> {
fn dim(&self) -> usize;
}
impl<N, IS, DS: Deref<Target = [N]>> VecDim<N> for CsVecBase<IS, DS> {
fn dim(&self) -> usize {
self.dim
}
}
impl<N, T: ?Sized> VecDim<N> for T
where
T: AsRef<[N]>,
{
fn dim(&self) -> usize {
self.as_ref().len()
}
}
pub struct VectorIterator<'a, N: 'a, I: 'a> {
ind_data: Zip<Iter<'a, I>, Iter<'a, N>>,
}
pub struct VectorIteratorPerm<'a, N: 'a, I: 'a> {
ind_data: Zip<Iter<'a, I>, Iter<'a, N>>,
perm: PermViewI<'a, I>,
}
pub struct VectorIteratorMut<'a, N: 'a, I: 'a> {
ind_data: Zip<Iter<'a, I>, IterMut<'a, N>>,
}
impl<'a, N: 'a, I: 'a + SpIndex> Iterator for VectorIterator<'a, N, I> {
type Item = (usize, &'a N);
fn next(&mut self) -> Option<<Self as Iterator>::Item> {
match self.ind_data.next() {
None => None,
Some((inner_ind, data)) => Some((inner_ind.index(), data)),
}
}
fn size_hint(&self) -> (usize, Option<usize>) {
self.ind_data.size_hint()
}
}
impl<'a, N: 'a, I: 'a + SpIndex> Iterator for VectorIteratorPerm<'a, N, I> {
type Item = (usize, &'a N);
fn next(&mut self) -> Option<<Self as Iterator>::Item> {
match self.ind_data.next() {
None => None,
Some((inner_ind, data)) => {
Some((self.perm.at(inner_ind.index()), data))
}
}
}
fn size_hint(&self) -> (usize, Option<usize>) {
self.ind_data.size_hint()
}
}
impl<'a, N: 'a, I: 'a + SpIndex> Iterator for VectorIteratorMut<'a, N, I> {
type Item = (usize, &'a mut N);
fn next(&mut self) -> Option<<Self as Iterator>::Item> {
match self.ind_data.next() {
None => None,
Some((inner_ind, data)) => Some((inner_ind.index(), data)),
}
}
fn size_hint(&self) -> (usize, Option<usize>) {
self.ind_data.size_hint()
}
}
pub trait SparseIterTools: Iterator {
fn nnz_or_zip<'a, I, N1, N2>(
self,
other: I,
) -> NnzOrZip<'a, Self, I::IntoIter, N1, N2>
where
Self: Iterator<Item = (usize, &'a N1)> + Sized,
I: IntoIterator<Item = (usize, &'a N2)>,
{
NnzOrZip {
left: self.peekable(),
right: other.into_iter().peekable(),
life: PhantomData,
}
}
fn nnz_zip<'a, I, N1, N2>(
self,
other: I,
) -> FilterMap<
NnzOrZip<'a, Self, I::IntoIter, N1, N2>,
fn(NnzEither<'a, N1, N2>) -> Option<(usize, &'a N1, &'a N2)>,
>
where
Self: Iterator<Item = (usize, &'a N1)> + Sized,
I: IntoIterator<Item = (usize, &'a N2)>,
{
let nnz_or_iter = NnzOrZip {
left: self.peekable(),
right: other.into_iter().peekable(),
life: PhantomData,
};
nnz_or_iter.filter_map(filter_both_nnz)
}
}
impl<T: Iterator> SparseIterTools for Enumerate<T> {}
impl<'a, N: 'a, I: 'a + SpIndex> SparseIterTools for VectorIterator<'a, N, I> {}
pub trait IntoSparseVecIter<'a, N: 'a> {
type IterType;
fn into_sparse_vec_iter(
self,
) -> <Self as IntoSparseVecIter<'a, N>>::IterType
where
<Self as IntoSparseVecIter<'a, N>>::IterType:
Iterator<Item = (usize, &'a N)>;
fn dim(&self) -> usize;
fn is_dense(&self) -> bool {
false
}
#[allow(unused_variables)]
fn index(self, idx: usize) -> &'a N
where
Self: Sized,
{
panic!("cannot be called on a vector that is not dense");
}
}
impl<'a, N: 'a, I: 'a> IntoSparseVecIter<'a, N> for CsVecViewI<'a, N, I>
where
I: SpIndex,
{
type IterType = VectorIterator<'a, N, I>;
fn dim(&self) -> usize {
self.dim()
}
fn into_sparse_vec_iter(self) -> VectorIterator<'a, N, I> {
self.iter_rbr()
}
}
impl<'a, N: 'a, I: 'a, IS, DS> IntoSparseVecIter<'a, N>
for &'a CsVecBase<IS, DS>
where
I: SpIndex,
IS: Deref<Target = [I]>,
DS: Deref<Target = [N]>,
{
type IterType = VectorIterator<'a, N, I>;
fn dim(&self) -> usize {
(*self).dim()
}
fn into_sparse_vec_iter(self) -> VectorIterator<'a, N, I> {
self.iter()
}
}
impl<'a, N: 'a> IntoSparseVecIter<'a, N> for &'a [N] {
type IterType = Enumerate<Iter<'a, N>>;
fn dim(&self) -> usize {
self.len()
}
fn into_sparse_vec_iter(self) -> Enumerate<Iter<'a, N>> {
self.into_iter().enumerate()
}
fn is_dense(&self) -> bool {
true
}
fn index(self, idx: usize) -> &'a N {
&self[idx]
}
}
impl<'a, N: 'a> IntoSparseVecIter<'a, N> for &'a Vec<N> {
type IterType = Enumerate<Iter<'a, N>>;
fn dim(&self) -> usize {
self.len()
}
fn into_sparse_vec_iter(self) -> Enumerate<Iter<'a, N>> {
self.into_iter().enumerate()
}
fn is_dense(&self) -> bool {
true
}
fn index(self, idx: usize) -> &'a N {
&self[idx]
}
}
impl<'a, N: 'a, S> IntoSparseVecIter<'a, N> for &'a ArrayBase<S, Ix1>
where
S: ndarray::Data<Elem = N>,
{
type IterType = Enumerate<ndarray::iter::Iter<'a, N, Ix1>>;
fn dim(&self) -> usize {
self.shape()[0]
}
fn into_sparse_vec_iter(
self,
) -> Enumerate<ndarray::iter::Iter<'a, N, Ix1>> {
self.iter().enumerate()
}
fn is_dense(&self) -> bool {
true
}
fn index(self, idx: usize) -> &'a N {
&self[[idx]]
}
}
pub trait DenseVector<N> {
fn dim(&self) -> usize;
fn index(&self, idx: usize) -> &N;
}
impl<'a, N: 'a> DenseVector<N> for &'a [N] {
fn dim(&self) -> usize {
self.len()
}
fn index(&self, idx: usize) -> &N {
&self[idx]
}
}
impl<N> DenseVector<N> for Vec<N> {
fn dim(&self) -> usize {
self.len()
}
fn index(&self, idx: usize) -> &N {
&self[idx]
}
}
impl<'a, N: 'a> DenseVector<N> for &'a Vec<N> {
fn dim(&self) -> usize {
self.len()
}
fn index(&self, idx: usize) -> &N {
&self[idx]
}
}
impl<N, S> DenseVector<N> for ArrayBase<S, Ix1>
where
S: ndarray::Data<Elem = N>,
{
fn dim(&self) -> usize {
self.shape()[0]
}
fn index(&self, idx: usize) -> &N {
&self[[idx]]
}
}
pub struct NnzOrZip<'a, Ite1, Ite2, N1: 'a, N2: 'a>
where
Ite1: Iterator<Item = (usize, &'a N1)>,
Ite2: Iterator<Item = (usize, &'a N2)>,
{
left: Peekable<Ite1>,
right: Peekable<Ite2>,
life: PhantomData<(&'a N1, &'a N2)>,
}
#[derive(PartialEq, Debug)]
pub enum NnzEither<'a, N1: 'a, N2: 'a> {
Both((usize, &'a N1, &'a N2)),
Left((usize, &'a N1)),
Right((usize, &'a N2)),
}
fn filter_both_nnz<'a, N: 'a, M: 'a>(
elem: NnzEither<'a, N, M>,
) -> Option<(usize, &'a N, &'a M)> {
match elem {
NnzEither::Both((ind, lval, rval)) => Some((ind, lval, rval)),
_ => None,
}
}
impl<'a, Ite1, Ite2, N1: 'a, N2: 'a> Iterator
for NnzOrZip<'a, Ite1, Ite2, N1, N2>
where
Ite1: Iterator<Item = (usize, &'a N1)>,
Ite2: Iterator<Item = (usize, &'a N2)>,
{
type Item = NnzEither<'a, N1, N2>;
fn next(&mut self) -> Option<(NnzEither<'a, N1, N2>)> {
match (self.left.peek(), self.right.peek()) {
(None, Some(&(_, _))) => {
let (rind, rval) = self.right.next().unwrap();
Some(NnzEither::Right((rind, rval)))
}
(Some(&(_, _)), None) => {
let (lind, lval) = self.left.next().unwrap();
Some(NnzEither::Left((lind, lval)))
}
(None, None) => None,
(Some(&(lind, _)), Some(&(rind, _))) => {
if lind < rind {
let (lind, lval) = self.left.next().unwrap();
Some(NnzEither::Left((lind, lval)))
} else if rind < lind {
let (rind, rval) = self.right.next().unwrap();
Some(NnzEither::Right((rind, rval)))
} else {
let (lind, lval) = self.left.next().unwrap();
let (_, rval) = self.right.next().unwrap();
Some(NnzEither::Both((lind, lval, rval)))
}
}
}
}
#[inline]
fn size_hint(&self) -> (usize, Option<usize>) {
let (left_lower, left_upper) = self.left.size_hint();
let (right_lower, right_upper) = self.right.size_hint();
let upper = match (left_upper, right_upper) {
(Some(x), Some(y)) => Some(x + y),
(Some(x), None) => Some(x),
(None, Some(y)) => Some(y),
(None, None) => None,
};
(cmp::max(left_lower, right_lower), upper)
}
}
impl<N, I: SpIndex> CsVecBase<Vec<I>, Vec<N>> {
pub fn new(n: usize, mut indices: Vec<I>, mut data: Vec<N>) -> CsVecI<N, I>
where
N: Copy,
{
let mut buf = Vec::with_capacity(indices.len());
utils::sort_indices_data_slices(
&mut indices[..],
&mut data[..],
&mut buf,
);
let v = CsVecI {
dim: n,
indices,
data,
};
v.check_structure().and(Ok(v)).unwrap()
}
pub fn empty(dim: usize) -> CsVecI<N, I> {
CsVecI {
dim,
indices: Vec::new(),
data: Vec::new(),
}
}
pub fn append(&mut self, ind: usize, val: N) {
match self.indices.last() {
None => (),
Some(&last_ind) => {
assert!(ind > last_ind.index(), "unsorted append")
}
}
assert!(ind <= self.dim, "out of bounds index");
self.indices.push(I::from_usize(ind));
self.data.push(val);
}
pub fn reserve(&mut self, size: usize) {
self.indices.reserve(size);
self.data.reserve(size);
}
pub fn reserve_exact(&mut self, exact_size: usize) {
self.indices.reserve_exact(exact_size);
self.data.reserve_exact(exact_size);
}
pub fn clear(&mut self) {
self.indices.clear();
self.data.clear();
}
}
impl<N, I, IStorage, DStorage> CsVecBase<IStorage, DStorage>
where
I: SpIndex,
IStorage: Deref<Target = [I]>,
DStorage: Deref<Target = [N]>,
{
pub fn view(&self) -> CsVecViewI<N, I> {
CsVecViewI {
dim: self.dim,
indices: &self.indices[..],
data: &self.data[..],
}
}
pub fn iter(&self) -> VectorIterator<N, I> {
VectorIterator {
ind_data: self.indices.iter().zip(self.data.iter()),
}
}
#[doc(hidden)]
pub fn iter_perm<'a, 'perm: 'a>(
&'a self,
perm: PermViewI<'perm, I>,
) -> VectorIteratorPerm<'a, N, I>
where
N: 'a,
{
VectorIteratorPerm {
ind_data: self.indices.iter().zip(self.data.iter()),
perm,
}
}
pub fn indices(&self) -> &[I] {
&self.indices
}
pub fn data(&self) -> &[N] {
&self.data
}
pub fn into_raw_storage(self) -> (IStorage, DStorage) {
let Self { indices, data, .. } = self;
(indices, data)
}
pub fn dim(&self) -> usize {
self.dim
}
pub fn nnz(&self) -> usize {
self.data.len()
}
pub fn check_structure(&self) -> Result<(), SprsError> {
if !self.indices.windows(2).all(|x| x[0] < x[1]) {
return Err(SprsError::NonSortedIndices);
}
if self.dim == 0 && self.indices.len() == 0 && self.data.len() == 0 {
return Ok(());
}
let max_ind = self.indices.iter().max().unwrap_or(&I::zero()).index();
if max_ind >= self.dim {
panic!("Out of bounds index");
}
Ok(())
}
pub fn to_owned(&self) -> CsVecI<N, I>
where
N: Clone,
{
CsVecI {
dim: self.dim,
indices: self.indices.to_vec(),
data: self.data.to_vec(),
}
}
pub fn to_other_types<I2>(&self) -> CsVecI<N, I2>
where
N: Clone,
I2: SpIndex,
{
let indices = self
.indices
.iter()
.map(|i| I2::from_usize(i.index()))
.collect();
let data = self.data.iter().cloned().collect();
CsVecI {
dim: self.dim,
indices,
data,
}
}
pub fn row_view(&self) -> CsMatVecView_<N, I> {
let indptr = Array2 {
data: [I::zero(), I::from_usize(self.indices.len())],
};
CsMatBase {
storage: CSR,
nrows: 1,
ncols: self.dim,
indptr,
indices: &self.indices[..],
data: &self.data[..],
}
}
pub fn col_view(&self) -> CsMatVecView_<N, I> {
let indptr = Array2 {
data: [I::zero(), I::from_usize(self.indices.len())],
};
CsMatBase {
storage: CSC,
nrows: self.dim,
ncols: 1,
indptr,
indices: &self.indices[..],
data: &self.data[..],
}
}
pub fn get<'a>(&'a self, index: usize) -> Option<&'a N>
where
I: 'a,
{
self.view().get_rbr(index)
}
pub fn nnz_index(&self, index: usize) -> Option<NnzIndex> {
self.indices
.binary_search(&I::from_usize(index))
.map(|i| NnzIndex(i.index()))
.ok()
}
pub fn dot<'b, T: IntoSparseVecIter<'b, N>>(&'b self, rhs: T) -> N
where
N: 'b + Num + Copy + Sum,
I: 'b,
<T as IntoSparseVecIter<'b, N>>::IterType:
Iterator<Item = (usize, &'b N)>,
T: Copy,
{
assert_eq!(self.dim(), rhs.dim());
if rhs.is_dense() {
self.iter()
.map(|(idx, val)| *val * *rhs.index(idx.index()))
.sum()
} else {
let mut lhs_iter = self.iter();
let mut rhs_iter = rhs.into_sparse_vec_iter().into_iter();
let mut sum = N::zero();
let mut left_nnz = lhs_iter.next();
let mut right_nnz = rhs_iter.next();
while left_nnz.is_some() && right_nnz.is_some() {
let (left_ind, left_val) = left_nnz.unwrap();
let (right_ind, right_val) = right_nnz.unwrap();
if left_ind == right_ind {
sum = sum + *left_val * *right_val;
}
if left_ind <= right_ind {
left_nnz = lhs_iter.next();
}
if left_ind >= right_ind {
right_nnz = rhs_iter.next();
}
}
sum
}
}
pub fn dot_dense<T>(&self, rhs: T) -> N
where
T: DenseVector<N>,
N: Num + Copy + Sum,
{
assert_eq!(self.dim(), rhs.dim());
self.iter()
.map(|(idx, val)| *val * *rhs.index(idx.index()))
.sum()
}
pub fn scatter(&self, out: &mut [N])
where
N: Clone,
{
for (ind, val) in self.iter() {
out[ind] = val.clone();
}
}
pub fn to_set(&self) -> HashSet<(usize, N)>
where
N: Hash + Eq + Clone,
{
self.indices()
.iter()
.map(|i| i.index())
.zip(self.data.iter().cloned())
.collect()
}
pub fn map<F>(&self, f: F) -> CsVecI<N, I>
where
F: FnMut(&N) -> N,
N: Clone,
{
let mut res = self.to_owned();
res.map_inplace(f);
res
}
}
impl<'a, N, I, IStorage, DStorage> CsVecBase<IStorage, DStorage>
where
N: 'a,
I: 'a + SpIndex,
IStorage: 'a + Deref<Target = [I]>,
DStorage: DerefMut<Target = [N]>,
{
fn data_mut(&mut self) -> &mut [N] {
&mut self.data[..]
}
pub fn view_mut(&mut self) -> CsVecViewMutI<N, I> {
CsVecBase {
dim: self.dim,
indices: &self.indices[..],
data: &mut self.data[..],
}
}
pub fn get_mut(&mut self, index: usize) -> Option<&mut N> {
if let Some(NnzIndex(position)) = self.nnz_index(index) {
Some(&mut self.data[position])
} else {
None
}
}
pub fn map_inplace<F>(&mut self, mut f: F)
where
F: FnMut(&N) -> N,
{
for val in &mut self.data[..] {
*val = f(val);
}
}
pub fn iter_mut(&mut self) -> VectorIteratorMut<N, I> {
VectorIteratorMut {
ind_data: self.indices.iter().zip(self.data.iter_mut()),
}
}
}
impl<'a, N: 'a, I: 'a + SpIndex> CsVecBase<&'a [I], &'a [N]> {
pub fn new_view(
n: usize,
indices: &'a [I],
data: &'a [N],
) -> Result<CsVecViewI<'a, N, I>, SprsError> {
let v = CsVecViewI {
dim: n,
indices,
data,
};
v.check_structure().and(Ok(v))
}
pub fn get_rbr(&self, index: usize) -> Option<&'a N> {
self.nnz_index(index)
.map(|NnzIndex(position)| &self.data[position])
}
fn iter_rbr(&self) -> VectorIterator<'a, N, I> {
VectorIterator {
ind_data: self.indices.iter().zip(self.data.iter()),
}
}
pub unsafe fn new_view_raw(
n: usize,
nnz: usize,
indices: *const I,
data: *const N,
) -> CsVecViewI<'a, N, I> {
CsVecViewI {
dim: n,
indices: slice::from_raw_parts(indices, nnz),
data: slice::from_raw_parts(data, nnz),
}
}
}
impl<'a, N, I> CsVecBase<&'a [I], &'a mut [N]>
where
N: 'a,
I: 'a + SpIndex,
{
pub unsafe fn new_view_mut_raw(
n: usize,
nnz: usize,
indices: *const I,
data: *mut N,
) -> CsVecViewMutI<'a, N, I> {
CsVecBase {
dim: n,
indices: slice::from_raw_parts(indices, nnz),
data: slice::from_raw_parts_mut(data, nnz),
}
}
}
impl<'a, 'b, N, I, IS1, DS1, IpS2, IS2, DS2>
Mul<&'b CsMatBase<N, I, IpS2, IS2, DS2>> for &'a CsVecBase<IS1, DS1>
where
N: 'a + Copy + Num + Default,
I: 'a + SpIndex,
IS1: 'a + Deref<Target = [I]>,
DS1: 'a + Deref<Target = [N]>,
IpS2: 'b + Deref<Target = [I]>,
IS2: 'b + Deref<Target = [I]>,
DS2: 'b + Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn mul(self, rhs: &CsMatBase<N, I, IpS2, IS2, DS2>) -> CsVecI<N, I> {
(&self.row_view() * rhs).outer_view(0).unwrap().to_owned()
}
}
impl<'a, 'b, N, I, IpS1, IS1, DS1, IS2, DS2> Mul<&'b CsVecBase<IS2, DS2>>
for &'a CsMatBase<N, I, IpS1, IS1, DS1>
where
N: Copy + Num + Default + Sum,
I: SpIndex,
IpS1: Deref<Target = [I]>,
IS1: Deref<Target = [I]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [I]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn mul(self, rhs: &CsVecBase<IS2, DS2>) -> CsVecI<N, I> {
if self.is_csr() {
prod::csr_mul_csvec(self.view(), rhs.view())
} else {
(self * &rhs.col_view()).outer_view(0).unwrap().to_owned()
}
}
}
impl<N, I, IS1, DS1, IS2, DS2> Add<CsVecBase<IS2, DS2>> for CsVecBase<IS1, DS1>
where
N: Copy + Num,
I: SpIndex,
IS1: Deref<Target = [I]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [I]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn add(self, rhs: CsVecBase<IS2, DS2>) -> CsVecI<N, I> {
&self + &rhs
}
}
impl<'a, N, I, IS1, DS1, IS2, DS2> Add<&'a CsVecBase<IS2, DS2>>
for CsVecBase<IS1, DS1>
where
N: Copy + Num,
I: SpIndex,
IS1: Deref<Target = [I]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [I]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn add(self, rhs: &CsVecBase<IS2, DS2>) -> CsVecI<N, I> {
&self + rhs
}
}
impl<'a, N, I, IS1, DS1, IS2, DS2> Add<CsVecBase<IS2, DS2>>
for &'a CsVecBase<IS1, DS1>
where
N: Copy + Num,
I: SpIndex,
IS1: Deref<Target = [I]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [I]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn add(self, rhs: CsVecBase<IS2, DS2>) -> CsVecI<N, I> {
self + &rhs
}
}
impl<'a, 'b, N, I, IS1, DS1, IS2, DS2> Add<&'b CsVecBase<IS2, DS2>>
for &'a CsVecBase<IS1, DS1>
where
N: Copy + Num,
I: SpIndex,
IS1: Deref<Target = [I]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [I]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVecI<N, I>;
fn add(self, rhs: &CsVecBase<IS2, DS2>) -> CsVecI<N, I> {
binop::csvec_binop(self.view(), rhs.view(), |&x, &y| x + y).unwrap()
}
}
impl<'a, 'b, N, IS1, DS1, IS2, DS2> Sub<&'b CsVecBase<IS2, DS2>>
for &'a CsVecBase<IS1, DS1>
where
N: Copy + Num,
IS1: Deref<Target = [usize]>,
DS1: Deref<Target = [N]>,
IS2: Deref<Target = [usize]>,
DS2: Deref<Target = [N]>,
{
type Output = CsVec<N>;
fn sub(self, rhs: &CsVecBase<IS2, DS2>) -> CsVec<N> {
binop::csvec_binop(self.view(), rhs.view(), |&x, &y| x - y).unwrap()
}
}
impl<N: Num + Copy + Neg<Output = N>, I: SpIndex> Neg for CsVecI<N, I> {
type Output = CsVecI<N, I>;
fn neg(mut self) -> CsVecI<N, I> {
for value in &mut self.data {
*value = -*value;
}
self
}
}
impl<N, IS, DS> Index<usize> for CsVecBase<IS, DS>
where
IS: Deref<Target = [usize]>,
DS: Deref<Target = [N]>,
{
type Output = N;
fn index(&self, index: usize) -> &N {
self.get(index).unwrap()
}
}
impl<N, IS, DS> IndexMut<usize> for CsVecBase<IS, DS>
where
IS: Deref<Target = [usize]>,
DS: DerefMut<Target = [N]>,
{
fn index_mut(&mut self, index: usize) -> &mut N {
self.get_mut(index).unwrap()
}
}
impl<N, IS, DS> Index<NnzIndex> for CsVecBase<IS, DS>
where
IS: Deref<Target = [usize]>,
DS: Deref<Target = [N]>,
{
type Output = N;
fn index(&self, index: NnzIndex) -> &N {
let NnzIndex(i) = index;
self.data().get(i).unwrap()
}
}
impl<N, IS, DS> IndexMut<NnzIndex> for CsVecBase<IS, DS>
where
IS: Deref<Target = [usize]>,
DS: DerefMut<Target = [N]>,
{
fn index_mut(&mut self, index: NnzIndex) -> &mut N {
let NnzIndex(i) = index;
self.data_mut().get_mut(i).unwrap()
}
}
impl<N: Num + Copy, I: SpIndex> Zero for CsVecI<N, I> {
fn zero() -> CsVecI<N, I> {
CsVecI::new(0, vec![], vec![])
}
fn is_zero(&self) -> bool {
self.data.iter().all(|x| x.is_zero())
}
}
#[cfg(feature = "alga")]
mod alga_impls {
use super::*;
use alga::general::*;
impl<N: Clone + Copy + Num, I: Clone + SpIndex> AbstractMagma<Additive>
for CsVecI<N, I>
{
fn operate(&self, right: &CsVecI<N, I>) -> CsVecI<N, I> {
self + right
}
}
impl<N: Copy + Num, I: SpIndex> Identity<Additive> for CsVecI<N, I> {
fn identity() -> CsVecI<N, I> {
CsVecI::zero()
}
}
impl<N: Copy + Num, I: SpIndex> AbstractSemigroup<Additive> for CsVecI<N, I> {}
impl<N: Copy + Num, I: SpIndex> AbstractMonoid<Additive> for CsVecI<N, I> {}
impl<N, I> Inverse<Additive> for CsVecI<N, I>
where
N: Clone + Neg<Output = N> + Copy + Num,
I: SpIndex,
{
fn inverse(&self) -> CsVecI<N, I> {
CsVecBase {
data: self.data.iter().map(|x| -*x).collect(),
indices: self.indices.clone(),
dim: self.dim,
}
}
}
impl<N: Copy + Num + Neg<Output = N>, I: SpIndex>
AbstractQuasigroup<Additive> for CsVecI<N, I>
{
}
impl<N: Copy + Num + Neg<Output = N>, I: SpIndex> AbstractLoop<Additive>
for CsVecI<N, I>
{
}
impl<N: Copy + Num + Neg<Output = N>, I: SpIndex> AbstractGroup<Additive>
for CsVecI<N, I>
{
}
impl<N: Copy + Num + Neg<Output = N>, I: SpIndex>
AbstractGroupAbelian<Additive> for CsVecI<N, I>
{
}
#[cfg(test)]
mod test {
use super::*;
#[test]
fn additive_operator_is_addition() {
let a = CsVec::new(2, vec![0], vec![2.]);
let b = CsVec::new(2, vec![0], vec![3.]);
assert_eq!(AbstractMagma::<Additive>::operate(&a, &b), &a + &b);
}
#[test]
fn additive_identity_is_zero() {
assert_eq!(CsVec::<f64>::zero(), Identity::<Additive>::identity());
}
#[test]
fn additive_inverse_is_negated() {
let vector = CsVec::new(2, vec![0], vec![2.]);
assert_eq!(-vector.clone(), Inverse::<Additive>::inverse(&vector));
}
}
}
#[cfg(test)]
mod test {
use super::SparseIterTools;
use ndarray::Array;
use num_traits::Zero;
use sparse::{CsVec, CsVecI};
fn test_vec1() -> CsVec<f64> {
let n = 8;
let indices = vec![0, 1, 4, 5, 7];
let data = vec![0., 1., 4., 5., 7.];
return CsVec::new(n, indices, data);
}
fn test_vec2() -> CsVecI<f64, usize> {
let n = 8;
let indices = vec![0, 2, 4, 6, 7];
let data = vec![0.5, 2.5, 4.5, 6.5, 7.5];
return CsVecI::new(n, indices, data);
}
#[test]
fn test_nnz_zip_iter() {
let vec1 = test_vec1();
let vec2 = test_vec2();
let mut iter = vec1.iter().nnz_zip(vec2.iter());
assert_eq!(iter.next().unwrap(), (0, &0., &0.5));
assert_eq!(iter.next().unwrap(), (4, &4., &4.5));
assert_eq!(iter.next().unwrap(), (7, &7., &7.5));
assert!(iter.next().is_none());
}
#[test]
fn test_nnz_or_zip_iter() {
use super::NnzEither::*;
let vec1 = test_vec1();
let vec2 = test_vec2();
let mut iter = vec1.iter().nnz_or_zip(vec2.iter());
assert_eq!(iter.next().unwrap(), Both((0, &0., &0.5)));
assert_eq!(iter.next().unwrap(), Left((1, &1.)));
assert_eq!(iter.next().unwrap(), Right((2, &2.5)));
assert_eq!(iter.next().unwrap(), Both((4, &4., &4.5)));
assert_eq!(iter.next().unwrap(), Left((5, &5.)));
assert_eq!(iter.next().unwrap(), Right((6, &6.5)));
assert_eq!(iter.next().unwrap(), Both((7, &7., &7.5)));
}
#[test]
fn dot_product() {
let vec1 = CsVec::new(8, vec![0, 2, 4, 6], vec![1.; 4]);
let vec2 = CsVec::new(8, vec![1, 3, 5, 7], vec![2.; 4]);
let vec3 = CsVec::new(8, vec![1, 2, 5, 6], vec![3.; 4]);
assert_eq!(0., vec1.dot(&vec2));
assert_eq!(4., vec1.dot(&vec1));
assert_eq!(16., vec2.dot(&vec2));
assert_eq!(6., vec1.dot(&vec3));
assert_eq!(12., vec2.dot(vec3.view()));
let dense_vec = vec![1., 2., 3., 4., 5., 6., 7., 8.];
{
let slice = &dense_vec[..];
assert_eq!(16., vec1.dot(&dense_vec));
assert_eq!(16., vec1.dot(slice));
assert_eq!(16., vec1.dot_dense(slice));
assert_eq!(16., vec1.dot_dense(&dense_vec));
}
assert_eq!(16., vec1.dot_dense(dense_vec));
let ndarray_vec = Array::linspace(1., 8., 8);
assert_eq!(16., vec1.dot(&ndarray_vec));
assert_eq!(16., vec1.dot_dense(ndarray_vec.view()));
}
#[test]
#[should_panic]
fn dot_product_panics() {
let vec1 = CsVec::new(8, vec![0, 2, 4, 6], vec![1.; 4]);
let vec2 = CsVec::new(9, vec![1, 3, 5, 7], vec![2.; 4]);
vec1.dot(&vec2);
}
#[test]
#[should_panic]
fn dot_product_panics2() {
let vec1 = CsVec::new(8, vec![0, 2, 4, 6], vec![1.; 4]);
let dense_vec = vec![0., 1., 2., 3., 4., 5., 6., 7., 8.];
vec1.dot(&dense_vec);
}
#[test]
fn nnz_index() {
let vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1.; 4]);
assert_eq!(vec.nnz_index(1), None);
assert_eq!(vec.nnz_index(9), None);
assert_eq!(vec.nnz_index(0), Some(super::NnzIndex(0)));
assert_eq!(vec.nnz_index(4), Some(super::NnzIndex(2)));
let index = vec.nnz_index(2).unwrap();
assert_eq!(vec[index], 1.);
let mut vec = vec;
vec[index] = 2.;
assert_eq!(vec[index], 2.);
}
#[test]
fn get_mut() {
let mut vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1.; 4]);
*vec.get_mut(4).unwrap() = 2.;
let expected = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 1., 2., 1.]);
assert_eq!(vec, expected);
vec[6] = 3.;
let expected = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 1., 2., 3.]);
assert_eq!(vec, expected);
}
#[test]
fn indexing() {
let vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 2., 3., 4.]);
assert_eq!(vec[0], 1.);
assert_eq!(vec[2], 2.);
assert_eq!(vec[4], 3.);
assert_eq!(vec[6], 4.);
}
#[test]
fn map_inplace() {
let mut vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 2., 3., 4.]);
vec.map_inplace(|&x| x + 1.);
let expected = CsVec::new(8, vec![0, 2, 4, 6], vec![2., 3., 4., 5.]);
assert_eq!(vec, expected);
}
#[test]
fn map() {
let vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 2., 3., 4.]);
let res = vec.map(|&x| x * 2.);
let expected = CsVec::new(8, vec![0, 2, 4, 6], vec![2., 4., 6., 8.]);
assert_eq!(res, expected);
}
#[test]
fn iter_mut() {
let mut vec = CsVec::new(8, vec![0, 2, 4, 6], vec![1., 2., 3., 4.]);
for (ind, val) in vec.iter_mut() {
if ind == 2 {
*val += 1.;
} else {
*val *= 2.;
}
}
let expected = CsVec::new(8, vec![0, 2, 4, 6], vec![2., 3., 6., 8.]);
assert_eq!(vec, expected);
}
#[test]
fn adds_vectors_by_value() {
let (a, b, expected_sum) = addition_sample();
assert_eq!(expected_sum, a + b);
}
#[test]
fn adds_vectors_by_left_value_and_right_reference() {
let (a, b, expected_sum) = addition_sample();
assert_eq!(expected_sum, a + &b);
}
#[test]
fn adds_vectors_by_left_reference_and_right_value() {
let (a, b, expected_sum) = addition_sample();
assert_eq!(expected_sum, &a + b);
}
#[test]
fn adds_vectors_by_reference() {
let (a, b, expected_sum) = addition_sample();
assert_eq!(expected_sum, &a + &b);
}
fn addition_sample() -> (CsVec<f64>, CsVec<f64>, CsVec<f64>) {
let dim = 8;
let a = CsVec::new(dim, vec![0, 3, 5, 7], vec![2., -3., 7., -1.]);
let b = CsVec::new(dim, vec![1, 3, 4, 5], vec![4., 2., -3., 1.]);
let expected_sum = CsVec::new(
dim,
vec![0, 1, 3, 4, 5, 7],
vec![2., 4., -1., -3., 8., -1.],
);
(a, b, expected_sum)
}
#[test]
fn negates_vectors() {
let vector = CsVec::new(4, vec![0, 3], vec![2., -3.]);
let negated = CsVec::new(4, vec![0, 3], vec![-2., 3.]);
assert_eq!(-vector, negated);
}
#[test]
fn can_construct_zero_sized_vectors() {
CsVec::<f64>::new(0, vec![], vec![]);
}
#[test]
fn zero_element_vanishes_when_added() {
let zero = CsVec::<f64>::zero();
let vector = CsVec::new(3, vec![0, 2], vec![1., 2.]);
assert_eq!(&vector + &zero, vector);
}
#[test]
fn zero_element_is_identified_as_zero() {
assert!(CsVec::<f32>::zero().is_zero());
}
#[test]
fn larger_zero_vector_is_identified_as_zero() {
let vector = CsVec::new(3, vec![1, 2], vec![0., 0.]);
assert!(vector.is_zero());
}
}