#include "csaps.h"

#define CATCH_CONFIG_MAIN

#include "catch.hpp"

TEST_CASE("Univariate auto smoothing", "[csaps]")

{

SECTION("21 points")

{

const size_t pcount = 21;

csaps::DoubleArray xdata(pcount);

xdata <<

0, 0.314159265358979, 0.628318530717959, 0.942477796076938, 1.25663706143592,

1.5707963267949, 1.88495559215388, 2.19911485751286, 2.51327412287183,

2.82743338823081, 3.14159265358979, 3.45575191894877, 3.76991118430775,

4.08407044966673, 4.39822971502571, 4.71238898038469, 5.02654824574367,

5.34070751110265, 5.65486677646163, 5.96902604182061, 6.28318530717959;

csaps::DoubleArray ydata(pcount);

ydata <<

-0.0308480549621108, 0.321227871478931, 0.581558026193185, 0.837552852746324,

0.979054097107034, 0.977259260528264, 0.928702941809463, 0.839204212128449,

0.633599187660844, 0.346610257849157, -0.0142182730042132, -0.308917481822601,

-0.569438958775259, -0.787746791676657, -0.964031440816114, -0.993880381393437,

-0.950748199764373, -0.857640149415879, -0.560502590131236, -0.270752875311336,

-0.013511401609863;

const size_t xidata_size = 120;

csaps::UnivariateCubicSmoothingSpline sp(xdata, ydata);

csaps::DoubleArray xidata;

csaps::DoubleArray yidata = sp(xidata_size, xidata);

csaps::DoubleArray desired_yidata(xidata_size);

desired_yidata <<

-0.0235609972734076, 0.0342554130011887, 0.0917604768962524, 0.148642848032251,

0.204591180029653, 0.259294126508924, 0.312440351240669, 0.363812477806949,

0.413516746584544, 0.461729017734914, 0.508625151419519, 0.554381007799819,

0.59917235322775, 0.643050468816528, 0.685696744725134, 0.726724183969348,

0.765745789564952, 0.802374564527724, 0.836223661222756, 0.866972935842014,

0.894473725242276, 0.918604542323229, 0.939243899984561, 0.956270311125961,

0.969562963467305, 0.979144106608908, 0.985355180101302, 0.988580811987222,

0.989205630309408, 0.987614263110598, 0.984190455154498, 0.979211074442117,

0.972745418402277, 0.964838935929978, 0.955537075920216, 0.944885287267991,

0.932927490881626, 0.919589386774692, 0.904596282338286, 0.88765407387356,

0.868468657681665, 0.846745930063754, 0.822194253635299, 0.794653654806466,

0.76415937272086, 0.73076244387651, 0.694513904771452, 0.655464791903716,

0.613667457460308, 0.569225401037061, 0.522308564522944, 0.473091330257213,

0.421748080579123, 0.368453197827928, 0.313386409949078, 0.256884901954999,

0.199465270307742, 0.141653871417018, 0.0839770616925367, 0.0269611975440111,

-0.0288790108903068, -0.0832967177157068, -0.136312941281018, -0.187960346206531,

-0.238271597112534, -0.287279358619317, -0.335016296758048, -0.381515103495269,

-0.426808493559216, -0.47092918245087, -0.513909885671217, -0.555783318721242,

-0.596577904130849, -0.636257828862702, -0.674737831212602, -0.711931377484918,

-0.74775193398402, -0.782112967014277, -0.814916147256374, -0.845917383266534,

-0.874774272419685, -0.901142570534573, -0.92467803342995, -0.945036416924564,

-0.961906593475871, -0.975319318974765, -0.985507042193372, -0.99270481875293,

-0.997147704274675, -0.999070754379844, -0.998693087710525, -0.996095112164348,

-0.991285814362011, -0.984273590665723, -0.975066837437697, -0.963673951040142,

-0.95006546257481, -0.933932055202886, -0.914838983410284, -0.892350910038224,

-0.86603249792793, -0.835448409920623, -0.800282146364173, -0.76096707154776,

-0.718228237424455, -0.672791349033563, -0.625382111414392, -0.576726229606248,

-0.527490657611685, -0.478024578050658, -0.428570228926154, -0.379369767649891,

-0.330665351633583, -0.282699138288946, -0.23565484652653, -0.189444552928148,

-0.143901464676976, -0.098858780436212, -0.0541496988690512, -0.00960741863869263;

REQUIRE(yidata.isApprox(desired_yidata));

}

SECTION("4 points")

{

const size_t pcount = 4;

csaps::DoubleArray xdata(pcount);

xdata << 1, 2, 4, 6;

csaps::DoubleArray ydata(pcount);

ydata << 2, 4, 5, 7;

csaps::UnivariateCubicSmoothingSpline sp(xdata, ydata);

const size_t xidata_size = 10;

csaps::DoubleArray xidata;

csaps::DoubleArray yidata = sp(xidata_size, xidata);

csaps::DoubleArray desired_yidata(xidata_size);

desired_yidata <<

2.2579392157892, 3.0231172855707, 3.6937304019483,

4.21971044584031, 4.65026761247821, 5.04804510368134,

5.47288175793241, 5.94265482897362, 6.44293945952166,

6.95847986982311;

REQUIRE(yidata.isApprox(desired_yidata));

}

SECTION("2 points")

{

const size_t pcount = 2;

csaps::DoubleArray xdata(pcount); xdata << 1., 2.;

csaps::DoubleArray ydata(pcount); ydata << 3., 4.;

csaps::DoubleArray xidata(pcount + 1); xidata << 1., 1.5, 2.;

csaps::UnivariateCubicSmoothingSpline sp(xdata, ydata);

csaps::DoubleArray yidata = sp(xidata);

csaps::DoubleArray desired_yidata(pcount + 1);

desired_yidata << 3., 3.5, 4.;

REQUIRE(yidata.isApprox(desired_yidata));

}

}

