From 7facca7fcc7ee23301bba4f0a05f5eefec3fd381 Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Wed, 23 Oct 2019 22:49:20 +0700
Subject: [PATCH 01/12] wip

---
 source/mir/interpolate/pchip.d   | 105 +---------
 source/mir/interpolate/spline.d  | 330 ++++++++++++++++++++++++++-----
 source/mir/interpolate/utility.d |  20 ++
 3 files changed, 312 insertions(+), 143 deletions(-)

diff --git a/source/mir/interpolate/pchip.d b/source/mir/interpolate/pchip.d
index 01ab3b54..2ebd377d 100644
--- a/source/mir/interpolate/pchip.d
+++ b/source/mir/interpolate/pchip.d
@@ -11,22 +11,25 @@ Macros:
 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolate, $1)$(NBSP)
 T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
 +/
+deprecated
 module mir.interpolate.pchip;
 
 public import mir.interpolate.spline: Spline;
-
-import std.traits;
-private alias AliasSeq(T...) = T;
+import mir.interpolate;
+import mir.interpolate.spline;
+import mir.math.common: fmamath;
 import mir.ndslice.slice;
 import mir.primitives;
-import mir.math.common: fmamath;
-import mir.ndslice.internal: ConstIfPointer;
+import std.traits;
+
+private alias AliasSeq(T...) = T;
 
 @fmamath:
 
 /++
 Constructs piecewise cubic hermite interpolating polynomial with nodes on rectilinear grid.
 +/
+deprecated
 template pchip(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIterators = Repeat!(N - 1, FirstGridIterator))
     if (isFloatingPoint!T && is(T == Unqual!T) && N <= 6)
 {
@@ -48,11 +51,7 @@ template pchip(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridItera
         GridVectors grid,
         Slice!(yIterator, N, ykind) values) @safe
     {
-        static if (__VERSION__ >= 2085) import core.lifetime: move; else import std.algorithm.mutation: move; 
-        auto ret = Spline!T(grid);
-        ret._values = values;
-        pchipSlopes(grid, values, typeof(return).pickDataSubslice(ret._data, 1));
-        return ret.move;
+        return spline!T(grid, values, SplineKind.monotone);
     }
 }
 
@@ -121,89 +120,3 @@ version(mir_test)
     version(X86_64)
     assert(map!interpolant(points[0 .. $ - 1] +  0.5) == map!interpolantR(pointsR.retro[0 .. $ - 1] - 0.5));
 }
-
-
-/++
-Computes slopes for piecewise spline hermite interpolating polynomial.
-Params:
-    points = `x` values for interpolant
-    values = `f(x)` values for interpolant
-    slopes = uninitialized ndslice to write slopes into
-Constraints:
-    `points`, `values`, and `slopes` must have the same length >= 3
-+/
-void pchipSlopes(IG, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind skind)(
-    Slice!(IG, 1, gkind) points,
-    Slice!(IV, 1, vkind) values,
-    Slice!(IS, 1, skind) slopes) @trusted
-{
-    if (points.length < 3)
-        assert(0);
-    if (points.length != values.length)
-        assert(0);
-    if (points.length != slopes.length)
-        assert(0);
-
-    auto step0 = cast()(points[1] - points[0]);
-    auto step1 = cast()(points[2] - points[1]);
-    auto diff0 = cast()(values[1] - values[0]);
-    auto diff1 = cast()(values[2] - values[1]);
-    diff0 /= step0;
-    diff1 /= step1;
-
-    slopes.front = pchipTail(step0, step1, diff0, diff1);
-    for(size_t i = 1;;)
-    {
-        import mir.math.common: copysign;
-        if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1))
-        {
-            auto w0 = step1 * 2 + step0;
-            auto w1 = step0 * 2 + step1;
-            slopes[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1);
-        }
-        else
-        {
-            slopes[i] = 0;
-        }
-        if (++i == slopes.length - 1)
-        {
-            break;
-        }
-        step0 = step1;
-        diff0 = diff1;
-        step1 = points[i + 1] - points[i + 0];
-        diff1 = values[i + 1] - values[i + 0];
-        diff1 /= step1;
-    }
-    slopes.back = pchipTail(step1, step0, diff1, diff0);
-}
-
-auto pchipTail(T)(in T step0, in T step1, in T diff0, in T diff1)
-{
-    import mir.math.common: copysign, fabs;
-    if (!diff0)
-    {
-        return 0;
-    }
-    auto slope = ((step0 * 2 + step1) * diff0 - step0 * diff1) / (step0 + step1);
-    if (copysign(1f, slope) != copysign(1f, diff0))
-    {
-        return 0;
-    }
-    if ((copysign(1f, diff0) != copysign(1f, diff1)) && (fabs(slope) > fabs(diff0 * 3)))
-    {
-        return diff0 * 3;
-    }
-    return slope;
-}
-
-template Repeat(ulong n, L...)
-{
-    static if (n)
-    {
-        import std.meta: Repeat;
-        alias Repeat = std.meta.Repeat!(n, L);
-    }
-    else
-        alias Repeat = L[0 .. 0];
-}
diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index c029deaa..7a405bc9 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -332,6 +332,73 @@ unittest
     // assert(appreq(d[1][1][1], y_x0x1x2));
 }
 
+
+/// Monotone PCHIP
+version(mir_test)
+@safe unittest
+{
+    import mir.math.common: approxEqual;
+    import mir.algorithm.iteration: all;
+    import mir.ndslice.allocation: slice;
+    import mir.ndslice.slice: sliced;
+    import mir.ndslice.topology: vmap;
+
+    auto x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
+    auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;
+    auto interpolant = spline!double(x, y, SplineKind.monotone);
+
+    auto xs = x[0 .. $ - 1] + 0.5;
+
+    () @trusted {
+        auto ys = xs.vmap(interpolant);
+
+        assert(ys.all!approxEqual([
+            5.333333333333334,
+            2.500000000000000,
+            10.000000000000000,
+            4.288971807628524,
+            11.202580845771145,
+            16.250000000000000,
+            17.962962962962962,
+            5.558593750000000,
+            17.604662698412699,
+            ]));
+    }();
+}
+
+// Check direction equality
+version(mir_test)
+@safe unittest
+{
+    import mir.math.common: approxEqual;
+    import mir.ndslice.slice: sliced;
+    import mir.ndslice.allocation: slice;
+    import mir.ndslice.topology: retro, map;
+
+    auto points = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
+    auto values = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;
+
+    auto results = [
+        5.333333333333334,
+        2.500000000000000,
+        10.000000000000000,
+        4.288971807628524,
+        11.202580845771145,
+        16.250000000000000,
+        17.962962962962962,
+        5.558593750000000,
+        17.604662698412699,
+        ];
+    auto interpolant = spline!double(points, values, SplineKind.monotone);
+
+    auto pointsR = slice(-points.retro);
+    auto valuesR = values.retro.slice;
+    auto interpolantR = spline!double(pointsR, valuesR, SplineKind.monotone);
+
+    version(X86_64)
+    assert(map!interpolant(points[0 .. $ - 1] +  0.5) == map!interpolantR(pointsR.retro[0 .. $ - 1] - 0.5));
+}
+
 /++
 Constructs multivariate cubic spline in symmetrical form with nodes on rectilinear grid.
 Result has continues second derivatives throughout the curve / nd-surface.
