/* statistics/test_float_source.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, Brian Gough
 * 
 * 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.
 */

void FUNCTION (test, func) (const size_t stridea, const size_t strideb);

void
FUNCTION (test, func) (const size_t stridea, const size_t strideb)
{
  /* sample sets of doubles */
  size_t i;
  const size_t na = 14, nb = 14;

  const double rawa[] =
  {.0421, .0941, .1064, .0242, .1331,
   .0773, .0243, .0815, .1186, .0356,
   .0728, .0999, .0614, .0479};

  const double rawb[] =
  {.1081, .0986, .1566, .1961, .1125,
   .1942, .1079, .1021, .1583, .1673,
   .1675, .1856, .1688, .1512};

  const double raww[] = 
  {.0000, .0000, .0000, 3.000, .0000,
   1.000, 1.000, 1.000, 0.000, .5000,
   7.000, 5.000, 4.000, 0.123};

  BASE * sorted ;

  BASE * groupa = (BASE *) malloc (stridea * na * sizeof(BASE));
  BASE * groupb = (BASE *) malloc (strideb * nb * sizeof(BASE));
  BASE * w = (BASE *) malloc (strideb * na * sizeof(BASE));

#ifdef BASE_FLOAT
  double rel = 1e-6;
#else
  double rel = 1e-10;
#endif

  for (i = 0 ; i < na ; i++)
    groupa[i * stridea] = (BASE) rawa[i] ;

  for (i = 0 ; i < na ; i++)
    w[i * strideb] = (BASE) raww[i] ;

  for (i = 0 ; i < nb ; i++)
    groupb[i * strideb] = (BASE) rawb[i] ;


  {
    double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
    double expected = 0.0728;
    gsl_test_rel (mean, expected, rel, NAME(gsl_stats) "_mean");
  }

  {
    double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
    double var = FUNCTION(gsl_stats,variance_with_fixed_mean) (groupa, stridea, na, mean);
    double expected = 0.00113837428571429;
    gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance_with_fixed_mean");
  }


  {
    double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
    double var = FUNCTION(gsl_stats,sd_with_fixed_mean) (groupa, stridea, na, mean);
    double expected = 0.0337398026922845;
    gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_sd_with_fixed_mean");
  }


  {
    double var = FUNCTION(gsl_stats,variance) (groupb, strideb, nb);
    double expected = 0.00124956615384615;
    gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance");
  }

  {
    double sd = FUNCTION(gsl_stats,sd) (groupa, stridea, na);
    double expected = 0.0350134479659107;
    gsl_test_rel (sd, expected, rel, NAME(gsl_stats) "_sd");
  }

  {
    double ss = FUNCTION(gsl_stats,tss) (groupb, strideb, nb);
    double expected = 0.01624436;
    gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss");
  }

  {
    double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
    double ss = FUNCTION(gsl_stats,tss_m) (groupa, stridea, na, mean);
    double expected = 1.59372400000000e-02;
    gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss_m");
  }

  {
    double absdev = FUNCTION(gsl_stats,absdev) (groupa, stridea, na);
    double expected = 0.0287571428571429;
    gsl_test_rel (absdev, expected, rel, NAME(gsl_stats) "_absdev");
  }

  {
    double skew = FUNCTION(gsl_stats,skew) (groupa, stridea, na);
    double expected = 0.0954642051479004;
    gsl_test_rel (skew, expected, rel, NAME(gsl_stats) "_skew");
  }

  {
    double kurt = FUNCTION(gsl_stats,kurtosis) (groupa, stridea, na);
    double expected = -1.38583851548909 ;
    gsl_test_rel (kurt, expected, rel, NAME(gsl_stats) "_kurtosis");
  }

  {
    double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
    double expected = 0.0678111523670601;
    gsl_test_rel (wmean, expected, rel, NAME(gsl_stats) "_wmean");
  }

  {
    double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
    double wvar = FUNCTION(gsl_stats,wvariance_with_fixed_mean) (w, strideb, groupa, stridea, na, wmean);
    double expected = 0.000615793060878654;
    gsl_test_rel (wvar, expected, rel, NAME(gsl_stats) "_wvariance_with_fixed_mean");
  }

  {
    double est_wvar = FUNCTION(gsl_stats,wvariance) (w, strideb, groupa, stridea, na);
    double expected = 0.000769562962860317;
    gsl_test_rel (est_wvar, expected, rel, NAME(gsl_stats) "_wvariance");
  }

  {
    double wsd = FUNCTION(gsl_stats,wsd) (w, strideb, groupa, stridea, na);
    double expected = 0.0277409978706664;
    gsl_test_rel (wsd, expected, rel, NAME(gsl_stats) "_wsd");
  }


  {
    double wtss = FUNCTION(gsl_stats,wtss) (w, strideb, groupa, stridea, na);
    double expected =  1.39310864162578e-02;
    gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss");
  }

  {
    double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
    double wtss = FUNCTION(gsl_stats,wtss_m) (w, strideb, groupa, stridea, na, wmean);
    double expected =  1.39310864162578e-02;
    gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss_m");
  }



  {
    double wabsdev = FUNCTION(gsl_stats,wabsdev) (w, strideb, groupa, stridea, na);
    double expected = 0.0193205027504008;
    gsl_test_rel (wabsdev, expected, rel, NAME(gsl_stats) "_wabsdev");
  }

  {
    double wskew = FUNCTION(gsl_stats,wskew) (w, strideb, groupa, stridea, na);
    double expected = -0.373631000307076;
    gsl_test_rel (wskew, expected, rel, NAME(gsl_stats) "_wskew");
  }

  {
    double wkurt = FUNCTION(gsl_stats,wkurtosis) (w, strideb, groupa, stridea, na);
    double expected = -1.48114233353963;
    gsl_test_rel (wkurt, expected, rel, NAME(gsl_stats) "_wkurtosis");
  }

  {
    double c = FUNCTION(gsl_stats,covariance) (groupa, stridea, groupb, strideb, nb);
    double expected = -0.000139021538461539;
    gsl_test_rel (c, expected, rel, NAME(gsl_stats) "_covariance");
  }

  {
    double r = FUNCTION(gsl_stats,correlation) (groupa, stridea, groupb, strideb, nb);
    double expected = -0.112322712666074171;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_correlation");
  }

  {
    double *work = malloc(2 * na * sizeof(double));
    double r = FUNCTION(gsl_stats,spearman) (groupa, stridea, groupb, strideb, nb, work);
    double expected = -0.1604395604395604396;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_spearman");
    free(work);
  }

  {
    double pv = FUNCTION(gsl_stats,pvariance) (groupa, stridea, na, groupb, strideb, nb);
    double expected = 0.00123775384615385;
    gsl_test_rel (pv, expected, rel, NAME(gsl_stats) "_pvariance");
  }

  {
    double t = FUNCTION(gsl_stats,ttest) (groupa, stridea, na, groupb, strideb, nb);
    double expected = -5.67026326985851;
    gsl_test_rel (t, expected, rel, NAME(gsl_stats) "_ttest");
  }

  {
    BASE expected = (BASE)0.1331;
    gsl_test  (FUNCTION(gsl_stats,max) (groupa, stridea, na) != expected,
               NAME(gsl_stats) "_max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               FUNCTION(gsl_stats,max) (groupa, stridea, na), expected);
  }

  {
    BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
    BASE expected = (BASE)0.0242;
    gsl_test (min != expected,
              NAME(gsl_stats) "_min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
              min, expected);
  }

  {
    BASE min, max;
    BASE expected_max = (BASE)0.1331;
    BASE expected_min = (BASE)0.0242;
    
    FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
 
    gsl_test  (max != expected_max,
               NAME(gsl_stats) "_minmax max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               max, expected_max);
    gsl_test  (min != expected_min,
               NAME(gsl_stats) "_minmax min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               min, expected_min);
  }

  {
    int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
    int expected = 4;
    gsl_test (max_index != expected,
              NAME(gsl_stats) "_max_index (%d observed vs %d expected)",
              max_index, expected);
  }

  {
    int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
    int expected = 3;
    gsl_test (min_index != expected,
              NAME(gsl_stats) "_min_index (%d observed vs %d expected)",
              min_index, expected);
  }

  {
    size_t min_index, max_index;
    size_t expected_max_index = 4;
    size_t expected_min_index = 3;

    FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);

    gsl_test  (max_index != expected_max_index,
               NAME(gsl_stats) "_minmax_index max (%u observed vs %u expected)", 
               max_index, expected_max_index);
    gsl_test  (min_index != expected_min_index,
               NAME(gsl_stats) "_minmax_index min (%u observed vs %u expected)", 
               min_index, expected_min_index);
  }


  sorted = (BASE *) malloc(stridea * na * sizeof(BASE)) ;
  
  for (i = 0 ; i < na ; i++)
    sorted[stridea * i] = groupa[stridea * i] ;
  
  TYPE(gsl_sort)(sorted, stridea, na) ;

  {
    double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na) ;
    double expected = 0.07505;
    gsl_test_rel  (median,expected, rel,
                   NAME(gsl_stats) "_median_from_sorted_data (even)");
  }

  {
    double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na - 1) ;
    double expected = 0.0728;
    gsl_test_rel  (median,expected, rel,
                   NAME(gsl_stats) "_median_from_sorted_data (odd)");
  }

  {
    BASE * work = (BASE *) malloc (stridea * na * sizeof(BASE));
    double expected = 0.07505;
    double median;

    for (i = 0; i < na; i++)
      work[i * stridea] = (BASE) rawa[i];

    median = FUNCTION(gsl_stats,median)(work, stridea, na) ;

    gsl_test_rel (median, expected, rel,
                  NAME(gsl_stats) "_median (even)");

    free(work);
  }

  {
    BASE * work = (BASE *) malloc (stridea * (na - 1) * sizeof(BASE));
    double expected = 0.0728;
    double median;

    for (i = 0; i < na - 1; i++)
      work[i * stridea] = (BASE) sorted[i * stridea];

    median = FUNCTION(gsl_stats,median)(work, stridea, na - 1) ;

    gsl_test_rel (median, expected, rel,
                  NAME(gsl_stats) "_median (odd)");

    free(work);
  }

  {
    const double trim = 0.2;
    double trmean = FUNCTION(gsl_stats,trmean_from_sorted_data)(trim, sorted, stridea, na) ;
    double expected = 0.0719;
    gsl_test_rel  (trmean, expected, rel,
                   NAME(gsl_stats) "_trmean_from_sorted_data (even)");
  }

  {
    const double trim = 0.2;
    double trmean = FUNCTION(gsl_stats,trmean_from_sorted_data)(trim, sorted, stridea, na - 1) ;
    double expected = 0.06806666666666666;
    gsl_test_rel  (trmean, expected, rel,
                   NAME(gsl_stats) "_trmean_from_sorted_data (odd)");
  }

  {
    double gastwirth = FUNCTION(gsl_stats,gastwirth_from_sorted_data)(sorted, stridea, na) ;
    double expected = 0.07271;
    gsl_test_rel  (gastwirth, expected, rel,
                   NAME(gsl_stats) "_gastwirth_from_sorted_data (even)");
  }

  {
    double gastwirth = FUNCTION(gsl_stats,gastwirth_from_sorted_data)(sorted, stridea, na - 1) ;
    double expected = 0.06794;
    gsl_test_rel  (gastwirth, expected, rel,
                   NAME(gsl_stats) "_gastwirth_from_sorted_data (odd)");
  }

  {
    double zeroth = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.0) ;
    double expected = 0.0242;
    gsl_test_rel  (zeroth,expected, rel,
                   NAME(gsl_stats) "_quantile_from_sorted_data (0)");
  }

  {
    double top = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 1.0) ;
    double expected = 0.1331;
    gsl_test_rel  (top,expected, rel,
                   NAME(gsl_stats) "_quantile_from_sorted_data (100)");
  }

  {
    double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.5) ;
    double expected = 0.07505;
    gsl_test_rel  (median,expected, rel,
                   NAME(gsl_stats) "_quantile_from_sorted_data (50even)");
  }

  {
    double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na - 1, 0.5);
    double expected = 0.0728;
    gsl_test_rel  (median,expected, rel,
                   NAME(gsl_stats) "_quantile_from_sorted_data (50odd)");

  }

  {
    size_t k;
    BASE * work = (BASE *) malloc (stridea * na * sizeof(BASE));

    for (k = 0; k < na; ++k)
      {
        double expected = sorted[k * stridea];
        double kselect;

        /* copy rawa[] for each k, since gsl_stats_select() changes input array */
        for (i = 0; i < na; i++)
          work[i * stridea] = (BASE) rawa[i];

        kselect = FUNCTION(gsl_stats,select)(work, stridea, na, k);

        gsl_test_rel(kselect, expected, rel, NAME(gsl_stats) "_select");
      }

    free(work);
  }

  /* Test for IEEE handling - set third element to NaN */

  groupa [3*stridea] = GSL_NAN;

  {
    BASE max = FUNCTION(gsl_stats,max) (groupa, stridea, na);
    BASE expected = GSL_NAN;
    gsl_test  (!isnan(max),
               NAME(gsl_stats) "_max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               max, expected);
  }

  {
    BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
    BASE expected = GSL_NAN;
    gsl_test (!isnan(min),
              NAME(gsl_stats) "_min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
              min, expected);
  }

  {
    BASE min, max;
    BASE expected_max = GSL_NAN;
    BASE expected_min = GSL_NAN;
    
    FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
 
    gsl_test  (!isnan(max),
               NAME(gsl_stats) "_minmax max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               max, expected_max);
    gsl_test  (!isnan(min),
               NAME(gsl_stats) "_minmax min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               min, expected_min);
  }

#ifdef FAST
  {
    BASE min, max;
    BASE expected_max = GSL_NAN;
    BASE expected_min = GSL_NAN;
    
    FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na-1);
 
    gsl_test  (!isnan(max),
               NAME(gsl_stats) "_minmax(-1) max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               max, expected_max);
    gsl_test  (!isnan(min),
               NAME(gsl_stats) "_minmax(-1) min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)", 
               min, expected_min);
  }
#endif


  {
    int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
    int expected = 3;
    gsl_test (max_index != expected,
              NAME(gsl_stats) "_max_index NaN (%d observed vs %d expected)",
              max_index, expected);
  }

  {
    int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
    int expected = 3;
    gsl_test (min_index != expected,
              NAME(gsl_stats) "_min_index NaN (%d observed vs %d expected)",
              min_index, expected);
  }

  {
    size_t min_index, max_index;
    size_t expected_max_index = 3;
    size_t expected_min_index = 3;

    FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);

    gsl_test  (max_index != expected_max_index,
               NAME(gsl_stats) "_minmax_index max NaN (%u observed vs %u expected)", 
               max_index, expected_max_index);
    gsl_test  (min_index != expected_min_index,
               NAME(gsl_stats) "_minmax_index min NaN (%u observed vs %u expected)", 
               min_index, expected_min_index);
  }

  /* restore */
  groupa [3*stridea] = (BASE) rawa[3];

  {
    double * work = (double *) malloc (na * sizeof(double));
    double expected = 0.02925;
    double mad0 = FUNCTION(gsl_stats,mad0)(groupa, stridea, na, work);
    gsl_test_rel (mad0, expected, rel, NAME(gsl_stats) "_mad0 (even)");
    free(work);
  }

  {
    double * work = (double *) malloc ((na - 1) * sizeof(double));
    double expected = 0.02910;
    double mad0 = FUNCTION(gsl_stats,mad0)(groupa, stridea, na - 1, work);
    gsl_test_rel (mad0, expected, rel, NAME(gsl_stats) "_mad0 (odd)");
    free(work);
  }

  {
    BASE *work = malloc(na * sizeof(BASE));
    double r = FUNCTION(gsl_stats,Sn_from_sorted_data) (sorted, stridea, na, work);
    double expected = 0.04007136;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Sn_from_sorted_data (even)");
    free(work);
  }

  {
    BASE *work = malloc((na - 1) * sizeof(BASE));
    double r = FUNCTION(gsl_stats,Sn_from_sorted_data) (sorted, stridea, na - 1, work);
    double expected = 0.03728599834710744;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Sn_from_sorted_data (odd)");
    free(work);
  }

  {
    BASE *work = malloc(3 * na * sizeof(BASE));
    int *work_int = malloc(5 * na * sizeof(int));
    double r = FUNCTION(gsl_stats,Qn_from_sorted_data) (sorted, stridea, na, work, work_int);
    double expected = 0.04113672759664409;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Qn_from_sorted_data (even)");
    free(work);
    free(work_int);
  }

  {
    BASE *work = malloc(3 * (na - 1) * sizeof(BASE));
    int *work_int = malloc(5 * (na - 1) * sizeof(int));
    double r = FUNCTION(gsl_stats,Qn_from_sorted_data) (sorted, stridea, na - 1, work, work_int);
    double expected = 0.03684305546303433;
    gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Qn_from_sorted_data (odd)");
    free(work);
    free(work_int);
  }

  free (sorted);
  free (groupa);
  free (groupb);
  free (w);

}