TEST_CASE("Zero based indices cause crash", "[csaps]") {

const size_t pcount = 83;

csaps::DoubleArray xdata(pcount);

csaps::DoubleArray ydata(pcount);

csaps::DoubleArray desired_ydata(pcount);

xdata <<

14700, 14716.7, 14733.3, 14750, 14766.7, 14783.3, 14733.3, 14750, 14766.7,

14783.3, 14800, 14816.7, 14766.7, 14783.3, 14800, 14816.7, 14833.3, 14850,

14866.7, 14883.3, 14900, 14866.7, 14883.3, 14900, 14933.3, 14950, 14966.7,

14983.3, 15000, 15016.7, 15033.3, 15050, 15066.7, 15083.3, 15100, 15116.7,

15133.3, 15150, 15166.7, 15183.3, 15200, 15216.7, 15233.3, 15250, 15266.7,

15283.3, 15300, 15316.7, 15333.3, 15350, 15366.7, 15383.3, 15400, 15433.3,

15450, 15466.7, 15483.3, 15500, 15516.7, 15533.3, 15550, 15566.7, 15583.3,

15600, 15616.7, 15633.3, 15650, 15666.7, 15683.3, 15700, 15716.7, 15733.3,

15750, 15766.7, 15833.3, 15850, 15866.7, 15883.3, 15900, 15916.7, 15933.3,

15950, 15966.7;

ydata <<

-0.0414725, -0.0326667, -0.0326667, -0.0326667, -0.0326667, -0.0282625,

-0.0305985, -0.0227109, -0.0227109, -0.0227109, -0.0227109, -0.0188219,

-0.0320486, -0.0258779, -0.0251352, -0.0235144, -0.0235144, -0.0110714,

0.00259269, 0.00259269, 0.00536645, -0.00531659, 0.00446568, 0.00446568,

0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998,

0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998,

0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998,

0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998,

0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.00853998, 0.0115007,

0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007,

0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007,

0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0115007, 0.0154065,

0.0148546, 0.0148546, 0.0162699, 0.0171716, 0.0171716, 0.0171716, 0.0180766,

0.0179968;

desired_ydata <<

-0.04141771374345, -0.03278908150104, -0.03258429432982, -0.03268043844157,

-0.03280513020903, -0.02274215961948, -0.03749490144114, -0.06895467068492,

-0.0320388447726, -0.02590728278666, -0.02504346191236, -0.02367221083675,

-0.1406232270049, -0.264993765197, -0.4569275259756, -0.02367221083675,

-0.02333221857753, -0.0110606663237, 0.00247470897021, 0.002266711647101,

0.004463335796486, -2.159191881555, -2.881886011007, 0.004463335796486,

0.008514325430511, 0.008557881993718, 0.008536063171583, 0.008540762037834,

0.008539828158979, 0.008540009712038, 0.008539974167637, 0.008539981129271,

0.008539979779126, 0.008539980043353, 0.008539979991606, 0.008539980001642,

0.008539979999678, 0.008539980000061, 0.008539979999997, 0.008539979999954,

0.008539980000247, 0.00853997999872, 0.008539980006523, 0.008539979966649,

0.008539980172252, 0.008539979122427, 0.008539984486767, 0.008539956828557,

0.008540097999324, 0.008539378355551, 0.008543035208958, 0.008526022436719,

0.008560702288716, 0.01147998229411, 0.0115146363266, 0.01149765547087,

0.01150130774154, 0.01150058200463, 0.01150072308903, 0.01150069546772,

0.01150070087752, 0.01150069982849, 0.01150070003302, 0.01150069999691,

0.01150069998353, 0.01150070009016, 0.01150069953778, 0.01150070238749,

0.01150068783734, 0.01150076215467, 0.01150037993163, 0.0115023014714,

0.01149343128828, 0.01150932664402, 0.01539522395144, 0.01485708306995,

0.01487359427, 0.01625739938059, 0.01716728084523, 0.01716788462739,

0.01719369011052, 0.01805479951794, 0.01800415289884;

csaps::UnivariateCubicSmoothingSpline sp(xdata, ydata, 0.05);

csaps::DoubleArray output_ydata = sp(xdata);

REQUIRE(output_ydata.isApprox(desired_ydata));

}