@@ -359,7 +426,19 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         in T valueOfBoundaryConditions = 0,
         )
     {
-        return spline(grid, values, SplineBoundaryCondition!T(typeOfBoundaries, valueOfBoundaryConditions));
+        return spline(grid, values, SplineKind.c2, 0, typeOfBoundaries, valueOfBoundaryConditions);
+    }
+
+    Spline!(T, N, GridIterators) spline(yIterator, SliceKind ykind)(
+        GridVectors grid,
+        Slice!(yIterator, N, ykind) values,
+        SplineKind kind,
+        in T param = 0,
+        SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot,
+        in T valueOfBoundaryConditions = 0,
+        )
+    {
+        return spline(grid, values, SplineBoundaryCondition!T(typeOfBoundaries, valueOfBoundaryConditions), kind, param);
     }
 
     /++
@@ -375,9 +454,11 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         GridVectors grid,
         Slice!(yIterator, N, ykind) values,
         SplineBoundaryCondition!T boundaries,
+        SplineKind kind = SplineKind.c2,
+        in T param = 0,
         )
     {
-        return spline(grid, values, boundaries, boundaries);
+        return spline(grid, values, boundaries, boundaries, kind, param);
     }
 
     /++
@@ -395,12 +476,14 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         Slice!(yIterator, N, ykind) values,
         SplineBoundaryCondition!T rBoundary,
         SplineBoundaryCondition!T lBoundary,
+        SplineKind kind = SplineKind.c2,
+        in T param = 0,
         )
     {
         static if (__VERSION__ >= 2085) import core.lifetime: move; else import std.algorithm.mutation: move; 
         auto ret = typeof(return)(grid);
         ret._values = values;
-        ret._computeDerivatives(rBoundary, lBoundary);
+        ret._computeDerivatives(kind, param, rBoundary, lBoundary);
         return ret.move;
     }
 }
@@ -412,16 +495,31 @@ See_also: $(LREF SplineBoundaryCondition)
 +/
 enum SplineBoundaryType
 {
-    /// not implemented
+    /++
+     not implemented
+    +/
     periodic = -1,
-    /// (default)
+    /++
+     (default)
+    +/
     notAKnot,
-    /// set the first derivative
+    /++
+     set the first derivative
+    +/
     firstDerivative,
-    /// set the second derivative
+    /++
+     set the second derivative
+    +/
     secondDerivative,
-    ///
+    /++
+    +/
     parabolic,
+    /++
+    +/
+    monotone,
+    /++
+    +/
+    akima,
 }
 
 /++
@@ -439,6 +537,27 @@ struct SplineBoundaryCondition(T)
 
 // private auto iter(alias s) = s.iterator;
 
+/++
++/
+enum SplineKind
+{
+    /++
+    +/
+    c2,
+    /++
+    +/
+    cardinal,
+    /++
+    +/
+    monotone,
+    /++
+    +/
+    doubleQuadratic,
+    /++
+    +/
+    akima,
+}
+
 /++
 Multivariate cubic spline with nodes on rectilinear grid.
 +/
@@ -465,7 +584,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
     {
         size_t length = 1;
         size_t[N] shape;
-        enum  msg =  "spline/pchip interpolant: minimal allowed length for the grid equals 2.";
+        enum  msg =  "spline interpolant: minimal allowed length for the grid equals 2.";
         version(D_Exceptions)
             static immutable exc = new Exception(msg);
         foreach(i, ref x; grid)
@@ -513,18 +632,18 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
 
     $(RED For internal use.)
     +/
-    void _computeDerivatives()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc
+    void _computeDerivatives(SplineKind kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc
     {
         import mir.algorithm.iteration: maxLength;
         auto ml = this._data.maxLength;
         auto temp = RCArray!F(ml);
         auto tempSlice = temp[].sliced;
-        _computeDerivativesTemp(lbc, rbc, tempSlice);
+        _computeDerivativesTemp(kind, param, lbc, rbc, tempSlice);
     }
 
     /// ditto
     pragma(inline, false)
-    void _computeDerivativesTemp()(SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc
+    void _computeDerivativesTemp(SplineKind kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc
     {
         import mir.algorithm.iteration: maxLength, each;
         import mir.ndslice.topology: map, byDim, evertPack;
@@ -533,7 +652,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
 
         static if (N == 1)
         {
-            splineSlopes!(F, F)(_grid.sliced(_data._lengths[0]), pickDataSubslice(_data, 0), pickDataSubslice(_data, 1), temp[0 .. _data._lengths[0]], lbc, rbc);
+            splineSlopes!(F, F)(_grid.sliced(_data._lengths[0]), pickDataSubslice(_data, 0), pickDataSubslice(_data, 1), temp[0 .. _data._lengths[0]], kind, param, lbc, rbc);
         }
         else
         foreach_reverse(i, ref x; _grid)
@@ -549,7 +668,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
                         auto y = pickDataSubslice(d, l);
                         auto s = pickDataSubslice(d, L + l);
                         // debug printf("ptr = %ld, stride = %ld, stride = %ld, d = %ld i = %ld l = %ld\n", d.iterator, d._stride!0, y._stride!0, d.length, i, l);
-                        splineSlopes!(F, F)(x.sliced(_data._lengths[i]), y, s, temp[0 .. _data._lengths[i]], lbc, rbc);
+                        splineSlopes!(F, F)(x.sliced(_data._lengths[i]), y, s, temp[0 .. _data._lengths[i]], kind, param, lbc, rbc);
                         // debug{
                         //     (cast(void delegate() @nogc)(){
                         //     writeln("y = ", y);
@@ -564,9 +683,9 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
 @trusted:
 
     ///
-    Spline lightConst()() const @property { return *cast(Spline*)&this; }
+    Spline lightConst() const @property { return *cast(Spline*)&this; }
     ///
-    Spline lightImmutable()() immutable @property { return *cast(Spline*)&this; }
+    Spline lightImmutable() immutable @property { return *cast(Spline*)&this; }
 
     ///
     GridVectors[dimension] grid(size_t dimension = 0)() scope return const @property
@@ -585,7 +704,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
     }
 
     ///
-    size_t[N] gridShape()() scope const @property
+    size_t[N] gridShape() scope const @property
     {
         return _data.shape;
     }
@@ -774,6 +893,8 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     Slice!(IV, 1, vkind) values,
     Slice!(IS, 1, skind) slopes,
     Slice!(T*) temp,
+    SplineKind kind,
+    F param,
     SplineBoundaryCondition!F lbc,
     SplineBoundaryCondition!F rbc,
     ) @trusted
@@ -793,27 +914,61 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     auto xd = points.diff;
     auto yd = values.diff;
 
-    /// special case
-    static assert(SplineBoundaryType.notAKnot == 0);
-    if (n <= 3
-        && (rbc.type == SplineBoundaryType.parabolic || rbc.type == SplineBoundaryType.notAKnot)
-        && (lbc.type == SplineBoundaryType.parabolic || lbc.type == SplineBoundaryType.notAKnot)
-    )
+    with(SplineKind) final switch(kind)
     {
-        import mir.interpolate.utility;
-        if (n == 3)
-        {
-            auto derivatives = parabolaDerivatives(points[0], points[1], points[2], values[0], values[1], values[2]);
-            slopes[0] = derivatives[0];
-            slopes[1] = derivatives[1];
-            slopes[2] = derivatives[2];
-        }
-        else
+        case c2:
+            break;
+        case cardinal:
+            if (lbc.type == SplineBoundaryType.notAKnot)
+                lbc.type = SplineBoundaryType.parabolic;
+            if (rbc.type == SplineBoundaryType.notAKnot)
+                rbc.type = SplineBoundaryType.parabolic;
+            break;
+        case monotone:
+            if (lbc.type == SplineBoundaryType.notAKnot)
+                lbc.type = SplineBoundaryType.monotone;
+            if (rbc.type == SplineBoundaryType.notAKnot)
+                rbc.type = SplineBoundaryType.monotone;
+            break;
+        case doubleQuadratic:
+            if (lbc.type == SplineBoundaryType.notAKnot)
+                lbc.type = SplineBoundaryType.parabolic;
+            if (rbc.type == SplineBoundaryType.notAKnot)
+                rbc.type = SplineBoundaryType.parabolic;
+            break;
+        case akima:
+            if (lbc.type == SplineBoundaryType.notAKnot)
+                lbc.type = SplineBoundaryType.akima;
+            if (rbc.type == SplineBoundaryType.notAKnot)
+                rbc.type = SplineBoundaryType.akima;
+            break;
+    }
+
+    if (n <= 3)
+    {
+        if (lbc.type == SplineBoundaryType.notAKnot)
+            lbc.type = SplineBoundaryType.parabolic;
+        if (rbc.type == SplineBoundaryType.notAKnot)
+            rbc.type = SplineBoundaryType.parabolic;
+        /// special case
+        if (rbc.type == SplineBoundaryType.parabolic
+         && lbc.type == SplineBoundaryType.parabolic)
         {
-            assert(slopes.length == 2);
-            slopes.back = slopes.front = yd.front / xd.front;
+            import mir.interpolate.utility;
+            if (n == 3)
+            {
+                auto derivatives = parabolaDerivatives(points[0], points[1], points[2], values[0], values[1], values[2]);
+                slopes[0] = derivatives[0];
+                slopes[1] = derivatives[1];
+                slopes[2] = derivatives[2];
+            }
+            else
+            {
+                assert(slopes.length == 2);
+                slopes.back = slopes.front = yd.front / xd.front;
+            }
+            return;
         }
-        return;
     }
 
     with(SplineBoundaryType) final switch(lbc.type)
@@ -843,7 +998,6 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             slopes.front = 1;
             first = 0;
             temp.front = yd.front / xd.front;
-
         }
         break;
     
@@ -867,6 +1021,19 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         first = 1;
         temp.front = 2 * (yd.front / xd.front);
         break;
+    
+    case monotone:
+        slopes.front = 1;
+        first = 0;
+        temp.front = pchipTail(xd[0], xd[1], yd[0]/xd[0], yd[1]/xd[1]);
+        break;
+
+    case akima:
+        // TODO
+        slopes.front = 2;
+        first = 1;
+        temp.front = 2 * (yd.front / xd.front);
+        break;
     }
 
     with(SplineBoundaryType) final switch(rbc.type)
@@ -915,22 +1082,72 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         last = 1;
         temp.back = 2 * (yd.back / xd.back);
         break;
+
+    case monotone:
+        slopes.back = 1;
+        last = 0;
+        temp.back = pchipTail(xd[$ - 1], xd[$ - 2], yd[$ - 1] / xd[$ - 1], yd[$ - 2] / xd[$ - 2]);
+        break;
+
+    case akima:
+        // TODO
+        slopes.back = 2;
+        last = 1;
+        temp.back = 2 * (yd.back / xd.back);
+        break;
     }
 
