# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1668896623 -7200
# Node ID 465fa2121ccb130324a7b514f06093738cc355fb
# Parent  3297d14c7bffd0a88ae9058580723e5981d4a3c5
Better Linear and Mapping structure that can provide consuming and reference `apply`.

diff -r 3297d14c7bff -r 465fa2121ccb src/bisection_tree/btfn.rs
--- a/src/bisection_tree/btfn.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/bisection_tree/btfn.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -4,8 +4,8 @@
 use std::marker::PhantomData;
 use std::sync::Arc;
 use crate::types::Float;
-use crate::mapping::Mapping;
-use crate::linops::Linear;
+use crate::mapping::{Apply, Mapping};
+//use crate::linops::{Apply, Linear};
 use crate::sets::Set;
 use crate::sets::Cube;
 use crate::loc::Loc;
@@ -392,33 +392,33 @@
 // Mapping
 //
 
-impl<'a, F : Float, G, BT, V, const N : usize> Mapping<&'a Loc<F,N>>
+impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a Loc<F, N>>
 for BTFN<F, G, BT, N>
 where BT : BTImpl<F, N>,
       G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<&'a Loc<F,N>, Codomain = V>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<&'a Loc<F, N>, Output = V>,
       V : Sum {
 
-    type Codomain = V;
+    type Output = V;
 
-    fn value(&self, x : &'a Loc<F,N>) -> Self::Codomain {
+    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
         self.bt.iter_at(x)
-            .map(|&d| self.generator.support_for(d).value(x)).sum()
+            .map(|&d| self.generator.support_for(d).apply(x)).sum()
     }
 }
 
-impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F,N>>
+impl<F : Float, G, BT, V, const N : usize> Apply<Loc<F, N>>
 for BTFN<F, G, BT, N>
 where BT : BTImpl<F, N>,
       G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F,N>, Codomain = V>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<Loc<F, N>, Output = V>,
       V : Sum {
 
-    type Codomain = V;
+    type Output = V;
 
-    fn value(&self, x : Loc<F,N>) -> Self::Codomain {
+    fn apply(&self, x : Loc<F, N>) -> Self::Output {
         self.bt.iter_at(&x)
-            .map(|&d| self.generator.support_for(d).value(x)).sum()
+            .map(|&d| self.generator.support_for(d).apply(x)).sum()
     }
 }
 
@@ -439,19 +439,42 @@
 // that are linear functionals over BTFN.
 //
 
+/*
+impl<'b, X, F : Float, G, BT, const N : usize> Apply<&'b X, F>
+for BTFN<F, G, BT, N>
+where BT : BTImpl<F, N>,
+      G : SupportGenerator<F, N, Id=BT::Data>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {
+
+    #[inline]
+    fn apply(&self, x : &'b X) -> F {
+        x.apply(&self)
+    }
+}
+
+impl<X, F : Float, G, BT, const N : usize> Apply<X, F>
+for BTFN<F, G, BT, N>
+where BT : BTImpl<F, N>,
+      G : SupportGenerator<F, N, Id=BT::Data>,
+      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {
+
+    #[inline]
+    fn apply(&self, x : X) -> F {
+        x.apply(&self)
+    }
+}
+
 impl<X, F : Float, G, BT, const N : usize> Linear<X>
 for BTFN<F, G, BT, N>
 where BT : BTImpl<F, N>,
       G : SupportGenerator<F, N, Id=BT::Data>,
       G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-      X : Linear<BTFN<F, G, BT, N>, Codomain=F> {
+      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {
     type Codomain = F;
-
-    #[inline]
-    fn apply(&self, x : &X) -> F {
-        x.apply(self)
-    }
 }
+*/
 
 /// Helper trait for performing approximate minimisation using P2 elements.
 ///
@@ -564,7 +587,7 @@
 for P2Refiner<F, RefineMax>
 where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
       G : SupportGenerator<F, N>,
-      G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F>
+      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
                        + LocalAnalysis<F, Bounds<F>, N> {
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = UpperBoundSorting<F>;
@@ -584,8 +607,8 @@
         }
 
         // g gives the negative of the value of the function presented by `data` and `generator`.
-        let g = move |x : &Loc<F,N>| {
-            let f = move |&d| generator.support_for(d).value(x);
+        let g = move |x : &Loc<F, N>| {
+            let f = move |&d| generator.support_for(d).apply(x);
             -data.iter().map(f).sum::<F>()
         };
         // … so the negative of the minimum is the maximm we want.
@@ -621,7 +644,7 @@
 for P2Refiner<F, RefineMin>
 where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
       G : SupportGenerator<F, N>,
-      G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F>
+      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
                        + LocalAnalysis<F, Bounds<F>, N> {
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = LowerBoundSorting<F>;
@@ -641,8 +664,8 @@
         }
 
         // g gives the value of the function presented by `data` and `generator`.
-        let g = move |x : &Loc<F,N>| {
-            let f = move |&d| generator.support_for(d).value(x);
+        let g = move |x : &Loc<F, N>| {
+            let f = move |&d| generator.support_for(d).apply(x);
             data.iter().map(f).sum::<F>()
         };
         // Minimise it.
@@ -783,7 +806,7 @@
 impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N>
 where BT : BTSearch<F, N, Agg=Bounds<F>>,
       G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F>
+      G::SupportType : Mapping<Loc<F, N>,Codomain=F>
                        + LocalAnalysis<F, Bounds<F>, N>,
       Cube<F, N> : P2Minimise<Loc<F, N>, F> {
 
diff -r 3297d14c7bff -r 465fa2121ccb src/bisection_tree/either.rs
--- a/src/bisection_tree/either.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/bisection_tree/either.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -3,7 +3,7 @@
 use std::sync::Arc;
 
 use crate::types::*;
-use crate::mapping::Mapping;
+use crate::mapping::Apply;
 use crate::iter::{Mappable,MapF,MapZ};
 use crate::sets::Cube;
 use crate::loc::Loc;
@@ -177,15 +177,15 @@
     }
 }
 
-impl<F, S1, S2, X> Mapping<X> for EitherSupport<S1, S2>
-where S1 : Mapping<X, Codomain=F>,
-      S2 : Mapping<X, Codomain=F> {
-    type Codomain = F;
+impl<F, S1, S2, X> Apply<X> for EitherSupport<S1, S2>
+where S1 : Apply<X, Output=F>,
+      S2 : Apply<X, Output=F> {
+    type Output = F;
     #[inline]
-    fn value(&self, x : X) -> F {
+    fn apply(&self, x : X) -> F {
         match self {
-            EitherSupport::Left(ref a) => a.value(x),
-            EitherSupport::Right(ref b) => b.value(x),
+            EitherSupport::Left(ref a) => a.apply(x),
+            EitherSupport::Right(ref b) => b.apply(x),
         }
     }
 }
diff -r 3297d14c7bff -r 465fa2121ccb src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/bisection_tree/support.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -6,7 +6,7 @@
 use std::ops::{MulAssign,DivAssign,Neg};
 use crate::types::{Float, Num};
 use crate::maputil::map2;
-use crate::mapping::Mapping;
+use crate::mapping::Apply;
 use crate::sets::Cube;
 use crate::loc::Loc;
 use super::aggregator::Bounds;
@@ -127,21 +127,21 @@
     base_fn : T,
 }
 
-impl<'a, T, V, F : Float, const N : usize> Mapping<&'a Loc<F, N>> for Shift<T,F,N>
-where T : for<'b> Mapping<&'b Loc<F,N>,Codomain=V> {
-    type Codomain = V;
+impl<'a, T, V, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Shift<T,F,N>
+where T : Apply<Loc<F, N>, Output=V> {
+    type Output = V;
     #[inline]
-    fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain {
-        self.base_fn.value(&(x - &self.shift))
+    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.base_fn.apply(x - &self.shift)
     }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Mapping<Loc<F, N>> for Shift<T,F,N>
-where T : for<'b> Mapping<Loc<F,N>,Codomain=V> {
-    type Codomain = V;
+impl<'a, T, V, F : Float, const N : usize> Apply<Loc<F, N>> for Shift<T,F,N>
+where T : Apply<Loc<F, N>, Output=V> {
+    type Output = V;
     #[inline]
-    fn value(&self, x : Loc<F, N>) -> Self::Codomain {
-        self.base_fn.value(x - &self.shift)
+    fn apply(&self, x : Loc<F, N>) -> Self::Output {
+        self.base_fn.apply(x - &self.shift)
     }
 }
 
@@ -210,25 +210,25 @@
     pub base_fn : T,
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Mapping<&'a Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Mapping<&'b Loc<F,N>,Codomain=V>,
+impl<'a, T, V, F : Float, C, const N : usize> Apply<&'a Loc<F, N>> for Weighted<T, C>
+where T : for<'b> Apply<&'b Loc<F, N>, Output=V>,
       V : std::ops::Mul<F,Output=V>,
       C : Constant<Type=F> {
-    type Codomain = V;
+    type Output = V;
     #[inline]
-    fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain {
-        self.base_fn.value(x) * self.weight.value()
+    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
+        self.base_fn.apply(x) * self.weight.value()
     }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Mapping<Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Mapping<Loc<F,N>,Codomain=V>,
+impl<'a, T, V, F : Float, C, const N : usize> Apply<Loc<F, N>> for Weighted<T, C>
+where T : Apply<Loc<F, N>, Output=V>,
       V : std::ops::Mul<F,Output=V>,
       C : Constant<Type=F> {
-    type Codomain = V;
+    type Output = V;
     #[inline]
-    fn value(&self, x : Loc<F, N>) -> Self::Codomain {
-        self.base_fn.value(x) * self.weight.value()
+    fn apply(&self, x : Loc<F, N>) -> Self::Output {
+        self.base_fn.apply(x) * self.weight.value()
     }
 }
 
@@ -340,23 +340,23 @@
     pub T
 );
 
-impl<'a, T, F : Float, const N : usize> Mapping<&'a Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + for<'b> Mapping<&'b Loc<F,N>, Codomain=F> {
-    type Codomain = F;
+impl<'a, T, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Normalised<T>
+where T : Norm<F, L1> + for<'b> Apply<&'b Loc<F, N>, Output=F> {
+    type Output = F;
     #[inline]
-    fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain {
+    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
         let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { self.0.value(x) / w }
+        if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w }
     }
 }
 
-impl<'a, T, F : Float, const N : usize> Mapping<Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + for<'b> Mapping<Loc<F,N>, Codomain=F> {
-    type Codomain = F;
+impl<'a, T, F : Float, const N : usize> Apply<Loc<F, N>> for Normalised<T>
+where T : Norm<F, L1> + Apply<Loc<F,N>, Output=F> {
+    type Output = F;
     #[inline]
-    fn value(&self, x : Loc<F, N>) -> Self::Codomain {
+    fn apply(&self, x : Loc<F, N>) -> Self::Output {
         let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { self.0.value(x) / w }
+        if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w }
     }
 }
 
diff -r 3297d14c7bff -r 465fa2121ccb src/linops.rs
--- a/src/linops.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/linops.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -6,13 +6,12 @@
 use std::marker::PhantomData;
 use crate::types::*;
 use serde::Serialize;
+pub use crate::mapping::Apply;
 
 /// Trait for linear operators on `X`.
-pub trait Linear<X> {
-    /// The range space of the operator.
+pub trait Linear<X> : Apply<X, Output=Self::Codomain>
+                      + for<'a> Apply<&'a X, Output=Self::Codomain> {
     type Codomain;
-    /// Apply the linear operator to `x`.
-    fn apply(&self, x : &X) -> Self::Codomain;
 }
 
 /// Efficient in-place summation.
@@ -103,19 +102,33 @@
 
 /// The identity operator
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct IdOp<X : Clone> (PhantomData<X>);
+pub struct IdOp<X> (PhantomData<X>);
+
+impl<X> IdOp<X> {
+    fn new() -> IdOp<X> { IdOp(PhantomData) }
+}
+
+impl<X> Apply<X> for IdOp<X> {
+    type Output = X;
 
-impl<X> IdOp<X> where X : Clone {
-    fn new() -> IdOp<X> { IdOp(PhantomData) }
+    fn apply(&self, x : X) -> X {
+        x
+    }
+}
+
+impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
+    type Output = X;
+    
+    fn apply(&self, x : &'a X) -> X {
+        x.clone()
+    }
 }
 
 impl<X> Linear<X> for IdOp<X> where X : Clone {
     type Codomain = X;
-    fn apply(&self, x : &X) -> X {
-        x.clone()
-    }
 }
 
+
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
diff -r 3297d14c7bff -r 465fa2121ccb src/mapping.rs
--- a/src/mapping.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/mapping.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -7,37 +7,75 @@
 use serde::Serialize;
 use crate::loc::Loc;
 
-/// A mapping from `Domain` to `Codomain`.
-pub trait Mapping<Domain> {
-    type Codomain;
+/// Trait for application of `Self` as a mathematical function or operator on `X`.
+pub trait Apply<X> {
+    type Output;
+
+    /// Compute the value of `self` at `x`.
+    fn apply(&self, x : X) -> Self::Output;
+}
 
-    /// Calculate the value of the mapping at `x`.
-    fn value(&self, x : Domain) -> Self::Codomain;
+/// This blanket implementation is a workaround helper to Rust trait system limitations.
+///
+/// It is introduced because the trait system does not allow blanket implementations of both
+/// [`Apply<X>`] and [`Apply<&'a X>`]. With this, the latter is implemented automatically for
+/// the reference, which can be sufficient to apply the operation in another blanket implementation.
+impl<'a, T, X> Apply<X> for &'a T where T : Apply<X> {
+    type Output = <T as Apply<X>>::Output;
+
+    #[inline]
+    fn apply(&self, x : X) -> Self::Output {
+        (*self).apply(x)
+    }
 }
 
-/// A helper trait alias for referring to `Mapping`s from references to floats.
-pub trait RealRefMapping<F : Float, const N : usize>
-: for<'a> Mapping<&'a Loc<F, N>, Codomain=F> {}
+/// A mapping from `Domain` to `Codomain`.
+///
+/// This is automatically implemented when the relevant [`Apply`] are implemented.
+pub trait Mapping<Domain> : Apply<Domain, Output=Self::Codomain>
+                            + for<'a> Apply<&'a Domain, Output=Self::Codomain> {
+    type Codomain;
+}
+
+impl<Domain, Codomain, T> Mapping<Domain> for T
+where T : Apply<Domain, Output=Codomain> + for<'a> Apply<&'a Domain, Output=Codomain> {
+    type Codomain = Codomain;
+}
+
 
-impl<F : Float, T, const N : usize> RealRefMapping<F, N> for T
-where T : for<'a> Mapping<&'a Loc<F, N>, Codomain=F> {}
+/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to `F` a [`Float`].
+pub trait RealMapping<F : Float, const N : usize> : Mapping<Loc<F, N>, Codomain = F> {}
+
+impl<F : Float, T, const N : usize> RealMapping<F, N> for T
+where T : Mapping<Loc<F, N>, Codomain = F> {}
+
+
+/// Trait for calculation the differential of `Self` as a mathematical function on `X`.
+pub trait Differentiate<X> {
+    type Output;
+
+    /// Compute the differential of `self` at `x`.
+    fn differential(&self, x : X) -> Self::Output;
+}
 
 
 /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials
 /// `Differential`.
-pub trait DifferentiableMapping<Domain> : Mapping<Domain> {
+///
+/// This is automatically implemented when the relevant [`Differentiate`] are implemented.
+pub trait DifferentiableMapping<Domain>
+: Mapping<Domain>
+  + Differentiate<Domain, Output=Self::Differential>
+  + for<'a> Differentiate<&'a Domain, Output=Self::Differential>{
     type Differential;
-
-    /// Calculate the differentialeof the mapping at `x`.
-    fn differential(&self, x : Domain) -> Self::Differential;
 }
 
-/// A `Mapping` whose minimum and maximum can be computed.
-pub trait RealMapping<Domain> : Mapping<Domain> where Self::Codomain : Float {
-    /// Calculate a minimum and a minimiser of the mapping.
-    fn minimise(&self, tolerance : Self::Codomain) -> (Domain, Self::Codomain);
-    /// Calculate a maximum and a maximiser of the mapping.
-    fn maximise(&self, tolerance : Self::Codomain) -> (Domain, Self::Codomain);
+
+impl<Domain, Differential, T> DifferentiableMapping<Domain> for T
+where T : Mapping<Domain>
+          + Differentiate<Domain, Output=Differential>
+          + for<'a> Differentiate<&'a Domain, Output=Differential> {
+    type Differential = Differential;
 }
 
 /// A sum of [`Mapping`]s.
@@ -55,27 +93,25 @@
 }
 
 
-impl<Domain, M> Mapping<Domain> for Sum<Domain,M>
+impl<Domain : Copy, M> Apply<Domain> for Sum<Domain, M>
 where M : Mapping<Domain>,
-      M :: Codomain : std::iter::Sum,
-      Domain : Copy {
+      M::Codomain : std::iter::Sum {
+    type Output = M::Codomain;
 
-    type Codomain = M::Codomain;
-
-    fn value(&self, x : Domain) -> Self::Codomain {
-        self.components.iter().map(|c| c.value(x)).sum()
+    fn apply(&self, x : Domain) -> Self::Output {
+        self.components.iter().map(|c| c.apply(x)).sum()
     }
 }
 
-impl<Domain, M> DifferentiableMapping<Domain> for Sum<Domain,M>
+impl<Domain, M> Differentiate<Domain> for Sum<Domain, M>
 where M : DifferentiableMapping<Domain>,
       M :: Codomain : std::iter::Sum,
       M :: Differential : std::iter::Sum,
       Domain : Copy {
 
-    type Differential = M::Differential;
+    type Output = M::Differential;
 
-    fn differential(&self, x : Domain) -> Self::Differential {
+    fn differential(&self, x : Domain) -> Self::Output {
         self.components.iter().map(|c| c.differential(x)).sum()
     }
 }
diff -r 3297d14c7bff -r 465fa2121ccb src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/nalgebra_support.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -26,7 +26,37 @@
 use crate::types::Float;
 use crate::norms::*;
 
-impl<SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
+impl<SM,SV,N,M,K,E> Apply<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
+where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
+        N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
+        DefaultAllocator : Allocator<E,N,K>,
+        DefaultAllocator : Allocator<E,M,K>,
+        DefaultAllocator : Allocator<E,N,M>,
+        DefaultAllocator : Allocator<E,M,N> {
+    type Output = OMatrix<E,N,K>;
+
+    #[inline]
+    fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output {
+        self.mul(x)
+    }
+}
+
+impl<'a, SM,SV,N,M,K,E> Apply<&'a Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
+where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
+        N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
+        DefaultAllocator : Allocator<E,N,K>,
+        DefaultAllocator : Allocator<E,M,K>,
+        DefaultAllocator : Allocator<E,N,M>,
+        DefaultAllocator : Allocator<E,M,N> {
+    type Output = OMatrix<E,N,K>;
+
+    #[inline]
+    fn apply(&self, x : &'a Matrix<E,M,K,SV>) -> Self::Output {
+        self.mul(x)
+    }
+}
+
+impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
 where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
         N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
         DefaultAllocator : Allocator<E,N,K>,
@@ -34,11 +64,6 @@
         DefaultAllocator : Allocator<E,N,M>,
         DefaultAllocator : Allocator<E,M,N> {
     type Codomain = OMatrix<E,N,K>;
-
-    #[inline]
-    fn apply(&self, x : &Matrix<E,M,K,SV>) -> Self::Codomain {
-        self.mul(x)
-    }
 }
 
 impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM>
diff -r 3297d14c7bff -r 465fa2121ccb src/types.rs
--- a/src/types.rs	Fri Nov 18 10:34:04 2022 +0200
+++ b/src/types.rs	Sun Nov 20 00:23:43 2022 +0200
@@ -159,3 +159,4 @@
     pub trait CompatibleSigned<F : Float> = Signed + Into<F>;
 }
 */
+