diff options
Diffstat (limited to 'gsl-1.9/wavelet/test.c')
| -rw-r--r-- | gsl-1.9/wavelet/test.c | 270 |
1 files changed, 270 insertions, 0 deletions
diff --git a/gsl-1.9/wavelet/test.c b/gsl-1.9/wavelet/test.c new file mode 100644 index 0000000..997e4aa --- /dev/null +++ b/gsl-1.9/wavelet/test.c @@ -0,0 +1,270 @@ +#include <config.h> +#include <stdio.h> +#include <math.h> + +#include <gsl/gsl_errno.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_test.h> + +#include <gsl/gsl_wavelet.h> +#include <gsl/gsl_wavelet2d.h> + +#define N_BS 11 + +double urand (void); + +double +urand (void) +{ + static unsigned long int x = 1; + x = (1103515245 * x + 12345) & 0x7fffffffUL; + return x / 2147483648.0; +} + +const size_t member[N_BS] = + { 309, 307, 305, 303, 301, 208, 206, 204, 202, 105, 103 }; + +void +test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member); + +void +test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type); + +int +main (int argc, char **argv) +{ + size_t i, N, stride, tda; + const int S = 1, NS = 2; /* Standard & Non-standard transforms */ + + /* One-dimensional tests */ + + for (N = 1; N <= 16384; N *= 2) + { + for (stride = 1; stride <= 5; stride++) + { + for (i = 0; i < N_BS; i++) + { + test_1d (N, stride, gsl_wavelet_bspline, member[i]); + test_1d (N, stride, gsl_wavelet_bspline_centered, member[i]); + } + + for (i = 4; i <= 20; i += 2) + { + test_1d (N, stride, gsl_wavelet_daubechies, i); + test_1d (N, stride, gsl_wavelet_daubechies_centered, i); + } + + test_1d (N, stride, gsl_wavelet_haar, 2); + test_1d (N, stride, gsl_wavelet_haar_centered, 2); + } + } + + /* Two-dimensional tests */ + + for (N = 1; N <= 64; N *= 2) + { + for (tda = N; tda <= N + 5; tda++) + { + for (i = 0; i < N_BS; i++) + { + test_2d (N, tda, gsl_wavelet_bspline, member[i], S); + test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], S); + + test_2d (N, tda, gsl_wavelet_bspline, member[i], NS); + test_2d (N, tda, gsl_wavelet_bspline_centered, member[i], NS); + } + + for (i = 4; i <= 20; i += 2) + { + test_2d (N, tda, gsl_wavelet_daubechies, i, S); + test_2d (N, tda, gsl_wavelet_daubechies_centered, i, S); + + test_2d (N, tda, gsl_wavelet_daubechies, i, NS); + test_2d (N, tda, gsl_wavelet_daubechies_centered, i, NS); + } + + test_2d (N, tda, gsl_wavelet_haar, 2, S); + test_2d (N, tda, gsl_wavelet_haar_centered, 2, S); + + test_2d (N, tda, gsl_wavelet_haar, 2, NS); + test_2d (N, tda, gsl_wavelet_haar_centered, 2, NS); + } + } + + exit (gsl_test_summary ()); +} + + +void +test_1d (size_t N, size_t stride, const gsl_wavelet_type * T, size_t member) +{ + gsl_wavelet_workspace *work; + gsl_vector *v1, *v2, *vdelta; + gsl_vector_view v; + gsl_wavelet *w; + + size_t i; + double *data = malloc (N * stride * sizeof (double)); + + for (i = 0; i < N * stride; i++) + data[i] = 12345.0 + i; + + v = gsl_vector_view_array_with_stride (data, stride, N); + v1 = &(v.vector); + + for (i = 0; i < N; i++) + { + gsl_vector_set (v1, i, urand ()); + } + + v2 = gsl_vector_alloc (N); + gsl_vector_memcpy (v2, v1); + + vdelta = gsl_vector_alloc (N); + + work = gsl_wavelet_workspace_alloc (N); + + w = gsl_wavelet_alloc (T, member); + + gsl_wavelet_transform_forward (w, v2->data, v2->stride, v2->size, work); + gsl_wavelet_transform_inverse (w, v2->data, v2->stride, v2->size, work); + + for (i = 0; i < N; i++) + { + double x1 = gsl_vector_get (v1, i); + double x2 = gsl_vector_get (v2, i); + gsl_vector_set (vdelta, i, fabs (x1 - x2)); + } + + { + double x1, x2; + i = gsl_vector_max_index (vdelta); + x1 = gsl_vector_get (v1, i); + x2 = gsl_vector_get (v2, i); + + gsl_test (fabs (x2 - x1) > N * 1e-15, + "%s(%d), n = %d, stride = %d, maxerr = %g", + gsl_wavelet_name (w), member, N, stride, fabs (x2 - x1)); + } + + if (stride > 1) + { + int status = 0; + + for (i = 0; i < N * stride; i++) + { + if (i % stride == 0) + continue; + + status |= (data[i] != (12345.0 + i)); + } + + gsl_test (status, "%s(%d) other data untouched, n = %d, stride = %d", + gsl_wavelet_name (w), member, N, stride); + } + + gsl_wavelet_workspace_free (work); + gsl_wavelet_free (w); + gsl_vector_free (vdelta); + gsl_vector_free (v2); + free (data); +} + + +void +test_2d (size_t N, size_t tda, const gsl_wavelet_type * T, size_t member, int type) +{ + gsl_wavelet_workspace *work; + gsl_matrix *m2; + gsl_wavelet *w; + gsl_matrix *m1; + gsl_matrix *mdelta; + gsl_matrix_view m; + size_t i; + size_t j; + + double *data = malloc (N * tda * sizeof (double)); + + const char * typename; + + typename = (type == 1) ? "standard" : "nonstd" ; + + for (i = 0; i < N * tda; i++) + data[i] = 12345.0 + i; + + m = gsl_matrix_view_array_with_tda (data, N, N, tda); + m1 = &(m.matrix); + + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + gsl_matrix_set (m1, i, j, urand()); + } + } + + m2 = gsl_matrix_alloc (N, N); + gsl_matrix_memcpy (m2, m1); + + mdelta = gsl_matrix_alloc (N, N); + + work = gsl_wavelet_workspace_alloc (N); + + w = gsl_wavelet_alloc (T, member); + + switch (type) + { + case 1: + gsl_wavelet2d_transform_matrix_forward (w, m2, work); + gsl_wavelet2d_transform_matrix_inverse (w, m2, work); + break; + case 2: + gsl_wavelet2d_nstransform_matrix_forward (w, m2, work); + gsl_wavelet2d_nstransform_matrix_inverse (w, m2, work); + break; + } + + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + double x1 = gsl_matrix_get (m1, i, j); + double x2 = gsl_matrix_get (m2, i, j ); + gsl_matrix_set (mdelta, i, j, fabs (x1 - x2)); + } + } + + { + double x1, x2; + gsl_matrix_max_index (mdelta, &i, &j); + x1 = gsl_matrix_get (m1, i, j); + x2 = gsl_matrix_get (m2, i, j); + + gsl_test (fabs (x2 - x1) > N * 1e-15, + "%s(%d)-2d %s, n = %d, tda = %d, maxerr = %g", + gsl_wavelet_name (w), member, typename, N, tda, fabs (x2 - x1)); + } + + if (tda > N) + { + int status = 0; + + for (i = 0; i < N ; i++) + { + for (j = N; j < tda; j++) + { + status |= (data[i*tda+j] != (12345.0 + (i*tda+j))); + } + } + + gsl_test (status, "%s(%d)-2d %s other data untouched, n = %d, tda = %d", + gsl_wavelet_name (w), member, typename, N, tda); + } + + free (data); + gsl_wavelet_workspace_free (work); + gsl_wavelet_free (w); + gsl_matrix_free (m2); + gsl_matrix_free (mdelta); +} |