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

   1 /* histogram/test2d.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 <stdio.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_machine.h>
  26 #include <gsl/gsl_histogram2d.h>
  27 #include <gsl/gsl_test.h>
  28 #include <gsl/gsl_ieee_utils.h>
  29
  30 #define M 107
  31 #define N 239
  32 #define M1 17
  33 #define N1 23
  34 #define MR 10
  35 #define NR 5
  36
  37 void
  38 test2d (void)
  39 {
  40   double xr[MR + 1] =
  41     { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
  42
  43   double yr[NR + 1] = { 90.0, 91.0, 92.0, 93.0, 94.0, 95.0 };
  44
  45   gsl_histogram2d *h, *h1, *g, *hr;
  46   size_t i, j, k;
  47
  48   gsl_ieee_env_setup ();
  49
  50   h = gsl_histogram2d_calloc (M, N);
  51   h1 = gsl_histogram2d_calloc (M, N);
  52   g = gsl_histogram2d_calloc (M, N);
  53
  54   gsl_test (h->xrange == 0,
  55             "gsl_histogram2d_calloc returns valid xrange pointer");
  56   gsl_test (h->yrange == 0,
  57             "gsl_histogram2d_calloc returns valid yrange pointer");
  58   gsl_test (h->bin == 0, "gsl_histogram2d_calloc returns valid bin pointer");
  59
  60   gsl_test (h->nx != M, "gsl_histogram2d_calloc returns valid nx");
  61   gsl_test (h->ny != N, "gsl_histogram2d_calloc returns valid ny");
  62
  63   hr = gsl_histogram2d_calloc_range (MR, NR, xr, yr);
  64
  65   gsl_test (hr->xrange == 0,
  66             "gsl_histogram2d_calloc_range returns valid xrange pointer");
  67   gsl_test (hr->yrange == 0,
  68             "gsl_histogram2d_calloc_range returns valid yrange pointer");
  69   gsl_test (hr->bin == 0,
  70             "gsl_histogram2d_calloc_range returns valid bin pointer");
  71
  72   gsl_test (hr->nx != MR, "gsl_histogram2d_calloc_range returns valid nx");
  73   gsl_test (hr->ny != NR, "gsl_histogram2d_calloc_range returns valid ny");
  74
  75   {
  76     int status = 0;
  77     for (i = 0; i <= MR; i++)
  78       {
  79         if (hr->xrange[i] != xr[i])
  80           {
  81             status = 1;
  82           }
  83       };
  84
  85     gsl_test (status,
  86               "gsl_histogram2d_calloc_range creates xrange");
  87   }
  88
  89   {
  90     int status = 0;
  91     for (i = 0; i <= NR; i++)
  92       {
  93         if (hr->yrange[i] != yr[i])
  94           {
  95             status = 1;
  96           }
  97       };
  98
  99     gsl_test (status,
 100               "gsl_histogram2d_calloc_range creates yrange");
 101   }
 102
 103   for (i = 0; i <= MR; i++)
 104     {
 105       hr->xrange[i] = 0.0;
 106     }
 107
 108   for (i = 0; i <= NR; i++)
 109     {
 110       hr->yrange[i] = 0.0;
 111     }
 112
 113   {
 114     int status = gsl_histogram2d_set_ranges (hr, xr, MR + 1, yr, NR + 1);
 115
 116     for (i = 0; i <= MR; i++)
 117       {
 118         if (hr->xrange[i] != xr[i])
 119           {
 120             status = 1;
 121           }
 122       };
 123
 124     gsl_test (status, "gsl_histogram2d_set_ranges sets xrange");
 125   }
 126
 127   {
 128     int status = 0;
 129     for (i = 0; i <= NR; i++)
 130       {
 131         if (hr->yrange[i] != yr[i])
 132           {
 133             status = 1;
 134           }
 135       };
 136
 137     gsl_test (status, "gsl_histogram2d_set_ranges sets yrange");
 138   }
 139
 140
 141   k = 0;
 142   for (i = 0; i < M; i++)
 143     {
 144       for (j = 0; j < N; j++)
 145         {
 146           k++;
 147           gsl_histogram2d_accumulate (h, (double) i, (double) j, (double) k);
 148         };
 149     }
 150
 151   {
 152     int status = 0;
 153     k = 0;
 154     for (i = 0; i < M; i++)
 155       {
 156         for (j = 0; j < N; j++)
 157           {
 158             k++;
 159             if (h->bin[i * N + j] != (double) k)
 160               {
 161                 status = 1;
 162               }
 163           }
 164       }
 165     gsl_test (status,
 166               "gsl_histogram2d_accumulate writes into array");
 167   }
 168
 169   {
 170     int status = 0;
 171     k = 0;
 172     for (i = 0; i < M; i++)
 173       {
 174         for (j = 0; j < N; j++)
 175           {
 176             k++;
 177             if (gsl_histogram2d_get (h, i, j) != (double) k)
 178               status = 1;
 179           };
 180       }
 181     gsl_test (status, "gsl_histogram2d_get reads from array");
 182   }
 183
 184   for (i = 0; i <= M; i++)
 185     {
 186       h1->xrange[i] = 100.0 + i;
 187     }
 188
 189   for (i = 0; i <= N; i++)
 190     {
 191       h1->yrange[i] = 900.0 + i * i;
 192     }
 193
 194   gsl_histogram2d_memcpy (h1, h);
 195
 196   {
 197     int status = 0;
 198     for (i = 0; i <= M; i++)
 199       {
 200         if (h1->xrange[i] != h->xrange[i])
 201           status = 1;
 202       };
 203     gsl_test (status, "gsl_histogram2d_memcpy copies bin xranges");
 204   }
 205
 206   {
 207     int status = 0;
 208     for (i = 0; i <= N; i++)
 209       {
 210         if (h1->yrange[i] != h->yrange[i])
 211           status = 1;
 212       };
 213     gsl_test (status, "gsl_histogram2d_memcpy copies bin yranges");
 214   }
 215
 216   {
 217     int status = 0;
 218     for (i = 0; i < M; i++)
 219       {
 220         for (j = 0; j < N; j++)
 221           {
 222             if (gsl_histogram2d_get (h1, i, j) !=
 223                 gsl_histogram2d_get (h, i, j))
 224               status = 1;
 225           }
 226       }
 227     gsl_test (status, "gsl_histogram2d_memcpy copies bin values");
 228   }
 229
 230   gsl_histogram2d_free (h1);
 231
 232   h1 = gsl_histogram2d_clone (h);
 233
 234   {
 235     int status = 0;
 236     for (i = 0; i <= M; i++)
 237       {
 238         if (h1->xrange[i] != h->xrange[i])
 239           status = 1;
 240       };
 241     gsl_test (status, "gsl_histogram2d_clone copies bin xranges");
 242   }
 243
 244   {
 245     int status = 0;
 246     for (i = 0; i <= N; i++)
 247       {
 248         if (h1->yrange[i] != h->yrange[i])
 249           status = 1;
 250       };
 251     gsl_test (status, "gsl_histogram2d_clone copies bin yranges");
 252   }
 253
 254   {
 255     int status = 0;
 256     for (i = 0; i < M; i++)
 257       {
 258         for (j = 0; j < N; j++)
 259           {
 260             if (gsl_histogram2d_get (h1, i, j) !=
 261                 gsl_histogram2d_get (h, i, j))
 262               status = 1;
 263           }
 264       }
 265     gsl_test (status, "gsl_histogram2d_clone copies bin values");
 266   }
 267
 268
 269   gsl_histogram2d_reset (h);
 270
 271   {
 272     int status = 0;
 273
 274     for (i = 0; i < M * N; i++)
 275       {
 276         if (h->bin[i] != 0)
 277           status = 1;
 278       }
 279     gsl_test (status, "gsl_histogram2d_reset zeros array");
 280   }
 281
 282   gsl_histogram2d_free (h);
 283   h = gsl_histogram2d_calloc (M1, N1);
 284
 285   {
 286
 287     int status = 0;
 288
 289     for (i = 0; i < M1; i++)
 290       {
 291         for (j = 0; j < N1; j++)
 292           {
 293             gsl_histogram2d_increment (h, (double) i, (double) j);
 294
 295             for (k = 0; k <= i * N1 + j; k++)
 296               {
 297                 if (h->bin[k] != 1)
 298                   {
 299                     status = 1;
 300                   }
 301               }
 302
 303             for (k = i * N1 + j + 1; k < M1 * N1; k++)
 304               {
 305                 if (h->bin[k] != 0)
 306                   {
 307                     status = 1;
 308                   }
 309               }
 310           }
 311       }
 312     gsl_test (status, "gsl_histogram2d_increment increases bin value");
 313   }
 314
 315   gsl_histogram2d_free (h);
 316   h = gsl_histogram2d_calloc (M, N);
 317
 318   {
 319     int status = 0;
 320     for (i = 0; i < M; i++)
 321       {
 322         double x0 = 0, x1 = 0;
 323         gsl_histogram2d_get_xrange (h, i, &x0, &x1);
 324
 325         if (x0 != i || x1 != i + 1)
 326           {
 327             status = 1;
 328           }
 329       }
 330     gsl_test (status,
 331               "gsl_histogram2d_get_xlowerlimit and xupperlimit");
 332   }
 333
 334
 335   {
 336     int status = 0;
 337     for (i = 0; i < N; i++)
 338       {
 339         double y0 = 0, y1 = 0;
 340         gsl_histogram2d_get_yrange (h, i, &y0, &y1);
 341
 342         if (y0 != i || y1 != i + 1)
 343           {
 344             status = 1;
 345           }
 346       }
 347     gsl_test (status,
 348               "gsl_histogram2d_get_ylowerlimit and yupperlimit");
 349   }
 350
 351
 352   {
 353     int status = 0;
 354     if (gsl_histogram2d_xmax (h) != M)
 355       status = 1;
 356     gsl_test (status, "gsl_histogram2d_xmax");
 357   }
 358
 359   {
 360     int status = 0;
 361     if (gsl_histogram2d_xmin (h) != 0)
 362       status = 1;
 363     gsl_test (status, "gsl_histogram2d_xmin");
 364   }
 365
 366   {
 367     int status = 0;
 368     if (gsl_histogram2d_nx (h) != M)
 369       status = 1;
 370     gsl_test (status, "gsl_histogram2d_nx");
 371   }
 372
 373   {
 374     int status = 0;
 375     if (gsl_histogram2d_ymax (h) != N)
 376       status = 1;
 377     gsl_test (status, "gsl_histogram2d_ymax");
 378   }
 379
 380   {
 381     int status = 0;
 382     if (gsl_histogram2d_ymin (h) != 0)
 383       status = 1;
 384     gsl_test (status, "gsl_histogram2d_ymin");
 385   }
 386
 387   {
 388     int status = 0;
 389     if (gsl_histogram2d_ny (h) != N)
 390       status = 1;
 391     gsl_test (status, "gsl_histogram2d_ny");
 392   }
 393
 394   h->bin[3 * N + 2] = 123456.0;
 395   h->bin[4 * N + 3] = -654321;
 396
 397   {
 398     double max = gsl_histogram2d_max_val (h);
 399     gsl_test (max != 123456.0, "gsl_histogram2d_max_val finds maximum value");
 400   }
 401
 402   {
 403     double min = gsl_histogram2d_min_val (h);
 404     gsl_test (min != -654321.0,
 405               "gsl_histogram2d_min_val finds minimum value");
 406   }
 407
 408   {
 409     size_t imax, jmax;
 410     gsl_histogram2d_max_bin (h, &imax, &jmax);
 411     gsl_test (imax != 3
 412               || jmax != 2,
 413               "gsl_histogram2d_max_bin finds maximum value bin");
 414   }
 415
 416   {
 417     size_t imin, jmin;
 418     gsl_histogram2d_min_bin (h, &imin, &jmin);
 419     gsl_test (imin != 4
 420               || jmin != 3, "gsl_histogram2d_min_bin find minimum value bin");
 421   }
 422
 423   for (i = 0; i < M * N; i++)
 424     {
 425       h->bin[i] = i + 27;
 426       g->bin[i] = (i + 27) * (i + 1);
 427     }
 428
 429   {
 430     double sum = gsl_histogram2d_sum (h);
 431     gsl_test (sum != N * M * 27 + ((N * M - 1) * N * M) / 2,
 432               "gsl_histogram2d_sum sums all bin values");
 433   }
 434
 435   {
 436     /* first test... */
 437     const double xpos = 0.6;
 438     const double ypos = 0.85;
 439     double xmean;
 440     double ymean;
 441     size_t xbin;
 442     size_t ybin;
 443     gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
 444     gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
 445     gsl_histogram2d_increment (h3, xpos, ypos);
 446     gsl_histogram2d_find (h3, xpos, ypos, &xbin, &ybin);
 447     xmean = gsl_histogram2d_xmean (h3);
 448     ymean = gsl_histogram2d_ymean (h3);
 449
 450     {
 451       double expected_xmean = (h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0;
 452       double expected_ymean = (h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0;
 453       gsl_test_abs (xmean, expected_xmean, 100.0 * GSL_DBL_EPSILON,
 454                     "gsl_histogram2d_xmean");
 455       gsl_test_abs (ymean, expected_ymean, 100.0 * GSL_DBL_EPSILON,
 456                     "gsl_histogram2d_ymean");
 457     };
 458     gsl_histogram2d_free (h3);
 459   }
 460
 461   {
 462     /* test it with bivariate normal distribution */
 463     const double xmean = 0.7;
 464     const double ymean = 0.7;
 465     const double xsigma = 0.1;
 466     const double ysigma = 0.1;
 467     const double correl = 0.5;
 468     const double norm =
 469       10.0 / M_PI / xsigma / ysigma / sqrt (1.0 - correl * correl);
 470     size_t xbin;
 471     size_t ybin;
 472     gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
 473     gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
 474     /* initialize with 2d gauss pdf in two directions */
 475     for (xbin = 0; xbin < M; xbin++)
 476       {
 477         double xi =
 478           ((h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0 - xmean) / xsigma;
 479         for (ybin = 0; ybin < N; ybin++)
 480           {
 481             double yi =
 482               ((h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0 -
 483                ymean) / ysigma;
 484             double prob =
 485               norm * exp (-(xi * xi - 2.0 * correl * xi * yi + yi * yi) /
 486                           2.0 / (1 - correl * correl));
 487             h3->bin[xbin * N + ybin] = prob;
 488           }
 489       }
 490     {
 491       double xs = gsl_histogram2d_xsigma (h3);
 492       double ys = gsl_histogram2d_ysigma (h3);
 493       /* evaluate results and compare with parameters */
 494
 495       gsl_test_abs (gsl_histogram2d_xmean (h3), xmean, 2.0/M,
 496                     "gsl_histogram2d_xmean histogram mean(x)");
 497       gsl_test_abs (gsl_histogram2d_ymean (h3), ymean, 2.0/N,
 498                     "gsl_histogram2d_ymean histogram mean(y)");
 499       gsl_test_abs (xs, xsigma, 2.0/M,
 500                     "gsl_histogram2d_xsigma histogram stdev(x)");
 501       gsl_test_abs (ys, ysigma, 2.0/N,
 502                     "gsl_histogram2d_ysigma histogram stdev(y)");
 503       gsl_test_abs (gsl_histogram2d_cov (h3) / xs / ys, correl,
 504                     2.0/((M < N) ? M : N),
 505                     "gsl_histogram2d_cov histogram covariance");
 506     }
 507     gsl_histogram2d_free (h3);
 508   }
 509
 510   gsl_histogram2d_memcpy (h1, g);
 511   gsl_histogram2d_add (h1, h);
 512
 513   {
 514     int status = 0;
 515     for (i = 0; i < M * N; i++)
 516       {
 517         if (h1->bin[i] != g->bin[i] + h->bin[i])
 518           status = 1;
 519       }
 520     gsl_test (status, "gsl_histogram2d_add histogram addition");
 521   }
 522
 523   gsl_histogram2d_memcpy (h1, g);
 524   gsl_histogram2d_sub (h1, h);
 525
 526   {
 527     int status = 0;
 528     for (i = 0; i < M * N; i++)
 529       {
 530         if (h1->bin[i] != g->bin[i] - h->bin[i])
 531           status = 1;
 532       }
 533     gsl_test (status, "gsl_histogram2d_sub histogram subtraction");
 534   }
 535
 536
 537   gsl_histogram2d_memcpy (h1, g);
 538   gsl_histogram2d_mul (h1, h);
 539
 540   {
 541     int status = 0;
 542     for (i = 0; i < M * N; i++)
 543       {
 544         if (h1->bin[i] != g->bin[i] * h->bin[i])
 545           status = 1;
 546       }
 547     gsl_test (status, "gsl_histogram2d_mul histogram multiplication");
 548   }
 549
 550   gsl_histogram2d_memcpy (h1, g);
 551   gsl_histogram2d_div (h1, h);
 552
 553   {
 554     int status = 0;
 555     for (i = 0; i < M * N; i++)
 556       {
 557         if (h1->bin[i] != g->bin[i] / h->bin[i])
 558           status = 1;
 559       }
 560     gsl_test (status, "gsl_histogram2d_div histogram division");
 561   }
 562
 563   gsl_histogram2d_memcpy (h1, g);
 564   gsl_histogram2d_scale (h1, 0.5);
 565
 566   {
 567     int status = 0;
 568     for (i = 0; i < M * N; i++)
 569       {
 570         if (h1->bin[i] != 0.5 * g->bin[i])
 571           status = 1;
 572       }
 573     gsl_test (status, "gsl_histogram2d_scale histogram scaling");
 574   }
 575
 576   gsl_histogram2d_memcpy (h1, g);
 577   gsl_histogram2d_shift (h1, 0.25);
 578
 579   {
 580     int status = 0;
 581     for (i = 0; i < M * N; i++)
 582       {
 583         if (h1->bin[i] != 0.25 + g->bin[i])
 584           status = 1;
 585       }
 586     gsl_test (status, "gsl_histogram2d_shift histogram shift");
 587   }
 588
 589   gsl_histogram2d_free (h);     /* free whatever is in h */
 590
 591   h = gsl_histogram2d_calloc_uniform (M1, N1, 0.0, 5.0, 0.0, 5.0);
 592
 593   gsl_test (h->xrange == 0,
 594             "gsl_histogram2d_calloc_uniform returns valid range pointer");
 595   gsl_test (h->yrange == 0,
 596             "gsl_histogram2d_calloc_uniform returns valid range pointer");
 597   gsl_test (h->bin == 0,
 598             "gsl_histogram2d_calloc_uniform returns valid bin pointer");
 599   gsl_test (h->nx != M1, "gsl_histogram2d_calloc_uniform returns valid nx");
 600   gsl_test (h->ny != N1, "gsl_histogram2d_calloc_uniform returns valid ny");
 601
 602   gsl_histogram2d_accumulate (h, 0.0, 3.01, 1.0);
 603   gsl_histogram2d_accumulate (h, 0.1, 2.01, 2.0);
 604   gsl_histogram2d_accumulate (h, 0.2, 1.01, 3.0);
 605   gsl_histogram2d_accumulate (h, 0.3, 0.01, 4.0);
 606
 607   {
 608     size_t i1, i2, i3, i4;
 609     size_t j1, j2, j3, j4;
 610     double expected;
 611     int status;
 612     status = gsl_histogram2d_find (h, 0.0, 3.01, &i1, &j1);
 613     status = gsl_histogram2d_find (h, 0.1, 2.01, &i2, &j2);
 614     status = gsl_histogram2d_find (h, 0.2, 1.01, &i3, &j3);
 615     status = gsl_histogram2d_find (h, 0.3, 0.01, &i4, &j4);
 616
 617     for (i = 0; i < M1; i++)
 618       {
 619         for (j = 0; j < N1; j++)
 620           {
 621             if (i == i1 && j == j1)
 622               {
 623                 expected = 1.0;
 624               }
 625             else if (i == i2 && j == j2)
 626               {
 627                 expected = 2.0;
 628               }
 629             else if (i == i3 && j == j3)
 630               {
 631                 expected = 3.0;
 632               }
 633             else if (i == i4 && j == j4)
 634               {
 635                 expected = 4.0;
 636               }
 637             else
 638               {
 639                 expected = 0.0;
 640               }
 641
 642             if (h->bin[i * N1 + j] != expected)
 643               {
 644                 status = 1;
 645               }
 646           }
 647       }
 648     gsl_test (status, "gsl_histogram2d_find returns index");
 649   }
 650
 651   {
 652     FILE *f = fopen ("test.txt", "w");
 653     gsl_histogram2d_fprintf (f, h, "%.19e", "%.19e");
 654     fclose (f);
 655   }
 656
 657   {
 658     FILE *f = fopen ("test.txt", "r");
 659     gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
 660     int status = 0;
 661
 662     gsl_histogram2d_fscanf (f, hh);
 663
 664     for (i = 0; i <= M1; i++)
 665       {
 666         if (h->xrange[i] != hh->xrange[i])
 667           {
 668             printf ("xrange[%d] : %g orig vs %g\n",
 669                     (int) i, h->xrange[i], hh->xrange[i]);
 670             status = 1;
 671           }
 672       }
 673
 674     for (j = 0; j <= N1; j++)
 675       {
 676         if (h->yrange[j] != hh->yrange[j])
 677           {
 678             printf ("yrange[%d] : %g orig vs %g\n",
 679                     (int) j, h->yrange[j], hh->yrange[j]);
 680             status = 1;
 681           }
 682       }
 683
 684     for (i = 0; i < M1 * N1; i++)
 685       {
 686         if (h->bin[i] != hh->bin[i])
 687           {
 688             printf ("bin[%d] : %g orig vs %g\n",
 689                     (int) i, h->bin[i], hh->bin[i]);
 690             status = 1;
 691           }
 692       }
 693
 694     gsl_test (status, "gsl_histogram2d_fprintf and fscanf");
 695
 696     gsl_histogram2d_free (hh);
 697     fclose (f);
 698   }
 699
 700   {
 701     FILE *f = fopen ("test.dat", "wb");
 702     gsl_histogram2d_fwrite (f, h);
 703     fclose (f);
 704   }
 705
 706   {
 707     FILE *f = fopen ("test.dat", "rb");
 708     gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
 709     int status = 0;
 710
 711     gsl_histogram2d_fread (f, hh);
 712
 713     for (i = 0; i <= M1; i++)
 714       {
 715         if (h->xrange[i] != hh->xrange[i])
 716           {
 717             printf ("xrange[%d] : %g orig vs %g\n",
 718                     (int) i, h->xrange[i], hh->xrange[i]);
 719             status = 1;
 720           }
 721       }
 722
 723     for (j = 0; j <= N1; j++)
 724       {
 725         if (h->yrange[j] != hh->yrange[j])
 726           {
 727             printf ("yrange[%d] : %g orig vs %g\n",
 728                     (int) j, h->yrange[j], hh->yrange[j]);
 729             status = 1;
 730           }
 731       }
 732
 733     for (i = 0; i < M1 * N1; i++)
 734       {
 735         if (h->bin[i] != hh->bin[i])
 736           {
 737             printf ("bin[%d] : %g orig vs %g\n",
 738                     (int) i, h->bin[i], hh->bin[i]);
 739             status = 1;
 740           }
 741       }
 742
 743     gsl_test (status, "gsl_histogram2d_fwrite and fread");
 744
 745     gsl_histogram2d_free (hh);
 746     fclose (f);
 747   }
 748
 749   gsl_histogram2d_free (h);
 750   gsl_histogram2d_free (h1);
 751   gsl_histogram2d_free (g);
 752   gsl_histogram2d_free (hr);
 753 }