htswanalysis/MACS/lib/gsl/gsl-1.11/wavelet/test.c

   1 #include <config.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4
   5 #include <gsl/gsl_errno.h>
   6 #include <gsl/gsl_vector.h>
   7 #include <gsl/gsl_blas.h>
   8 #include <gsl/gsl_test.h>
   9
  10 #include <gsl/gsl_wavelet.h>
  11 #include <gsl/gsl_wavelet2d.h>
  12
  13 #define N_BS 11
  14
  15 double urand (void);
  16
  17 double
  18 urand (void)
  19 {
  20   static unsigned long int x = 1;
  21   x = (1103515245 * x + 12345) & 0x7fffffffUL;
  22   return x / 2147483648.0;
  23 }
  24
  25 const size_t member[N_BS] =
  26   { 309, 307, 305, 303, 301, 208, 206, 204, 202, 105, 103 };
  27
  28 void
  29 test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member);
  30
  31 void
  32 test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type);
  33
  34 int
  35 main (int argc, char **argv)
  36 {
  37   size_t i, N, stride, tda;
  38   const int S = 1, NS = 2;  /* Standard & Non-standard transforms */
  39
  40   /* One-dimensional tests */
  41
  42   for (N = 1; N <= 16384; N *= 2)
  43     {
  44       for (stride = 1; stride <= 5; stride++)
  45         {
  46           for (i = 0; i < N_BS; i++)
  47             {
  48               test_1d (N, stride, gsl_wavelet_bspline, member[i]);
  49               test_1d (N, stride, gsl_wavelet_bspline_centered, member[i]);
  50             }
  51
  52           for (i = 4; i <= 20; i += 2)
  53             {
  54               test_1d (N, stride, gsl_wavelet_daubechies, i);
  55               test_1d (N, stride, gsl_wavelet_daubechies_centered, i);
  56             }
  57
  58           test_1d (N, stride, gsl_wavelet_haar, 2);
  59           test_1d (N, stride, gsl_wavelet_haar_centered, 2);
  60         }
  61     }
  62
  63   /* Two-dimensional tests */
  64
  65   for (N = 1; N <= 64; N *= 2)
  66     {
  67       for (tda = N; tda <= N + 5; tda++)
  68         {
  69           for (i = 0; i < N_BS; i++)
  70             {
  71               test_2d (N, tda, gsl_wavelet_bspline, member[i], S);
  72               test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], S);
  73
  74               test_2d (N, tda, gsl_wavelet_bspline, member[i], NS);
  75               test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], NS);
  76             }
  77
  78           for (i = 4; i <= 20; i += 2)
  79             {
  80               test_2d (N, tda, gsl_wavelet_daubechies, i, S);
  81               test_2d (N, tda, gsl_wavelet_daubechies_centered, i, S);
  82
  83               test_2d (N, tda, gsl_wavelet_daubechies, i, NS);
  84               test_2d (N, tda, gsl_wavelet_daubechies_centered, i, NS);
  85             }
  86
  87           test_2d (N, tda, gsl_wavelet_haar, 2, S);
  88           test_2d (N, tda, gsl_wavelet_haar_centered, 2, S);
  89
  90           test_2d (N, tda, gsl_wavelet_haar, 2, NS);
  91           test_2d (N, tda, gsl_wavelet_haar_centered, 2, NS);
  92         }
  93     }
  94
  95   exit (gsl_test_summary ());
  96 }
  97
  98
  99 void
 100 test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member)
 101 {
 102   gsl_wavelet_workspace *work;
 103   gsl_vector *v1, *v2, *vdelta;
 104   gsl_vector_view v;
 105   gsl_wavelet *w;
 106
 107   size_t i;
 108   double *data = malloc (N * stride * sizeof (double));
 109
 110   for (i = 0; i < N * stride; i++)
 111     data[i] = 12345.0 + i;
 112
 113   v = gsl_vector_view_array_with_stride (data, stride, N);
 114   v1 = &(v.vector);
 115
 116   for (i = 0; i < N; i++)
 117     {
 118       gsl_vector_set (v1, i, urand ());
 119     }
 120
 121   v2 = gsl_vector_alloc (N);
 122   gsl_vector_memcpy (v2, v1);
 123
 124   vdelta = gsl_vector_alloc (N);
 125
 126   work = gsl_wavelet_workspace_alloc (N);
 127
 128   w = gsl_wavelet_alloc (T, member);
 129
 130   gsl_wavelet_transform_forward (w, v2->data, v2->stride, v2->size, work);
 131   gsl_wavelet_transform_inverse (w, v2->data, v2->stride, v2->size, work);
 132
 133   for (i = 0; i < N; i++)
 134     {
 135       double x1 = gsl_vector_get (v1, i);
 136       double x2 = gsl_vector_get (v2, i);
 137       gsl_vector_set (vdelta, i, fabs (x1 - x2));
 138     }
 139
 140   {
 141     double x1, x2;
 142     i = gsl_vector_max_index (vdelta);
 143     x1 = gsl_vector_get (v1, i);
 144     x2 = gsl_vector_get (v2, i);
 145
 146     gsl_test (fabs (x2 - x1) > N * 1e-15,
 147               "%s(%d), n = %d, stride = %d, maxerr = %g",
 148               gsl_wavelet_name (w), member, N, stride, fabs (x2 - x1));
 149   }
 150
 151   if (stride > 1)
 152     {
 153       int status = 0;
 154
 155       for (i = 0; i < N * stride; i++)
 156         {
 157           if (i % stride == 0)
 158             continue;
 159
 160           status |= (data[i] != (12345.0 + i));
 161         }
 162
 163       gsl_test (status, "%s(%d) other data untouched, n = %d, stride = %d",
 164                 gsl_wavelet_name (w), member, N, stride);
 165     }
 166
 167   gsl_wavelet_workspace_free (work);
 168   gsl_wavelet_free (w);
 169   gsl_vector_free (vdelta);
 170   gsl_vector_free (v2);
 171   free (data);
 172 }
 173
 174
 175 void
 176 test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type)
 177 {
 178   gsl_wavelet_workspace *work;
 179   gsl_matrix *m2;
 180   gsl_wavelet *w;
 181   gsl_matrix *m1;
 182   gsl_matrix *mdelta;
 183   gsl_matrix_view m;
 184   size_t i;
 185   size_t j;
 186
 187   double *data = malloc (N * tda * sizeof (double));
 188
 189   const char * typename;
 190
 191   typename = (type == 1) ? "standard" : "nonstd" ;
 192
 193   for (i = 0; i < N * tda; i++)
 194     data[i] = 12345.0 + i;
 195
 196   m = gsl_matrix_view_array_with_tda (data, N, N, tda);
 197   m1 = &(m.matrix);
 198
 199   for (i = 0; i < N; i++)
 200     {
 201       for (j = 0; j < N; j++)
 202         {
 203             gsl_matrix_set (m1, i, j, urand());
 204         }
 205     }
 206
 207   m2 = gsl_matrix_alloc (N, N);
 208   gsl_matrix_memcpy (m2, m1);
 209
 210   mdelta = gsl_matrix_alloc (N, N);
 211
 212   work = gsl_wavelet_workspace_alloc (N);
 213
 214   w = gsl_wavelet_alloc (T, member);
 215
 216   switch (type)
 217     {
 218     case 1:
 219       gsl_wavelet2d_transform_matrix_forward (w, m2, work);
 220       gsl_wavelet2d_transform_matrix_inverse (w, m2, work);
 221       break;
 222     case 2:
 223       gsl_wavelet2d_nstransform_matrix_forward (w, m2, work);
 224       gsl_wavelet2d_nstransform_matrix_inverse (w, m2, work);
 225       break;
 226     }
 227
 228   for (i = 0; i < N; i++)
 229     {
 230       for (j = 0; j < N; j++)
 231         {
 232           double x1 = gsl_matrix_get (m1, i, j);
 233           double x2 = gsl_matrix_get (m2, i, j );
 234           gsl_matrix_set (mdelta, i, j, fabs (x1 - x2));
 235         }
 236     }
 237
 238   {
 239     double x1, x2;
 240     gsl_matrix_max_index (mdelta, &i, &j);
 241     x1 = gsl_matrix_get (m1, i, j);
 242     x2 = gsl_matrix_get (m2, i, j);
 243
 244     gsl_test (fabs (x2 - x1) > N * 1e-15,
 245               "%s(%d)-2d %s, n = %d, tda = %d, maxerr = %g",
 246               gsl_wavelet_name (w), member, typename, N, tda, fabs (x2 - x1));
 247   }
 248
 249   if (tda > N)
 250     {
 251       int status = 0;
 252
 253       for (i = 0; i < N ; i++)
 254         {
 255           for (j = N; j < tda; j++)
 256             {
 257               status |= (data[i*tda+j] != (12345.0 + (i*tda+j)));
 258             }
 259         }
 260
 261       gsl_test (status, "%s(%d)-2d %s other data untouched, n = %d, tda = %d",
 262                 gsl_wavelet_name (w), member, typename, N, tda);
 263     }
 264
 265   free (data);
 266   gsl_wavelet_workspace_free (work);
 267   gsl_wavelet_free (w);
 268   gsl_matrix_free (m2);
 269   gsl_matrix_free (mdelta);
 270 }