testMath.cc

Go to the documentation of this file.

#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/BinomialCoefficients.h>
#include <utl/Polynomial.h>

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>

using namespace std;
using namespace utl;
using namespace tst;


#define EQUAL(x, y) CPPUNIT_ASSERT(Verify<Equal>(x, y))
#define CLOSE(x, y) CPPUNIT_ASSERT(Verify<CloseTo>(x, y))
#define CLOSEEPS(x, y, eps) CPPUNIT_ASSERT(Verify<CloseTo>(x, y, eps))


namespace std {
  template<typename T>
  ostream&
  operator<<(ostream& os, const vector<T>& v)
  {
    copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
    return os;
  }
}


class MathTest : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(MathTest);
  CPPUNIT_TEST(TestKolmogorovProbability);
  CPPUNIT_TEST(TestInverseNormalCDF);
  CPPUNIT_TEST(TestBinomialCoefficients);
  CPPUNIT_TEST(TestLambertW);
  CPPUNIT_TEST(TestChi2Probability);
  CPPUNIT_TEST(TestChi2Probability2);
  CPPUNIT_TEST(TestSign);
  CPPUNIT_TEST(TestPolynomial);
  CPPUNIT_TEST(TestEvalPoly);
  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() { }

  void tearDown() { }

  void
  TestKolmogorovProbability()
  {
    const double kolmDist = 0.75;
    for (int i = 1; i <= 10; ++i)
      cout << "KolmogorovProbability(" << i << ", " << kolmDist << ") = "
           << KolmogorovProbability(i, kolmDist) << endl;
    for (int i = 15; i <= 25; i += 5)
      cout << "KolmogorovProbability(" << i << ", " << kolmDist << ") = "
           << KolmogorovProbability(i, kolmDist) << endl;
    for (int i = 30; i <= 50; i += 10)
      cout << "KolmogorovProbability(" << i << ", " << kolmDist << ") = "
           << KolmogorovProbability(i, kolmDist) << endl;
  }

  void
  TestInverseNormalCDF()
  {
    const double deltap = 0.1;
    for (double p = deltap; p <= 0.5; p += deltap) {
      cout << "NormalDistribution (" << p << ") = "
           << InverseNormalCDF(p) << endl;
      const double sum = InverseNormalCDF(p) + InverseNormalCDF(1 - p);
      const double identity = 0.5*(1 + erf(InverseNormalCDF(p)/kSqrt2));
      cout << "InverseNormalDistribution (" << p << ") + InverseNormalDistribution (" << 1 - p << ") = "
           << sum << endl;
      cout << "erf(InverseNormalDistribution (" <<  p << ")) = "
           << identity << endl;
      CPPUNIT_ASSERT(Verify<CloseTo>(sum, 0.));
      CPPUNIT_ASSERT(Verify<CloseTo>(identity, p));
    }
  }

  void
  TestBinomialCoefficients()
  {
    BinomialCoefficients& binomial = BinomialCoefficients::GetInstance();

    EQUAL(binomial(0, 0), 1.);
    for (int n = 1; n <= 100; ++n) {
      EQUAL(binomial(n, -10), 0.);
      EQUAL(binomial(n, -1), 0.);
      EQUAL(binomial(n, 0), 1.);
      EQUAL(binomial(n, 1), double(n));
      double sum = binomial(n, 0) + binomial(n, n);
      double sum2 = 0*binomial(n, 0) + n*binomial(n, n);
      for (int k = 1; k < n; ++k) {
        sum += binomial(n, k);
        sum2 += k * binomial(n, k);
        // Pascal's rule
        EQUAL(binomial(n, k), binomial(n-1, k-1) + binomial(n-1, k));
        // symmetry
        EQUAL(binomial(n, k), binomial(n, n - k));
        // other relations
        CLOSE(binomial(n, k), double(n)/k * binomial(n-1, k-1));
        CLOSE((n-2*k)*binomial(n, k)/n, binomial(n-1, k) - binomial(n-1, k-1));
        // Vandermonde's relation
        for (int m = 0; m <= n; ++m) {
          double sum3 = 0;
          for (int j = 0; j <= n; ++j)
            sum3 += binomial(m, j) * binomial(n-m, k-j);
          CLOSE(sum3, binomial(n, k));
        }
      }
      CLOSE(sum, pow(2., n));
      CLOSE(sum2, n*pow(2., n-1));
      EQUAL(binomial(n, n-1), double(n));
      EQUAL(binomial(n, n), 1.);
      EQUAL(binomial(n, n+1), 0.);
      EQUAL(binomial(n, n+100), 0.);
    }
  }

  void
  TestLambertW()
  {
    EQUAL(LambertW<0>(0), 0.);
    // second order expansion around (0, 0)
    for (double p = 10; p >= 3; p -= 0.1) {
      const double x = pow(10., -p);
      CLOSE(LambertW<0>(x), 0.5*(-1 + sqrt(1 + 4*x)));
      CLOSE(LambertW<0>(-x), 0.5*(-1 + sqrt(1 - 4*x)));
    }
    // branch point
    const double bp = -1 / kE;
    CLOSE(LambertW<0>(bp), -1.);
    CLOSE(LambertW<-1>(bp), -1.);
    // expansion around branch point
    for (double p = 10; p >= 4.5; p -= 0.1) {
      const double eps = pow(10., -p);
      const double x = bp + eps;
      CLOSEEPS(LambertW<0>(x), -1 + sqrt(2*kE*eps), 1e-4);
      CLOSEEPS(LambertW<-1>(x), -1 - sqrt(2*kE*eps), 1e-4);
    }
    // overall
    for (double x = bp; x < 3; x += 0.01) {
      const double y = LambertW<0>(x);
      CLOSE(y*exp(y), x);
    }
    for (double t = 0; t < 20; t += 0.1) {
      const double x = bp * exp(-t);
      const double y = LambertW<-1>(x);
      CLOSE(y*exp(y), x);
    }
  }

  void
  TestChi2Probability()
  {
    // this crappy table is from http://en.wikipedia.org/wiki/Chi-squared_distribution
    const unsigned int nNdof = 10;
    const unsigned int nChi2 = 11;
    const double chi2[][nChi2] = {
      { },
      { 0.004, 0.02, 0.06, 0.15, 0.46,  1.07,  1.64,  2.71,  3.84,  6.64, 10.83 },  // 1
      { 0.10,  0.21, 0.45, 0.71, 1.39,  2.41,  3.22,  4.60,  5.99,  9.21, 13.82 },  // 2
      { 0.35,  0.58, 1.01, 1.42, 2.37,  3.66,  4.64,  6.25,  7.82, 11.34, 16.27 },  // 3
      { 0.71,  1.06, 1.65, 2.20, 3.36,  4.88,  5.99,  7.78,  9.49, 13.28, 18.47 },  // 4
      { 1.14,  1.61, 2.34, 3.00, 4.35,  6.06,  7.29,  9.24, 11.07, 15.09, 20.52 },  // 5
      { 1.63,  2.20, 3.07, 3.83, 5.35,  7.23,  8.56, 10.64, 12.59, 16.81, 22.46 },  // 6
      { 2.17,  2.83, 3.82, 4.67, 6.35,  8.38,  9.80, 12.02, 14.07, 18.48, 24.32 },  // 7
      { 2.73,  3.49, 4.59, 5.53, 7.34,  9.52, 11.03, 13.36, 15.51, 20.09, 26.12 },  // 8
      { 3.32,  4.17, 5.38, 6.39, 8.34, 10.66, 12.24, 14.68, 16.92, 21.67, 27.88 },  // 9
      { 3.94,  4.86, 6.18, 7.27, 9.34, 11.78, 13.44, 15.99, 18.31, 23.21, 29.59 }   // 10
    };
    const double probability[nChi2] = {
      0.95, 0.90, 0.80, 0.70, 0.50, 0.30, 0.20, 0.10, 0.05, 0.01, 0.001
    };

    for (unsigned int ndof = 1; ndof <= nNdof; ++ndof)
      for (unsigned int i = 0; i < nChi2; ++i)
        CLOSEEPS(Chi2Probability(chi2[ndof][i], ndof), probability[i], 0.02);
  }

  void
  TestChi2Probability2()
  {
    // this table was made with Mathematica function
    // 2 InverseGammaRegularized[ndof/2, 0, 1 - probability]
    const unsigned int nNdof = 10;
    const unsigned int nChi2 = 11;
    const double chi2[][nChi2] = {
      { },
      { 0.00393, 0.01579, 0.06418, 0.14847, 0.45494, 1.07419, 1.64237, 2.70554,
        3.84146, 6.6349, 10.82757 },
      { 0.10259, 0.21072, 0.44629, 0.71335, 1.38629, 2.40795, 3.21888, 4.60517,
        5.99146, 9.21034, 13.81551 },
      { 0.35185, 0.58437, 1.00517, 1.42365, 2.36597, 3.66487, 4.64163, 6.25139,
        7.81473, 11.34487, 16.26624},
      { 0.71072, 1.06362, 1.64878, 2.1947, 3.35669, 4.87843, 5.98862, 7.77944,
        9.48773, 13.2767, 18.46683 },
      { 1.14548, 1.61031, 2.34253, 2.99991, 4.35146, 6.06443, 7.28928, 9.23636,
        11.0705, 15.08627, 20.51501},
      { 1.63538, 2.20413, 3.07009, 3.82755, 5.34812, 7.23114, 8.55806, 10.64464,
        12.59159, 16.81189, 22.45774 },
      { 2.16735, 2.83311, 3.82232, 4.67133, 6.34581, 8.38343, 9.80325, 12.01704,
        14.06714, 18.47531, 24.32189},
      { 2.73264, 3.48954, 4.59357, 5.52742, 7.34412, 9.52446, 11.03009, 13.36157,
        15.50731, 20.09024, 26.12448},
      { 3.32511, 4.16816, 5.38005, 6.39331, 8.34283, 10.65637, 12.24215, 14.68366,
        16.91898, 21.66599, 27.87716 },
      { 3.9403, 4.86518, 6.17908, 7.26722, 9.34182, 11.78072, 13.44196, 15.98718,
        18.30704, 23.20925, 29.5883 }
    };
    const double probability[nChi2] = {
      0.95, 0.90, 0.80, 0.70, 0.50, 0.30, 0.20, 0.10, 0.05, 0.01, 0.001
    };

    for (unsigned int ndof = 1; ndof <= nNdof; ++ndof)
      for (unsigned int i = 0; i < nChi2; ++i)
        CLOSEEPS(Chi2Probability(chi2[ndof][i], ndof), probability[i], 1e-4);
  }

  void
  TestSign()
  {
    unsigned int ui = 0;
    EQUAL(Sign(ui), 0);
    ui = 1;
    EQUAL(Sign(ui), 1);
    int i = -1;
    EQUAL(Sign(i), -1);
    i = 0;
    EQUAL(Sign(i), 0);
    i = 1;
    EQUAL(Sign(i), 1);
    i = 100;
    EQUAL(Sign(i), 1);
    i = -100;
    EQUAL(Sign(i), -1);
    double d = -2;
    EQUAL(Sign(d), -1);
    d = 0;
    EQUAL(Sign(d), 0);
    d = 1e19;
    EQUAL(Sign(d), 1);
  }

  void
  TestPolynomial()
  {
    const double coeff[] = { 1, 2, 3, 4, 5 };
    const Polynomial p(coeff);
    CPPUNIT_ASSERT(p);
    for (unsigned int i = 0; i < Length(coeff); ++i)
      EQUAL(p.GetCoefficients()[i], coeff[i]);
    const double res[] = { 1, 15, 129, 547, 1593, 3711, 7465, 13539, 22737, 35983 };
    for (int i = 0; i < 10; ++i)
      CLOSE(p(i), res[i]);
    const int intCoeff[] = { 1, 2, 3, 4, 5 };
    const Polynomial q(intCoeff);
    for (int i = 0; i < 10; ++i)
      CLOSE(q(i), res[i]);
    cout << q << endl;
    stringstream ss;
    ss << q;
    Polynomial r;
    CPPUNIT_ASSERT(!r);
    cout << ss.str() << endl;
    ss >> r;
    CPPUNIT_ASSERT(r);
    EQUAL(q.GetCoefficients(), r.GetCoefficients());
  }

  void
  TestEvalPoly()
  {
    const double res[] = { 1, 15, 129, 547, 1593, 3711, 7465, 13539, 22737, 35983 };
    // test initializer_list<int>
    for (int i = 0; i < 10; ++i)
      CLOSE(EvalPoly({1, 2, 3, 4, 5}, i), res[i]);
    // test initializer_list<double>
    for (int i = 0; i < 10; ++i)
      CLOSE(EvalPoly({1., 2., 3., 4., 5.}, i), res[i]);
    {
      // test vector<double>
      const vector<double> coeff({1, 2, 3, 4, 5});
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(coeff, i), res[i]);
    }
    {
      // test vector<int>
      const vector<int> coeff({1, 2, 3, 4, 5});
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(coeff, i), res[i]);
    }
    {
      // test temporary vector<double>
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(vector<double>({1, 2, 3, 4, 5}), i), res[i]);
    }
    {
      // test temporary vector<int>
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(vector<int>({1, 2, 3, 4, 5}), i), res[i]);
    }
    {
      // test double c-array
      const double coeff[] = { 1, 2, 3, 4, 5 };
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(coeff, i), res[i]);
    }
    {
      // test int c-array
      const int coeff[] = { 1, 2, 3, 4, 5 };
      for (int i = 0; i < 10; ++i)
        CLOSE(EvalPoly(coeff, i), res[i]);
    }
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(MathTest);