-    foreach (i; 1 .. n - 1)
+    if (kind == SplineKind.monotone)
+    {
+        auto step0 = cast()(points[1] - points[0]);
+        auto step1 = cast()(points[2] - points[1]);
+        auto diff0 = cast()(values[1] - values[0]);
+        auto diff1 = cast()(values[2] - values[1]);
+        diff0 /= step0;
+        diff1 /= step1;
+
+        for(size_t i = 1;;)
+        {
+            import mir.math.common: copysign;
+            slopes[i] = 1;
+            if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1))
+            {
+                auto w0 = step1 * 2 + step0;
+                auto w1 = step0 * 2 + step1;
+                temp[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1);
+            }
+            else
+            {
+                temp[i] = 0;
+            }
+            if (++i == n - 1)
+            {
+                break;
+            }
+            step0 = step1;
+            diff0 = diff1;
+            step1 = points[i + 1] - points[i + 0];
+            diff1 = values[i + 1] - values[i + 0];
+            diff1 /= step1;
+        }
+    }
+    else
     {
-        auto dx0 = xd[i - 1];
-        auto dx1 = xd[i - 0];
-        auto dy0 = yd[i - 1];
-        auto dy1 = yd[i - 0];
-        slopes[i] = 2 * (dx0 + dx1);
-        temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
+        foreach (i; 1 .. n - 1)
+        {
+            auto dx0 = xd[i - 1];
+            auto dx1 = xd[i - 0];
+            auto dy0 = yd[i - 1];
+            auto dy1 = yd[i - 0];
+            slopes[i] = 2 * (dx0 + dx1);
+            temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
+        }
     }
 
     foreach (i; 0 .. n - 1)
     {
-        auto c = i ==     0 ? first : xd[i - 1];
-        auto a = i == n - 2 ?  last : xd[i + 1];
+        auto c = i ==     0 ? first : kind == SplineKind.c2 ? xd[i - 1] : 0;
+        auto a = i == n - 2 ?  last : kind == SplineKind.c2 ? xd[i + 1] : 0;
         auto w = a / slopes[i];
         slopes[i + 1] -= w * c;
         temp[i + 1] -= w * temp[i];
@@ -940,7 +1157,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     foreach_reverse (i; 0 .. n - 1)
     {
-        auto c = i ==     0 ? first : xd[i - 1];
+        auto c = i ==     0 ? first : kind == SplineKind.c2 ? xd[i - 1] : 0;
         slopes[i] = (temp[i] - c * slopes[i + 1]) / slopes[i];
     }
 }
@@ -954,7 +1171,7 @@ struct SplineKernel(X)
     X wq = 0;
 
     ///
-    this()(X x0, X x1, X x)
+    this(X x0, X x1, X x)
     {
         step = x1 - x0;
         auto c0 = x - x0;
@@ -1009,3 +1226,22 @@ struct SplineKernel(X)
     ///
     alias withTwoDerivatives = opCall!2;
 }
+
+package T pchipTail(T)(in T step0, in T step1, in T diff0, in T diff1)
+{
+    import mir.math.common: copysign, fabs;
+    if (!diff0)
+    {
+        return 0;
+    }
+    auto slope = ((step0 * 2 + step1) * diff0 - step0 * diff1) / (step0 + step1);
+    if (copysign(1f, slope) != copysign(1f, diff0))
+    {
+        return 0;
+    }
+    if ((copysign(1f, diff0) != copysign(1f, diff1)) && (fabs(slope) > fabs(diff0 * 3)))
+    {
+        return diff0 * 3;
+    }
+    return slope;
+}
diff --git a/source/mir/interpolate/utility.d b/source/mir/interpolate/utility.d
index b0ed917f..5c3ffb5a 100644
--- a/source/mir/interpolate/utility.d
+++ b/source/mir/interpolate/utility.d
@@ -230,3 +230,23 @@ unittest
     assert(p.opCall!2(13)[2].approxEqual(s(13)));
     assert(p.opCall!3(13)[3].approxEqual(18));
 }
+
+
+package auto pchipTail(T)(in T step0, in T step1, in T diff0, in T diff1)
+{
+    import mir.math.common: copysign, fabs;
+    if (!diff0)
+    {
+        return 0;
+    }
+    auto slope = ((step0 * 2 + step1) * diff0 - step0 * diff1) / (step0 + step1);
+    if (copysign(1f, slope) != copysign(1f, diff0))
+    {
+        return 0;
+    }
+    if ((copysign(1f, diff0) != copysign(1f, diff1)) && (fabs(slope) > fabs(diff0 * 3)))
+    {
+        return diff0 * 3;
+    }
+    return slope;
+}

From e69cef36baa2d026ceb55e6164fa0d622d9a03d7 Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Wed, 23 Oct 2019 23:54:23 +0700
Subject: [PATCH 02/12] wip

