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

assert(skewness([1.0, 2, 3]).approxEqual(0.0));

assert(skewness([1.0, 2, 4]).approxEqual((2.222222 / 3) / pow(4.666667 / 2, 1.5) * (3.0 ^^ 2) / (2.0 * 1.0)));
assert(skewness([1.0, 2, 4], true).approxEqual((2.222222 / 3) / pow(4.666667 / 3, 1.5)));

assert(skewness!float([0, 1, 2, 3, 4, 6].sliced(3, 2)).approxEqual(0.462910));

static assert(is(typeof(skewness!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.skewness.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.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.skewness.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.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]
].fuse;
auto result = [-1.090291, 1.732051, 1.732051, 1.229809];

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

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

import mir.math.common: approxEqual, pow, sqrt;
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.skewness!"online";
assert(!y.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));

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

// However, the two-pass and three-pass algorithms are numerically stable in this case
auto z1 = x.skewness!"twoPass";
assert(z1.approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));
assert(!z1.approxEqual(y));
auto z2 = x.skewness!"threePass";
assert(z2.approxEqual((12.000999 / 12) * sqrt(12.0 * 11.0) / 10.0));
assert(!z2.approxEqual(y));

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

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

//Set population skewness, skewness algorithm, sum algorithm or output type

auto a = [1.0, 1e98, 1, -1e98].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 squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.skewness(false).approxEqual(0.0));
assert(x.skewness(true).approxEqual(0.0));

assert(x.skewness!("online").approxEqual(0.0));
assert(x.skewness!("online", "kbn").approxEqual(0.0));
assert(x.skewness!("online", "kb2").approxEqual(0.0));
assert(x.skewness!("online", "precise").approxEqual(0.0));
assert(x.skewness!(double, "online", "precise").approxEqual(0.0));
assert(x.skewness!(double, "online", "precise")(true).approxEqual(0.0));

auto y = [uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.skewness!(ulong, "threePass");
assert(z == 0.0);
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.skewness;
assert(y.approxEqual(0.925493));
static assert(is(typeof(y) == double));

assert(x.skewness!float.approxEqual(0.925493));

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

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f)];
assert(foo.skewness == 0f);

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],
    [3.0, 4, 5],
    [6.0, 7, 7],
].fuse;

assert(x.skewness.approxEqual(-0.308571));

auto m0 = [0, 0.0, 0.0];
assert(x.alongDim!0.map!skewness.all!approxEqual(m0));
assert(x.alongDim!(-2).map!skewness.all!approxEqual(m0));

auto m1 = [0.935220, 0.0, -1.732051];
assert(x.alongDim!1.map!skewness.all!approxEqual(m1));
assert(x.alongDim!(-1).map!skewness.all!approxEqual(m1));
assert(iota(3, 4, 5, 6).as!double.alongDim!0.map!skewness.all!approxEqual(repeat(0.0, 4, 5, 6)));

assert(skewness(1.0, 2, 3) == 0.0);
assert(skewness!float(1, 2, 3) == 0f);