diff -r 9f2214c961cb -r cf8ef9463664 src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Fri Apr 28 14:02:18 2023 +0300
+++ b/src/bisection_tree/support.rs	Tue Jul 18 15:44:10 2023 +0300
@@ -145,22 +145,48 @@
     }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Differentiable<&'a Loc<F, N>> for Shift<T,F,N>
-where T : Differentiable<Loc<F, N>, Output=V> {
+impl<'a, T, V, W, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
+for Shift<T,F,N>
+where T : Differentiable<Loc<F, N>, Output=V> + Apply<Loc<F, N>, Output=W> {
     type Output = V;
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+    fn differential(&self, x : &'a Loc<F, N>) -> V {
         self.base_fn.differential(x - &self.shift)
     }
+
+    #[inline]
+    fn linearisation_error(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>) -> W {
+        self.base_fn
+            .linearisation_error(x - &self.shift, y - &self.shift)
+    }
+
+    #[inline]
+    fn linearisation_error_gen(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>, z : &'a Loc<F, N>) -> W {
+        self.base_fn
+            .linearisation_error_gen(x - &self.shift, y - &self.shift, z - &self.shift)
+    }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Differentiable<Loc<F, N>> for Shift<T,F,N>
-where T : Differentiable<Loc<F, N>, Output=V> {
+impl<'a, T, V, W, F : Float, const N : usize> Differentiable<Loc<F, N>>
+for Shift<T,F,N>
+where T : Differentiable<Loc<F, N>, Output=V> + Apply<Loc<F, N>, Output=W> {
     type Output = V;
     #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+    fn differential(&self, x : Loc<F, N>) -> V {
         self.base_fn.differential(x - &self.shift)
     }
+
+    #[inline]
+    fn linearisation_error(&self, x : Loc<F, N>, y : Loc<F, N>) -> W {
+        self.base_fn
+            .linearisation_error(x - &self.shift, y - &self.shift)
+    }
+
+    #[inline]
+    fn linearisation_error_gen(&self, x : Loc<F, N>, y : Loc<F, N>, z : Loc<F, N>) -> W {
+        self.base_fn
+            .linearisation_error_gen(x - &self.shift, y - &self.shift, z - &self.shift)
+    }
 }
 
 impl<'a, T, F : Float, const N : usize> Support<F,N> for Shift<T,F,N>
@@ -250,26 +276,52 @@
     }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<&'a Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Differentiable<&'b Loc<F, N>, Output=V>,
+impl<'a, T, V, W, F : Float, C, const N : usize> Differentiable<&'a Loc<F, N>> for Weighted<T, C>
+where T : for<'b> Differentiable<&'b Loc<F, N>, Output=V>
+          + for<'b> Apply<&'b Loc<F, N>, Output=W>,
       V : std::ops::Mul<F, Output=V>,
+      W : std::ops::Mul<F, Output=W>,
       C : Constant<Type=F> {
     type Output = V;
+
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
+    fn differential(&self, x : &'a Loc<F, N>) -> V {
         self.base_fn.differential(x) * self.weight.value()
     }
+
+    #[inline]
+    fn linearisation_error(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>) -> W {
+        self.base_fn.linearisation_error(x, y) * self.weight.value()
+    }
+
+    #[inline]
+    fn linearisation_error_gen(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>, z : &'a Loc<F, N>) -> W {
+        self.base_fn.linearisation_error_gen(x, y, z) * self.weight.value()
+    }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<Loc<F, N>> for Weighted<T, C>
-where T : Differentiable<Loc<F, N>, Output=V>,
+impl<'a, T, V, W, F : Float, C, const N : usize> Differentiable<Loc<F, N>>
+for Weighted<T, C>
+where T : Differentiable<Loc<F, N>, Output=V> + Apply<Loc<F, N>, Output=W>,
       V : std::ops::Mul<F, Output=V>,
+      W : std::ops::Mul<F, Output=W>,
       C : Constant<Type=F> {
     type Output = V;
+
     #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Output {
+    fn differential(&self, x : Loc<F, N>) -> V {
         self.base_fn.differential(x) * self.weight.value()
     }
+
+    #[inline]
+    fn linearisation_error(&self, x : Loc<F, N>, y : Loc<F, N>) -> W {
+        self.base_fn.linearisation_error(x, y) * self.weight.value()
+    }
+
+    #[inline]
+    fn linearisation_error_gen(&self, x : Loc<F, N>, y : Loc<F, N>, z : Loc<F, N>) -> W {
+        self.base_fn.linearisation_error_gen(x, y, z) * self.weight.value()
+    }
 }
 
 impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C>