---
 source/mir/interpolate/spline.d | 217 ++++++++++++++++++--------------
 1 file changed, 124 insertions(+), 93 deletions(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index 7a405bc9..bf438842 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -26,7 +26,7 @@ import std.meta;
 import mir.primitives;
 import mir.functional;
 import mir.internal.utility: Iota;
-import mir.math.common: fmamath;
+import mir.math.common;
 import mir.ndslice.internal;
 import mir.ndslice.slice;
 import mir.ndslice.traits;
@@ -913,6 +913,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     auto xd = points.diff;
     auto yd = values.diff;
+    auto dd = yd / xd;
 
     with(SplineKind) final switch(kind)
     {
@@ -950,6 +951,16 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             lbc.type = SplineBoundaryType.parabolic;
         if (rbc.type == SplineBoundaryType.notAKnot)
             rbc.type = SplineBoundaryType.parabolic;
+
+        if (n == 2)
+        {
+            if (lbc.type == SplineBoundaryType.monotone
+             || lbc.type == SplineBoundaryType.akima)
+                lbc.type = SplineBoundaryType.parabolic;
+            if (rbc.type == SplineBoundaryType.monotone
+             || rbc.type == SplineBoundaryType.akima)
+                rbc.type = SplineBoundaryType.parabolic;
+        }
         /// special case
         if (rbc.type == SplineBoundaryType.parabolic
          && lbc.type == SplineBoundaryType.parabolic)
@@ -979,26 +990,16 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     case notAKnot:
 
-        if (n > 2)
-        {
-            auto dx0 = xd[0];
-            auto dx1 = xd[1];
-            auto dy0 = yd[0];
-            auto dy1 = yd[1];
-            auto dd0 = dy0 / dx0;
-            auto dd1 = dy1 / dx1;
-
-            slopes.front = dx1;
-            first = dx0 + dx1;
-            temp.front = ((dx0 + 2 * first) * dx1 * dd0 + dx0 ^^ 2 * dd1) / first;
-            break;
-        }
-        else
-        {
-            slopes.front = 1;
-            first = 0;
-            temp.front = yd.front / xd.front;
-        }
+        auto dx0 = xd[0];
+        auto dx1 = xd[1];
+        auto dy0 = yd[0];
+        auto dy1 = yd[1];
+        auto dd0 = dy0 / dx0;
+        auto dd1 = dy1 / dx1;
+
+        slopes.front = dx1;
+        first = dx0 + dx1;
+        temp.front = ((dx0 + 2 * first) * dx1 * dd0 + dx0 ^^ 2 * dd1) / first;
         break;
     
     case firstDerivative:
@@ -1012,27 +1013,28 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
         slopes.front = 2;
         first = 1;
-        temp.front = 3 * (yd.front / xd.front) - 0.5 * lbc.value * xd.front;
+        temp.front = 3 * dd.front - 0.5 * lbc.value * xd.front;
         break;
 
     case parabolic:
 
         slopes.front = 2;
         first = 1;
-        temp.front = 2 * (yd.front / xd.front);
+        temp.front = 2 * dd.front;
         break;
     
     case monotone:
+
         slopes.front = 1;
         first = 0;
-        temp.front = pchipTail(xd[0], xd[1], yd[0]/xd[0], yd[1]/xd[1]);
+        temp.front = pchipTail(xd[0], xd[1], dd[0], dd[1]);
         break;
 
     case akima:
-        // TODO
-        slopes.front = 2;
-        first = 1;
-        temp.front = 2 * (yd.front / xd.front);
+
+        slopes.front = 1;
+        first = 0;
+        temp.back = akimaTail(dd[0], dd[1]);
         break;
     }
 
@@ -1042,24 +1044,16 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         assert(0);
 
     case notAKnot:
-        if (n > 2)
-        {
-            auto dx0 = xd[$ - 1];
-            auto dx1 = xd[$ - 2];
-            auto dy0 = yd[$ - 1];
-            auto dy1 = yd[$ - 2];
-            auto dd0 = dy0 / dx0;
-            auto dd1 = dy1 / dx1;
-            slopes.back = dx1;
-            last = dx0 + dx1;
-            temp.back = ((dx0 + 2 * last) * dx1 * dd0 + dx0 ^^ 2 * dd1) / last;
-        }
-        else
-        {
-            slopes.back = 1;
-            last = 0;
-            temp.back = yd.back / xd.back;
-        }
+
+        auto dx0 = xd[$ - 1];
+        auto dx1 = xd[$ - 2];
+        auto dy0 = yd[$ - 1];
+        auto dy1 = yd[$ - 2];
+        auto dd0 = dy0 / dx0;
+        auto dd1 = dy1 / dx1;
+        slopes.back = dx1;
+        last = dx0 + dx1;
+        temp.back = ((dx0 + 2 * last) * dx1 * dd0 + dx0 ^^ 2 * dd1) / last;
         break;
     
     case firstDerivative:
@@ -1073,75 +1067,87 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
         slopes.back = 2;
         last = 1;
-        temp.back = 3 * (yd.back / xd.back) + 0.5 * rbc.value * xd.back;
+        temp.back = 3 * dd.back + 0.5 * rbc.value * xd.back;
         break;
 
     case parabolic:
 
         slopes.back = 2;
         last = 1;
-        temp.back = 2 * (yd.back / xd.back);
+        temp.back = 2 * dd.back;
         break;
 
     case monotone:
+
         slopes.back = 1;
         last = 0;
-        temp.back = pchipTail(xd[$ - 1], xd[$ - 2], yd[$ - 1] / xd[$ - 1], yd[$ - 2] / xd[$ - 2]);
+        temp.back = pchipTail(xd[$ - 1], xd[$ - 2], dd[$ - 1], dd[$ - 2]);
         break;
 
     case akima:
-        // TODO
-        slopes.back = 2;
-        last = 1;
-        temp.back = 2 * (yd.back / xd.back);
+
+        slopes.back = 1;
+        last = 0;
+        temp.back = akimaTail(dd[$ - 1], dd[$ - 2]);
         break;
     }
 
-    if (kind == SplineKind.monotone)
+    with(SplineKind) final switch(kind)
     {
-        auto step0 = cast()(points[1] - points[0]);
-        auto step1 = cast()(points[2] - points[1]);
-        auto diff0 = cast()(values[1] - values[0]);
-        auto diff1 = cast()(values[2] - values[1]);
-        diff0 /= step0;
-        diff1 /= step1;
-
-        for(size_t i = 1;;)
-        {
-            import mir.math.common: copysign;
-            slopes[i] = 1;
-            if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1))
-            {
-                auto w0 = step1 * 2 + step0;
-                auto w1 = step0 * 2 + step1;
-                temp[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1);
-            }
-            else
+        case c2:
+            foreach (i; 1 .. n - 1)
             {
-                temp[i] = 0;
+                auto dx0 = xd[i - 1];
+                auto dx1 = xd[i - 0];
+                auto dy0 = yd[i - 1];
+                auto dy1 = yd[i - 0];
+                slopes[i] = 2 * (dx0 + dx1);
+                temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
             }
-            if (++i == n - 1)
+            break;
+        case cardinal:
+            // TODO
+            break;
+        case monotone:
             {
-                break;
+                auto step0 = cast()xd[0];
+                auto step1 = cast()xd[1];
+                auto diff0 = cast()yd[0];
+                auto diff1 = cast()yd[1];
+                diff0 /= step0;
+                diff1 /= step1;
+
+                for(size_t i = 1;;)
+                {
+                    slopes[i] = 1;
+                    if (diff0 && diff1 && copysign(1f, diff0) == copysign(1f, diff1))
+                    {
+                        auto w0 = step1 * 2 + step0;
+                        auto w1 = step0 * 2 + step1;
+                        temp[i] = (w0 + w1) / (w0 / diff0 + w1 / diff1);
+                    }
+                    else
+                    {
+                        temp[i] = 0;
+                    }
+                    if (++i == n - 1)
+                    {
+                        break;
+                    }
+                    step0 = step1;
+                    diff0 = diff1;
+                    step1 = xd[i];
+                    diff1 = yd[i];
+                    diff1 /= step1;
+                }
             }
