/* filter/test_impulse.c

 * 

 * Copyright (C) 2018 Patrick Alken

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * 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 <gsl/gsl_math.h>

#include <gsl/gsl_vector.h>

#include <gsl/gsl_filter.h>

#include <gsl/gsl_movstat.h>

#include <gsl/gsl_test.h>

#include <gsl/gsl_rng.h>

#include <gsl/gsl_randist.h>

static int

vector_sum(const gsl_vector_int * v)

{

  size_t i;

  int sum = 0;

  for (i = 0; i < v->size; ++i)

    {

      int vi = gsl_vector_int_get(v, i);

      sum += vi;

    }

  return sum;

}

static void

test_impulse_proc(const double tol, const size_t n, const size_t K, const double nsigma,

                  const gsl_filter_end_t etype, const gsl_filter_scale_t stype,

                  const double outlier_percentage, gsl_rng * r)

{

  gsl_vector * x = gsl_vector_alloc(n);

  gsl_vector * y = gsl_vector_alloc(n);

  gsl_vector * z = gsl_vector_alloc(n);

  gsl_vector * y_med = gsl_vector_alloc(n);

  gsl_vector * xmedian = gsl_vector_alloc(n);

  gsl_vector * xsigma = gsl_vector_alloc(n);

  size_t noutlier;

  gsl_vector_int * ioutlier = gsl_vector_int_alloc(n);

  gsl_vector_int * ioutlier_exact = gsl_vector_int_alloc(n);

  size_t i;

  gsl_filter_impulse_workspace *impulse_p = gsl_filter_impulse_alloc(K);

  gsl_filter_median_workspace *median_p = gsl_filter_median_alloc(K);

  size_t noutlier_exact = 0;

  char buf[1024];

  gsl_vector_int_set_zero(ioutlier_exact);

  for (i = 0; i < n; ++i)

    {

      double xi = gsl_ran_gaussian(r, 1.0);

      double vi = gsl_rng_uniform(r);

      if (vi <= outlier_percentage)

        {

          xi += 15.0 * GSL_SIGN(xi);

          ++noutlier_exact;

          gsl_vector_int_set(ioutlier_exact, i, 1);

        }

      gsl_vector_set(x, i, xi);

    }

  /* first test that median filter is equal to impulse filter with nsigma = 0 */

  gsl_filter_median(etype, x, y_med, median_p);

  gsl_filter_impulse(etype, stype, 0.0, x, y, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  sprintf(buf, "impulse nsigma=0 smf comparison, etype=%u stype=%u", etype, stype);

  compare_vectors(tol, y, y_med, buf);

  /* second test: filter y = impulse(x) with given nsigma */

  gsl_filter_impulse(etype, stype, nsigma, x, y, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  /* test correct number of outliers detected */

  gsl_test(noutlier != noutlier_exact, "impulse [n=%zu,K=%zu,nsigma=%g,outlier_percentage=%g] noutlier=%zu exact=%zu",

           n, K, nsigma, outlier_percentage, noutlier, noutlier_exact);

#if 0

  {

    for (i = 0; i < n; ++i)

      {

        printf("%.12e %.12e %d %.12e %.12e\n",

               gsl_vector_get(x, i),

               gsl_vector_get(y, i),

               gsl_vector_int_get(ioutlier, i),

               gsl_vector_get(xmedian, i) + nsigma * gsl_vector_get(xsigma, i),

               gsl_vector_get(xmedian, i) - nsigma * gsl_vector_get(xsigma, i));

      }

  }

#endif

  /* test outliers found in correct locations */

  for (i = 0; i < n; ++i)

    {

      int val = gsl_vector_int_get(ioutlier, i);

      int val_exact = gsl_vector_int_get(ioutlier_exact, i);

      gsl_test(val != val_exact, "test_impulse: outlier mismatch [i=%zu,K=%zu,nsigma=%g,outlier_percentage=%g] ioutlier=%d ioutlier_exact=%d",

               i, K, nsigma, outlier_percentage, val, val_exact);

    }

  /* test noutlier = sum(ioutlier) */

  {

    size_t iout_sum = vector_sum(ioutlier);

    gsl_test(noutlier != iout_sum, "impulse [K=%zu,nsigma=%g,outlier_percentage=%g] noutlier=%zu sum(ioutlier)=%zu",

             K, nsigma, outlier_percentage, noutlier, iout_sum);

  }

  /* third test: test in-place filter z = impulse(z) */

  gsl_vector_memcpy(z, x);

  gsl_filter_impulse(etype, stype, nsigma, z, z, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  sprintf(buf, "impulse in-place nsigma=%g,n=%zu,K=%zu,etype=%u stype=%u", nsigma, n, K, etype, stype);

  compare_vectors(GSL_DBL_EPSILON, z, y, buf);

  gsl_vector_free(x);

  gsl_vector_free(y);

  gsl_vector_free(z);

  gsl_vector_free(y_med);

  gsl_vector_free(xmedian);

  gsl_vector_free(xsigma);

  gsl_vector_int_free(ioutlier);

  gsl_vector_int_free(ioutlier_exact);

  gsl_filter_impulse_free(impulse_p);

  gsl_filter_median_free(median_p);

}

static void

test_impulse(gsl_rng * r)

{

  const double tol = 1.0e-10;

  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_QN, 0.05, r);

  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_SN, 0.05, r);

  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_MAD, 0.05, r);

  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_IQR, 0.05, r);

}