1 /* movstat/test_minmax.c
2  *
3  * Copyright (C) 2018 Patrick Alken
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 <gsl/gsl_math.h>
21 #include <gsl/gsl_vector.h>
22 #include <gsl/gsl_test.h>
23 #include <gsl/gsl_rng.h>
24 #include <gsl/gsl_movstat.h>
25 
26 /* compute filtered data by explicitely constructing window and finding min/max */
27 int
slow_minmax(const gsl_movstat_end_t etype,const gsl_vector * x,gsl_vector * y_min,gsl_vector * y_max,const int H,const int J)28 slow_minmax(const gsl_movstat_end_t etype, const gsl_vector * x, gsl_vector * y_min, gsl_vector * y_max,
29             const int H, const int J)
30 {
31   const size_t n = x->size;
32   const int K = H + J + 1;
33   double *window = malloc(K * sizeof(double));
34   size_t i;
35 
36   for (i = 0; i < n; ++i)
37     {
38       size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
39       gsl_vector_view v = gsl_vector_view_array(window, wsize);
40       double min, max;
41 
42       gsl_vector_minmax(&v.vector, &min, &max);
43       gsl_vector_set(y_min, i, min);
44       gsl_vector_set(y_max, i, max);
45     }
46 
47   free(window);
48 
49   return GSL_SUCCESS;
50 }
51 
52 static double
func_min(const size_t n,double x[],void * params)53 func_min(const size_t n, double x[], void * params)
54 {
55   gsl_vector_view v = gsl_vector_view_array(x, n);
56   (void) params;
57   return gsl_vector_min(&v.vector);
58 }
59 
60 static double
func_max(const size_t n,double x[],void * params)61 func_max(const size_t n, double x[], void * params)
62 {
63   gsl_vector_view v = gsl_vector_view_array(x, n);
64   (void) params;
65   return gsl_vector_max(&v.vector);
66 }
67 
68 static void
test_minmax_x(const double tol,const gsl_vector * x,const int H,const int J,const gsl_movstat_end_t endtype,const char * desc)69 test_minmax_x(const double tol, const gsl_vector * x, const int H, const int J,
70               const gsl_movstat_end_t endtype, const char * desc)
71 {
72   const size_t n = x->size;
73   gsl_vector * u_min = gsl_vector_alloc(n);
74   gsl_vector * y_min = gsl_vector_alloc(n);
75   gsl_vector * z_min = gsl_vector_alloc(n);
76   gsl_vector * u_max = gsl_vector_alloc(n);
77   gsl_vector * y_max = gsl_vector_alloc(n);
78   gsl_vector * z_max = gsl_vector_alloc(n);
79   gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
80   gsl_movstat_function F1, F2;
81   char buf[2048];
82 
83   F1.function = func_min;
84   F1.params = NULL;
85 
86   F2.function = func_max;
87   F2.params = NULL;
88 
89   /* compute moving min/max */
90   gsl_movstat_min(endtype, x, u_min, w);
91   gsl_movstat_max(endtype, x, u_max, w);
92   gsl_movstat_minmax(endtype, x, y_min, y_max, w);
93 
94   /* compute moving min/max with slow brute force method */
95   slow_minmax(endtype, x, z_min, z_max, H, J);
96 
97   sprintf(buf, "test_minmax: %s min endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
98   compare_vectors(tol, u_min, z_min, buf);
99 
100   sprintf(buf, "test_minmax: %s max endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
101   compare_vectors(tol, u_max, z_max, buf);
102 
103   sprintf(buf, "test_minmax: %s minmax(minimum) endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
104   compare_vectors(tol, y_min, z_min, buf);
105 
106   sprintf(buf, "test_minmax: %s minmax(maximum) endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
107   compare_vectors(tol, y_max, z_max, buf);
108 
109   /* in-place tests */
110 
111   gsl_vector_memcpy(u_min, x);
112   gsl_vector_memcpy(u_max, x);
113 
114   gsl_movstat_min(endtype, u_min, u_min, w);
115   gsl_movstat_max(endtype, u_max, u_max, w);
116 
117   sprintf(buf, "test_minmax: %s in-place min endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
118   compare_vectors(tol, u_min, z_min, buf);
119 
120   sprintf(buf, "test_minmax: %s in-place max endtype=%d n=%zu H=%d J=%d", desc, endtype, n, H, J);
121   compare_vectors(tol, u_max, z_max, buf);
122 
123   /* user-defined function tests */
124 
125   gsl_movstat_apply(endtype, &F1, x, z_min, w);
126   sprintf(buf, "n=%zu H=%d J=%d endtype=%u min user", n, H, J, endtype);
127   compare_vectors(tol, z_min, y_min, buf);
128 
129   gsl_movstat_apply(endtype, &F2, x, z_max, w);
130   sprintf(buf, "n=%zu H=%d J=%d endtype=%u max user", n, H, J, endtype);
131   compare_vectors(tol, z_max, y_max, buf);
132 
133   gsl_vector_free(u_min);
134   gsl_vector_free(y_min);
135   gsl_vector_free(z_min);
136   gsl_vector_free(u_max);
137   gsl_vector_free(y_max);
138   gsl_vector_free(z_max);
139   gsl_movstat_free(w);
140 }
141 
142 /* test alternating sequence [a,b,a,b,...] input */
143 static void
test_minmax_alt(const double tol,const size_t n,const int H,const int J,const gsl_movstat_end_t endtype)144 test_minmax_alt(const double tol, const size_t n, const int H, const int J,
145                 const gsl_movstat_end_t endtype)
146 {
147   const double a = 5.0;
148   const double b = -5.0;
149   gsl_vector * x = gsl_vector_alloc(n);
150   size_t i;
151 
152   for (i = 0; i < n; ++i)
153     {
154       if (i % 2 == 0)
155         gsl_vector_set(x, i, a);
156       else
157         gsl_vector_set(x, i, b);
158     }
159 
160   test_minmax_x(tol, x, H, J, endtype, "alternating");
161 
162   gsl_vector_free(x);
163 }
164 
165 /* test noisy sine wave input */
166 static void
test_minmax_sine(const double tol,const size_t n,const int H,const int J,const gsl_movstat_end_t endtype,gsl_rng * rng_p)167 test_minmax_sine(const double tol, const size_t n, const int H, const int J,
168                  const gsl_movstat_end_t endtype, gsl_rng * rng_p)
169 {
170   gsl_vector * x = gsl_vector_alloc(n);
171 
172   /* construct noisy sine signal */
173   test_noisy_sine(0.5, x, rng_p);
174 
175   test_minmax_x(tol, x, H, J, endtype, "noisy_sine");
176 
177   gsl_vector_free(x);
178 }
179 
180 /* test random input */
181 static void
test_minmax_random(const double tol,const size_t n,const int H,const int J,const gsl_movstat_end_t endtype,gsl_rng * rng_p)182 test_minmax_random(const double tol, const size_t n, const int H, const int J,
183                    const gsl_movstat_end_t endtype, gsl_rng * rng_p)
184 {
185   gsl_vector * x = gsl_vector_alloc(n);
186 
187   /* construct random input signal */
188   random_vector(x, rng_p);
189 
190   test_minmax_x(tol, x, H, J, endtype, "random");
191 
192   gsl_vector_free(x);
193 }
194 
195 static void
test_minmax(gsl_rng * rng_p)196 test_minmax(gsl_rng * rng_p)
197 {
198   /* alternating input */
199 
200   test_minmax_alt(GSL_DBL_EPSILON, 1000, 7, 7, GSL_MOVSTAT_END_PADZERO);
201   test_minmax_alt(GSL_DBL_EPSILON, 1000, 5, 2, GSL_MOVSTAT_END_PADZERO);
202   test_minmax_alt(GSL_DBL_EPSILON, 500, 1, 3, GSL_MOVSTAT_END_PADZERO);
203   test_minmax_alt(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO);
204   test_minmax_alt(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO);
205 
206   /* noisy sine wave input */
207 
208   test_minmax_sine(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADZERO, rng_p);
209   test_minmax_sine(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADZERO, rng_p);
210   test_minmax_sine(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
211   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
212   test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
213   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
214 
215   test_minmax_sine(GSL_DBL_EPSILON, 500, 5, 7, GSL_MOVSTAT_END_PADVALUE, rng_p);
216   test_minmax_sine(GSL_DBL_EPSILON, 1000, 10, 20, GSL_MOVSTAT_END_PADVALUE, rng_p);
217   test_minmax_sine(GSL_DBL_EPSILON, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
218   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
219   test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
220   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
221 
222   test_minmax_sine(GSL_DBL_EPSILON, 500, 5, 7, GSL_MOVSTAT_END_TRUNCATE, rng_p);
223   test_minmax_sine(GSL_DBL_EPSILON, 1000, 10, 20, GSL_MOVSTAT_END_TRUNCATE, rng_p);
224   test_minmax_sine(GSL_DBL_EPSILON, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
225   test_minmax_sine(GSL_DBL_EPSILON, 1000, 30, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
226   test_minmax_sine(GSL_DBL_EPSILON, 1000, 5, 30, GSL_MOVSTAT_END_TRUNCATE, rng_p);
227   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
228   test_minmax_sine(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
229   test_minmax_sine(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
230 
231   /* random input */
232 
233   test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
234   test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADZERO, rng_p);
235   test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADZERO, rng_p);
236   test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
237   test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
238   test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
239   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
240   test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
241   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
242 
243   test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
244   test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_PADVALUE, rng_p);
245   test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_PADVALUE, rng_p);
246   test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
247   test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
248   test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
249   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
250   test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
251   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
252 
253   test_minmax_random(GSL_DBL_EPSILON, 1000, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
254   test_minmax_random(GSL_DBL_EPSILON, 1000, 5, 7, GSL_MOVSTAT_END_TRUNCATE, rng_p);
255   test_minmax_random(GSL_DBL_EPSILON, 2000, 0, 2, GSL_MOVSTAT_END_TRUNCATE, rng_p);
256   test_minmax_random(GSL_DBL_EPSILON, 500, 3, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
257   test_minmax_random(GSL_DBL_EPSILON, 500, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
258   test_minmax_random(GSL_DBL_EPSILON, 500, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
259   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
260   test_minmax_random(GSL_DBL_EPSILON, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
261   test_minmax_random(GSL_DBL_EPSILON, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
262 }
263