-            step0 = step1;
-            diff0 = diff1;
-            step1 = points[i + 1] - points[i + 0];
-            diff1 = values[i + 1] - values[i + 0];
-            diff1 /= step1;
-        }
-    }
-    else
-    {
-        foreach (i; 1 .. n - 1)
-        {
-            auto dx0 = xd[i - 1];
-            auto dx1 = xd[i - 0];
-            auto dy0 = yd[i - 1];
-            auto dy1 = yd[i - 0];
-            slopes[i] = 2 * (dx0 + dx1);
-            temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
-        }
+            break;
+        case doubleQuadratic:
+            // TODO
+            break;
+        case akima:
+            // TODO
+            break;
     }
 
     foreach (i; 0 .. n - 1)
@@ -1162,6 +1168,31 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     }
 }
 
+private F akimaTail(F)(in F d2, in F d3)
+{
+    auto d1 = 2 * d2 - d3;
+    auto d0 = 2 * d1 - d2;
+    return akimaSlope(d0, d1, d2, d3);
+}
+
+private F akimaSlope(F)(in F d0, in F d1, in F d2, in F d3)
+{
+    if (d1 == d2)
+        return d1;
+    if (d0 == d1 && d2 == d3)
+        return (d1 + d2) * 0.5f;
+    if (d0 == d1)
+        return d1;
+    if (d2 == d3)
+        return d2;
+    auto w0 = fabs(d1 - d0);
+    auto w1 = fabs(d3 - d2);
+    auto ws = w0 + w1;
+    w0 /= ws;
+    w1 /= ws;
+    return w0 * d2 + w1 * d1;
+}
+
 ///
 struct SplineKernel(X)
 {

From 4c336f1c360a4be8254fd4f3dd4e7413890f9f20 Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Wed, 23 Oct 2019 23:58:53 +0700
Subject: [PATCH 03/12] wip

---
 source/mir/interpolate/spline.d | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index bf438842..dcc84466 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -1175,6 +1175,12 @@ private F akimaTail(F)(in F d2, in F d3)
     return akimaSlope(d0, d1, d2, d3);
 }
 
+private F akimaTail(F)(in F d1, in F d2, in F d3)
+{
+    auto d0 = 2 * d1 - d2;
+    return akimaSlope(d0, d1, d2, d3);
+}
+
 private F akimaSlope(F)(in F d0, in F d1, in F d2, in F d3)
 {
     if (d1 == d2)

From 8eb164a2b5db33f9213111d414b0145055a60e3d Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 00:17:11 +0700
Subject: [PATCH 04/12] wip

---
 source/mir/interpolate/spline.d | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index dcc84466..e069ca3c 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -1095,6 +1095,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     with(SplineKind) final switch(kind)
     {
         case c2:
+
             foreach (i; 1 .. n - 1)
             {
                 auto dx0 = xd[i - 1];
@@ -1105,9 +1106,16 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
                 temp[i] = 3 * (dy0 / dx0 * dx1 + dy1 / dx1 * dx0);
             }
             break;
+
         case cardinal:
-            // TODO
+
+            foreach (i; 1 .. n - 1)
+            {
+                slopes[i] = 1;
+                temp[i] = (1 - param) * ((values[i + 1] - values[i - 1]) / (points[i + 1] - points[i - 1]));
+            }
             break;
+
         case monotone:
             {
                 auto step0 = cast()xd[0];
@@ -1142,11 +1150,24 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
                 }
             }
             break;
+
         case doubleQuadratic:
             // TODO
             break;
         case akima:
-            // TODO
+
+            // if (n > 3)
+            // {
+            //     F d3 = dd[2];
+            //     F d2 = dd[1];
+            //     F d1 = dd[0];
+            //     F d0 = 2 * d1 - d2;
+            //     foreach (i; 1 .. n - 1)
+            //     {
+            //         slopes[i] = 1;
+            //         temp[i] = akimaSlope(dd[i - 1], dd[i - 0], dd[i + 1], dd[i + 2]);
+            //     }
+            // }
             break;
     }
 

From 35925b11368b84c67be2b54c77e9d25844d4553b Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 12:36:40 +0700
Subject: [PATCH 05/12] rename enum

---
 source/mir/interpolate/spline.d | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index e069ca3c..65327ca2 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -345,7 +345,7 @@ version(mir_test)
 
     auto x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
     auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;
-    auto interpolant = spline!double(x, y, SplineKind.monotone);
+    auto interpolant = spline!double(x, y, SplineType.monotone);
 
     auto xs = x[0 .. $ - 1] + 0.5;
 
@@ -389,11 +389,11 @@ version(mir_test)
         5.558593750000000,
         17.604662698412699,
         ];
-    auto interpolant = spline!double(points, values, SplineKind.monotone);
+    auto interpolant = spline!double(points, values, SplineType.monotone);
 
     auto pointsR = slice(-points.retro);
     auto valuesR = values.retro.slice;
-    auto interpolantR = spline!double(pointsR, valuesR, SplineKind.monotone);
+    auto interpolantR = spline!double(pointsR, valuesR, SplineType.monotone);
 
     version(X86_64)
     assert(map!interpolant(points[0 .. $ - 1] +  0.5) == map!interpolantR(pointsR.retro[0 .. $ - 1] - 0.5));