TEST_CASE("Diff 4 elements", "[diff]")

{

csaps::DoubleArray arr(4);

arr << 1, 2, 4, 6;

csaps::DoubleArray d = csaps::Diff(arr);

csaps::DoubleArray dd(3);

dd << 1, 2, 2;

REQUIRE(d.isApprox(dd));

}

TEST_CASE("Make diagonal sparse matrix", "[spdiag][hide]")

{

csaps::DoubleArray2D diags(3, 3);

diags <<

1, 2, 3,

4, 5, 6,

7, 8, 9;

csaps::IndexArray offsets(3);

SECTION("Make 3x3 matrix with offsets -1 0 1")

{

offsets << -1, 0, 1;

csaps::DoubleSparseMatrix m = csaps::MakeSparseDiagMatrix(diags, offsets, 3, 3);

csaps::DoubleSparseMatrix dm(3, 3);

dm.coeffRef(0, 0) = 4;

dm.coeffRef(1, 0) = 1;

dm.coeffRef(0, 1) = 8;

dm.coeffRef(1, 1) = 5;

dm.coeffRef(2, 1) = 2;

dm.coeffRef(1, 2) = 9;

dm.coeffRef(2, 2) = 6;

/*

REQUIRE(m.isApprox(dm));

}

SECTION("Make 3x5 matrix with offsets -1 0 1")

{

offsets << -1, 0, 1;

csaps::DoubleSparseMatrix m = csaps::MakeSparseDiagMatrix(diags, offsets, 3, 5);

csaps::DoubleSparseMatrix dm(3, 5);

dm.coeffRef(0, 0) = 4;

dm.coeffRef(1, 0) = 2;

dm.coeffRef(0, 1) = 7;

dm.coeffRef(1, 1) = 5;

dm.coeffRef(2, 1) = 3;

dm.coeffRef(1, 2) = 8;

dm.coeffRef(2, 2) = 6;

dm.coeffRef(2, 3) = 9;

/*

REQUIRE(m.isApprox(dm));

}

SECTION("Make 5x3 matrix with offsets -1 0 1")

{

offsets << -1, 0, 1;

csaps::DoubleSparseMatrix m = csaps::MakeSparseDiagMatrix(diags, offsets, 5, 3);

csaps::DoubleSparseMatrix dm(5, 3);

dm.coeffRef(0, 0) = 4;

dm.coeffRef(1, 0) = 1;

dm.coeffRef(0, 1) = 8;

dm.coeffRef(1, 1) = 5;

dm.coeffRef(2, 1) = 2;

dm.coeffRef(1, 2) = 9;

dm.coeffRef(2, 2) = 6;

dm.coeffRef(3, 2) = 3;

/*

REQUIRE(m.isApprox(dm));

}

SECTION("Make 5x3 matrix with offsets -2 -1 0")

{

offsets << -2, -1, 0;

csaps::DoubleSparseMatrix m = csaps::MakeSparseDiagMatrix(diags, offsets, 5, 3);

csaps::DoubleSparseMatrix dm(5, 3);

dm.coeffRef(0, 0) = 7;

dm.coeffRef(1, 0) = 4;

dm.coeffRef(2, 0) = 1;

dm.coeffRef(1, 1) = 8;

dm.coeffRef(2, 1) = 5;

dm.coeffRef(3, 1) = 2;

dm.coeffRef(2, 2) = 9;

dm.coeffRef(3, 2) = 6;

dm.coeffRef(4, 2) = 3;

/*

REQUIRE(m.isApprox(dm));

}

SECTION("Make 3x3 matrix with offset 0")

{

csaps::DoubleArray diag(3);

diag << 1, 2, 3;

csaps::IndexArray offset(1);

offset << 0;

csaps::DoubleSparseMatrix m = csaps::MakeSparseDiagMatrix(diag.transpose(), offset, 3, 3);

std::cout << m << std::endl;

csaps::DoubleSparseMatrix dm(3, 3);

dm.coeffRef(0, 0) = 1;

dm.coeffRef(1, 1) = 2;

dm.coeffRef(2, 2) = 3;

/*

REQUIRE(m.isApprox(dm));

}

}