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