@@ -426,13 +426,13 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         in T valueOfBoundaryConditions = 0,
         )
     {
-        return spline(grid, values, SplineKind.c2, 0, typeOfBoundaries, valueOfBoundaryConditions);
+        return spline(grid, values, SplineType.c2, 0, typeOfBoundaries, valueOfBoundaryConditions);
     }
 
     Spline!(T, N, GridIterators) spline(yIterator, SliceKind ykind)(
         GridVectors grid,
         Slice!(yIterator, N, ykind) values,
-        SplineKind kind,
+        SplineType kind,
         in T param = 0,
         SplineBoundaryType typeOfBoundaries = SplineBoundaryType.notAKnot,
         in T valueOfBoundaryConditions = 0,
@@ -454,7 +454,7 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         GridVectors grid,
         Slice!(yIterator, N, ykind) values,
         SplineBoundaryCondition!T boundaries,
-        SplineKind kind = SplineKind.c2,
+        SplineType kind = SplineType.c2,
         in T param = 0,
         )
     {
@@ -476,7 +476,7 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         Slice!(yIterator, N, ykind) values,
         SplineBoundaryCondition!T rBoundary,
         SplineBoundaryCondition!T lBoundary,
-        SplineKind kind = SplineKind.c2,
+        SplineType kind = SplineType.c2,
         in T param = 0,
         )
     {
@@ -539,7 +539,7 @@ struct SplineBoundaryCondition(T)
 
 /++
 +/
-enum SplineKind
+enum SplineType
 {
     /++
     +/
@@ -632,7 +632,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
 
     $(RED For internal use.)
     +/
-    void _computeDerivatives(SplineKind kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc
+    void _computeDerivatives(SplineType kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc) scope @trusted nothrow @nogc
     {
         import mir.algorithm.iteration: maxLength;
         auto ml = this._data.maxLength;
@@ -643,7 +643,7 @@ struct Spline(F, size_t N = 1, FirstGridIterator = immutable(F)*, NextGridIterat
 
     /// ditto
     pragma(inline, false)
-    void _computeDerivativesTemp(SplineKind kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc
+    void _computeDerivativesTemp(SplineType kind, F param, SplineBoundaryCondition!F lbc, SplineBoundaryCondition!F rbc, Slice!(F*) temp) scope @system nothrow @nogc
     {
         import mir.algorithm.iteration: maxLength, each;
         import mir.ndslice.topology: map, byDim, evertPack;
@@ -893,7 +893,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     Slice!(IV, 1, vkind) values,
     Slice!(IS, 1, skind) slopes,
     Slice!(T*) temp,
-    SplineKind kind,
+    SplineType kind,
     F param,
     SplineBoundaryCondition!F lbc,
     SplineBoundaryCondition!F rbc,
@@ -915,7 +915,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     auto yd = values.diff;
     auto dd = yd / xd;
 
-    with(SplineKind) final switch(kind)
+    with(SplineType) final switch(kind)
     {
         case c2:
             break;
@@ -1092,7 +1092,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         break;
     }
 
-    with(SplineKind) final switch(kind)
+    with(SplineType) final switch(kind)
     {
         case c2:
 
@@ -1173,8 +1173,8 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     foreach (i; 0 .. n - 1)
     {
-        auto c = i ==     0 ? first : kind == SplineKind.c2 ? xd[i - 1] : 0;
-        auto a = i == n - 2 ?  last : kind == SplineKind.c2 ? xd[i + 1] : 0;
+        auto c = i ==     0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0;
+        auto a = i == n - 2 ?  last : kind == SplineType.c2 ? xd[i + 1] : 0;
         auto w = a / slopes[i];
         slopes[i + 1] -= w * c;
         temp[i + 1] -= w * temp[i];
@@ -1184,7 +1184,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     foreach_reverse (i; 0 .. n - 1)
     {
-        auto c = i ==     0 ? first : kind == SplineKind.c2 ? xd[i - 1] : 0;
+        auto c = i ==     0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0;
         slopes[i] = (temp[i] - c * slopes[i + 1]) / slopes[i];
     }
 }

From 6581234b01a3148ad5712eacb78f98d998f40853 Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 18:48:03 +0700
Subject: [PATCH 06/12] implement akima inner splines

---
 source/mir/interpolate/pchip.d  |  2 +-
 source/mir/interpolate/spline.d | 38 ++++++++++++++++-----------------
 2 files changed, 19 insertions(+), 21 deletions(-)

diff --git a/source/mir/interpolate/pchip.d b/source/mir/interpolate/pchip.d
index 2ebd377d..5e369dc0 100644
--- a/source/mir/interpolate/pchip.d
+++ b/source/mir/interpolate/pchip.d
@@ -51,7 +51,7 @@ template pchip(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridItera
         GridVectors grid,
         Slice!(yIterator, N, ykind) values) @safe
     {
-        return spline!T(grid, values, SplineKind.monotone);
+        return spline!T(grid, values, SplineType.monotone);
     }
 }
 
diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index 65327ca2..4ebfbcb7 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -1036,6 +1036,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         first = 0;
         temp.back = akimaTail(dd[0], dd[1]);
         break;
+
     }
 
     with(SplineBoundaryType) final switch(rbc.type)
@@ -1090,6 +1091,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
         last = 0;
         temp.back = akimaTail(dd[$ - 1], dd[$ - 2]);
         break;
+
     }
 
     with(SplineType) final switch(kind)
@@ -1155,20 +1157,22 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             // TODO
             break;
         case akima:
-
-            // if (n > 3)
-            // {
-            //     F d3 = dd[2];
-            //     F d2 = dd[1];
-            //     F d1 = dd[0];
-            //     F d0 = 2 * d1 - d2;
-            //     foreach (i; 1 .. n - 1)
-            //     {
-            //         slopes[i] = 1;
-            //         temp[i] = akimaSlope(dd[i - 1], dd[i - 0], dd[i + 1], dd[i + 2]);
-            //     }
-            // }
-            break;
+            {
+                auto d3 = dd[1];
+                auto d2 = dd[0];
+                auto d1 = 2 * d2 - d3;
+                auto d0 = d1;
+                foreach (i; 1 .. n - 1)
+                {
+                    d0 = d1;
+                    d1 = d2;
+                    d2 = d3;
+                    d3 = i == n - 2 ? 2 * d2 - d1 : dd[i + 1];
+                    slopes[i] = 1;
+                    temp[i] = akimaSlope(d0, d1, d2, d3);
+                }
+                break;
+            }
     }
 
     foreach (i; 0 .. n - 1)
@@ -1196,12 +1200,6 @@ private F akimaTail(F)(in F d2, in F d3)
     return akimaSlope(d0, d1, d2, d3);
 }
 
-private F akimaTail(F)(in F d1, in F d2, in F d3)
-{
-    auto d0 = 2 * d1 - d2;
-    return akimaSlope(d0, d1, d2, d3);
-}
-
 private F akimaSlope(F)(in F d0, in F d1, in F d2, in F d3)
 {
     if (d1 == d2)

From f18b4cf8126bca290f415baebdd8b448b0c693ac Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 19:42:29 +0700
Subject: [PATCH 07/12] add double quadratic inner slopes computation

---
 source/mir/interpolate/spline.d | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index 4ebfbcb7..697aa9ec 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -899,7 +899,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     SplineBoundaryCondition!F rbc,
     ) @trusted
 {
-    import mir.ndslice.topology: diff, ipack;
+    import mir.ndslice.topology: diff, zip, slide;
 
     assert (points.length >= 2);
     assert (points.length == values.length);
@@ -914,6 +914,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     auto xd = points.diff;
     auto yd = values.diff;
     auto dd = yd / xd;
+    auto dd2 = points.zip(values).slide!(3, "(c[1] - a[1]) / (c[0] - a[0])");
 
     with(SplineType) final switch(kind)
     {
@@ -1114,7 +1115,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             foreach (i; 1 .. n - 1)
             {
                 slopes[i] = 1;
-                temp[i] = (1 - param) * ((values[i + 1] - values[i - 1]) / (points[i + 1] - points[i - 1]));
+                temp[i] = (1 - param) * dd2[i - 1];
             }
             break;
 
@@ -1154,7 +1155,12 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             break;
 
         case doubleQuadratic:
-            // TODO
+            foreach (i; 1 .. n - 1)
+            {
+                slopes[i] = 1;
+                temp[i] = dd[i - 1] + dd[i + 1] - dd2[i - 1];
+            }
+
             break;
         case akima:
             {
@@ -1179,7 +1185,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     {
         auto c = i ==     0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0;
         auto a = i == n - 2 ?  last : kind == SplineType.c2 ? xd[i + 1] : 0;
-        auto w = a / slopes[i];
+        auto w = slopes[i] == 1 ? a : a / slopes[i];
         slopes[i + 1] -= w * c;
         temp[i + 1] -= w * temp[i];
     }
@@ -1189,7 +1195,8 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
     foreach_reverse (i; 0 .. n - 1)
     {
         auto c = i ==     0 ? first : kind == SplineType.c2 ? xd[i - 1] : 0;
-        slopes[i] = (temp[i] - c * slopes[i + 1]) / slopes[i];
+        auto v = temp[i] - c * slopes[i + 1];
+        slopes[i] = slopes[i]  == 1 ? v : v / slopes[i];
     }
 }
 

From f46382871a596041d14c94d93795b38f3e170009 Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 20:51:03 +0700
Subject: [PATCH 08/12] add test for akima spline

---
 source/mir/interpolate/spline.d | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index 697aa9ec..5577fb90 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -89,6 +89,12 @@ import mir.ndslice.traits;
             12.37135558,   4.99638066,  10.74362441,   0.16008641,
             11.94073593,  16.47908148,  17.49841853,   5.26600921,
             16.21796051,   8.96102894]));
+    
+    /// Akima spline
+    interpolant = spline!double(x, y, SplineType.akima);
+    assert(xs.vmap(interpolant).all!approxEqual(
+        [11.40871379,  2.64278898,  9.55774317,  4.84791141, 11.24842121,
+         16.16794267, 18.58060557,  5.2531411 , 17.45509005,  1.86992521]));
 }
 
 ///
@@ -1035,7 +1041,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
         slopes.front = 1;
         first = 0;
-        temp.back = akimaTail(dd[0], dd[1]);
+        temp.front = akimaTail(dd[0], dd[1]);
         break;
 
     }
