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