/*
 **********************************************************************
 * Copyright (C) Miroslav Lichvar  2017
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of version 2 of the GNU General Public License as
 * published by the Free Software Foundation.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 * 
 **********************************************************************
 */
#include <regress.c>
#include "test.h"

#define POINTS 64

void
test_unit(void)
{
  double x[POINTS], x2[POINTS], y[POINTS], w[POINTS];
  double b0, b1, b2, s2, sb0, sb1, slope, slope2, intercept, sd, median;
  double xrange, yrange, wrange, x2range;
  int i, j, n, m, c1, c2, c3, runs, best_start, dof;

  for (n = 3; n <= POINTS; n++) {
    for (i = 0; i < 200; i++) {
      slope = TST_GetRandomDouble(-0.1, 0.1);
      intercept = TST_GetRandomDouble(-1.0, 1.0);
      sd = TST_GetRandomDouble(1e-6, 1e-4);
      slope2 = (random() % 2 ? 1 : -1) * TST_GetRandomDouble(0.1, 0.5);

      DEBUG_LOG("iteration %d n=%d intercept=%e slope=%e sd=%e",
                i, n, intercept, slope, sd);

      for (j = 0; j < n; j++) {
        x[j] = -j;
        y[j] = intercept + slope * x[j] + (j % 2 ? 1 : -1) * TST_GetRandomDouble(1e-6, sd);
        w[j] = TST_GetRandomDouble(1.0, 2.0);
        x2[j] = (y[j] - intercept - slope * x[j]) / slope2;
      }

      RGR_WeightedRegression(x, y, w, n, &b0, &b1, &s2, &sb0, &sb1);
      DEBUG_LOG("WR b0=%e b1=%e s2=%e sb0=%e sb1=%e", b0, b1, s2, sb0, sb1);
      TEST_CHECK(fabs(b0 - intercept) < sd + 1e-3);
      TEST_CHECK(fabs(b1 - slope) < sd);

      if (RGR_FindBestRegression(x, y, w, n, 0, 3, &b0, &b1, &s2, &sb0, &sb1,
                                 &best_start, &runs, &dof)) {
        DEBUG_LOG("BR b0=%e b1=%e s2=%e sb0=%e sb1=%e runs=%d bs=%d dof=%d",
                  b0, b1, s2, sb0, sb1, runs, best_start, dof);

        TEST_CHECK(fabs(b0 - intercept) < sd + 1e-3);
        TEST_CHECK(fabs(b1 - slope) < sd);
      }

      if (RGR_MultipleRegress(x, x2, y, n, &b2)) {
        DEBUG_LOG("MR b2=%e", b2);
        TEST_CHECK(fabs(b2 - slope2) < 1e-6);
      }

      for (j = 0; j < n / 7; j++)
        y[random() % n] += 100 * sd;

      if (RGR_FindBestRobustRegression(x, y, n, 1e-8, &b0, &b1, &runs, &best_start)) {
        DEBUG_LOG("BRR b0=%e b1=%e runs=%d bs=%d", b0, b1, runs, best_start);

        TEST_CHECK(fabs(b0 - intercept) < sd + 1e-2);
        TEST_CHECK(fabs(b1 - slope) < 5.0 * sd);
      }

      for (j = 0; j < n; j++)
        x[j] = random() % 4 * TST_GetRandomDouble(-1000, 1000);

      median = RGR_FindMedian(x, n);

      for (j = c1 = c2 = c3 = 0; j < n; j++) {
        if (x[j] < median)
          c1++;
        if (x[j] > median)
          c3++;
        else
          c2++;
      }

      TEST_CHECK(c1 + c2 >= c3 && c1 <= c2 + c3);

      xrange = TST_GetRandomDouble(1e-6, pow(10.0, random() % 10));
      yrange = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10));
      wrange = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10));
      x2range = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10));
      m = random() % n;

      for (j = 0; j < n; j++) {
        x[j] = (j ? x[j - 1] : 0.0) + TST_GetRandomDouble(1e-6, xrange);
        y[j] = TST_GetRandomDouble(-yrange, yrange);
        w[j] = 1.0 + TST_GetRandomDouble(0.0, wrange);
        x2[j] = TST_GetRandomDouble(-x2range, x2range);
      }

      RGR_WeightedRegression(x, y, w, n, &b0, &b1, &s2, &sb0, &sb1);

      if (RGR_FindBestRegression(x + m, y + m, w, n - m, m, 3, &b0, &b1, &s2, &sb0, &sb1,
                                 &best_start, &runs, &dof))
        ;
      if (RGR_MultipleRegress(x, x2, y, n, &b2))
        ;
      if (RGR_FindBestRobustRegression(x, y, n, 1e-8, &b0, &b1, &runs, &best_start))
        ;
    }
  }
}