import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

assert(kurtosis([1.0, 2, 3, 4]).approxEqual(-1.2));

assert(kurtosis([1.0, 2, 4, 5]).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) * (3.0 * 5.0) / (2.0 * 1.0) - 3.0 * (3.0 * 3.0) / (2.0 * 1.0)));
// population excess kurtosis
assert(kurtosis([1.0, 2, 4, 5], true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) - 3.0));
// sample raw kurtosis
assert(kurtosis([1.0, 2, 4, 5], false, true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) * (3.0 * 5.0) / (2.0 * 1.0) - 3.0 * (3.0 * 3.0) / (2.0 * 1.0) + 3.0));
// population raw kurtosis
assert(kurtosis([1.0, 2, 4, 5], true, true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0)));

assert(kurtosis!float([0, 1, 2, 3, 4, 6].sliced(3, 2)).approxEqual(-0.2999999));

static assert(is(typeof(kurtosis!float([1, 2, 3])) == float));

import mir.math.common: approxEqual, pow;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0];

assert(x.kurtosis.approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));

import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.kurtosis.approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: alongDim, byDim, map;

auto x = [
    [0.0,  1.0,  1.5, 2.0], 
    [3.5, 4.25,  2.0, 7.5],
    [5.0,  1.0,  1.5, 0.0],
    [1.5,  4.5, 4.75, 0.5]
].fuse;
auto result = [-2.067182, -5.918089, 3.504056, 2.690240];

// Use byDim or alongDim with map to compute kurtosis of row/column.
assert(x.byDim!1.map!kurtosis.all!approxEqual(result));
assert(x.alongDim!0.map!kurtosis.all!approxEqual(result));

// FIXME
// Without using map, computes the kurtosis of the whole slice
// assert(x.byDim!1.kurtosis == x.sliced.kurtosis);
// assert(x.alongDim!0.kurtosis == x.sliced.kurtosis);

import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto x = a + 100_000_000_000;

// The online algorithm is numerically unstable in this case
auto y = x.kurtosis!"online";
assert(!y.approxEqual((792.78411865 / 12) / pow(54.76562500 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));

// The naive algorithm has an assert error in this case because standard
// deviation is calculated naively as zero. The kurtosis formula would then
// be dividing by zero. 
//auto z0 = x.kurtosis!(real, "naive");

// The two-pass algorithm is also numerically unstable in this case
auto z1 = x.kurtosis!"twoPass";
assert(!z1.approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(!z1.approxEqual(y));

// However, the three-pass algorithm is numerically stable in this case
auto z2 = x.kurtosis!"threePass";
assert(z2.approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(!z2.approxEqual(y));

// And the assumeZeroMean algorithm provides the incorrect answer, as expected
auto z3 = x.kurtosis!"assumeZeroMean";
assert(!z3.approxEqual(y));

import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

// Set population/sample kurtosis, excess/raw kurtosis, kurtosis algorithm,
// sum algorithm or output type

auto a = [1.0, 1e72, 1, -1e72].sliced;
auto x = a * 10_000;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and taking `x` to the fourth power, the first and third numbers in the slice
have precision too low to be included in the centered sum of cubes. 
+/
assert(x.kurtosis.approxEqual(1.5));
assert(x.kurtosis(false).approxEqual(1.5));
assert(x.kurtosis(true).approxEqual(-1.0));
assert(x.kurtosis(true, true).approxEqual(2.0));
assert(x.kurtosis(false, true).approxEqual(4.5));

assert(x.kurtosis!("online").approxEqual(1.5));
assert(x.kurtosis!("online", "kbn").approxEqual(1.5));
assert(x.kurtosis!("online", "kb2").approxEqual(1.5));
assert(x.kurtosis!("online", "precise").approxEqual(1.5));
assert(x.kurtosis!(double, "online", "precise").approxEqual(1.5));
assert(x.kurtosis!(double, "online", "precise")(true).approxEqual(-1.0));
assert(x.kurtosis!(double, "online", "precise")(true, true).approxEqual(2.0));

auto y = [uint.max - 3, uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.kurtosis!(ulong, "threePass");
assert(z.approxEqual(-1.2));
static assert(is(typeof(z) == double));

import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.kurtosis;
assert(y.approxEqual(0.223394));
static assert(is(typeof(y) == double));

assert(x.kurtosis!float.approxEqual(0.223394));

import mir.math.common: approxEqual;

static struct Foo {
    float x;
    alias x this;
}

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f), Foo(4f)];
assert(foo.kurtosis.approxEqual(-1.2f));

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: as, iota, alongDim, map, repeat;

auto x = [
    [0.0,  1,  3,  5],
    [3.0,  4,  5,  7],
    [6.0,  7, 10, 11],
    [9.0, 12, 15, 12]
].fuse;

assert(x.kurtosis.approxEqual(-0.770040));

auto m0 = [-1.200000, -0.152893, -1.713859, -3.869005];
assert(x.alongDim!0.map!kurtosis.all!approxEqual(m0));
assert(x.alongDim!(-2).map!kurtosis.all!approxEqual(m0));

auto m1 = [-1.699512, 0.342857, -4.339100, 1.500000];
assert(x.alongDim!1.map!kurtosis.all!approxEqual(m1));
assert(x.alongDim!(-1).map!kurtosis.all!approxEqual(m1));

assert(iota(4, 5, 6, 7).as!double.alongDim!0.map!kurtosis.all!approxEqual(repeat(-1.2, 5, 6, 7)));

import mir.math.common: approxEqual;

assert(kurtosis(1.0, 2, 3, 4).approxEqual(-1.2));
assert(kurtosis!float(1, 2, 3, 4).approxEqual(-1.2f));