#include <config.h>

#include <stdlib.h>

#include <string.h>

#include <gsl/gsl_ntuple.h>

#include <gsl/gsl_test.h>

#include <gsl/gsl_ieee_utils.h>

struct data

{

  int num;

  double x;

  double y;

  double z;

};

int sel_func (void *ntuple_data, void * params);

double val_func (void *ntuple_data, void * params);

int

main (void)

{

  struct data ntuple_row;

  int i;

  double x[1000], y[1000], z[1000], f[100];

  gsl_ntuple_select_fn S;

  gsl_ntuple_value_fn V;

  double scale = 1.5;

  gsl_ieee_env_setup ();

  /* zero struct including padding bytes to avoid valgrind errors */

  memset(&ntuple_row, 0, sizeof(struct data));

  S.function = &sel_func;

  S.params = &scal;;

  V.function = &val_func;

  V.params = &scal;;

  {

    gsl_ntuple *ntuple = gsl_ntuple_create ("test.dat", &ntuple_row, 

                                            sizeof (ntuple_row));

    int status = 0;

    for (i = 0; i < 100; i++) f[i] = 0;

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

      {

        double xi = 1.0 / (i + 1.5);

        double yi = xi * xi ;

        double zi = xi * xi * xi;

        ntuple_row.x = xi;

        ntuple_row.y = yi;

        ntuple_row.z = zi;

        ntuple_row.num = i;

        x[i] = xi; y[i] = yi; z[i] = zi;

        if (xi * scale < 0.1)

          {

            double v = xi + yi + zi;

            int k = (int)(100.0*v*scale);

            f[k]++;

          }

        /* printf ("x,y,z = %f,%f,%f; n=%x \n", ntuple_row.x,

           ntuple_row.y, ntuple_row.z, ntuple_row.num); */

        {

          int s = gsl_ntuple_bookdata (ntuple);

          if (s != GSL_SUCCESS)

            {

              status = 1;

            }

        }

      }

    gsl_ntuple_close (ntuple);

    gsl_test (status, "writing ntuples");

  }

  {

    gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row, 

                                          sizeof (ntuple_row));

    int status = 0;

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

      {

        gsl_ntuple_read (ntuple);

        status = (ntuple_row.num != i);

        status |= (ntuple_row.x != x[i]);

        status |= (ntuple_row.y != y[i]);

        status |= (ntuple_row.z != z[i]);

        /* printf ("x,y,z = %f,%f,%f; n=%d\n", ntuple_row.x,

                ntuple_row.y, ntuple_row.z, ntuple_row.num); */

      }

    gsl_ntuple_close (ntuple);

    gsl_test (status, "reading ntuples");

  }    

  {

    int status = 0;

    gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row, 

                                          sizeof (ntuple_row));

    gsl_histogram *h = gsl_histogram_calloc_uniform (100, 0., 1.);

    gsl_ntuple_project (h, ntuple, &V, &S);

    gsl_ntuple_close (ntuple);

    /* gsl_histogram_fprintf (stdout, h, "%f", "%f"); */

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

      {

        /* printf ("h  %g f  %g\n", h->bin[i], f[i]); */

        if (h->bin[i] != f[i])

          {

            status = 1;

          }

      }

    gsl_test (status, "histogramming ntuples");

    gsl_histogram_free (h);

  }

  exit (gsl_test_summary());

}

int

sel_func (void *ntuple_data, void * params)

{

  double x, y, z, scale;

  scale = *(double *)params;

  x = ((struct data *) ntuple_data)->x;

  y = ((struct data *) ntuple_data)->y;

  z = ((struct data *) ntuple_data)->z;

  return (x*scale < 0.1);

}

double

val_func (void *ntuple_data, void * params)

{

  double x, y, z, scale;

  scale = *(double *)params;

  x = ((struct data *) ntuple_data)->x;

  y = ((struct data *) ntuple_data)->y;

  z = ((struct data *) ntuple_data)->z;

  return (x + y + z) * scale;

}