@@ -1155,13 +1161,14 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             break;
 
         case doubleQuadratic:
+
             foreach (i; 1 .. n - 1)
             {
                 slopes[i] = 1;
                 temp[i] = dd[i - 1] + dd[i + 1] - dd2[i - 1];
             }
-
             break;
+
         case akima:
             {
                 auto d3 = dd[1];

From 51447bbe76bfeef670804be25886bbb9105890ce Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 22:43:27 +0700
Subject: [PATCH 09/12] update docs

---
 doc/Makefile                     |   2 +-
 source/mir/interpolate/package.d |   9 ++-
 source/mir/interpolate/pchip.d   |  72 +---------------------
 source/mir/interpolate/spline.d  | 101 +++++++++++++++++++++----------
 4 files changed, 77 insertions(+), 107 deletions(-)

diff --git a/doc/Makefile b/doc/Makefile
index e59cd183..441db156 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -36,7 +36,7 @@ PACKAGE_mir_cpp_export = numeric
 PACKAGE_mir_combinatorics = package
 PACKAGE_mir_container = binaryheap
 PACKAGE_mir_graph = tarjan package
-PACKAGE_mir_interpolate = package constant linear spline pchip utility polynomial
+PACKAGE_mir_interpolate = package constant linear spline utility polynomial
 PACKAGE_mir_math = sum numeric
 PACKAGE_mir_math_func = expdigamma
 PACKAGE_mir_ndslice_connect = cpython
diff --git a/source/mir/interpolate/package.d b/source/mir/interpolate/package.d
index dbd6ece8..8ebc8756 100644
--- a/source/mir/interpolate/package.d
+++ b/source/mir/interpolate/package.d
@@ -5,12 +5,11 @@ $(BOOKTABLE $(H2 Interpolation modules),
 $(TR $(TH Module) $(TH Interpolation kind))
 $(T2M constant, Constant Interpolant)
 $(T2M linear, Linear Interpolant)
-$(T2M pchip, Piecewise Cubic Hermite Interpolating Polynomial)
-$(T2M polynomial, Lagrange Barycentric Interpolation)
-$(T2M spline, Cubic Spline Interpolant)
+$(T2M polynomial, Lagrange Barycentric Interpolant)
+$(T2M spline, Piecewise Cubic Hermite Interpolant Spline: C2 (with contiguous second derivative), cardinal, monotone (aka PCHIP), double-quadratic, Akima)
 )
-
-Copyright: Copyright © 2017, Kaleidic Associates Advisory Limited
+]
+Copyright: Copyright © 2019, Symmetry Investments & Kaleidic Associates Advisory Limited
 Authors:   Ilya Yaroshenko
 
 Macros:
diff --git a/source/mir/interpolate/pchip.d b/source/mir/interpolate/pchip.d
index 5e369dc0..7717cee0 100644
--- a/source/mir/interpolate/pchip.d
+++ b/source/mir/interpolate/pchip.d
@@ -1,4 +1,4 @@
-/++
+/+
 $(H2 Piecewise Cubic Hermite Interpolating Polynomial)
 
 See_also: $(REF_ALTTEXT $(TT interp1), interp1, mir, interpolate)
@@ -11,7 +11,7 @@ Macros:
 SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, interpolate, $1)$(NBSP)
 T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
 +/
-deprecated
+deprecated("The PCHIP API was changed, the new one looks like spline!double(xSlice, ySlice, SplineType.monotone)")
 module mir.interpolate.pchip;
 
 public import mir.interpolate.spline: Spline;
@@ -29,7 +29,7 @@ private alias AliasSeq(T...) = T;
 /++
 Constructs piecewise cubic hermite interpolating polynomial with nodes on rectilinear grid.
 +/
-deprecated
+
 template pchip(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIterators = Repeat!(N - 1, FirstGridIterator))
     if (isFloatingPoint!T && is(T == Unqual!T) && N <= 6)
 {
@@ -54,69 +54,3 @@ template pchip(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridItera
         return spline!T(grid, values, SplineType.monotone);
     }
 }
-
-///
-version(mir_test)
-@safe unittest
-{
-    import mir.math.common: approxEqual;
-    import mir.algorithm.iteration: all;
-    import mir.ndslice.allocation: slice;
-    import mir.ndslice.slice: sliced;
-    import mir.ndslice.topology: vmap;
-
-    auto x = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
-    auto y = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;
-    auto interpolant = pchip!double(x, y);
-
-    auto xs = x[0 .. $ - 1] + 0.5;
-
-    () @trusted {
-        auto ys = xs.vmap(interpolant);
-
-        assert(ys.all!approxEqual([
-            5.333333333333334,
-            2.500000000000000,
-            10.000000000000000,
-            4.288971807628524,
-            11.202580845771145,
-            16.250000000000000,
-            17.962962962962962,
-            5.558593750000000,
-            17.604662698412699,
-            ]));
-    }();
-}
-
-// Check direction equality
-version(mir_test)
-@safe unittest
-{
-    import mir.math.common: approxEqual;
-    import mir.ndslice.slice: sliced;
-    import mir.ndslice.allocation: slice;
-    import mir.ndslice.topology: retro, map;
-
-    auto points = [1.0, 2, 4, 5, 8, 10, 12, 15, 19, 22].idup.sliced;
-    auto values = [17.0, 0, 16, 4, 10, 15, 19, 5, 18, 6].idup.sliced;
-
-    auto results = [
-        5.333333333333334,
-        2.500000000000000,
-        10.000000000000000,
-        4.288971807628524,
-        11.202580845771145,
-        16.250000000000000,
-        17.962962962962962,
-        5.558593750000000,
-        17.604662698412699,
-        ];
-    auto interpolant = pchip!double(points, values);
-
-    auto pointsR = slice(-points.retro);
-    auto valuesR = values.retro.slice;
-    auto interpolantR = pchip!double(pointsR, valuesR);
-
-    version(X86_64)
-    assert(map!interpolant(points[0 .. $ - 1] +  0.5) == map!interpolantR(pointsR.retro[0 .. $ - 1] - 0.5));
-}
diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index 5577fb90..fe27f8b7 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -95,6 +95,22 @@ import mir.ndslice.traits;
     assert(xs.vmap(interpolant).all!approxEqual(
         [11.40871379,  2.64278898,  9.55774317,  4.84791141, 11.24842121,
          16.16794267, 18.58060557,  5.2531411 , 17.45509005,  1.86992521]));
+
+    /// Double Quadratic spline
+    interpolant = spline!double(x, y, SplineType.doubleQuadratic);
+    import mir.interpolate.utility: ParabolaKernel;
+    auto kernel1 = ParabolaKernel!double(x[2], x[3], x[4],        y[2], y[3], y[4]);
+    auto kernel2 = ParabolaKernel!double(      x[3], x[4], x[5],        y[3], y[4], y[5]);
+    // weighted sum of quadratic functions
+    auto c = 0.35; // from [0 .. 1]
+    auto xp = c * x[3] + (1 - c) * x[4];
+    auto yp = c * kernel1(xp) + (1 - c) * kernel2(xp);
+    assert(interpolant(xp).approxEqual(yp));
+    // check parabolic extrapolation of the boundary intervals
+    kernel1 = ParabolaKernel!double(x[0], x[1], x[2], y[0], y[1], y[2]);
+    kernel2 = ParabolaKernel!double(x[$ - 3], x[$ - 2], x[$ - 1], y[$ - 3], y[$ - 2], y[$ - 1]);
+    assert(interpolant(x[0] - 23.421).approxEqual(kernel1(x[0] - 23.421)));
+    assert(interpolant(x[$ - 1] + 23.421).approxEqual(kernel2(x[$ - 1] + 23.421)));
 }
 
 ///
