htswanalysis/MACS/lib/gsl/gsl-1.11/histogram/test2d_resample.c

   1 /* histogram/test2d_resample.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <gsl/gsl_histogram2d.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_ieee_utils.h>
  26
  27 #include "urand.c"
  28
  29 void
  30 test2d_resample (void)
  31 {
  32   size_t i, j;
  33   int status = 0;
  34   double total = 0;
  35   size_t N = 200000;
  36
  37   gsl_histogram2d *h;
  38
  39   gsl_ieee_env_setup ();
  40
  41   h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 0.0, 1.0);
  42
  43   for (i = 0; i < 10; i++)
  44     {
  45       for (j = 0; j < 10; j++)
  46         {
  47           double w = 10.0 * i + j;
  48           total += w;
  49           gsl_histogram2d_accumulate (h, 0.1 * i, 0.1 * i, w);
  50         }
  51     }
  52
  53   {
  54     gsl_histogram2d_pdf *p = gsl_histogram2d_pdf_alloc (10,10);
  55
  56     gsl_histogram2d *hh = gsl_histogram2d_calloc_uniform (20, 20,
  57                                                           0.0, 1.0,
  58                                                           0.0, 1.0);
  59
  60     gsl_histogram2d_pdf_init (p, h);
  61
  62     for (i = 0; i < N; i++)
  63       {
  64         double u = urand();
  65         double v = urand();
  66         double x, y;
  67         status = gsl_histogram2d_pdf_sample (p, u, v, &x, &y);
  68         status = gsl_histogram2d_increment (hh, x, y);
  69       }
  70
  71     status = 0;
  72     for (i = 0; i < 20; i++)
  73       {
  74         for (j = 0; j < 20; j++)
  75           {
  76             double z = 4 * total * gsl_histogram2d_get (hh, i, j) / (double) N;
  77             size_t k1, k2;
  78             double ya;
  79             double x, xmax, y, ymax;
  80
  81             gsl_histogram2d_get_xrange (hh, i, &x, &xmax);
  82             gsl_histogram2d_get_yrange (hh, j, &y, &ymax);
  83
  84             gsl_histogram2d_find (h, x, y, &k1, &k2);
  85             ya = gsl_histogram2d_get (h, k1, k2);
  86
  87             if (ya == 0)
  88               {
  89                 if (z != 0)
  90                   {
  91                     status = 1;
  92                     printf ("(%d,%d): %g vs %g\n", (int)i, (int)j, z, ya);
  93                   }
  94               }
  95             else
  96               {
  97                 double err = 1 / sqrt (gsl_histogram2d_get (hh, i, j));
  98                 double sigma = fabs ((z - ya) / (ya * err));
  99                 if (sigma > 3)
 100                   {
 101                     status = 1;
 102                     printf ("%g vs %g err=%g sigma=%g\n", z, ya, err, sigma);
 103                   }
 104               }
 105           }
 106       }
 107
 108     gsl_histogram2d_pdf_free (p) ;
 109     gsl_histogram2d_free (hh) ;
 110
 111     gsl_test (status, "gsl_histogram2d_pdf_sample within statistical errors");
 112   }
 113
 114   gsl_histogram2d_free (h) ;
 115 }