@@ -405,6 +421,38 @@ version(mir_test)
     assert(map!interpolant(points[0 .. $ - 1] +  0.5) == map!interpolantR(pointsR.retro[0 .. $ - 1] - 0.5));
 }
 
+/++
+Cubic Spline types.
+
+The first derivatives are guaranteed to be continuous for all cubic splines.
++/
+enum SplineType
+{
+    /++
+    Spline with contiguous second derivative.
+    +/
+    c2,
+    /++
+    $(HTTPS en.wikipedia.org/wiki/Cubic_Hermite_spline#Cardinal_spline, Cardinal) and Catmull–Rom splines.
+    +/
+    cardinal,
+    /++
+    The interpolant preserves monotonicity in the interpolation data and does not overshoot if the data is not smooth.
+    It is also known as   $(HTTPS docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.PchipInterpolator.html, PCHIP)
+    in numpy and Matlab.
+    +/
+    monotone,
+    /++
+    Weighted sum of two nearbor quadratic functions.
+    It is used in $(HTTPS s3-eu-west-1.amazonaws.com/og-public-downloads/smile-interpolation-extrapolation.pdf, financial analysis).
+    +/
+    doubleQuadratic,
+    /++
+    $(HTTPS en.wikipedia.org/wiki/Akima_spline, Akima spline).
+    +/
+    akima,
+}
+
 /++
 Constructs multivariate cubic spline in symmetrical form with nodes on rectilinear grid.
 Result has continues second derivatives throughout the curve / nd-surface.
@@ -452,6 +500,9 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         grid = immutable `x` values for interpolant
         values = `f(x)` values for interpolant
         boundaries = $(LREF SplineBoundaryCondition) for both tails.
+        kind = $(LREF SplineType) type of cubic spline.
+        param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types).
+            Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline.
     Constraints:
         `grid` and `values` must have the same length >= 3
     Returns: $(LREF Spline)
@@ -473,6 +524,9 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
         values = `f(x)` values for interpolant
         rBoundary = $(LREF SplineBoundaryCondition) for left tail.
         lBoundary = $(LREF SplineBoundaryCondition) for right tail.
+        kind = $(LREF SplineType) type of cubic spline.
+        param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types).
+            Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline.
     Constraints:
         `grid` and `values` must have the same length >= 3
     Returns: $(LREF Spline)
@@ -495,35 +549,41 @@ template spline(T, size_t N = 1, FirstGridIterator = immutable(T)*, NextGridIter
 }
 
 /++
-Cubic Spline  Boundary Condition Type
+Cubic Spline Boundary Condition Type.
 
-See_also: $(LREF SplineBoundaryCondition)
+See_also: $(LREF SplineBoundaryCondition) $(LREF SplineType)
 +/
 enum SplineBoundaryType
 {
     /++
-     not implemented
+    Not implemented.
     +/
     periodic = -1,
     /++
-     (default)
+    Not-a-knot (or cubic) boundary condition.
+    It is an aggresive boundary condition that is used only for C2 splines and is default for all API calls.
+    For other then C2 splines `notAKnot` is changed internally to
+    a default boundary type for used $(LREF SplineType).
     +/
     notAKnot,
     /++
-     set the first derivative
+    Set the first derivative.
     +/
     firstDerivative,
     /++
-     set the second derivative
+    Set the second derivative.
     +/
     secondDerivative,
     /++
+    Default for Cardinal and Double-Quadratic splines.
     +/
     parabolic,
     /++
+    Default for monotone (aka PHCIP ) splines.
     +/
     monotone,
     /++
+    Default for Akima splines.
     +/
     akima,
 }
@@ -541,29 +601,6 @@ struct SplineBoundaryCondition(T)
     T value = 0;
 }
 
-// private auto iter(alias s) = s.iterator;
-
-/++
-+/
-enum SplineType
-{
-    /++
-    +/
-    c2,
-    /++
-    +/
-    cardinal,
-    /++
-    +/
-    monotone,
-    /++
-    +/
-    doubleQuadratic,
-    /++
-    +/
-    akima,
-}
-
 /++
 Multivariate cubic spline with nodes on rectilinear grid.
 +/
@@ -1025,7 +1062,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     case parabolic:
 
-        slopes.front = 2;
+        slopes.front = 1;
         first = 1;
         temp.front = 2 * dd.front;
         break;
@@ -1080,7 +1117,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
 
     case parabolic:
 
-        slopes.back = 2;
+        slopes.back = 1;
         last = 1;
         temp.back = 2 * dd.back;
         break;
@@ -1165,7 +1202,7 @@ void splineSlopes(F, T, IP, IV, IS, SliceKind gkind, SliceKind vkind, SliceKind
             foreach (i; 1 .. n - 1)
             {
                 slopes[i] = 1;
-                temp[i] = dd[i - 1] + dd[i + 1] - dd2[i - 1];
+                temp[i] = dd[i - 1] + dd[i] - dd2[i - 1];
             }
             break;
 

From 379192285ce8ee2f71b52c6e0d615e2abf319c2d Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 22:55:39 +0700
Subject: [PATCH 10/12] fixup

---
 source/mir/interpolate/spline.d | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index fe27f8b7..fdfb3d38 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -562,7 +562,7 @@ enum SplineBoundaryType
     /++
     Not-a-knot (or cubic) boundary condition.
     It is an aggresive boundary condition that is used only for C2 splines and is default for all API calls.
-    For other then C2 splines `notAKnot` is changed internally to
+    For other then C2 splines, `notAKnot` is changed internally to
     a default boundary type for used $(LREF SplineType).
     +/
     notAKnot,

From d1298344824e4c6eb1062391ff305fbaccec7c4e Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 22:58:52 +0700
Subject: [PATCH 11/12] ditto

---
 source/mir/interpolate/spline.d | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/source/mir/interpolate/spline.d b/source/mir/interpolate/spline.d
index fdfb3d38..bfff19ae 100644
--- a/source/mir/interpolate/spline.d
+++ b/source/mir/interpolate/spline.d
@@ -925,6 +925,9 @@ Params:
     values = `f(x)` values for interpolant
     slopes = uninitialized ndslice to write slopes into
     temp = uninitialized temporary ndslice
+    kind = $(LREF SplineType) type of cubic spline.
+    param = tangent power parameter for cardinal $(LREF SplineType) (ignored by other spline types).
+        Use `1` for zero derivatives at knots and `0` for Catmull–Rom spline.
     lbc = left boundary condition
     rbc = right boundary condition
 Constraints:

From 376083001fb1e71bcdaa495bb3af2d91bb08d4be Mon Sep 17 00:00:00 2001
From: Ilya Yaroshenko <ilyayaroshenko@gmail.com>
Date: Thu, 24 Oct 2019 23:01:54 +0700
Subject: [PATCH 12/12] fix travis build

---
 test_examples.sh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test_examples.sh b/test_examples.sh
index 4f0cb3df..941fdb11 100755
--- a/test_examples.sh
+++ b/test_examples.sh
@@ -12,5 +12,5 @@ git clone https://github.com/libmir/mir-core
 # compile the examples
 for file in $(find out -name "*.d") ; do
     echo "Testing: $file"
-    $DMD -de -unittest -Isource -Imir-core/source -c -o- "$file"
+    $DMD -unittest -Isource -Imir-core/source -c -o- "$file"
 done