1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3
4 // This file is part of the Compressed Continuous Computation (C3) Library
5 // Author: Alex A. Gorodetsky
6 // Contact: alex@alexgorodetsky.com
7
8 // All rights reserved.
9
10 // Redistribution and use in source and binary forms, with or without modification,
11 // are permitted provided that the following conditions are met:
12
13 // 1. Redistributions of source code must retain the above copyright notice,
14 // this list of conditions and the following disclaimer.
15
16 // 2. Redistributions in binary form must reproduce the above copyright notice,
17 // this list of conditions and the following disclaimer in the documentation
18 // and/or other materials provided with the distribution.
19
20 // 3. Neither the name of the copyright holder nor the names of its contributors
21 // may be used to endorse or promote products derived from this software
22 // without specific prior written permission.
23
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
25 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
27 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
28 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
29 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
30 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
32 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
33 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34
35 //Code
36
37 #include <stdlib.h>
38 #include <stdio.h>
39 #include <math.h>
40 #include <string.h>
41 #include <assert.h>
42 #include <float.h>
43
44 #include "CuTest.h"
45 #include "testfunctions.h"
46
47 #include "array.h"
48 #include "lib_linalg.h"
49 #include "lib_optimization.h"
50 #include "lib_funcs.h"
51
52 typedef struct OrthPolyExpansion* opoly_t;
53 #define POLY_EVAL orth_poly_expansion_eval
54 #define POLY_FREE orth_poly_expansion_free
55
56 static void
compute_error(double lb,double ub,size_t N,opoly_t cpoly,double (* func)(double,void *),void * arg,double * abs_err,double * func_norm)57 compute_error(double lb,double ub, size_t N, opoly_t cpoly,
58 double (*func)(double,void*), void* arg,
59 double * abs_err, double * func_norm)
60 {
61 double * xtest = linspace(lb,ub,N);
62 size_t ii;
63 *abs_err = 0.0;
64 *func_norm = 0.0;
65 for (ii = 0; ii < N; ii++){
66 *abs_err += pow(POLY_EVAL(cpoly,xtest[ii]) - func(xtest[ii],arg),2);
67 *func_norm += pow(func(xtest[ii],arg),2);
68 }
69 free(xtest); xtest = NULL;
70 }
71
72 static void
compute_error_vec(double lb,double ub,size_t N,opoly_t cpoly,int (* func)(size_t,const double *,double *,void *),void * arg,double * abs_err,double * func_norm)73 compute_error_vec(double lb,double ub, size_t N, opoly_t cpoly,
74 int (*func)(size_t, const double *,double *,void*),
75 void* arg,
76 double * abs_err, double * func_norm)
77 {
78 double * xtest = linspace(lb,ub,N);
79 size_t ii;
80 *abs_err = 0.0;
81 *func_norm = 0.0;
82 double val;
83 for (ii = 0; ii < N; ii++){
84 func(1,xtest+ii,&val,arg);
85 *abs_err += pow(POLY_EVAL(cpoly,xtest[ii]) - val,2);
86 *func_norm += pow(val,2);
87 }
88 free(xtest); xtest = NULL;
89 }
90
91
Test_cheb_approx(CuTest * tc)92 void Test_cheb_approx(CuTest * tc){
93
94 printf("Testing function: cheb_approx\n");
95
96 // function
97 struct Fwrap * fw = fwrap_create(1,"general-vec");
98 fwrap_set_fvec(fw,Sin3xTx2,NULL);
99
100 // approximation
101 size_t N = 50;
102 double lb=-1.0,ub=1.0;
103 opoly_t cpoly = orth_poly_expansion_init(CHEBYSHEV,N,lb,ub);
104 int res = orth_poly_expansion_approx_vec(cpoly,fw,NULL);
105 CuAssertIntEquals(tc,0,res);
106
107 // compute error
108 double abs_err;
109 double func_norm;
110 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
111 double err = abs_err / func_norm;
112 CuAssertDblEquals(tc, 0.0, err, 1e-13);
113
114 fwrap_destroy(fw);
115 POLY_FREE(cpoly);
116
117 }
118
Test_cheb_approx_nonnormal(CuTest * tc)119 void Test_cheb_approx_nonnormal(CuTest * tc){
120
121 printf("Testing function: cheb_approx on (a,b)\n");
122
123 // function
124 struct Fwrap * fw = fwrap_create(1,"general-vec");
125 fwrap_set_fvec(fw,Sin3xTx2,NULL);
126
127 size_t N = 50;
128 double lb = -2,ub = 3;
129 opoly_t cpoly = orth_poly_expansion_init(CHEBYSHEV,N,lb,ub);
130 orth_poly_expansion_approx_vec(cpoly,fw,NULL);
131
132 double abs_err;
133 double func_norm;
134 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
135 double err = abs_err / func_norm;
136 CuAssertDblEquals(tc, 0.0, err, 1e-13);
137
138 POLY_FREE(cpoly);
139 fwrap_destroy(fw);
140 }
141
Test_cheb_approx_adapt(CuTest * tc)142 void Test_cheb_approx_adapt(CuTest * tc){
143
144 printf("Testing function: cheb_approx_adapt\n");
145
146 // function
147 struct Fwrap * fw = fwrap_create(1,"general-vec");
148 fwrap_set_fvec(fw,Sin3xTx2,NULL);
149
150 // approximation
151 double lb = -1.0, ub = 1.0;
152 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
153 ope_opts_set_start(opts,10);
154 ope_opts_set_coeffs_check(opts,4);
155 ope_opts_set_tol(opts,1e-8);
156 ope_opts_set_lb(opts,lb);
157 ope_opts_set_ub(opts,ub);
158 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
159
160 // error
161 double abs_err;
162 double func_norm;
163 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
164 double err = abs_err / func_norm;
165 CuAssertDblEquals(tc, 0.0, err, 1e-13);
166
167 POLY_FREE(cpoly);
168 ope_opts_free(opts);
169 fwrap_destroy(fw);
170 }
171
Test_cheb_approx_adapt_weird(CuTest * tc)172 void Test_cheb_approx_adapt_weird(CuTest * tc){
173
174 printf("Testing function: cheb_approx_adapt on (a,b)\n");
175 // function
176 struct Fwrap * fw = fwrap_create(1,"general-vec");
177 fwrap_set_fvec(fw,Sin3xTx2,NULL);
178
179 // approximation
180 double lb = -2.0, ub = -1.0;
181 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
182 ope_opts_set_start(opts,10);
183 ope_opts_set_coeffs_check(opts,4);
184 ope_opts_set_tol(opts,1e-8);
185 ope_opts_set_lb(opts,lb);
186 ope_opts_set_ub(opts,ub);
187 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
188
189 // error
190 double abs_err;
191 double func_norm;
192 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
193 double err = abs_err / func_norm;
194 CuAssertDblEquals(tc, 0.0, err, 1e-13);
195
196 POLY_FREE(cpoly);
197 ope_opts_free(opts);
198 fwrap_destroy(fw);
199
200 }
201
Test_cheb_derivative(CuTest * tc)202 void Test_cheb_derivative(CuTest * tc){
203
204 printf("Testing function: orth_poly_expansion_deriv with chebyshev poly on (a,b) (1/2) \n");
205
206 // function
207 struct Fwrap * fw = fwrap_create(1,"general-vec");
208 fwrap_set_fvec(fw,Sin3xTx2,NULL);
209
210 // approximation
211 double lb = -2.0, ub = -1.0;
212 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
213 ope_opts_set_start(opts,10);
214 ope_opts_set_coeffs_check(opts,4);
215 ope_opts_set_tol(opts,1e-8);
216 ope_opts_set_lb(opts,lb);
217 ope_opts_set_ub(opts,ub);
218 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
219 opoly_t der = orth_poly_expansion_deriv(cpoly);
220
221 // error
222 double abs_err;
223 double func_norm;
224 compute_error(lb,ub,1000,der,funcderiv,NULL,&abs_err,&func_norm);
225 double err = abs_err / func_norm;
226 CuAssertDblEquals(tc, 0.0, err, 1e-20);
227
228 POLY_FREE(cpoly);
229 POLY_FREE(der);
230 ope_opts_free(opts);
231 fwrap_destroy(fw);
232
233 }
234
Test_cheb_derivative2(CuTest * tc)235 void Test_cheb_derivative2(CuTest * tc){
236
237 printf("Testing function: orth_poly_expansion_deriv with chebyshev poly on (a,b) (2/2) \n");
238
239 // function
240 struct Fwrap * fw = fwrap_create(1,"general-vec");
241 fwrap_set_fvec(fw,x3minusx,NULL);
242
243 // approximation
244 double lb = -2.0, ub = -1.0;
245 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
246 ope_opts_set_start(opts,10);
247 ope_opts_set_coeffs_check(opts,4);
248 ope_opts_set_tol(opts,1e-20);
249 ope_opts_set_lb(opts,lb);
250 ope_opts_set_ub(opts,ub);
251 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
252 opoly_t der = orth_poly_expansion_deriv(cpoly);
253
254 // error
255 double abs_err;
256 double func_norm;
257 compute_error(lb,ub,1000,der,x3minusxd,NULL,&abs_err,&func_norm);
258 double err = abs_err / func_norm;
259 CuAssertDblEquals(tc, 0.0, err, 1e-20);
260
261 POLY_FREE(cpoly);
262 POLY_FREE(der);
263 ope_opts_free(opts);
264 fwrap_destroy(fw);
265
266 }
267
Test_cheb_dderivative(CuTest * tc)268 void Test_cheb_dderivative(CuTest * tc){
269
270 printf("Testing function: orth_poly_expansion_dderiv with chebyshev poly on (a,b) (1/2) \n");
271
272 // function
273 struct Fwrap * fw = fwrap_create(1,"general-vec");
274 fwrap_set_fvec(fw,x3minusx,NULL);
275
276 // approximation
277 double lb = -10.0, ub = -5.0;
278 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
279 ope_opts_set_start(opts,10);
280 ope_opts_set_coeffs_check(opts,4);
281 ope_opts_set_tol(opts,1e-40);
282 ope_opts_set_lb(opts,lb);
283 ope_opts_set_ub(opts,ub);
284 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
285 opoly_t der = orth_poly_expansion_dderiv(cpoly);
286
287 // error
288 double abs_err;
289 double func_norm;
290 compute_error(lb,ub,1000,der,x3minusxdd,NULL,&abs_err,&func_norm);
291 double err = abs_err / func_norm;
292 CuAssertDblEquals(tc, 0.0, err, 1e-14);
293
294 POLY_FREE(cpoly);
295 POLY_FREE(der);
296 ope_opts_free(opts);
297 fwrap_destroy(fw);
298
299 }
300
301
Test_cheb_dderivative2(CuTest * tc)302 void Test_cheb_dderivative2(CuTest * tc){
303
304 printf("Testing function: orth_poly_expansion_dderiv with chebyshev poly on (a,b) (2/2) \n");
305
306 // function
307 struct Fwrap * fw = fwrap_create(1,"general-vec");
308 fwrap_set_fvec(fw,Sin3xTx2,NULL);
309
310 // approximation
311 double lb = -1.0, ub = 1.0;
312 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
313 ope_opts_set_start(opts,10);
314 ope_opts_set_coeffs_check(opts,4);
315 ope_opts_set_tol(opts,1e-20);
316 ope_opts_set_lb(opts,lb);
317 ope_opts_set_ub(opts,ub);
318 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
319 opoly_t der = orth_poly_expansion_dderiv(cpoly);
320
321 // error
322 double abs_err;
323 double func_norm;
324 compute_error(lb,ub,1000,der,funcdderiv,NULL,&abs_err,&func_norm);
325 double err = abs_err / func_norm;
326 CuAssertDblEquals(tc, 0.0, err, 1e-15);
327
328 POLY_FREE(cpoly);
329 POLY_FREE(der);
330 ope_opts_free(opts);
331 fwrap_destroy(fw);
332
333 }
334
Test_cheb_integrate(CuTest * tc)335 void Test_cheb_integrate(CuTest * tc){
336
337 printf("Testing function: cheb_integrate2\n");
338
339 // function
340 struct Fwrap * fw = fwrap_create(1,"general-vec");
341 fwrap_set_fvec(fw,powX2,NULL);
342
343 // approximation
344 double lb = -2.0, ub = 3.0;
345 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
346 ope_opts_set_start(opts,10);
347 ope_opts_set_coeffs_check(opts,4);
348 ope_opts_set_tol(opts,1e-8);
349 ope_opts_set_lb(opts,lb);
350 ope_opts_set_ub(opts,ub);
351 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
352
353 // error
354 double intshould = (pow(ub,3) - pow(lb,3))/3;
355 double intis = cheb_integrate2(cpoly);
356 CuAssertDblEquals(tc, intshould, intis, 1e-13);
357
358 POLY_FREE(cpoly);
359 ope_opts_free(opts);
360 fwrap_destroy(fw);
361 }
362
Test_cheb_inner(CuTest * tc)363 void Test_cheb_inner(CuTest * tc){
364
365 printf("Testing function: orth_poly_expansion_inner with chebyshev poly \n");
366
367 // function
368 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
369 fwrap_set_fvec(fw1,powX2,NULL);
370
371 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
372 fwrap_set_fvec(fw2,TwoPowX3,NULL);
373
374 double lb = -2.0, ub = 3.0;
375 /* double lb = -1.0, ub = 1.0; */
376 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
377 ope_opts_set_start(opts,10);
378 ope_opts_set_coeffs_check(opts,4);
379 ope_opts_set_tol(opts,1e-10);
380 ope_opts_set_lb(opts,lb);
381 ope_opts_set_ub(opts,ub);
382
383 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
384 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
385
386 double intshould = (pow(ub,6) - pow(lb,6))/3;
387 double intis = orth_poly_expansion_inner(cpoly,cpoly2);
388 CuAssertDblEquals(tc, intshould, intis, 1e-10);
389
390 POLY_FREE(cpoly);
391 POLY_FREE(cpoly2);
392 ope_opts_free(opts);
393 fwrap_destroy(fw1);
394 fwrap_destroy(fw2);
395 }
396
Test_cheb_norm(CuTest * tc)397 void Test_cheb_norm(CuTest * tc){
398
399 printf("Testing function: orth_poly_expansion_norm with chebyshev poly\n");
400
401 // function
402 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
403 fwrap_set_fvec(fw1,powX2,NULL);
404
405 double lb = -2, ub = 3.0;
406 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
407 ope_opts_set_start(opts,10);
408 ope_opts_set_coeffs_check(opts,4);
409 ope_opts_set_tol(opts,1e-8);
410 ope_opts_set_lb(opts,lb);
411 ope_opts_set_ub(opts,ub);
412
413 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
414
415 double intshould = (pow(ub,5) - pow(lb,5))/5;
416 double intis = orth_poly_expansion_norm(cpoly);
417 CuAssertDblEquals(tc, sqrt(intshould), intis, 1e-10);
418
419 POLY_FREE(cpoly);
420 ope_opts_free(opts);
421 fwrap_destroy(fw1);
422 }
423
Test_cheb_product(CuTest * tc)424 void Test_cheb_product(CuTest * tc){
425
426 printf("Testing function: orth_poly_expansion_product with chebyshev poly \n");
427
428 // function
429 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
430 fwrap_set_fvec(fw1,powX2,NULL);
431
432 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
433 fwrap_set_fvec(fw2,TwoPowX3,NULL);
434
435 // approximation
436 double lb = -3.0, ub = 2.0;
437 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
438 ope_opts_set_start(opts,10);
439 ope_opts_set_coeffs_check(opts,4);
440 ope_opts_set_tol(opts,1e-10);
441 ope_opts_set_lb(opts,lb);
442 ope_opts_set_ub(opts,ub);
443
444 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
445 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
446 opoly_t cpoly3 = orth_poly_expansion_prod(cpoly,cpoly2);
447
448 size_t N = 100;
449 double * pts = linspace(lb,ub,N);
450 size_t ii;
451 for (ii = 0; ii < N; ii++){
452 double eval1 = POLY_EVAL(cpoly3,pts[ii]);
453 double eval2 = POLY_EVAL(cpoly,pts[ii]) *
454 POLY_EVAL(cpoly2,pts[ii]);
455 double diff= fabs(eval1-eval2);
456 CuAssertDblEquals(tc, 0.0, diff, 1e-10);
457 }
458 free(pts); pts = NULL;
459
460 POLY_FREE(cpoly);
461 POLY_FREE(cpoly2);
462 POLY_FREE(cpoly3);
463 ope_opts_free(opts);
464 fwrap_destroy(fw1);
465 fwrap_destroy(fw2);
466 }
467
Test_cheb_orth_poly_expansion_create_with_params_and_grad(CuTest * tc)468 void Test_cheb_orth_poly_expansion_create_with_params_and_grad(CuTest * tc){
469
470 printf("Testing functions: orth_poly_expansion_create_with_params and grad with chebyshev poly\n");
471
472 size_t nparams = 10;
473 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
474 0.6,-0.9,-.2,0.04,1.2};
475 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
476 struct OrthPolyExpansion * ope = NULL;
477 ope = orth_poly_expansion_create_with_params(opts,nparams,params);
478
479
480 double grad[10];
481 double xloc = 0.4;
482 int res =
483 orth_poly_expansion_param_grad_eval(ope,1,&xloc,grad);
484 CuAssertIntEquals(tc,0,res);
485
486
487 // numerical derivative
488 struct OrthPolyExpansion * ope1 = NULL;
489 struct OrthPolyExpansion * ope2 = NULL;
490
491 size_t dim = nparams;
492 double * x1 = calloc_double(dim);
493 double * x2 = calloc_double(dim);
494 for (size_t ii = 0; ii < dim; ii++){
495 x1[ii] = params[ii];
496 x2[ii] = params[ii];
497 }
498
499 double diff = 0.0;
500 double v1,v2;
501 double norm = 0.0;
502 double eps = 1e-8;
503 for (size_t ii = 0; ii < dim; ii++){
504 x1[ii] += eps;
505 x2[ii] -= eps;
506 ope1 = orth_poly_expansion_create_with_params(opts,nparams,x1);
507 v1 = orth_poly_expansion_eval(ope1,xloc);
508
509 ope2 = orth_poly_expansion_create_with_params(opts,nparams,x2);
510 v2 = orth_poly_expansion_eval(ope2,xloc);
511
512 double diff_iter = pow( (v1-v2)/2.0/eps - grad[ii], 2 );
513 /* printf("current diff = %G\n",diff_iter); */
514 /* printf("\t norm = %G\n",grad[ii]); */
515 diff += diff_iter;
516 norm += pow( (v1-v2)/2.0/eps,2);
517
518 x1[ii] -= eps;
519 x2[ii] += eps;
520
521 orth_poly_expansion_free(ope1); ope1 = NULL;
522 orth_poly_expansion_free(ope2); ope2 = NULL;
523 }
524 if (norm > 1){
525 diff /= norm;
526 }
527 free(x1); x1 = NULL;
528 free(x2); x2 = NULL;
529 CuAssertDblEquals(tc,0.0,diff,1e-7);
530
531 ope_opts_free(opts); opts = NULL;
532 orth_poly_expansion_free(ope); ope = NULL;
533 }
534
Test_orth_poly_expansion_real_roots_chebyshev(CuTest * tc)535 void Test_orth_poly_expansion_real_roots_chebyshev(CuTest * tc){
536
537 printf("Testing function: orth_poly_expansion_real_roots with chebyshev\n");
538
539
540 struct Fwrap * fw = fwrap_create(1,"general-vec");
541 fwrap_set_fvec(fw,polyroots,NULL);
542
543 // approximation
544 double lb = -3.0, ub = 2.0;
545 /* double lb = -1.0, ub = 1.0; */
546 struct OpeOpts * opts = ope_opts_alloc(CHEBYSHEV);
547 ope_opts_set_start(opts,6);
548 ope_opts_set_coeffs_check(opts,2);
549 ope_opts_set_tol(opts,1e-13);
550 ope_opts_set_lb(opts,lb);
551 ope_opts_set_ub(opts,ub);
552 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
553
554 size_t nroots;
555 double * roots = orth_poly_expansion_real_roots(pl, &nroots);
556
557 /* printf("roots are: "); */
558 /* dprint(nroots, roots); */
559
560 CuAssertIntEquals(tc, 5, nroots);
561 CuAssertDblEquals(tc, -3.0, roots[0], 1e-9);
562 CuAssertDblEquals(tc, 0.0, roots[1], 1e-9);
563 CuAssertDblEquals(tc, 1.0, roots[2], 1e-5);
564 CuAssertDblEquals(tc, 1.0, roots[3], 1e-5);
565 CuAssertDblEquals(tc, 2.0, roots[4], 1e-9);
566
567 free(roots);
568 POLY_FREE(pl);
569 ope_opts_free(opts);
570 fwrap_destroy(fw);
571 }
572
573
ChebGetSuite()574 CuSuite * ChebGetSuite(){
575
576 CuSuite * suite = CuSuiteNew();
577 SUITE_ADD_TEST(suite, Test_cheb_approx);
578 SUITE_ADD_TEST(suite, Test_cheb_approx_nonnormal);
579 SUITE_ADD_TEST(suite, Test_cheb_approx_adapt);
580 SUITE_ADD_TEST(suite, Test_cheb_approx_adapt_weird);
581 SUITE_ADD_TEST(suite, Test_cheb_derivative);
582 SUITE_ADD_TEST(suite, Test_cheb_derivative2);
583 SUITE_ADD_TEST(suite, Test_cheb_dderivative);
584 SUITE_ADD_TEST(suite, Test_cheb_dderivative2);
585 SUITE_ADD_TEST(suite, Test_cheb_integrate);
586 SUITE_ADD_TEST(suite, Test_cheb_inner);
587 SUITE_ADD_TEST(suite, Test_cheb_norm);
588 SUITE_ADD_TEST(suite, Test_cheb_orth_poly_expansion_create_with_params_and_grad);
589 SUITE_ADD_TEST(suite, Test_orth_poly_expansion_real_roots_chebyshev);
590 SUITE_ADD_TEST(suite, Test_cheb_product);
591
592 return suite;
593 }
594
595 ///////////////////////////////////////////////////////////////////////////
596
Test_legendre_approx(CuTest * tc)597 void Test_legendre_approx(CuTest * tc){
598
599 printf("Testing function: legendre_approx\n");
600
601 // function
602 struct Fwrap * fw = fwrap_create(1,"general-vec");
603 fwrap_set_fvec(fw,Sin3xTx2,NULL);
604
605 // approximation
606 size_t N = 50;
607 double lb=-1.0,ub=1.0;
608 opoly_t cpoly = orth_poly_expansion_init(LEGENDRE,N,lb,ub);
609 int res = orth_poly_expansion_approx_vec(cpoly,fw,NULL);
610 CuAssertIntEquals(tc,0,res);
611
612 // compute error
613 double abs_err;
614 double func_norm;
615 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
616 double err = abs_err / func_norm;
617 CuAssertDblEquals(tc, 0.0, err, 1e-13);
618
619 fwrap_destroy(fw);
620 POLY_FREE(cpoly);
621 }
622
623
624
Test_legendre_approx_gaussbump(CuTest * tc)625 void Test_legendre_approx_gaussbump(CuTest * tc){
626
627 printf("Testing function: legendre_approx_gaussbump\n");
628
629 // function
630 struct Fwrap * fw = fwrap_create(1,"general-vec");
631 fwrap_set_fvec(fw,gaussbump,NULL);
632
633 // approximation
634 /* size_t N = 50; */
635 double lb=-1.0,ub=1.0;
636 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
637 ope_opts_set_start(opts,10);
638 ope_opts_set_coeffs_check(opts,4);
639 ope_opts_set_tol(opts,1e-8);
640 ope_opts_set_lb(opts,lb);
641 ope_opts_set_ub(opts,ub);
642 /* opoly_t cpoly = orth_poly_expansion_init(LEGENDRE,N,lb,ub); */
643 /* int res = orth_poly_expansion_approx_vec(cpoly,fw,NULL); */
644 /* CuAssertIntEquals(tc,0,res); */
645 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
646
647 // compute error
648 double abs_err;
649 double func_norm;
650 compute_error_vec(lb,ub,1000,cpoly,gaussbump,NULL,&abs_err,&func_norm);
651 double err = abs_err / func_norm;
652 CuAssertDblEquals(tc, 0.0, err, 1e-15);
653
654 ope_opts_free(opts);
655 fwrap_destroy(fw);
656 POLY_FREE(cpoly);
657 }
658
Test_legendre_approx_nonnormal(CuTest * tc)659 void Test_legendre_approx_nonnormal(CuTest * tc){
660
661 printf("Testing function: legendre_approx on (a,b)\n");
662
663 // function
664 struct Fwrap * fw = fwrap_create(1,"general-vec");
665 fwrap_set_fvec(fw,Sin3xTx2,NULL);
666
667 // approximation
668 size_t N = 50;
669 double lb = -2,ub = 3;
670 opoly_t cpoly = orth_poly_expansion_init(LEGENDRE,N,lb,ub);
671 int res = orth_poly_expansion_approx_vec(cpoly,fw,NULL);
672 CuAssertIntEquals(tc,0,res);
673
674 // compute error
675 double abs_err;
676 double func_norm;
677 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
678 double err = abs_err / func_norm;
679 CuAssertDblEquals(tc, 0.0, err, 1e-13);
680
681 fwrap_destroy(fw);
682 POLY_FREE(cpoly);
683 }
684
Test_legendre_approx_adapt(CuTest * tc)685 void Test_legendre_approx_adapt(CuTest * tc){
686
687 printf("Testing function: legendre_approx_adapt\n");
688
689 // function
690 struct Fwrap * fw = fwrap_create(1,"general-vec");
691 fwrap_set_fvec(fw,Sin3xTx2,NULL);
692
693 // approximation
694 double lb = -1.0, ub = 1.0;
695 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
696 ope_opts_set_start(opts,10);
697 ope_opts_set_coeffs_check(opts,4);
698 ope_opts_set_tol(opts,1e-8);
699 ope_opts_set_lb(opts,lb);
700 ope_opts_set_ub(opts,ub);
701 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
702
703 // error
704 double abs_err;
705 double func_norm;
706 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
707 double err = abs_err / func_norm;
708 CuAssertDblEquals(tc, 0.0, err, 1e-13);
709
710 POLY_FREE(cpoly);
711 ope_opts_free(opts);
712 fwrap_destroy(fw);
713 }
714
Test_legendre_approx_adapt_weird(CuTest * tc)715 void Test_legendre_approx_adapt_weird(CuTest * tc){
716
717 printf("Testing function: legendre_approx_adapt on (a,b)\n");
718
719 // function
720 struct Fwrap * fw = fwrap_create(1,"general-vec");
721 fwrap_set_fvec(fw,Sin3xTx2,NULL);
722
723 // approximation
724 double lb = -2.0, ub = -1.0;
725 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
726 ope_opts_set_start(opts,10);
727 ope_opts_set_coeffs_check(opts,4);
728 ope_opts_set_tol(opts,1e-8);
729 ope_opts_set_lb(opts,lb);
730 ope_opts_set_ub(opts,ub);
731 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
732
733 // error
734 double abs_err;
735 double func_norm;
736 compute_error_vec(lb,ub,1000,cpoly,Sin3xTx2,NULL,&abs_err,&func_norm);
737 double err = abs_err / func_norm;
738 CuAssertDblEquals(tc, 0.0, err, 1e-13);
739
740 POLY_FREE(cpoly);
741 ope_opts_free(opts);
742 fwrap_destroy(fw);
743 }
744
Test_legendre_derivative(CuTest * tc)745 void Test_legendre_derivative(CuTest * tc){
746
747 printf("Testing function: orth_poly_expansion_deriv on (a,b)\n");
748
749 // function
750 struct Fwrap * fw = fwrap_create(1,"general-vec");
751 fwrap_set_fvec(fw,Sin3xTx2,NULL);
752
753 // approximation
754 double lb = -2.0, ub = -1.0;
755 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
756 ope_opts_set_start(opts,10);
757 ope_opts_set_coeffs_check(opts,4);
758 ope_opts_set_tol(opts,1e-8);
759 ope_opts_set_lb(opts,lb);
760 ope_opts_set_ub(opts,ub);
761 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
762 opoly_t der = orth_poly_expansion_deriv(cpoly);
763
764 // error
765 double abs_err;
766 double func_norm;
767 compute_error(lb,ub,1000,der,funcderiv,NULL,&abs_err,&func_norm);
768 double err = abs_err / func_norm;
769 CuAssertDblEquals(tc, 0.0, err, 1e-13);
770
771 POLY_FREE(cpoly);
772 POLY_FREE(der);
773 ope_opts_free(opts);
774 fwrap_destroy(fw);
775
776 }
777
Test_legendre_linear(CuTest * tc)778 void Test_legendre_linear(CuTest * tc){
779
780 printf("Testing function: orth_poly_expansion_linear with legendre poly \n");
781
782 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
783 double a = 2.0, offset=3.0;
784 opoly_t poly = orth_poly_expansion_linear(a,offset,opts);
785
786 size_t N = 100;
787 double * pts = linspace(-1,1,N);
788 size_t ii;
789 for (ii = 0; ii < N; ii++){
790 double eval1 = POLY_EVAL(poly,pts[ii]);
791 double eval2 = a*pts[ii] + offset;
792 double diff= fabs(eval1-eval2);
793 CuAssertDblEquals(tc, 0.0, diff, 1e-7);
794 }
795 free(pts); pts = NULL;
796
797 POLY_FREE(poly);
798 ope_opts_free(opts);
799 }
800
Test_legendre_quadratic(CuTest * tc)801 void Test_legendre_quadratic(CuTest * tc){
802
803 printf("Testing function: orth_poly_expansion_quadratic with legendre poly \n");
804
805 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
806 double a = 2.0, offset=3.0;
807 opoly_t poly = orth_poly_expansion_quadratic(a,offset,opts);
808
809 size_t N = 100;
810 double * pts = linspace(-1,1,N);
811 size_t ii;
812 for (ii = 0; ii < N; ii++){
813 double eval1 = POLY_EVAL(poly,pts[ii]);
814 double eval2 = a*pow(pts[ii] - offset,2);
815 double diff= fabs(eval1-eval2);
816 CuAssertDblEquals(tc, 0.0, diff, 1e-7);
817 }
818 free(pts); pts = NULL;
819
820 POLY_FREE(poly);
821 ope_opts_free(opts);
822 }
823
824
Test_legendre_genorder(CuTest * tc)825 void Test_legendre_genorder(CuTest * tc){
826
827 printf("Testing function: legendre_genorder\n");
828 // function
829
830 // approximation
831 double lb = -2.0, ub = 3.0;
832 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
833 ope_opts_set_start(opts,10);
834 ope_opts_set_coeffs_check(opts,4);
835 ope_opts_set_tol(opts,1e-8);
836 ope_opts_set_lb(opts,lb);
837 ope_opts_set_ub(opts,ub);
838 for (size_t ii = 0; ii < 10; ii++){
839 opoly_t cpoly = orth_poly_expansion_genorder(ii,opts);
840
841 for (size_t jj = 0; jj < 10; jj++){
842 opoly_t cpoly2 = orth_poly_expansion_genorder(jj,opts);
843
844 double val = orth_poly_expansion_inner(cpoly,cpoly2);
845 if (ii == jj){
846 CuAssertDblEquals(tc,1.0,val,1e-15);
847 }
848 else{
849 CuAssertDblEquals(tc,0.0,val,1e-15);
850 }
851
852 POLY_FREE(cpoly2);
853 }
854 POLY_FREE(cpoly);
855 }
856
857 ope_opts_free(opts);
858 }
859
860
Test_legendre_integrate(CuTest * tc)861 void Test_legendre_integrate(CuTest * tc){
862
863 printf("Testing function: legendre_integrate\n");
864 // function
865 struct Fwrap * fw = fwrap_create(1,"general-vec");
866 fwrap_set_fvec(fw,powX2,NULL);
867
868 // approximation
869 double lb = -2.0, ub = 3.0;
870 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
871 ope_opts_set_start(opts,10);
872 ope_opts_set_coeffs_check(opts,4);
873 ope_opts_set_tol(opts,1e-8);
874 ope_opts_set_lb(opts,lb);
875 ope_opts_set_ub(opts,ub);
876 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
877
878 double intshould = (pow(ub,3) - pow(lb,3))/3;
879 double intis = legendre_integrate(cpoly);
880 CuAssertDblEquals(tc, intshould, intis, 1e-13);
881
882 POLY_FREE(cpoly);
883 ope_opts_free(opts);
884 fwrap_destroy(fw);
885 }
886
Test_legendre_integrate_weighted(CuTest * tc)887 void Test_legendre_integrate_weighted(CuTest * tc){
888
889 printf("Testing function: legendre_integrate_weighted\n");
890 // function
891 struct Fwrap * fw = fwrap_create(1,"general-vec");
892 fwrap_set_fvec(fw,x3minusx,NULL);
893
894 // approximation
895 double lb = -3.0, ub = 9.0;
896 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
897 ope_opts_set_start(opts,10);
898 ope_opts_set_coeffs_check(opts,4);
899 ope_opts_set_tol(opts,1e-8);
900 ope_opts_set_lb(opts,lb);
901 ope_opts_set_ub(opts,ub);
902 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
903
904 double ubterm = 3.0/4.0 * pow(ub,4) - 1.0/2.0 * pow(ub,2);
905 double lbterm = 3.0/4.0 * pow(lb,4) - 1.0/2.0 * pow(lb,2);
906 double intshould = (ubterm-lbterm)/(ub-lb);
907 double intis = orth_poly_expansion_integrate_weighted(cpoly);
908 CuAssertDblEquals(tc, intshould, intis, 1e-13);
909
910 POLY_FREE(cpoly);
911 ope_opts_free(opts);
912 fwrap_destroy(fw);
913 }
914
915
Test_legendre_inner(CuTest * tc)916 void Test_legendre_inner(CuTest * tc){
917
918 printf("Testing function: orth_poly_expansion_inner with legendre poly \n");
919
920 // function
921 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
922 fwrap_set_fvec(fw1,powX2,NULL);
923
924 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
925 fwrap_set_fvec(fw2,TwoPowX3,NULL);
926
927 double lb = -2.0, ub = 3.0;
928 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
929 ope_opts_set_start(opts,10);
930 ope_opts_set_coeffs_check(opts,4);
931 ope_opts_set_tol(opts,1e-10);
932 ope_opts_set_lb(opts,lb);
933 ope_opts_set_ub(opts,ub);
934
935 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
936 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
937
938 double intshould = (pow(ub,6) - pow(lb,6))/3;
939 double intis = orth_poly_expansion_inner(cpoly,cpoly2);
940 CuAssertDblEquals(tc, intshould, intis, 1e-10);
941
942 POLY_FREE(cpoly);
943 POLY_FREE(cpoly2);
944 ope_opts_free(opts);
945 fwrap_destroy(fw1);
946 fwrap_destroy(fw2);
947 }
948
Test_legendre_inner_w(CuTest * tc)949 void Test_legendre_inner_w(CuTest * tc){
950
951 printf("Testing function: orth_poly_expansion_inner with legendre poly \n");
952
953 // function
954 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
955 fwrap_set_fvec(fw1,powX2,NULL);
956
957 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
958 fwrap_set_fvec(fw2,TwoPowX3,NULL);
959
960 double lb = -2.0, ub = 3.0;
961 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
962 ope_opts_set_start(opts,10);
963 ope_opts_set_coeffs_check(opts,4);
964 ope_opts_set_tol(opts,1e-10);
965 ope_opts_set_lb(opts,lb);
966 ope_opts_set_ub(opts,ub);
967
968 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
969 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
970
971 double intshould = (pow(ub,6) - pow(lb,6))/3/5;
972 double intis = orth_poly_expansion_inner_w(cpoly,cpoly2);
973 CuAssertDblEquals(tc, intshould, intis, 1e-10);
974
975 POLY_FREE(cpoly);
976 POLY_FREE(cpoly2);
977 ope_opts_free(opts);
978 fwrap_destroy(fw1);
979 fwrap_destroy(fw2);
980 }
981
Test_legendre_norm_w(CuTest * tc)982 void Test_legendre_norm_w(CuTest * tc){
983
984 printf("Testing function: orth_poly_expansion_norm_w with legendre poly\n");
985
986 // function
987 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
988 fwrap_set_fvec(fw1,powX2,NULL);
989
990 double lb = -2, ub = 3.0;
991 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
992 ope_opts_set_start(opts,10);
993 ope_opts_set_coeffs_check(opts,4);
994 ope_opts_set_tol(opts,1e-8);
995 ope_opts_set_lb(opts,lb);
996 ope_opts_set_ub(opts,ub);
997
998 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
999
1000 double intshould = sqrt((pow(ub,5) - pow(lb,5))/5/5);
1001 double intis = orth_poly_expansion_norm_w(cpoly);
1002 CuAssertDblEquals(tc, sqrt(intshould), intis, 1e-10);
1003
1004 POLY_FREE(cpoly);
1005 ope_opts_free(opts);
1006 fwrap_destroy(fw1);
1007 }
1008
Test_legendre_norm(CuTest * tc)1009 void Test_legendre_norm(CuTest * tc){
1010
1011 printf("Testing function: orth_poly_expansion_norm with legendre poly\n");
1012
1013 // function
1014 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1015 fwrap_set_fvec(fw1,powX2,NULL);
1016
1017 double lb = -2, ub = 3.0;
1018 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
1019 ope_opts_set_start(opts,10);
1020 ope_opts_set_coeffs_check(opts,4);
1021 ope_opts_set_tol(opts,1e-8);
1022 ope_opts_set_lb(opts,lb);
1023 ope_opts_set_ub(opts,ub);
1024
1025 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
1026
1027 double intshould = (pow(ub,5) - pow(lb,5))/5.0;
1028 double intis = orth_poly_expansion_norm(cpoly);
1029 CuAssertDblEquals(tc, sqrt(intshould), intis, 1e-10);
1030
1031 POLY_FREE(cpoly);
1032 ope_opts_free(opts);
1033 fwrap_destroy(fw1);
1034 }
1035
Test_legendre_product(CuTest * tc)1036 void Test_legendre_product(CuTest * tc){
1037
1038 printf("Testing function: orth_poly_expansion_product with legendre poly \n");
1039
1040 // function
1041 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1042 fwrap_set_fvec(fw1,powX2,NULL);
1043
1044 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1045 fwrap_set_fvec(fw2,TwoPowX3,NULL);
1046
1047 // approximation
1048 double lb = -3.0, ub = 2.0;
1049 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
1050 ope_opts_set_start(opts,10);
1051 ope_opts_set_coeffs_check(opts,4);
1052 ope_opts_set_tol(opts,1e-10);
1053 ope_opts_set_lb(opts,lb);
1054 ope_opts_set_ub(opts,ub);
1055
1056 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
1057 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
1058 opoly_t cpoly3 = orth_poly_expansion_prod(cpoly,cpoly2);
1059
1060 size_t N = 100;
1061 double * pts = linspace(lb,ub,N);
1062 size_t ii;
1063 for (ii = 0; ii < N; ii++){
1064 double eval1 = POLY_EVAL(cpoly3,pts[ii]);
1065 double eval2 = POLY_EVAL(cpoly,pts[ii]) *
1066 POLY_EVAL(cpoly2,pts[ii]);
1067 double diff= fabs(eval1-eval2);
1068 CuAssertDblEquals(tc, 0.0, diff, 1e-10);
1069 }
1070 free(pts); pts = NULL;
1071
1072 POLY_FREE(cpoly);
1073 POLY_FREE(cpoly2);
1074 POLY_FREE(cpoly3);
1075 ope_opts_free(opts);
1076 fwrap_destroy(fw1);
1077 fwrap_destroy(fw2);
1078 }
1079
Test_legendre_axpy(CuTest * tc)1080 void Test_legendre_axpy(CuTest * tc){
1081
1082 printf("Testing function: orth_poly_expansion_axpy with legendre poly \n");
1083
1084 // function
1085 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1086 fwrap_set_fvec(fw1,powX2,NULL);
1087
1088 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1089 fwrap_set_fvec(fw2,TwoPowX3,NULL);
1090
1091 // approximation
1092 double lb = -3.0, ub = 2.0;
1093 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
1094 ope_opts_set_start(opts,10);
1095 ope_opts_set_coeffs_check(opts,4);
1096 ope_opts_set_tol(opts,1e-10);
1097 ope_opts_set_lb(opts,lb);
1098 ope_opts_set_ub(opts,ub);
1099
1100 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw1);
1101 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
1102
1103
1104 int success = orth_poly_expansion_axpy(2.0,cpoly2,cpoly);
1105 CuAssertIntEquals(tc,0,success);
1106
1107 size_t N = 100;
1108 double * pts = linspace(lb,ub,N);
1109 size_t ii;
1110 for (ii = 0; ii < N; ii++){
1111 double eval1 = POLY_EVAL(cpoly,pts[ii]);
1112 double eval2a,eval2b;
1113 TwoPowX3(1,pts+ii,&eval2a,NULL);
1114 powX2(1,pts+ii,&eval2b,NULL);
1115 double eval2 = 2.0 * eval2a + eval2b;
1116 double diff= fabs(eval1-eval2);
1117 CuAssertDblEquals(tc, 0.0, diff, 1e-10);
1118 }
1119 free(pts); pts = NULL;
1120
1121 POLY_FREE(cpoly);
1122 POLY_FREE(cpoly2);
1123 ope_opts_free(opts);
1124 fwrap_destroy(fw1);
1125 fwrap_destroy(fw2);
1126 }
1127
Test_legendre_orth_poly_expansion_create_with_params_and_grad(CuTest * tc)1128 void Test_legendre_orth_poly_expansion_create_with_params_and_grad(CuTest * tc){
1129
1130 printf("Testing function: orth_poly_expansion_create_with_params and grad with legendre poly\n");
1131
1132 size_t nparams = 10;
1133 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
1134 0.6,-0.9,-.2,0.04,1.2};
1135 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
1136 struct OrthPolyExpansion * ope = NULL;
1137 ope = orth_poly_expansion_create_with_params(opts,nparams,params);
1138
1139
1140 double grad[10];
1141 double xloc = 0.4;
1142 int res =
1143 orth_poly_expansion_param_grad_eval(ope,1,&xloc,grad);
1144 CuAssertIntEquals(tc,0,res);
1145
1146
1147 // numerical derivative
1148 struct OrthPolyExpansion * ope1 = NULL;
1149 struct OrthPolyExpansion * ope2 = NULL;
1150
1151 size_t dim = nparams;
1152 double * x1 = calloc_double(dim);
1153 double * x2 = calloc_double(dim);
1154 for (size_t ii = 0; ii < dim; ii++){
1155 x1[ii] = params[ii];
1156 x2[ii] = params[ii];
1157 }
1158
1159 double diff = 0.0;
1160 double v1,v2;
1161 double norm = 0.0;
1162 double eps = 1e-8;
1163 for (size_t ii = 0; ii < dim; ii++){
1164 x1[ii] += eps;
1165 x2[ii] -= eps;
1166 ope1 = orth_poly_expansion_create_with_params(opts,nparams,x1);
1167 v1 = orth_poly_expansion_eval(ope1,xloc);
1168
1169 ope2 = orth_poly_expansion_create_with_params(opts,nparams,x2);
1170 v2 = orth_poly_expansion_eval(ope2,xloc);
1171
1172 double diff_iter = pow( (v1-v2)/2.0/eps - grad[ii], 2 );
1173 /* printf("current diff = %G\n",diff_iter); */
1174 /* printf("\t norm = %G\n",grad[ii]); */
1175 diff += diff_iter;
1176 norm += pow( (v1-v2)/2.0/eps,2);
1177
1178 x1[ii] -= eps;
1179 x2[ii] += eps;
1180
1181 orth_poly_expansion_free(ope1); ope1 = NULL;
1182 orth_poly_expansion_free(ope2); ope2 = NULL;
1183 }
1184 if (norm > 1){
1185 diff /= norm;
1186 }
1187 free(x1); x1 = NULL;
1188 free(x2); x2 = NULL;
1189 CuAssertDblEquals(tc,0.0,diff,1e-7);
1190
1191 ope_opts_free(opts); opts = NULL;
1192 orth_poly_expansion_free(ope); ope = NULL;
1193 }
1194
1195
Test_legendre_orth_poly_expansion_squared_norm_param_grad(CuTest * tc)1196 void Test_legendre_orth_poly_expansion_squared_norm_param_grad(CuTest * tc){
1197
1198 printf("Testing function: orth_poly_expansion_squared_norm_param_grad with legendre poly\n");
1199
1200 size_t nparams = 10;
1201 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
1202 0.6,-0.9,-.2,0.04,1.2};
1203 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
1204 ope_opts_set_lb(opts,-3.0);
1205 ope_opts_set_ub(opts,2.0);
1206 struct OrthPolyExpansion * ope = NULL;
1207 ope = orth_poly_expansion_create_with_params(opts,nparams,params);
1208
1209
1210 double grad[10] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
1211 double w = 0.5;
1212 int res =
1213 orth_poly_expansion_squared_norm_param_grad(ope,w,grad);
1214 CuAssertIntEquals(tc,0,res);
1215
1216
1217 // numerical derivative
1218 struct OrthPolyExpansion * ope1 = NULL;
1219 struct OrthPolyExpansion * ope2 = NULL;
1220
1221 size_t dim = nparams;
1222 double * x1 = calloc_double(dim);
1223 double * x2 = calloc_double(dim);
1224 for (size_t ii = 0; ii < dim; ii++){
1225 x1[ii] = params[ii];
1226 x2[ii] = params[ii];
1227 }
1228
1229 double diff = 0.0;
1230 double v1,v2;
1231 double norm = 0.0;
1232 double eps = 1e-8;
1233 for (size_t ii = 0; ii < dim; ii++){
1234 x1[ii] += eps;
1235 x2[ii] -= eps;
1236 ope1 = orth_poly_expansion_create_with_params(opts,nparams,x1);
1237 v1 = orth_poly_expansion_inner(ope1,ope1)*w;
1238
1239 ope2 = orth_poly_expansion_create_with_params(opts,nparams,x2);
1240 v2 = orth_poly_expansion_inner(ope2,ope2)*w;
1241
1242 double diff_iter = pow( (v1-v2)/2.0/eps - grad[ii], 2 );
1243 /* printf("current diff = %G\n",diff_iter); */
1244 /* printf("\t norm = %G\n",grad[ii]); */
1245 diff += diff_iter;
1246 norm += pow( (v1-v2)/2.0/eps,2);
1247
1248 x1[ii] -= eps;
1249 x2[ii] += eps;
1250
1251 orth_poly_expansion_free(ope1); ope1 = NULL;
1252 orth_poly_expansion_free(ope2); ope2 = NULL;
1253 }
1254 if (norm > 1){
1255 diff /= norm;
1256 }
1257 free(x1); x1 = NULL;
1258 free(x2); x2 = NULL;
1259 CuAssertDblEquals(tc,0.0,diff,1e-7);
1260
1261 ope_opts_free(opts); opts = NULL;
1262 orth_poly_expansion_free(ope); ope = NULL;
1263 }
1264
1265
LegGetSuite()1266 CuSuite * LegGetSuite(){
1267
1268 CuSuite * suite = CuSuiteNew();
1269 SUITE_ADD_TEST(suite, Test_legendre_approx);
1270 SUITE_ADD_TEST(suite, Test_legendre_approx_gaussbump);
1271 SUITE_ADD_TEST(suite, Test_legendre_approx_nonnormal);
1272 SUITE_ADD_TEST(suite, Test_legendre_approx_adapt);
1273 SUITE_ADD_TEST(suite, Test_legendre_approx_adapt_weird);
1274 SUITE_ADD_TEST(suite, Test_legendre_linear);
1275 SUITE_ADD_TEST(suite, Test_legendre_quadratic);
1276 SUITE_ADD_TEST(suite, Test_legendre_genorder);
1277 SUITE_ADD_TEST(suite, Test_legendre_derivative);
1278 SUITE_ADD_TEST(suite, Test_legendre_integrate);
1279 SUITE_ADD_TEST(suite, Test_legendre_integrate_weighted);
1280 SUITE_ADD_TEST(suite, Test_legendre_inner);
1281 SUITE_ADD_TEST(suite, Test_legendre_inner_w);
1282 SUITE_ADD_TEST(suite, Test_legendre_norm_w);
1283 SUITE_ADD_TEST(suite, Test_legendre_norm);
1284 SUITE_ADD_TEST(suite, Test_legendre_product);
1285 SUITE_ADD_TEST(suite, Test_legendre_axpy);
1286 SUITE_ADD_TEST(suite, Test_legendre_orth_poly_expansion_create_with_params_and_grad);
1287 SUITE_ADD_TEST(suite, Test_legendre_orth_poly_expansion_squared_norm_param_grad);
1288
1289 return suite;
1290 }
1291
1292
Test_fourier_approx(CuTest * tc)1293 void Test_fourier_approx(CuTest * tc){
1294
1295 printf("Testing function: fourier_approx\n");
1296
1297 // function
1298 struct Fwrap * fw = fwrap_create(1,"general-vec");
1299 fwrap_set_fvec(fw,gaussbump,NULL);
1300
1301 // approximation
1302 size_t N = 2*2*2*2*2*2*2*2 + 1;
1303 double lb=-5.0,ub=5.0;
1304 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1305 ope_opts_set_lb(opts, lb);
1306 ope_opts_set_ub(opts, ub);
1307
1308 struct OrthPolyExpansion * fourier =
1309 orth_poly_expansion_init_from_opts(opts, N);
1310
1311 int res = orth_poly_expansion_approx_vec(fourier,fw,NULL);
1312 CuAssertIntEquals(tc,0,res);
1313
1314 // compute error
1315 double abs_err;
1316 double func_norm;
1317 compute_error_vec(lb,ub,1000,fourier,gaussbump,NULL,&abs_err,&func_norm);
1318 double err = abs_err / func_norm;
1319 CuAssertDblEquals(tc, 0.0, err, 1e-13);
1320
1321 fwrap_destroy(fw);
1322 ope_opts_free(opts); opts = NULL;
1323 POLY_FREE(fourier);
1324 }
1325
Test_fourier_orth_basis(CuTest * tc)1326 void Test_fourier_orth_basis(CuTest * tc){
1327
1328 printf("Testing function: fourier_orth_basis\n");
1329 // function
1330
1331 // approximation
1332 double lb = -2.0, ub = 3.0;
1333 size_t n = 32;
1334 size_t N = n+1;
1335 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1336 ope_opts_set_start(opts,N);
1337 ope_opts_set_lb(opts,lb);
1338 ope_opts_set_ub(opts,ub);
1339 struct OrthPolyExpansion *f[32];
1340
1341 orth_poly_expansion_orth_basis(32, f, opts);
1342
1343 for (size_t ii = 0; ii < n; ii++){
1344 for (size_t jj = 0; jj < n; jj++){
1345 double inner = orth_poly_expansion_inner(f[ii], f[jj]);
1346 if (ii == jj){
1347 CuAssertDblEquals(tc, 1, inner, 1e-10);
1348 }
1349 else{
1350 CuAssertDblEquals(tc, 0, inner, 1e-10);
1351 }
1352 /* printf("ii = %zu, jj = %zu\n", ii, jj); */
1353 /* printf("\t inner=%3.5G\n", inner); */
1354 }
1355 }
1356 for (size_t ii = 0; ii < n; ii++){
1357 orth_poly_expansion_free(f[ii]); f[ii] = NULL;
1358 }
1359 ope_opts_free(opts);
1360 }
1361
1362
Test_fourier_derivative(CuTest * tc)1363 void Test_fourier_derivative(CuTest * tc){
1364
1365 printf("Testing function: orth_poly_expansion_deriv: fourier\n");
1366
1367 // function
1368 struct Fwrap * fw = fwrap_create(1,"general-vec");
1369 fwrap_set_fvec(fw,gaussbump,NULL);
1370
1371 double lb = -8.0, ub = 8.0;
1372 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1373 ope_opts_set_lb(opts, lb);
1374 ope_opts_set_ub(opts, ub);
1375 size_t N = 512+1;
1376 struct OrthPolyExpansion * cpoly =
1377 orth_poly_expansion_init_from_opts(opts, N);
1378 int res = orth_poly_expansion_approx_vec(cpoly, fw, NULL);
1379 CuAssertIntEquals(tc,0,res);
1380
1381 opoly_t der = orth_poly_expansion_deriv(cpoly);
1382
1383 // error
1384 double abs_err;
1385 double func_norm;
1386 compute_error_vec(lb,ub,1000,der,gaussbump_deriv,NULL,&abs_err,&func_norm);
1387 double err = abs_err / func_norm;
1388 CuAssertDblEquals(tc, 0.0, err, 1e-13);
1389
1390 POLY_FREE(cpoly);
1391 POLY_FREE(der);
1392 ope_opts_free(opts);
1393 fwrap_destroy(fw);
1394
1395 }
1396
Test_fourier_dderivative(CuTest * tc)1397 void Test_fourier_dderivative(CuTest * tc){
1398
1399 printf("Testing function: orth_poly_expansion_dderiv: fourier\n");
1400
1401 // function
1402 struct Fwrap * fw = fwrap_create(1,"general-vec");
1403 fwrap_set_fvec(fw,gaussbump,NULL);
1404
1405 double lb = -8.0, ub = 8.0;
1406 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1407 ope_opts_set_lb(opts, lb);
1408 ope_opts_set_ub(opts, ub);
1409 size_t N = 512+1;
1410 struct OrthPolyExpansion * cpoly =
1411 orth_poly_expansion_init_from_opts(opts, N);
1412 int res = orth_poly_expansion_approx_vec(cpoly, fw, NULL);
1413 CuAssertIntEquals(tc,0,res);
1414
1415 opoly_t der = orth_poly_expansion_dderiv(cpoly);
1416
1417 // error
1418 double abs_err;
1419 double func_norm;
1420 compute_error_vec(lb,ub,1000,der,gaussbump_dderiv,NULL,&abs_err,&func_norm);
1421 double err = abs_err / func_norm;
1422 CuAssertDblEquals(tc, 0.0, err, 1e-13);
1423
1424 POLY_FREE(cpoly);
1425 POLY_FREE(der);
1426 ope_opts_free(opts);
1427 fwrap_destroy(fw);
1428
1429 }
1430
1431
Test_fourier_integrate(CuTest * tc)1432 void Test_fourier_integrate(CuTest * tc){
1433
1434 printf("Testing function: fourier_integrate\n");
1435 // function
1436 struct Fwrap * fw = fwrap_create(1,"general-vec");
1437 fwrap_set_fvec(fw,powX2,NULL);
1438
1439 // approximation
1440 double lb = -8.0, ub = 8.0;
1441 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1442 ope_opts_set_lb(opts, lb);
1443 ope_opts_set_ub(opts, ub);
1444 size_t N = 512+1;
1445 struct OrthPolyExpansion * cpoly =
1446 orth_poly_expansion_init_from_opts(opts, N);
1447 int res = orth_poly_expansion_approx_vec(cpoly, fw, NULL);
1448 CuAssertIntEquals(tc,0,res);
1449
1450 double intshould = (pow(ub,3) - pow(lb,3))/3;
1451 double intis = orth_poly_expansion_integrate(cpoly);
1452 CuAssertDblEquals(tc, intshould, intis, 1e-3);
1453
1454 POLY_FREE(cpoly);
1455 ope_opts_free(opts);
1456 fwrap_destroy(fw);
1457 }
1458
1459
Test_fourier_inner(CuTest * tc)1460 void Test_fourier_inner(CuTest * tc){
1461
1462 printf("Testing function: orth_poly_expansion_inner with fourier poly \n");
1463
1464 // function
1465 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1466 fwrap_set_fvec(fw1,powX2,NULL);
1467
1468 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1469 fwrap_set_fvec(fw2,powX2,NULL);
1470
1471 double lb = -8.0, ub = 8.0;
1472 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1473 ope_opts_set_lb(opts, lb);
1474 ope_opts_set_ub(opts, ub);
1475 size_t N = 512+1;
1476
1477 struct OrthPolyExpansion * cpoly = orth_poly_expansion_init_from_opts(opts, N);
1478 orth_poly_expansion_approx_vec(cpoly, fw1, NULL);
1479
1480 struct OrthPolyExpansion * cpoly2 = orth_poly_expansion_init_from_opts(opts, N);
1481 orth_poly_expansion_approx_vec(cpoly2, fw2, NULL);
1482
1483 double intshould = (pow(ub,5) - pow(lb,5))/5;
1484 double intis = orth_poly_expansion_inner(cpoly,cpoly2);
1485 double relerr = fabs(intshould - intis) / fabs(intshould);
1486 CuAssertDblEquals(tc, 0, relerr, 1e-5);
1487
1488 POLY_FREE(cpoly);
1489 POLY_FREE(cpoly2);
1490 ope_opts_free(opts);
1491 fwrap_destroy(fw1);
1492 fwrap_destroy(fw2);
1493 }
1494
Test_fourier_product(CuTest * tc)1495 void Test_fourier_product(CuTest * tc){
1496
1497 printf("Testing function: orth_poly_expansion_product with fourier poly \n");
1498
1499 // function
1500 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1501 fwrap_set_fvec(fw1,powX2,NULL);
1502
1503 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1504 fwrap_set_fvec(fw2,TwoPowX3,NULL);
1505
1506 // approximation
1507 double lb = -3.0, ub = 2.0;
1508 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1509 ope_opts_set_lb(opts,lb);
1510 ope_opts_set_ub(opts,ub);
1511
1512 size_t N = 512+1;
1513 struct OrthPolyExpansion * cpoly = orth_poly_expansion_init_from_opts(opts, N);
1514 orth_poly_expansion_approx_vec(cpoly, fw1, NULL);
1515
1516 struct OrthPolyExpansion * cpoly2 = orth_poly_expansion_init_from_opts(opts, N);
1517 orth_poly_expansion_approx_vec(cpoly, fw2, NULL);
1518
1519 struct OrthPolyExpansion * cpoly3 = orth_poly_expansion_prod(cpoly, cpoly2);
1520
1521 N = 100;
1522 double * pts = linspace(lb,ub,N);
1523 size_t ii;
1524 for (ii = 0; ii < N; ii++){
1525 double eval1 = POLY_EVAL(cpoly3,pts[ii]);
1526 double eval2 = POLY_EVAL(cpoly,pts[ii]) *
1527 POLY_EVAL(cpoly2,pts[ii]);
1528 double diff= fabs(eval1-eval2);
1529 CuAssertDblEquals(tc, 0.0, diff, 1e-10);
1530 }
1531 free(pts); pts = NULL;
1532
1533 POLY_FREE(cpoly);
1534 POLY_FREE(cpoly2);
1535 POLY_FREE(cpoly3);
1536 ope_opts_free(opts);
1537 fwrap_destroy(fw1);
1538 fwrap_destroy(fw2);
1539 }
1540
1541
Test_fourier_axpy(CuTest * tc)1542 void Test_fourier_axpy(CuTest * tc){
1543
1544 printf("Testing function: orth_poly_expansion_axpy with fourier poly \n");
1545
1546 // function
1547 struct Fwrap * fw1 = fwrap_create(1,"general-vec");
1548 fwrap_set_fvec(fw1,powX2,NULL);
1549
1550 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1551 /* fwrap_set_fvec(fw2,TwoPowX3,NULL); */
1552 fwrap_set_fvec(fw2,powX2,NULL);
1553
1554 // approximation
1555 double lb = -8.0, ub = 8.0;
1556 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
1557 ope_opts_set_lb(opts,lb);
1558 ope_opts_set_ub(opts,ub);
1559
1560 size_t N = 4096+1;
1561
1562 struct OrthPolyExpansion * cpoly = orth_poly_expansion_init_from_opts(opts, N);
1563 orth_poly_expansion_approx_vec(cpoly, fw1, NULL);
1564
1565 struct OrthPolyExpansion * cpoly2 = orth_poly_expansion_init_from_opts(opts, N);
1566 orth_poly_expansion_approx_vec(cpoly2, fw2, NULL);
1567
1568 int success = orth_poly_expansion_axpy(2.0,cpoly2,cpoly);
1569 CuAssertIntEquals(tc,0,success);
1570
1571 N = 100;
1572 double * pts = linspace(lb,ub,N);
1573 size_t ii;
1574 for (ii = 0; ii < N; ii++){
1575 double eval1 = POLY_EVAL(cpoly,pts[ii]);
1576 double eval2a,eval2b;
1577 /* TwoPowX3(1,pts+ii,&eval2a,NULL); */
1578 powX2(1,pts+ii,&eval2a,NULL);
1579 powX2(1,pts+ii,&eval2b,NULL);
1580 double eval2 = 2.0 * eval2a + eval2b;
1581 double diff= fabs(eval1-eval2);
1582 /* printf("eval2 = %3.15G %3.15G\n", eval2, diff); */
1583 CuAssertDblEquals(tc, 0.0, diff/fabs(eval2), 1e-1);
1584 }
1585 free(pts); pts = NULL;
1586
1587 POLY_FREE(cpoly);
1588 POLY_FREE(cpoly2);
1589 ope_opts_free(opts);
1590 fwrap_destroy(fw1);
1591 fwrap_destroy(fw2);
1592 }
1593
FourierGetSuite()1594 CuSuite * FourierGetSuite(){
1595 CuSuite * suite = CuSuiteNew();
1596 SUITE_ADD_TEST(suite, Test_fourier_approx);
1597 SUITE_ADD_TEST(suite, Test_fourier_orth_basis);
1598 SUITE_ADD_TEST(suite, Test_fourier_derivative);
1599 SUITE_ADD_TEST(suite, Test_fourier_dderivative);
1600 SUITE_ADD_TEST(suite, Test_fourier_integrate);
1601 SUITE_ADD_TEST(suite, Test_fourier_inner);
1602 SUITE_ADD_TEST(suite, Test_fourier_product);
1603 SUITE_ADD_TEST(suite, Test_fourier_axpy);
1604
1605 return suite;
1606 }
1607
fherm1(size_t N,const double * x,double * out,void * arg)1608 int fherm1(size_t N, const double * x, double * out, void * arg)
1609 {
1610 (void)(arg);
1611 for (size_t ii = 0; ii < N; ii++){
1612 out[ii]= x[ii] + x[ii]*x[ii];
1613 }
1614 return 0;
1615 }
1616
Test_hermite_approx(CuTest * tc)1617 void Test_hermite_approx(CuTest * tc){
1618
1619 printf("Testing function: hermite_approx\n");
1620
1621 // function
1622 struct Fwrap * fw = fwrap_create(1,"general-vec");
1623 fwrap_set_fvec(fw,fherm1,NULL);
1624
1625 // approximation
1626 size_t N = 20;
1627 opoly_t cpoly = orth_poly_expansion_init(HERMITE,N,-DBL_MAX,DBL_MAX);
1628 int res = orth_poly_expansion_approx_vec(cpoly,fw,NULL);
1629 CuAssertIntEquals(tc,0,res);
1630
1631 // compute error
1632 double abs_err;
1633 double func_norm;
1634 compute_error_vec(-1.0,1.0,1000,cpoly,fherm1,NULL,&abs_err,&func_norm);
1635 double err = abs_err / func_norm;
1636 CuAssertDblEquals(tc, 0.0, err, 1e-13);
1637
1638 fwrap_destroy(fw);
1639 POLY_FREE(cpoly);
1640 }
1641
fherm2(size_t N,const double * x,double * out,void * arg)1642 int fherm2(size_t N, const double * x, double * out, void * arg)
1643 {
1644 (void)(arg);
1645 for (size_t ii = 0; ii < N; ii++){
1646 out[ii]= sin(2.0*x[ii]);
1647 }
1648 return 0;
1649 }
1650
1651
Test_hermite_approx_adapt(CuTest * tc)1652 void Test_hermite_approx_adapt(CuTest * tc)
1653 {
1654 printf("Testing function: hermite_approx_adapt\n");
1655
1656 // function
1657 struct Fwrap * fw = fwrap_create(1,"general-vec");
1658 fwrap_set_fvec(fw,fherm2,NULL);
1659
1660 // approximation
1661 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1662 ope_opts_set_start(opts,10);
1663 ope_opts_set_coeffs_check(opts,4);
1664 ope_opts_set_tol(opts,1e-15);
1665 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1666
1667 // error
1668 double abs_err;
1669 double func_norm;
1670 compute_error_vec(-1.0,1.0,1000,cpoly,fherm2,NULL,&abs_err,&func_norm);
1671 double err = abs_err / func_norm;
1672 CuAssertDblEquals(tc, 0.0, err, 1e-10);
1673
1674 POLY_FREE(cpoly);
1675 ope_opts_free(opts);
1676 fwrap_destroy(fw);
1677 }
1678
Test_hermite_derivative(CuTest * tc)1679 void Test_hermite_derivative(CuTest * tc){
1680
1681 printf("Testing function: orth_poly_expansion_deriv with hermite poly\n");
1682
1683 // function
1684 struct Fwrap * fw = fwrap_create(1,"general-vec");
1685 fwrap_set_fvec(fw,Sin3xTx2,NULL);
1686
1687 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1688 ope_opts_set_start(opts,5);
1689 ope_opts_set_coeffs_check(opts,4);
1690 ope_opts_set_tol(opts,1e-15);
1691 ope_opts_set_maxnum(opts,50);
1692 ope_opts_set_mean_and_std(opts,0.0,1.0);
1693 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1694 opoly_t der = orth_poly_expansion_deriv(cpoly);
1695
1696 // error
1697 double abs_err;
1698 double func_norm;
1699 compute_error(-2,2,1000,der,funcderiv,NULL,&abs_err,&func_norm);
1700 double err = abs_err / func_norm;
1701 CuAssertDblEquals(tc, 0.0, err, 1e-4);
1702
1703 POLY_FREE(cpoly);
1704 POLY_FREE(der);
1705 ope_opts_free(opts);
1706 fwrap_destroy(fw);
1707 }
1708
fherm3(size_t N,const double * x,double * out,void * arg)1709 int fherm3(size_t N, const double * x, double * out, void * arg)
1710 {
1711 (void)(arg);
1712 for (size_t ii = 0; ii < N; ii++){
1713 out[ii]= sin(2*x[ii]+3) + 3*pow(x[ii],3);
1714 }
1715 return 0;
1716 }
1717
Test_hermite_integrate(CuTest * tc)1718 void Test_hermite_integrate(CuTest * tc){
1719
1720 printf("Testing function: hermite_integrate\n");
1721
1722 // function
1723 struct Fwrap * fw = fwrap_create(1,"general-vec");
1724 fwrap_set_fvec(fw,fherm3,NULL);
1725
1726 // approximation
1727 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1728 ope_opts_set_start(opts,10);
1729 ope_opts_set_coeffs_check(opts,4);
1730 ope_opts_set_tol(opts,1e-8);
1731 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1732
1733 double intshould = sin(3)/exp(2);
1734 double intis = hermite_integrate(cpoly);
1735 CuAssertDblEquals(tc, intshould, intis, 1e-13);
1736
1737 POLY_FREE(cpoly);
1738 ope_opts_free(opts);
1739 fwrap_destroy(fw);
1740 }
1741
fherm_lin(size_t N,const double * x,double * out,void * arg)1742 int fherm_lin(size_t N, const double * x, double * out, void * arg)
1743 {
1744 (void)(arg);
1745 for (size_t ii = 0; ii < N; ii++){
1746 out[ii] = 3.0*x[ii] - 0.2;
1747 }
1748 return 0;
1749 }
1750
Test_hermite_non_std_normal(CuTest * tc)1751 void Test_hermite_non_std_normal(CuTest * tc)
1752 {
1753
1754 // function
1755 struct Fwrap * fw = fwrap_create(1,"general-vec");
1756 fwrap_set_fvec(fw,fherm_lin,NULL);
1757
1758 // approximation
1759 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1760 ope_opts_set_start(opts,10);
1761 ope_opts_set_coeffs_check(opts,4);
1762 ope_opts_set_tol(opts,1e-8);
1763 ope_opts_set_mean_and_std(opts,2.0,0.05);
1764 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1765
1766
1767 double mean_should = -0.2 + 3.0*2.0;
1768 double mean_is = hermite_integrate(cpoly);
1769 double second_moment = orth_poly_expansion_inner(cpoly,cpoly);
1770
1771 opoly_t cpoly_prod = orth_poly_expansion_prod(cpoly,cpoly);
1772 double second_moment2 = hermite_integrate(cpoly_prod);
1773
1774 double var = second_moment - mean_is * mean_is;
1775 double var2 = second_moment2 - mean_is * mean_is;
1776
1777 double var_should = 9.0*0.05*0.05;
1778
1779 CuAssertDblEquals(tc, mean_should, mean_is, 1e-13);
1780 CuAssertDblEquals(tc, var_should, var, 1e-13);
1781 CuAssertDblEquals(tc, var_should, var2, 1e-13);
1782
1783 POLY_FREE(cpoly_prod);
1784 POLY_FREE(cpoly);
1785 ope_opts_free(opts);
1786 fwrap_destroy(fw);
1787 }
1788
1789
fherm4(size_t N,const double * x,double * out,void * arg)1790 int fherm4(size_t N, const double * x, double * out, void * arg)
1791 {
1792 (void)(arg);
1793 for (size_t ii = 0; ii < N; ii++){
1794 out[ii] = sin(2.0*x[ii]+3.0);
1795 }
1796 return 0;
1797 }
1798
fherm5(size_t N,const double * x,double * out,void * arg)1799 int fherm5(size_t N, const double * x, double * out, void * arg)
1800 {
1801 (void)(arg);
1802 for (size_t ii = 0; ii < N; ii++){
1803 out[ii] = 3*pow(x[ii],3);
1804 }
1805 return 0;
1806 }
1807
Test_hermite_inner(CuTest * tc)1808 void Test_hermite_inner(CuTest * tc){
1809
1810 printf("Testing function: orth_poly_expansion_inner with hermite poly \n");
1811
1812 // function
1813 struct Fwrap * fw = fwrap_create(1,"general-vec");
1814 fwrap_set_fvec(fw,fherm4,NULL);
1815 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1816 fwrap_set_fvec(fw2,fherm5,NULL);
1817
1818 // approximation
1819 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1820 ope_opts_set_start(opts,10);
1821 ope_opts_set_coeffs_check(opts,4);
1822 ope_opts_set_tol(opts,1e-15);
1823 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1824 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
1825
1826 double intshould = -6.0*cos(3)/exp(2.0);
1827 double intis = orth_poly_expansion_inner(cpoly,cpoly2);
1828 CuAssertDblEquals(tc, intshould, intis, 1e-10);
1829
1830 POLY_FREE(cpoly);
1831 POLY_FREE(cpoly2);
1832 ope_opts_free(opts);
1833 fwrap_destroy(fw);
1834 fwrap_destroy(fw2);
1835 }
1836
fherm6(size_t N,const double * x,double * out,void * arg)1837 int fherm6(size_t N, const double * x, double * out, void * arg)
1838 {
1839 (void)(arg);
1840 for (size_t ii = 0; ii < N; ii++){
1841 out[ii] = pow(x[ii],2)*sin(x[ii]+0.5);
1842 }
1843 return 0;
1844 }
1845
Test_hermite_norm_w(CuTest * tc)1846 void Test_hermite_norm_w(CuTest * tc){
1847
1848 printf("Testing function: orth_poly_expansion_norm_w with hermite poly\n");
1849
1850 struct Fwrap * fw = fwrap_create(1,"general-vec");
1851 fwrap_set_fvec(fw,fherm6,NULL);
1852
1853 // approximation
1854 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1855 ope_opts_set_start(opts,10);
1856 ope_opts_set_coeffs_check(opts,4);
1857 ope_opts_set_tol(opts,1e-15);
1858 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1859
1860 double intshould = sqrt(5*cos(1)/2/exp(2) + 3.0/2.0);
1861 double intis = orth_poly_expansion_norm_w(cpoly);
1862 CuAssertDblEquals(tc, sqrt(intshould), intis, 1e-13);
1863
1864 POLY_FREE(cpoly);
1865 ope_opts_free(opts);
1866 fwrap_destroy(fw);
1867 }
1868
Test_hermite_norm(CuTest * tc)1869 void Test_hermite_norm(CuTest * tc){
1870
1871 printf("Testing function: orth_poly_expansion_norm with hermite poly\n");
1872
1873 struct Fwrap * fw = fwrap_create(1,"general-vec");
1874 fwrap_set_fvec(fw,fherm6,NULL);
1875
1876 // approximation
1877 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1878 ope_opts_set_start(opts,10);
1879 ope_opts_set_coeffs_check(opts,4);
1880 ope_opts_set_tol(opts,1e-15);
1881 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1882
1883 double intshould = sqrt(5*cos(1)/2/exp(2) + 3.0/2.0);
1884 double intis = orth_poly_expansion_norm_w(cpoly);
1885 CuAssertDblEquals(tc, sqrt(intshould), intis, 1e-13);
1886
1887 POLY_FREE(cpoly);
1888 ope_opts_free(opts);
1889 fwrap_destroy(fw);
1890 }
1891
1892
fherm7(size_t N,const double * x,double * out,void * arg)1893 int fherm7(size_t N, const double * x, double * out, void * arg)
1894 {
1895 (void)(arg);
1896 for (size_t ii = 0; ii < N; ii++){
1897 out[ii] = pow(x[ii],2);
1898 }
1899 return 0;
1900 }
1901
fherm8(size_t N,const double * x,double * out,void * arg)1902 int fherm8(size_t N, const double * x, double * out, void * arg)
1903 {
1904 (void)(arg);
1905 for (size_t ii = 0; ii < N; ii++){
1906 out[ii] = 2.0 + 3.0*pow(x[ii],5);
1907 }
1908 return 0;
1909 }
1910
Test_hermite_product(CuTest * tc)1911 void Test_hermite_product(CuTest * tc){
1912
1913 printf("Testing function: orth_poly_expansion_product with hermite poly \n");
1914
1915 struct Fwrap * fw = fwrap_create(1,"general-vec");
1916 fwrap_set_fvec(fw,fherm7,NULL);
1917
1918 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1919 fwrap_set_fvec(fw2,fherm8,NULL);
1920
1921 // approximation
1922 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1923 ope_opts_set_start(opts,10);
1924 ope_opts_set_coeffs_check(opts,4);
1925 ope_opts_set_tol(opts,1e-15);
1926 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1927 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
1928 opoly_t cpoly3 = orth_poly_expansion_prod(cpoly,cpoly2);
1929
1930 size_t N = 100;
1931 double * pts = linspace(-1,1,N);
1932 size_t ii;
1933 for (ii = 0; ii < N; ii++){
1934 double eval1 = POLY_EVAL(cpoly3,pts[ii]);
1935 double eval2 = POLY_EVAL(cpoly,pts[ii]) *
1936 POLY_EVAL(cpoly2,pts[ii]);
1937 double diff= fabs(eval1-eval2);
1938 CuAssertDblEquals(tc, 0.0, diff, 1e-10);
1939 }
1940
1941 free(pts); pts = NULL;
1942 POLY_FREE(cpoly);
1943 POLY_FREE(cpoly2);
1944 POLY_FREE(cpoly3);
1945 ope_opts_free(opts);
1946 fwrap_destroy(fw);
1947 fwrap_destroy(fw2);
1948 }
1949
1950
Test_hermite_axpy(CuTest * tc)1951 void Test_hermite_axpy(CuTest * tc){
1952
1953 printf("Testing function: orth_poly_expansion_axpy with hermite poly \n");
1954
1955 struct Fwrap * fw = fwrap_create(1,"general-vec");
1956 fwrap_set_fvec(fw,fherm6,NULL);
1957
1958 struct Fwrap * fw2 = fwrap_create(1,"general-vec");
1959 fwrap_set_fvec(fw2,fherm7,NULL);
1960
1961 // approximation
1962 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1963 ope_opts_set_start(opts,10);
1964 ope_opts_set_coeffs_check(opts,4);
1965 ope_opts_set_tol(opts,1e-15);
1966 opoly_t cpoly = orth_poly_expansion_approx_adapt(opts,fw);
1967 opoly_t cpoly2 = orth_poly_expansion_approx_adapt(opts,fw2);
1968 int success = orth_poly_expansion_axpy(2.0,cpoly2,cpoly);
1969 CuAssertIntEquals(tc,0,success);
1970
1971 size_t N = 100;
1972 double * pts = linspace(-1,1,N);
1973 size_t ii;
1974 for (ii = 0; ii < N; ii++){
1975 double eval1 = POLY_EVAL(cpoly,pts[ii]);
1976 double eval2a,eval2b;
1977 fherm7(1,pts+ii,&eval2a,NULL);
1978 fherm6(1,pts+ii,&eval2b,NULL);
1979 double eval2 = 2.0 * eval2a + eval2b;
1980 double diff= fabs(eval1-eval2);
1981 CuAssertDblEquals(tc, 0.0, diff, 1e-7);
1982 }
1983 free(pts); pts = NULL;
1984
1985 POLY_FREE(cpoly);
1986 POLY_FREE(cpoly2);
1987 ope_opts_free(opts);
1988 fwrap_destroy(fw);
1989 fwrap_destroy(fw2);
1990 }
1991
Test_hermite_linear(CuTest * tc)1992 void Test_hermite_linear(CuTest * tc){
1993
1994 printf("Testing function: orth_poly_expansion_linear with hermite poly \n");
1995
1996 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
1997 double a = 2.0, offset=3.0;
1998 opoly_t poly = orth_poly_expansion_linear(a,offset,opts);
1999
2000 size_t N = 100;
2001 double * pts = linspace(-1,1,N);
2002 size_t ii;
2003 for (ii = 0; ii < N; ii++){
2004 double eval1 = POLY_EVAL(poly,pts[ii]);
2005 double eval2 = a*pts[ii] + offset;
2006 double diff= fabs(eval1-eval2);
2007 CuAssertDblEquals(tc, 0.0, diff, 1e-7);
2008 }
2009 free(pts); pts = NULL;
2010
2011 POLY_FREE(poly);
2012 ope_opts_free(opts);
2013 }
2014
Test_hermite_quadratic(CuTest * tc)2015 void Test_hermite_quadratic(CuTest * tc){
2016
2017 printf("Testing function: orth_poly_expansion_quadratic with hermite poly \n");
2018
2019 struct OpeOpts * opts = ope_opts_alloc(HERMITE);
2020 double a = 2.0, offset=3.0;
2021 opoly_t poly = orth_poly_expansion_quadratic(a,offset,opts);
2022
2023 size_t N = 100;
2024 double * pts = linspace(-1,1,N);
2025 size_t ii;
2026 for (ii = 0; ii < N; ii++){
2027 double eval1 = POLY_EVAL(poly,pts[ii]);
2028 double eval2 = a*pow(pts[ii] - offset,2);
2029 double diff= fabs(eval1-eval2);
2030 CuAssertDblEquals(tc, 0.0, diff, 1e-7);
2031 }
2032 free(pts); pts = NULL;
2033
2034 POLY_FREE(poly);
2035 ope_opts_free(opts);
2036 }
2037
Test_hermite_orth_poly_expansion_create_with_params_and_grad(CuTest * tc)2038 void Test_hermite_orth_poly_expansion_create_with_params_and_grad(CuTest * tc){
2039
2040 printf("Testing functions: orth_poly_expansion_create_with_params and grad with hermite poly\n");
2041
2042 size_t nparams = 10;
2043 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2044 0.6,-0.9,-.2,0.04,1.2};
2045 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2046 struct OrthPolyExpansion * ope = NULL;
2047 ope = orth_poly_expansion_create_with_params(opts,nparams,params);
2048
2049
2050 double grad[10];
2051 double xloc = 0.4;
2052 int res =
2053 orth_poly_expansion_param_grad_eval(ope,1,&xloc,grad);
2054 CuAssertIntEquals(tc,0,res);
2055
2056
2057 // numerical derivative
2058 struct OrthPolyExpansion * ope1 = NULL;
2059 struct OrthPolyExpansion * ope2 = NULL;
2060
2061 size_t dim = nparams;
2062 double * x1 = calloc_double(dim);
2063 double * x2 = calloc_double(dim);
2064 for (size_t ii = 0; ii < dim; ii++){
2065 x1[ii] = params[ii];
2066 x2[ii] = params[ii];
2067 }
2068
2069 double diff = 0.0;
2070 double v1,v2;
2071 double norm = 0.0;
2072 double eps = 1e-8;
2073 for (size_t ii = 0; ii < dim; ii++){
2074 x1[ii] += eps;
2075 x2[ii] -= eps;
2076 ope1 = orth_poly_expansion_create_with_params(opts,nparams,x1);
2077 v1 = orth_poly_expansion_eval(ope1,xloc);
2078
2079 ope2 = orth_poly_expansion_create_with_params(opts,nparams,x2);
2080 v2 = orth_poly_expansion_eval(ope2,xloc);
2081
2082 double diff_iter = pow( (v1-v2)/2.0/eps - grad[ii], 2 );
2083 /* printf("current diff = %G\n",diff_iter); */
2084 /* printf("\t norm = %G\n",grad[ii]); */
2085 diff += diff_iter;
2086 norm += pow( (v1-v2)/2.0/eps,2);
2087
2088 x1[ii] -= eps;
2089 x2[ii] += eps;
2090
2091 orth_poly_expansion_free(ope1); ope1 = NULL;
2092 orth_poly_expansion_free(ope2); ope2 = NULL;
2093 }
2094 if (norm > 1){
2095 diff /= norm;
2096 }
2097 free(x1); x1 = NULL;
2098 free(x2); x2 = NULL;
2099 CuAssertDblEquals(tc,0.0,diff,1e-7);
2100
2101 ope_opts_free(opts); opts = NULL;
2102 orth_poly_expansion_free(ope); ope = NULL;
2103 }
2104
2105
HermGetSuite()2106 CuSuite * HermGetSuite(){
2107
2108 CuSuite * suite = CuSuiteNew();
2109 SUITE_ADD_TEST(suite, Test_hermite_approx);
2110 SUITE_ADD_TEST(suite, Test_hermite_approx_adapt);
2111 SUITE_ADD_TEST(suite, Test_hermite_derivative);
2112 SUITE_ADD_TEST(suite, Test_hermite_integrate);
2113 SUITE_ADD_TEST(suite, Test_hermite_non_std_normal);
2114 SUITE_ADD_TEST(suite, Test_hermite_inner);
2115 SUITE_ADD_TEST(suite, Test_hermite_norm_w);
2116 SUITE_ADD_TEST(suite, Test_hermite_norm);
2117 SUITE_ADD_TEST(suite, Test_hermite_product);
2118 SUITE_ADD_TEST(suite, Test_hermite_axpy);
2119 SUITE_ADD_TEST(suite, Test_hermite_linear);
2120 SUITE_ADD_TEST(suite, Test_hermite_quadratic);
2121 SUITE_ADD_TEST(suite, Test_hermite_orth_poly_expansion_create_with_params_and_grad);
2122 return suite;
2123 }
2124
2125
Test_orth_to_standard_poly(CuTest * tc)2126 void Test_orth_to_standard_poly(CuTest * tc){
2127
2128 printf("Testing function: orth_to_standard_poly \n");
2129 struct OrthPoly * leg = init_leg_poly();
2130 struct OrthPoly * cheb = init_cheb_poly();
2131
2132 struct StandardPoly * p = orth_to_standard_poly(leg,0);
2133 CuAssertDblEquals(tc, 1.0, p->coeff[0], 1e-13);
2134 standard_poly_free(p);
2135
2136 p = orth_to_standard_poly(leg,1);
2137 CuAssertDblEquals(tc, 0.0, p->coeff[0], 1e-13);
2138 CuAssertDblEquals(tc, 1.0*sqrt(2 * 1 + 1), p->coeff[1], 1e-13);
2139 standard_poly_free(p);
2140
2141 p = orth_to_standard_poly(leg,5);
2142 CuAssertDblEquals(tc, 0.0, p->coeff[0], 1e-13);
2143 CuAssertDblEquals(tc, 15.0/8.0*sqrt(2 * 5 + 1), p->coeff[1], 1e-13);
2144 CuAssertDblEquals(tc, 0.0, p->coeff[2], 1e-13);
2145 CuAssertDblEquals(tc, -70.0/8.0*sqrt(2 * 5 + 1), p->coeff[3], 1e-13);
2146 CuAssertDblEquals(tc, 0.0, p->coeff[4], 1e-13);
2147 CuAssertDblEquals(tc, 63.0/8.0*sqrt(2 * 5 + 1), p->coeff[5], 1e-13);
2148 standard_poly_free(p);
2149
2150 p = orth_to_standard_poly(cheb,5);
2151 CuAssertDblEquals(tc, 0.0, p->coeff[0], 1e-13);
2152 CuAssertDblEquals(tc, 5.0, p->coeff[1], 1e-13);
2153 CuAssertDblEquals(tc, 0.0, p->coeff[2], 1e-13);
2154 CuAssertDblEquals(tc, -20.0, p->coeff[3], 1e-13);
2155 CuAssertDblEquals(tc, 0.0, p->coeff[4], 1e-13);
2156 CuAssertDblEquals(tc, 16.0, p->coeff[5], 1e-13);
2157 standard_poly_free(p);
2158
2159
2160 free_orth_poly(leg);
2161 free_orth_poly(cheb);
2162 }
2163
Test_orth_poly_expansion_to_standard_poly(CuTest * tc)2164 void Test_orth_poly_expansion_to_standard_poly(CuTest * tc){
2165
2166 printf("Testing function: orth_poly_expansion_to_standard_poly \n");
2167
2168 struct OrthPolyExpansion * pl =
2169 orth_poly_expansion_init(LEGENDRE,10,-1.0,1.0);
2170 pl->coeff[0] = 5.0;
2171 pl->coeff[4] = 2.0 / sqrt(2.0*4+1);
2172 pl->coeff[7] = 3.0 / sqrt(2.0*7+1);
2173
2174 struct StandardPoly * p = orth_poly_expansion_to_standard_poly(pl);
2175
2176 CuAssertDblEquals(tc, 5.0 + 2.0*3.0/8.0, p->coeff[0], 1e-13);
2177 CuAssertDblEquals(tc, 3.0 * -35.0/16.0, p->coeff[1], 1e-13);
2178 CuAssertDblEquals(tc, 2.0 * -30.0/8.0, p->coeff[2], 1e-13);
2179 CuAssertDblEquals(tc, 3.0 * 315.0/16.0, p->coeff[3], 1e-13);
2180 CuAssertDblEquals(tc, 2.0 * 35.0/8.0, p->coeff[4], 1e-13);
2181 CuAssertDblEquals(tc, 3.0 * -693.0/16.0, p->coeff[5], 1e-13);
2182 CuAssertDblEquals(tc, 0.0, p->coeff[6], 1e-13);
2183 CuAssertDblEquals(tc, 3.0 * 429/16.0, p->coeff[7], 1e-13);
2184 CuAssertDblEquals(tc, 0.0, p->coeff[8], 1e-13);
2185 CuAssertDblEquals(tc, 0.0, p->coeff[9], 1e-13);
2186
2187 standard_poly_free(p);
2188 orth_poly_expansion_free(pl);
2189 }
2190
opeTosp(size_t N,const double * x,double * out,void * args)2191 int opeTosp(size_t N,const double * x, double * out, void * args)
2192 {
2193 (void)(args);
2194 for (size_t ii = 0; ii < N; ii++){
2195 out[ii] = 1.0 + 2.0 * x[ii] +
2196 5.0 * pow(x[ii],3) +
2197 2.0 * pow(x[ii],5) +
2198 1.5 * pow(x[ii],6);
2199 }
2200 return 0;
2201 }
2202
Test_orth_poly_expansion_to_standard_poly_ab(CuTest * tc)2203 void Test_orth_poly_expansion_to_standard_poly_ab(CuTest * tc){
2204
2205 printf("Testing function: orth_expansion_to_standard on (-3.0, 2.0) \n");
2206
2207 struct Fwrap * fw = fwrap_create(1,"general-vec");
2208 fwrap_set_fvec(fw,opeTosp,NULL);
2209
2210 // approximation
2211 double lb = -3.0, ub = 2.0;
2212 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2213 ope_opts_set_start(opts,10);
2214 ope_opts_set_coeffs_check(opts,4);
2215 ope_opts_set_tol(opts,1e-15);
2216 ope_opts_set_lb(opts,lb);
2217 ope_opts_set_ub(opts,ub);
2218 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
2219
2220 struct StandardPoly * p = orth_poly_expansion_to_standard_poly(pl);
2221
2222 size_t ii;
2223 for (ii = 0; ii < p->num_poly; ii++){
2224 if (ii == 0){
2225 CuAssertDblEquals(tc, 1.0, p->coeff[ii], 1e-10);
2226 }
2227 else if (ii == 1){
2228 CuAssertDblEquals(tc, 2.0, p->coeff[ii], 1e-10);
2229 }
2230 else if (ii == 3){
2231 CuAssertDblEquals(tc, 5.0, p->coeff[ii], 1e-10);
2232 }
2233 else if (ii == 5){
2234 CuAssertDblEquals(tc, 2.0, p->coeff[ii], 1e-8);
2235 }
2236 else if (ii == 6){
2237 CuAssertDblEquals(tc, 1.5, p->coeff[ii], 1e-8);
2238 }
2239 else{
2240 CuAssertDblEquals(tc, 0.0, p->coeff[ii], 1e-8);
2241 }
2242 }
2243
2244 standard_poly_free(p);
2245 ope_opts_free(opts);
2246 POLY_FREE(pl);
2247 fwrap_destroy(fw);
2248 }
2249
StandardPolyGetSuite()2250 CuSuite * StandardPolyGetSuite(){
2251
2252 CuSuite * suite = CuSuiteNew();
2253 SUITE_ADD_TEST(suite, Test_orth_to_standard_poly);
2254 SUITE_ADD_TEST(suite, Test_orth_poly_expansion_to_standard_poly);
2255 SUITE_ADD_TEST(suite, Test_orth_poly_expansion_to_standard_poly_ab);
2256
2257 return suite;
2258 }
2259
Test_orth_poly_expansion_real_roots(CuTest * tc)2260 void Test_orth_poly_expansion_real_roots(CuTest * tc){
2261
2262 printf("Testing function: orth_poly_expansion_real_roots \n");
2263
2264
2265 struct Fwrap * fw = fwrap_create(1,"general-vec");
2266 fwrap_set_fvec(fw,polyroots,NULL);
2267
2268 // approximation
2269 double lb = -3.0, ub = 2.0;
2270 /* double lb = -1.0, ub = 1.0; */
2271 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2272 ope_opts_set_start(opts,10);
2273 ope_opts_set_coeffs_check(opts,4);
2274 ope_opts_set_tol(opts,1e-15);
2275 ope_opts_set_lb(opts,lb);
2276 ope_opts_set_ub(opts,ub);
2277 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
2278
2279 size_t nroots;
2280 double * roots = orth_poly_expansion_real_roots(pl, &nroots);
2281
2282 /* printf("roots are: "); */
2283 /* dprint(nroots, roots); */
2284
2285 CuAssertIntEquals(tc, 5, nroots);
2286 CuAssertDblEquals(tc, -3.0, roots[0], 1e-9);
2287 CuAssertDblEquals(tc, 0.0, roots[1], 1e-9);
2288 CuAssertDblEquals(tc, 1.0, roots[2], 1e-5);
2289 CuAssertDblEquals(tc, 1.0, roots[3], 1e-5);
2290 CuAssertDblEquals(tc, 2.0, roots[4], 1e-9);
2291
2292 free(roots);
2293 POLY_FREE(pl);
2294 ope_opts_free(opts);
2295 fwrap_destroy(fw);
2296 }
2297
2298
Test_maxmin_poly_expansion(CuTest * tc)2299 void Test_maxmin_poly_expansion(CuTest * tc){
2300
2301 printf("Testing functions: absmax, max and min of orth_poly_expansion \n");
2302
2303 struct Fwrap * fw = fwrap_create(1,"general-vec");
2304 fwrap_set_fvec(fw,maxminpoly,NULL);
2305
2306 // approximation
2307 double lb = -1.0, ub = 2.0;
2308 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2309 ope_opts_set_start(opts,10);
2310 ope_opts_set_coeffs_check(opts,4);
2311 ope_opts_set_tol(opts,1e-15);
2312 ope_opts_set_lb(opts,lb);
2313 ope_opts_set_ub(opts,ub);
2314 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
2315
2316 double loc;
2317 double max = orth_poly_expansion_max(pl, &loc);
2318 double min = orth_poly_expansion_min(pl, &loc);
2319 double absmax = orth_poly_expansion_absmax(pl, &loc,NULL);
2320
2321 double diff;
2322
2323 diff = fabs(1.0-max);
2324 //printf("diff =%G\n",diff);
2325 CuAssertDblEquals(tc,0.0, diff, 1e-9);
2326 CuAssertDblEquals(tc, -1.0, min, 1e-9);
2327 CuAssertDblEquals(tc, 1.0, absmax, 1e-9);
2328
2329 POLY_FREE(pl);
2330 ope_opts_free(opts);
2331 fwrap_destroy(fw);
2332 }
2333
PolyAlgorithmsGetSuite()2334 CuSuite * PolyAlgorithmsGetSuite(){
2335
2336 CuSuite * suite = CuSuiteNew();
2337 SUITE_ADD_TEST(suite, Test_orth_poly_expansion_real_roots);
2338 SUITE_ADD_TEST(suite, Test_maxmin_poly_expansion);
2339
2340 return suite;
2341 }
2342
Test_serialize_orth_poly(CuTest * tc)2343 void Test_serialize_orth_poly(CuTest * tc){
2344
2345 printf("Testing functions: (de)serialize_orth_poly \n");
2346
2347 struct OrthPoly * poly = init_leg_poly();
2348
2349 unsigned char * text = serialize_orth_poly(poly);
2350
2351 struct OrthPoly * pt = deserialize_orth_poly(text);
2352 CuAssertIntEquals(tc,0,pt->ptype);
2353 free(text);
2354 free_orth_poly(poly);
2355 free_orth_poly(pt);
2356
2357 poly = init_cheb_poly();
2358 text = serialize_orth_poly(poly);
2359 pt = deserialize_orth_poly(text);
2360 CuAssertIntEquals(tc,1,pt->ptype);
2361 free(text);
2362 free_orth_poly(poly);
2363 free_orth_poly(pt);
2364
2365 poly = init_cheb_poly();
2366 text = serialize_orth_poly(poly);
2367 pt = deserialize_orth_poly(text);
2368 CuAssertIntEquals(tc,1,pt->ptype);
2369 free(text);
2370 free_orth_poly(poly);
2371 free_orth_poly(pt);
2372
2373 poly = init_hermite_poly();
2374 text = serialize_orth_poly(poly);
2375 pt = deserialize_orth_poly(text);
2376 CuAssertIntEquals(tc,2,pt->ptype);
2377 free(text);
2378 free_orth_poly(poly);
2379 free_orth_poly(pt);
2380
2381 poly = init_fourier_poly();
2382 text = serialize_orth_poly(poly);
2383 pt = deserialize_orth_poly(text);
2384 CuAssertIntEquals(tc,4,pt->ptype);
2385 free(text);
2386 free_orth_poly(poly);
2387 free_orth_poly(pt);
2388
2389 }
2390
Test_serialize_orth_poly_expansion(CuTest * tc)2391 void Test_serialize_orth_poly_expansion(CuTest * tc){
2392
2393 printf("Testing functions: (de)serializing orth_poly_expansion \n");
2394
2395 struct Fwrap * fw = fwrap_create(1,"general-vec");
2396 fwrap_set_fvec(fw,maxminpoly,NULL);
2397
2398 // approximation
2399 double lb = -1.0, ub = 2.0;
2400 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2401 ope_opts_set_start(opts,10);
2402 ope_opts_set_coeffs_check(opts,4);
2403 ope_opts_set_tol(opts,1e-15);
2404 ope_opts_set_lb(opts,lb);
2405 ope_opts_set_ub(opts,ub);
2406 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
2407
2408 unsigned char * text = NULL;
2409 size_t size_to_be;
2410 serialize_orth_poly_expansion(text, pl, &size_to_be);
2411 text = malloc(size_to_be * sizeof(char));
2412
2413 serialize_orth_poly_expansion(text, pl, NULL);
2414
2415 opoly_t pt = NULL;
2416 deserialize_orth_poly_expansion(text, &pt);
2417
2418 double * xtest = linspace(lb,ub,1000);
2419 size_t ii;
2420 double err = 0.0;
2421 for (ii = 0; ii < 1000; ii++){
2422 err += pow(POLY_EVAL(pl,xtest[ii]) - POLY_EVAL(pt,xtest[ii]),2);
2423 }
2424 err = sqrt(err);
2425 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2426 free(xtest);
2427 free(text);
2428
2429 POLY_FREE(pl);
2430 POLY_FREE(pt);
2431 ope_opts_free(opts);
2432 fwrap_destroy(fw);
2433 }
2434
Test_serialize_orth_poly_expansion_fourier(CuTest * tc)2435 void Test_serialize_orth_poly_expansion_fourier(CuTest * tc){
2436
2437 printf("Testing functions: (de)serializing orth_poly_expansion_fourier \n");
2438
2439 struct Fwrap * fw = fwrap_create(1,"general-vec");
2440 fwrap_set_fvec(fw,gaussbump,NULL);
2441
2442 // approximation
2443 size_t n = 16;
2444 double lb = -4.0, ub = 2.0;
2445 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
2446 ope_opts_set_lb(opts,lb);
2447 ope_opts_set_ub(opts,ub);
2448 ope_opts_set_start(opts, n+1);
2449
2450
2451 struct OrthPolyExpansion * pl = orth_poly_expansion_init_from_opts(opts, n+1);
2452 orth_poly_expansion_approx_vec(pl, fw, NULL);
2453
2454 unsigned char * text = NULL;
2455 size_t size_to_be;
2456 serialize_orth_poly_expansion(text, pl, &size_to_be);
2457 text = malloc(size_to_be * sizeof(char));
2458
2459 serialize_orth_poly_expansion(text, pl, NULL);
2460
2461 opoly_t pt = NULL;
2462 deserialize_orth_poly_expansion(text, &pt);
2463
2464 double * xtest = linspace(lb,ub,1000);
2465 size_t ii;
2466 double err = 0.0;
2467 for (ii = 0; ii < 1000; ii++){
2468 err += pow(POLY_EVAL(pl,xtest[ii]) - POLY_EVAL(pt,xtest[ii]),2);
2469 }
2470 err = sqrt(err);
2471 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2472 free(xtest);
2473 free(text);
2474
2475 POLY_FREE(pl);
2476 POLY_FREE(pt);
2477 ope_opts_free(opts);
2478 fwrap_destroy(fw);
2479 }
2480
Test_orth_poly_expansion_savetxt(CuTest * tc)2481 void Test_orth_poly_expansion_savetxt(CuTest * tc){
2482
2483 printf("Testing functions: orth_poly_expansion_savetxt and _loadtxt \n");
2484
2485 struct Fwrap * fw = fwrap_create(1,"general-vec");
2486 fwrap_set_fvec(fw,maxminpoly,NULL);
2487
2488 // approximation
2489 double lb = -1.0, ub = 2.0;
2490 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2491 ope_opts_set_start(opts,10);
2492 ope_opts_set_coeffs_check(opts,4);
2493 ope_opts_set_tol(opts,1e-15);
2494 ope_opts_set_lb(opts,lb);
2495 ope_opts_set_ub(opts,ub);
2496 opoly_t pl = orth_poly_expansion_approx_adapt(opts,fw);
2497
2498
2499 FILE * fp = fopen("testorthpoly.txt","w+");
2500 size_t prec = 21;
2501 orth_poly_expansion_savetxt(pl,fp,prec);
2502
2503 opoly_t pt = NULL;
2504 rewind(fp);
2505 pt = orth_poly_expansion_loadtxt(fp);
2506
2507 double * xtest = linspace(lb,ub,1000);
2508 size_t ii;
2509 double err = 0.0;
2510 for (ii = 0; ii < 1000; ii++){
2511 err += pow(POLY_EVAL(pl,xtest[ii]) - POLY_EVAL(pt,xtest[ii]),2);
2512 }
2513 err = sqrt(err);
2514 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2515 free(xtest);
2516
2517 fclose(fp);
2518 POLY_FREE(pl);
2519 POLY_FREE(pt);
2520 ope_opts_free(opts);
2521 fwrap_destroy(fw);
2522 }
2523
Test_serialize_generic_function(CuTest * tc)2524 void Test_serialize_generic_function(CuTest * tc){
2525
2526 printf("Testing functions: (de)serializing generic_function \n");
2527
2528 struct Fwrap * fw = fwrap_create(1,"general-vec");
2529 fwrap_set_fvec(fw,maxminpoly,NULL);
2530
2531 // approximation
2532 double lb = -1.0, ub = 2.0;
2533 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2534 ope_opts_set_start(opts,10);
2535 ope_opts_set_coeffs_check(opts,4);
2536 ope_opts_set_tol(opts,1e-15);
2537 ope_opts_set_lb(opts,lb);
2538 ope_opts_set_ub(opts,ub);
2539
2540 struct GenericFunction * pl =
2541 generic_function_approximate1d(POLYNOMIAL,opts,fw);
2542
2543 unsigned char * text = NULL;
2544 size_t size_to_be;
2545 serialize_generic_function(text, pl, &size_to_be);
2546 text = malloc(size_to_be * sizeof(char));
2547
2548 serialize_generic_function(text, pl, NULL);
2549
2550 struct GenericFunction * pt = NULL;
2551 deserialize_generic_function(text, &pt);
2552
2553 double * xtest = linspace(lb,ub,1000);
2554 size_t ii;
2555 double err = 0.0;
2556 for (ii = 0; ii < 1000; ii++){
2557 err += pow(generic_function_1d_eval(pl,xtest[ii]) -
2558 generic_function_1d_eval(pt,xtest[ii]),2);
2559 }
2560 err = sqrt(err);
2561 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2562 free(xtest);
2563 free(text);
2564
2565 ope_opts_free(opts);
2566 generic_function_free(pl);
2567 generic_function_free(pt);
2568 fwrap_destroy(fw);
2569 }
2570
Test_serialize_generic_function_fourier(CuTest * tc)2571 void Test_serialize_generic_function_fourier(CuTest * tc){
2572
2573 printf("Testing functions: (de)serializing generic_function_fourier \n");
2574
2575 struct Fwrap * fw = fwrap_create(1,"general-vec");
2576 fwrap_set_fvec(fw,maxminpoly,NULL);
2577
2578 size_t n = 16;
2579 double lb = -4.0, ub = 2.0;
2580 struct OpeOpts * opts = ope_opts_alloc(FOURIER);
2581 ope_opts_set_lb(opts,lb);
2582 ope_opts_set_ub(opts,ub);
2583 ope_opts_set_start(opts, n+1);
2584
2585 struct GenericFunction * pl =
2586 generic_function_approximate1d(POLYNOMIAL,opts,fw);
2587
2588 unsigned char * text = NULL;
2589 size_t size_to_be;
2590 serialize_generic_function(text, pl, &size_to_be);
2591 text = malloc(size_to_be * sizeof(char));
2592
2593 serialize_generic_function(text, pl, NULL);
2594
2595 struct GenericFunction * pt = NULL;
2596 deserialize_generic_function(text, &pt);
2597
2598 double * xtest = linspace(lb,ub,1000);
2599 size_t ii;
2600 double err = 0.0;
2601 for (ii = 0; ii < 1000; ii++){
2602 err += pow(generic_function_1d_eval(pl,xtest[ii]) -
2603 generic_function_1d_eval(pt,xtest[ii]),2);
2604 }
2605 err = sqrt(err);
2606 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2607 free(xtest);
2608 free(text);
2609
2610 ope_opts_free(opts);
2611 generic_function_free(pl);
2612 generic_function_free(pt);
2613 fwrap_destroy(fw);
2614 }
2615
Test_generic_function_savetxt(CuTest * tc)2616 void Test_generic_function_savetxt(CuTest * tc){
2617
2618 printf("Testing functions: generic_function_savetxt and _loadtxt \n");
2619
2620 struct Fwrap * fw = fwrap_create(1,"general-vec");
2621 fwrap_set_fvec(fw,maxminpoly,NULL);
2622
2623 // approximation
2624 double lb = -1.0, ub = 2.0;
2625 struct OpeOpts * opts = ope_opts_alloc(LEGENDRE);
2626 ope_opts_set_start(opts,10);
2627 ope_opts_set_coeffs_check(opts,4);
2628 ope_opts_set_tol(opts,1e-15);
2629 ope_opts_set_lb(opts,lb);
2630 ope_opts_set_ub(opts,ub);
2631
2632 struct GenericFunction * pl =
2633 generic_function_approximate1d(POLYNOMIAL,opts,fw);
2634
2635 FILE * fp = fopen("genfunctest.txt","w+");
2636 size_t prec = 21;
2637 generic_function_savetxt(pl,fp,prec);
2638
2639 struct GenericFunction * pt = NULL;
2640 rewind(fp);
2641 pt = generic_function_loadtxt(fp);
2642
2643 double * xtest = linspace(lb,ub,1000);
2644 size_t ii;
2645 double err = 0.0;
2646 for (ii = 0; ii < 1000; ii++){
2647 err += pow(generic_function_1d_eval(pl,xtest[ii]) -
2648 generic_function_1d_eval(pt,xtest[ii]),2);
2649 }
2650 err = sqrt(err);
2651 CuAssertDblEquals(tc, 0.0, err, 1e-15);
2652 free(xtest);
2653
2654
2655 fclose(fp);
2656 ope_opts_free(opts);
2657 generic_function_free(pl);
2658 generic_function_free(pt);
2659 fwrap_destroy(fw);
2660 }
2661
2662
2663
PolySerializationGetSuite()2664 CuSuite * PolySerializationGetSuite(){
2665
2666 CuSuite * suite = CuSuiteNew();
2667 SUITE_ADD_TEST(suite, Test_serialize_orth_poly);
2668 SUITE_ADD_TEST(suite, Test_serialize_orth_poly_expansion);
2669 SUITE_ADD_TEST(suite, Test_serialize_orth_poly_expansion_fourier);
2670 SUITE_ADD_TEST(suite, Test_orth_poly_expansion_savetxt);
2671 SUITE_ADD_TEST(suite, Test_serialize_generic_function);
2672 SUITE_ADD_TEST(suite, Test_serialize_generic_function_fourier);
2673 SUITE_ADD_TEST(suite, Test_generic_function_savetxt);
2674
2675 return suite;
2676 }
2677
regress_func(size_t N,const double * x,double * out)2678 static void regress_func(size_t N, const double * x, double * out)
2679 {
2680 for (size_t ii = 0; ii < N; ii++){
2681 out[ii] = -1.0 * pow(x[ii],7) + 2.0*pow(x[ii],2) + 0.2 * 0.5*x[ii];
2682 }
2683 }
2684
Test_LS_cheb_regress(CuTest * tc)2685 void Test_LS_cheb_regress(CuTest * tc){
2686
2687 printf("Testing functions: least squares regression with cheb\n");
2688
2689 size_t nparams = 10;
2690 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2691 0.6,-0.9,-.2,0.04,1.2};
2692 struct OpeOpts * aopts = ope_opts_alloc(CHEBYSHEV);
2693 ope_opts_set_nparams(aopts,nparams);
2694
2695 // create data
2696 size_t ndata = 5000;
2697 double * x = linspace(-1,1,ndata);
2698 double * y = calloc_double(ndata);
2699 regress_func(ndata,x,y);
2700 // // add noise
2701 for (size_t ii =0 ; ii < ndata; ii++){
2702 y[ii] += randn()*0.01;
2703 }
2704
2705 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
2706 c3opt_set_verbose(optimizer,0);
2707
2708 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
2709 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
2710 regress_1d_opts_set_initial_parameters(regopts,params);
2711
2712 // check derivative
2713 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
2714 double * deriv_diff = calloc_double(nparams);
2715 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
2716 /* for (size_t ii = 0; ii < nparams; ii++){ */
2717 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
2718 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
2719 /* } */
2720 CuAssertDblEquals(tc,0.0,gerr,1e-3);
2721 free(deriv_diff); deriv_diff = NULL;
2722
2723 int info;
2724 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
2725 CuAssertIntEquals(tc,1,info>-1);
2726
2727 double * xtest = linspace(-1.0,1.0,1000);
2728 double * vals = calloc_double(1000);
2729 regress_func(1000,xtest,vals);
2730 size_t ii;
2731 double err = 0.0;
2732 double norm = 0.0;
2733 for (ii = 0; ii < 1000; ii++){
2734 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
2735 norm += vals[ii]*vals[ii];
2736 }
2737 err = sqrt(err);
2738 norm = sqrt(norm);
2739 double rat = err/norm;
2740 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
2741 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
2742 free(xtest); xtest = NULL;
2743 free(vals); vals = NULL;
2744
2745 free(x); x = NULL;
2746 free(y); y = NULL;
2747 ope_opts_free(aopts); aopts = NULL;
2748 regress_1d_opts_destroy(regopts); regopts = NULL;
2749 c3opt_free(optimizer); optimizer = NULL;
2750 generic_function_free(gf); gf = NULL;;
2751 }
2752
Test_LS_leg_regress(CuTest * tc)2753 void Test_LS_leg_regress(CuTest * tc){
2754
2755 printf("Testing functions: least squares regression with legendre\n");
2756
2757 size_t nparams = 10;
2758 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2759 0.6,-0.9,-.2,0.04,1.2};
2760 struct OpeOpts * aopts = ope_opts_alloc(LEGENDRE);
2761 ope_opts_set_nparams(aopts,nparams);
2762
2763 // create data
2764 size_t ndata = 5000;
2765 double * x = linspace(-1,1,ndata);
2766 double * y = calloc_double(ndata);
2767 regress_func(ndata,x,y);
2768 // // add noise
2769 for (size_t ii =0 ; ii < ndata; ii++){
2770 y[ii] += randn()*0.01;
2771 }
2772
2773 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
2774 c3opt_set_verbose(optimizer,0);
2775
2776 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
2777 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
2778 regress_1d_opts_set_initial_parameters(regopts,params);
2779
2780 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
2781 double * deriv_diff = calloc_double(nparams);
2782 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
2783 /* for (size_t ii = 0; ii < nparams; ii++){ */
2784 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
2785 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
2786 /* } */
2787 CuAssertDblEquals(tc,0.0,gerr,1e-3);
2788 free(deriv_diff); deriv_diff = NULL;
2789
2790
2791 int info;
2792 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
2793 CuAssertIntEquals(tc,1,info>-1);
2794
2795 double * xtest = linspace(-1.0,1.0,1000);
2796 double * vals = calloc_double(1000);
2797 regress_func(1000,xtest,vals);
2798 size_t ii;
2799 double err = 0.0;
2800 double norm = 0.0;
2801 for (ii = 0; ii < 1000; ii++){
2802 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
2803 norm += vals[ii]*vals[ii];
2804 }
2805 err = sqrt(err);
2806 norm = sqrt(norm);
2807 double rat = err/norm;
2808 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
2809 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
2810 free(xtest); xtest = NULL;
2811 free(vals); vals = NULL;
2812
2813 free(x); x = NULL;
2814 free(y); y = NULL;
2815 ope_opts_free(aopts); aopts = NULL;
2816 regress_1d_opts_destroy(regopts); regopts = NULL;
2817 c3opt_free(optimizer); optimizer = NULL;
2818 generic_function_free(gf); gf = NULL;;
2819 }
2820
Test_LS_herm_regress(CuTest * tc)2821 void Test_LS_herm_regress(CuTest * tc){
2822
2823 printf("Testing functions: least squares regression with hermite\n");
2824
2825 size_t nparams = 10;
2826 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2827 0.6,-0.9,-.2,0.04,1.2};
2828 struct OpeOpts * aopts = ope_opts_alloc(HERMITE);
2829 ope_opts_set_nparams(aopts,nparams);
2830
2831 // create data
2832 size_t ndata = 5000;
2833 double * x = linspace(-1,1,ndata);
2834 double * y = calloc_double(ndata);
2835 regress_func(ndata,x,y);
2836 // // add noise
2837 for (size_t ii =0 ; ii < ndata; ii++){
2838 y[ii] += randn()*0.01;
2839 }
2840
2841 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
2842 c3opt_set_verbose(optimizer,0);
2843
2844 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
2845 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
2846 regress_1d_opts_set_initial_parameters(regopts,params);
2847
2848 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
2849 double * deriv_diff = calloc_double(nparams);
2850 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
2851 /* for (size_t ii = 0; ii < nparams; ii++){ */
2852 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
2853 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
2854 /* } */
2855 /* printf("gerr = %G\n",gerr); */
2856 CuAssertDblEquals(tc,0.0,gerr,1e-3);
2857 free(deriv_diff); deriv_diff = NULL;
2858
2859
2860 int info;
2861 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
2862 CuAssertIntEquals(tc,1,info>-1);
2863
2864 double * xtest = linspace(-1.0,1.0,1000);
2865 double * vals = calloc_double(1000);
2866 regress_func(1000,xtest,vals);
2867 size_t ii;
2868 double err = 0.0;
2869 double norm = 0.0;
2870 for (ii = 0; ii < 1000; ii++){
2871 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
2872 norm += vals[ii]*vals[ii];
2873 }
2874 err = sqrt(err);
2875 norm = sqrt(norm);
2876 double rat = err/norm;
2877 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
2878 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
2879 free(xtest); xtest = NULL;
2880 free(vals); vals = NULL;
2881
2882 free(x); x = NULL;
2883 free(y); y = NULL;
2884 ope_opts_free(aopts); aopts = NULL;
2885 regress_1d_opts_destroy(regopts); regopts = NULL;
2886 c3opt_free(optimizer); optimizer = NULL;
2887 generic_function_free(gf); gf = NULL;;
2888 }
2889
Test_RLS2_cheb_regress(CuTest * tc)2890 void Test_RLS2_cheb_regress(CuTest * tc){
2891
2892 printf("Testing functions: ridge regression with cheb\n");
2893
2894 size_t nparams = 10;
2895 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2896 0.6,-0.9,-.2,0.04,1.2};
2897 struct OpeOpts * aopts = ope_opts_alloc(CHEBYSHEV);
2898 ope_opts_set_nparams(aopts,nparams);
2899
2900 // create data
2901 size_t ndata = 5000;
2902 double * x = linspace(-1,1,ndata);
2903 double * y = calloc_double(ndata);
2904 regress_func(ndata,x,y);
2905 // // add noise
2906 for (size_t ii =0 ; ii < ndata; ii++){
2907 y[ii] += randn()*0.01;
2908 }
2909
2910 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
2911 c3opt_set_verbose(optimizer,0);
2912
2913 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLS2,ndata,x,y);
2914 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
2915 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
2916 regress_1d_opts_set_initial_parameters(regopts,params);
2917
2918 c3opt_add_objective(optimizer,param_RLS2regress_cost,regopts);
2919 double * deriv_diff = calloc_double(nparams);
2920 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
2921 /* for (size_t ii = 0; ii < nparams; ii++){ */
2922 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
2923 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
2924 /* } */
2925 /* printf("gerr = %G\n",gerr); */
2926 CuAssertDblEquals(tc,0.0,gerr,1e-3);
2927 free(deriv_diff); deriv_diff = NULL;
2928
2929 int info;
2930 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
2931 CuAssertIntEquals(tc,1,info>-1);
2932
2933 double * xtest = linspace(-1.0,1.0,1000);
2934 double * vals = calloc_double(1000);
2935 regress_func(1000,xtest,vals);
2936 size_t ii;
2937 double err = 0.0;
2938 double norm = 0.0;
2939 for (ii = 0; ii < 1000; ii++){
2940 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
2941 norm += vals[ii]*vals[ii];
2942 }
2943 err = sqrt(err);
2944 norm = sqrt(norm);
2945 double rat = err/norm;
2946 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
2947 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
2948 free(xtest); xtest = NULL;
2949 free(vals); vals = NULL;
2950
2951 free(x); x = NULL;
2952 free(y); y = NULL;
2953 ope_opts_free(aopts); aopts = NULL;
2954 regress_1d_opts_destroy(regopts); regopts = NULL;
2955 c3opt_free(optimizer); optimizer = NULL;
2956 generic_function_free(gf); gf = NULL;;
2957 }
2958
Test_RLS2_leg_regress(CuTest * tc)2959 void Test_RLS2_leg_regress(CuTest * tc){
2960
2961 printf("Testing functions: ridge regression with legendre\n");
2962
2963 size_t nparams = 10;
2964 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
2965 0.6,-0.9,-.2,0.04,1.2};
2966 struct OpeOpts * aopts = ope_opts_alloc(LEGENDRE);
2967 ope_opts_set_nparams(aopts,nparams);
2968
2969 // create data
2970 size_t ndata = 5000;
2971 double * x = linspace(-1,1,ndata);
2972 double * y = calloc_double(ndata);
2973 regress_func(ndata,x,y);
2974 // // add noise
2975 for (size_t ii =0 ; ii < ndata; ii++){
2976 y[ii] += randn()*0.01;
2977 }
2978
2979 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
2980 c3opt_set_verbose(optimizer,0);
2981
2982 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLS2,ndata,x,y);
2983 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
2984 regress_1d_opts_set_regularization_penalty(regopts,1e2/sqrt(ndata));
2985 regress_1d_opts_set_initial_parameters(regopts,params);
2986
2987 c3opt_add_objective(optimizer,param_RLS2regress_cost,regopts);
2988 double * deriv_diff = calloc_double(nparams);
2989 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
2990 /* for (size_t ii = 0; ii < nparams; ii++){ */
2991 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
2992 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
2993 /* } */
2994 /* printf("gerr = %G\n",gerr); */
2995 CuAssertDblEquals(tc,0.0,gerr,1e-3);
2996 free(deriv_diff); deriv_diff = NULL;
2997
2998 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
2999 int info;
3000 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3001 CuAssertIntEquals(tc,1,info>-1);
3002
3003 double * xtest = linspace(-1.0,1.0,1000);
3004 double * vals = calloc_double(1000);
3005 regress_func(1000,xtest,vals);
3006 size_t ii;
3007 double err = 0.0;
3008 double norm = 0.0;
3009 for (ii = 0; ii < 1000; ii++){
3010 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3011 norm += vals[ii]*vals[ii];
3012 }
3013 err = sqrt(err);
3014 norm = sqrt(norm);
3015 double rat = err/norm;
3016 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3017 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3018 free(xtest); xtest = NULL;
3019 free(vals); vals = NULL;
3020
3021 free(x); x = NULL;
3022 free(y); y = NULL;
3023 ope_opts_free(aopts); aopts = NULL;
3024 regress_1d_opts_destroy(regopts); regopts = NULL;
3025 c3opt_free(optimizer); optimizer = NULL;
3026 generic_function_free(gf); gf = NULL;;
3027 }
3028
Test_RLS2_herm_regress(CuTest * tc)3029 void Test_RLS2_herm_regress(CuTest * tc){
3030
3031 printf("Testing functions: ridge regression with hermite\n");
3032
3033 size_t nparams = 10;
3034 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3035 0.6,-0.9,-.2,0.04,1.2};
3036 struct OpeOpts * aopts = ope_opts_alloc(HERMITE);
3037 ope_opts_set_nparams(aopts,nparams);
3038
3039 // create data
3040 size_t ndata = 5000;
3041 double * x = linspace(-1,1,ndata);
3042 double * y = calloc_double(ndata);
3043 regress_func(ndata,x,y);
3044 // // add noise
3045 for (size_t ii =0 ; ii < ndata; ii++){
3046 y[ii] += randn()*0.01;
3047 }
3048
3049 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3050 c3opt_set_verbose(optimizer,0);
3051
3052 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLS2,ndata,x,y);
3053 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3054 regress_1d_opts_set_regularization_penalty(regopts,1e-4/sqrt(ndata));
3055 regress_1d_opts_set_initial_parameters(regopts,params);
3056
3057 c3opt_add_objective(optimizer,param_RLS2regress_cost,regopts);
3058 double * deriv_diff = calloc_double(nparams);
3059 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3060 /* for (size_t ii = 0; ii < nparams; ii++){ */
3061 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3062 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3063 /* } */
3064 /* printf("gerr = %G\n",gerr); */
3065 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3066 free(deriv_diff); deriv_diff = NULL;
3067
3068
3069 int info;
3070 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3071 CuAssertIntEquals(tc,1,info>-1);
3072
3073 double * xtest = linspace(-1.0,1.0,1000);
3074 double * vals = calloc_double(1000);
3075 regress_func(1000,xtest,vals);
3076 size_t ii;
3077 double err = 0.0;
3078 double norm = 0.0;
3079 for (ii = 0; ii < 1000; ii++){
3080 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3081 norm += vals[ii]*vals[ii];
3082 }
3083 err = sqrt(err);
3084 norm = sqrt(norm);
3085 double rat = err/norm;
3086 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3087 printf("\t note: must have smaller regularization\n");
3088 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3089 free(xtest); xtest = NULL;
3090 free(vals); vals = NULL;
3091
3092 free(x); x = NULL;
3093 free(y); y = NULL;
3094 ope_opts_free(aopts); aopts = NULL;
3095 regress_1d_opts_destroy(regopts); regopts = NULL;
3096 c3opt_free(optimizer); optimizer = NULL;
3097 generic_function_free(gf); gf = NULL;;
3098 }
3099
Test_RLSD2_cheb_regress(CuTest * tc)3100 void Test_RLSD2_cheb_regress(CuTest * tc){
3101
3102 printf("Testing functions: second deriv reg with cheb\n");
3103
3104 size_t nparams = 10;
3105 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3106 0.6,-0.9,-.2,0.04,1.2};
3107 struct OpeOpts * aopts = ope_opts_alloc(CHEBYSHEV);
3108 ope_opts_set_nparams(aopts,nparams);
3109
3110 // create data
3111 size_t ndata = 5000;
3112 double * x = linspace(-1,1,ndata);
3113 double * y = calloc_double(ndata);
3114 regress_func(ndata,x,y);
3115 // // add noise
3116 for (size_t ii =0 ; ii < ndata; ii++){
3117 y[ii] += randn()*0.01;
3118 }
3119
3120 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3121 c3opt_set_verbose(optimizer,0);
3122
3123 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSD2,ndata,x,y);
3124 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3125 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3126 regress_1d_opts_set_initial_parameters(regopts,params);
3127
3128 c3opt_add_objective(optimizer,param_RLSD2regress_cost,regopts);
3129 double * deriv_diff = calloc_double(nparams);
3130 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3131 /* for (size_t ii = 0; ii < nparams; ii++){ */
3132 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3133 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3134 /* } */
3135 /* printf("gerr = %G\n",gerr); */
3136 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3137 free(deriv_diff); deriv_diff = NULL;
3138
3139 int info;
3140 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3141 CuAssertIntEquals(tc,1,info>-1);
3142
3143 double * xtest = linspace(-1.0,1.0,1000);
3144 double * vals = calloc_double(1000);
3145 regress_func(1000,xtest,vals);
3146 size_t ii;
3147 double err = 0.0;
3148 double norm = 0.0;
3149 for (ii = 0; ii < 1000; ii++){
3150 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3151 norm += vals[ii]*vals[ii];
3152 }
3153 err = sqrt(err);
3154 norm = sqrt(norm);
3155 double rat = err/norm;
3156 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3157 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3158 free(xtest); xtest = NULL;
3159 free(vals); vals = NULL;
3160
3161 free(x); x = NULL;
3162 free(y); y = NULL;
3163 ope_opts_free(aopts); aopts = NULL;
3164 regress_1d_opts_destroy(regopts); regopts = NULL;
3165 c3opt_free(optimizer); optimizer = NULL;
3166 generic_function_free(gf); gf = NULL;;
3167 }
3168
Test_RLSD2_leg_regress(CuTest * tc)3169 void Test_RLSD2_leg_regress(CuTest * tc){
3170
3171 printf("Testing functions: second deriv reg with legendre\n");
3172
3173 size_t nparams = 10;
3174 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3175 0.6,-0.9,-.2,0.04,1.2};
3176 struct OpeOpts * aopts = ope_opts_alloc(LEGENDRE);
3177 ope_opts_set_nparams(aopts,nparams);
3178
3179 // create data
3180 size_t ndata = 5000;
3181 double * x = linspace(-1,1,ndata);
3182 double * y = calloc_double(ndata);
3183 regress_func(ndata,x,y);
3184 // // add noise
3185 for (size_t ii =0 ; ii < ndata; ii++){
3186 y[ii] += randn()*0.01;
3187 }
3188
3189 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3190 c3opt_set_verbose(optimizer,0);
3191
3192 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSD2,ndata,x,y);
3193 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3194 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3195 regress_1d_opts_set_initial_parameters(regopts,params);
3196
3197 c3opt_add_objective(optimizer,param_RLSD2regress_cost,regopts);
3198 double * deriv_diff = calloc_double(nparams);
3199 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3200 /* for (size_t ii = 0; ii < nparams; ii++){ */
3201 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3202 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3203 /* } */
3204 /* printf("gerr = %G\n",gerr); */
3205 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3206 free(deriv_diff); deriv_diff = NULL;
3207
3208 int info;
3209 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3210 CuAssertIntEquals(tc,1,info>-1);
3211
3212 double * xtest = linspace(-1.0,1.0,1000);
3213 double * vals = calloc_double(1000);
3214 regress_func(1000,xtest,vals);
3215 size_t ii;
3216 double err = 0.0;
3217 double norm = 0.0;
3218 for (ii = 0; ii < 1000; ii++){
3219 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3220 norm += vals[ii]*vals[ii];
3221 }
3222 err = sqrt(err);
3223 norm = sqrt(norm);
3224 double rat = err/norm;
3225 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3226 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3227 free(xtest); xtest = NULL;
3228 free(vals); vals = NULL;
3229
3230 free(x); x = NULL;
3231 free(y); y = NULL;
3232 ope_opts_free(aopts); aopts = NULL;
3233 regress_1d_opts_destroy(regopts); regopts = NULL;
3234 c3opt_free(optimizer); optimizer = NULL;
3235 generic_function_free(gf); gf = NULL;;
3236 }
3237
Test_RLSD2_herm_regress(CuTest * tc)3238 void Test_RLSD2_herm_regress(CuTest * tc){
3239
3240 printf("Testing functions: second deriv reg with hermite\n");
3241
3242 size_t nparams = 10;
3243 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3244 0.6,-0.9,-.2,0.04,1.2};
3245 struct OpeOpts * aopts = ope_opts_alloc(HERMITE);
3246 ope_opts_set_nparams(aopts,nparams);
3247
3248 // create data
3249 size_t ndata = 5000;
3250 double * x = linspace(-1,1,ndata);
3251 double * y = calloc_double(ndata);
3252 regress_func(ndata,x,y);
3253 // // add noise
3254 for (size_t ii =0 ; ii < ndata; ii++){
3255 y[ii] += randn()*0.01;
3256 }
3257
3258 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3259 c3opt_set_verbose(optimizer,0);
3260
3261 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSD2,ndata,x,y);
3262 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3263 regress_1d_opts_set_regularization_penalty(regopts,1e-4/sqrt(ndata));
3264 regress_1d_opts_set_initial_parameters(regopts,params);
3265
3266 c3opt_add_objective(optimizer,param_RLSD2regress_cost,regopts);
3267 double * deriv_diff = calloc_double(nparams);
3268 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3269 /* for (size_t ii = 0; ii < nparams; ii++){ */
3270 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3271 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3272 /* } */
3273 /* printf("gerr = %G\n",gerr); */
3274 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3275 free(deriv_diff); deriv_diff = NULL;
3276
3277 int info;
3278 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3279 CuAssertIntEquals(tc,1,info>-1);
3280
3281 double * xtest = linspace(-1.0,1.0,1000);
3282 double * vals = calloc_double(1000);
3283 regress_func(1000,xtest,vals);
3284 size_t ii;
3285 double err = 0.0;
3286 double norm = 0.0;
3287 for (ii = 0; ii < 1000; ii++){
3288 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3289 norm += vals[ii]*vals[ii];
3290 }
3291 err = sqrt(err);
3292 norm = sqrt(norm);
3293 double rat = err/norm;
3294 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3295 printf("\t note: must have smaller regularization\n");
3296 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3297 free(xtest); xtest = NULL;
3298 free(vals); vals = NULL;
3299
3300 free(x); x = NULL;
3301 free(y); y = NULL;
3302 ope_opts_free(aopts); aopts = NULL;
3303 regress_1d_opts_destroy(regopts); regopts = NULL;
3304 c3opt_free(optimizer); optimizer = NULL;
3305 generic_function_free(gf); gf = NULL;;
3306 }
3307
Test_RLSRKHS_alg_cheb_regress(CuTest * tc)3308 void Test_RLSRKHS_alg_cheb_regress(CuTest * tc){
3309
3310 printf("Testing functions: RKHS regularization with chebyshev, algebraic decay\n");
3311
3312 size_t nparams = 10;
3313 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3314 0.6,-0.9,-.2,0.04,1.2};
3315 struct OpeOpts * aopts = ope_opts_alloc(CHEBYSHEV);
3316 ope_opts_set_nparams(aopts,nparams);
3317
3318 // create data
3319 size_t ndata = 5000;
3320 double * x = linspace(-1,1,ndata);
3321 double * y = calloc_double(ndata);
3322 regress_func(ndata,x,y);
3323 // // add noise
3324 for (size_t ii =0 ; ii < ndata; ii++){
3325 y[ii] += randn()*0.01;
3326 }
3327
3328 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3329 c3opt_set_verbose(optimizer,0);
3330
3331 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3332 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3333 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3334 regress_1d_opts_set_RKHS_decay_rate(regopts,ALGEBRAIC,0.9);
3335 regress_1d_opts_set_initial_parameters(regopts,params);
3336
3337 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3338 double * deriv_diff = calloc_double(nparams);
3339 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3340 /* for (size_t ii = 0; ii < nparams; ii++){ */
3341 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3342 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3343 /* } */
3344 /* printf("gerr = %G\n",gerr); */
3345 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3346 free(deriv_diff); deriv_diff = NULL;
3347
3348
3349 int info;
3350 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3351 CuAssertIntEquals(tc,1,info>-1);
3352
3353 double * xtest = linspace(-1.0,1.0,1000);
3354 double * vals = calloc_double(1000);
3355 regress_func(1000,xtest,vals);
3356 size_t ii;
3357 double err = 0.0;
3358 double norm = 0.0;
3359 for (ii = 0; ii < 1000; ii++){
3360 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3361 norm += vals[ii]*vals[ii];
3362 }
3363 err = sqrt(err);
3364 norm = sqrt(norm);
3365 double rat = err/norm;
3366 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3367 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3368 free(xtest); xtest = NULL;
3369 free(vals); vals = NULL;
3370
3371 free(x); x = NULL;
3372 free(y); y = NULL;
3373 ope_opts_free(aopts); aopts = NULL;
3374 regress_1d_opts_destroy(regopts); regopts = NULL;
3375 c3opt_free(optimizer); optimizer = NULL;
3376 generic_function_free(gf); gf = NULL;;
3377 }
3378
Test_RLSRKHS_exp_cheb_regress(CuTest * tc)3379 void Test_RLSRKHS_exp_cheb_regress(CuTest * tc){
3380
3381 printf("Testing functions: RKHS regularization with chebyshev, exponential decay\n");
3382
3383 size_t nparams = 10;
3384 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3385 0.6,-0.9,-.2,0.04,1.2};
3386 struct OpeOpts * aopts = ope_opts_alloc(CHEBYSHEV);
3387 ope_opts_set_nparams(aopts,nparams);
3388
3389 // create data
3390 size_t ndata = 5000;
3391 double * x = linspace(-1,1,ndata);
3392 double * y = calloc_double(ndata);
3393 regress_func(ndata,x,y);
3394 // // add noise
3395 for (size_t ii =0 ; ii < ndata; ii++){
3396 y[ii] += randn()*0.01;
3397 }
3398
3399 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3400 c3opt_set_verbose(optimizer,0);
3401
3402 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3403 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3404 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3405 regress_1d_opts_set_RKHS_decay_rate(regopts,EXPONENTIAL,2.0);
3406 regress_1d_opts_set_initial_parameters(regopts,params);
3407
3408 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3409 double * deriv_diff = calloc_double(nparams);
3410 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3411 /* for (size_t ii = 0; ii < nparams; ii++){ */
3412 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3413 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3414 /* } */
3415 /* printf("gerr = %G\n",gerr); */
3416 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3417 free(deriv_diff); deriv_diff = NULL;
3418
3419 int info;
3420 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3421 CuAssertIntEquals(tc,1,info>-1);
3422
3423 double * xtest = linspace(-1.0,1.0,1000);
3424 double * vals = calloc_double(1000);
3425 regress_func(1000,xtest,vals);
3426 size_t ii;
3427 double err = 0.0;
3428 double norm = 0.0;
3429 for (ii = 0; ii < 1000; ii++){
3430 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3431 norm += vals[ii]*vals[ii];
3432 }
3433 err = sqrt(err);
3434 norm = sqrt(norm);
3435 double rat = err/norm;
3436 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3437 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3438 free(xtest); xtest = NULL;
3439 free(vals); vals = NULL;
3440
3441 free(x); x = NULL;
3442 free(y); y = NULL;
3443 ope_opts_free(aopts); aopts = NULL;
3444 regress_1d_opts_destroy(regopts); regopts = NULL;
3445 c3opt_free(optimizer); optimizer = NULL;
3446 generic_function_free(gf); gf = NULL;;
3447 }
3448
3449
Test_RLSRKHS_alg_leg_regress(CuTest * tc)3450 void Test_RLSRKHS_alg_leg_regress(CuTest * tc){
3451
3452 printf("Testing functions: RKHS regularization with legendre, algebraic decay\n");
3453
3454 size_t nparams = 10;
3455 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3456 0.6,-0.9,-.2,0.04,1.2};
3457 struct OpeOpts * aopts = ope_opts_alloc(LEGENDRE);
3458 ope_opts_set_nparams(aopts,nparams);
3459
3460 // create data
3461 size_t ndata = 5000;
3462 double * x = linspace(-1,1,ndata);
3463 double * y = calloc_double(ndata);
3464 regress_func(ndata,x,y);
3465 // // add noise
3466 for (size_t ii =0 ; ii < ndata; ii++){
3467 y[ii] += randn()*0.01;
3468 }
3469
3470 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3471 c3opt_set_verbose(optimizer,0);
3472
3473 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3474 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3475 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3476 regress_1d_opts_set_RKHS_decay_rate(regopts,ALGEBRAIC,0.9);
3477 regress_1d_opts_set_initial_parameters(regopts,params);
3478
3479 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3480 double * deriv_diff = calloc_double(nparams);
3481 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3482 /* for (size_t ii = 0; ii < nparams; ii++){ */
3483 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3484 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3485 /* } */
3486 /* printf("gerr = %G\n",gerr); */
3487 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3488 free(deriv_diff); deriv_diff = NULL;
3489
3490 int info;
3491 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3492 CuAssertIntEquals(tc,1,info>-1);
3493
3494 double * xtest = linspace(-1.0,1.0,1000);
3495 double * vals = calloc_double(1000);
3496 regress_func(1000,xtest,vals);
3497 size_t ii;
3498 double err = 0.0;
3499 double norm = 0.0;
3500 for (ii = 0; ii < 1000; ii++){
3501 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3502 norm += vals[ii]*vals[ii];
3503 }
3504 err = sqrt(err);
3505 norm = sqrt(norm);
3506 double rat = err/norm;
3507 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3508 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3509 free(xtest); xtest = NULL;
3510 free(vals); vals = NULL;
3511
3512 free(x); x = NULL;
3513 free(y); y = NULL;
3514 ope_opts_free(aopts); aopts = NULL;
3515 regress_1d_opts_destroy(regopts); regopts = NULL;
3516 c3opt_free(optimizer); optimizer = NULL;
3517 generic_function_free(gf); gf = NULL;;
3518 }
3519
Test_RLSRKHS_exp_leg_regress(CuTest * tc)3520 void Test_RLSRKHS_exp_leg_regress(CuTest * tc){
3521
3522 printf("Testing functions: RKHS regularization with legendre, exponential decay\n");
3523
3524 size_t nparams = 10;
3525 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3526 0.6,-0.9,-.2,0.04,1.2};
3527 struct OpeOpts * aopts = ope_opts_alloc(LEGENDRE);
3528 ope_opts_set_nparams(aopts,nparams);
3529
3530 // create data
3531 size_t ndata = 5000;
3532 double * x = linspace(-1,1,ndata);
3533 double * y = calloc_double(ndata);
3534 regress_func(ndata,x,y);
3535 // // add noise
3536 for (size_t ii =0 ; ii < ndata; ii++){
3537 y[ii] += randn()*0.01;
3538 }
3539
3540 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3541 c3opt_set_verbose(optimizer,0);
3542
3543 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3544 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3545 regress_1d_opts_set_regularization_penalty(regopts,1e-2/sqrt(ndata));
3546 regress_1d_opts_set_RKHS_decay_rate(regopts,EXPONENTIAL,2.0);
3547 regress_1d_opts_set_initial_parameters(regopts,params);
3548
3549 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3550 double * deriv_diff = calloc_double(nparams);
3551 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3552 /* for (size_t ii = 0; ii < nparams; ii++){ */
3553 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3554 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3555 /* } */
3556 /* printf("gerr = %G\n",gerr); */
3557 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3558 free(deriv_diff); deriv_diff = NULL;
3559
3560 int info;
3561 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3562 CuAssertIntEquals(tc,1,info>-1);
3563
3564 double * xtest = linspace(-1.0,1.0,1000);
3565 double * vals = calloc_double(1000);
3566 regress_func(1000,xtest,vals);
3567 size_t ii;
3568 double err = 0.0;
3569 double norm = 0.0;
3570 for (ii = 0; ii < 1000; ii++){
3571 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3572 norm += vals[ii]*vals[ii];
3573 }
3574 err = sqrt(err);
3575 norm = sqrt(norm);
3576 double rat = err/norm;
3577 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3578 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3579 free(xtest); xtest = NULL;
3580 free(vals); vals = NULL;
3581
3582 free(x); x = NULL;
3583 free(y); y = NULL;
3584 ope_opts_free(aopts); aopts = NULL;
3585 regress_1d_opts_destroy(regopts); regopts = NULL;
3586 c3opt_free(optimizer); optimizer = NULL;
3587 generic_function_free(gf); gf = NULL;;
3588 }
3589
3590
Test_RLSRKHS_alg_herm_regress(CuTest * tc)3591 void Test_RLSRKHS_alg_herm_regress(CuTest * tc){
3592
3593 printf("Testing functions: RKHS regularization with hermite, algebraic decay\n");
3594
3595 size_t nparams = 10;
3596 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3597 0.6,-0.9,-.2,0.04,1.2};
3598 struct OpeOpts * aopts = ope_opts_alloc(HERMITE);
3599 ope_opts_set_nparams(aopts,nparams);
3600
3601 // create data
3602 size_t ndata = 5000;
3603 double * x = linspace(-1,1,ndata);
3604 double * y = calloc_double(ndata);
3605 regress_func(ndata,x,y);
3606 // // add noise
3607 for (size_t ii =0 ; ii < ndata; ii++){
3608 y[ii] += randn()*0.01;
3609 }
3610
3611 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3612 c3opt_set_verbose(optimizer,0);
3613
3614 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3615 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3616 regress_1d_opts_set_regularization_penalty(regopts,1e-4/sqrt(ndata));
3617 regress_1d_opts_set_RKHS_decay_rate(regopts,ALGEBRAIC,0.9);
3618 regress_1d_opts_set_initial_parameters(regopts,params);
3619
3620 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3621 double * deriv_diff = calloc_double(nparams);
3622 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3623 /* for (size_t ii = 0; ii < nparams; ii++){ */
3624 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3625 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3626 /* } */
3627 /* printf("gerr = %G\n",gerr); */
3628 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3629 free(deriv_diff); deriv_diff = NULL;
3630
3631 int info;
3632 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3633 CuAssertIntEquals(tc,1,info>-1);
3634
3635 double * xtest = linspace(-1.0,1.0,1000);
3636 double * vals = calloc_double(1000);
3637 regress_func(1000,xtest,vals);
3638 size_t ii;
3639 double err = 0.0;
3640 double norm = 0.0;
3641 for (ii = 0; ii < 1000; ii++){
3642 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3643 norm += vals[ii]*vals[ii];
3644 }
3645 err = sqrt(err);
3646 norm = sqrt(norm);
3647 double rat = err/norm;
3648 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3649 printf("\t note: must have smaller regularization\n");
3650 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3651 free(xtest); xtest = NULL;
3652 free(vals); vals = NULL;
3653
3654 free(x); x = NULL;
3655 free(y); y = NULL;
3656 ope_opts_free(aopts); aopts = NULL;
3657 regress_1d_opts_destroy(regopts); regopts = NULL;
3658 c3opt_free(optimizer); optimizer = NULL;
3659 generic_function_free(gf); gf = NULL;;
3660 }
3661
Test_RLSRKHS_exp_herm_regress(CuTest * tc)3662 void Test_RLSRKHS_exp_herm_regress(CuTest * tc){
3663
3664 printf("Testing functions: RKHS regularization with hermite, exponential decay\n");
3665
3666 size_t nparams = 10;
3667 double params[10] = {0.12,2.04,9.0,-0.2,0.4,
3668 0.6,-0.9,-.2,0.04,1.2};
3669 struct OpeOpts * aopts = ope_opts_alloc(HERMITE);
3670 ope_opts_set_nparams(aopts,nparams);
3671
3672 // create data
3673 size_t ndata = 5000;
3674 double * x = linspace(-1,1,ndata);
3675 double * y = calloc_double(ndata);
3676 regress_func(ndata,x,y);
3677 // // add noise
3678 for (size_t ii =0 ; ii < ndata; ii++){
3679 y[ii] += randn()*0.01;
3680 }
3681
3682 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
3683 c3opt_set_verbose(optimizer,0);
3684
3685 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,RLSRKHS,ndata,x,y);
3686 regress_1d_opts_set_parametric_form(regopts,POLYNOMIAL,aopts);
3687 regress_1d_opts_set_regularization_penalty(regopts,1e-4/sqrt(ndata));
3688 regress_1d_opts_set_RKHS_decay_rate(regopts,EXPONENTIAL,1.1);
3689 regress_1d_opts_set_initial_parameters(regopts,params);
3690
3691 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,regopts);
3692 double * deriv_diff = calloc_double(nparams);
3693 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
3694 /* for (size_t ii = 0; ii < nparams; ii++){ */
3695 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
3696 /* /\* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); *\/ */
3697 /* } */
3698 /* printf("gerr = %G\n",gerr); */
3699 CuAssertDblEquals(tc,0.0,gerr,1e-3);
3700 free(deriv_diff); deriv_diff = NULL;
3701
3702 int info;
3703 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
3704 CuAssertIntEquals(tc,1,info>-1);
3705
3706 double * xtest = linspace(-1.0,1.0,1000);
3707 double * vals = calloc_double(1000);
3708 regress_func(1000,xtest,vals);
3709 size_t ii;
3710 double err = 0.0;
3711 double norm = 0.0;
3712 for (ii = 0; ii < 1000; ii++){
3713 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
3714 norm += vals[ii]*vals[ii];
3715 }
3716 err = sqrt(err);
3717 norm = sqrt(norm);
3718 double rat = err/norm;
3719 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
3720 CuAssertDblEquals(tc, 0.0, rat, 1e-3);
3721 free(xtest); xtest = NULL;
3722 free(vals); vals = NULL;
3723
3724 free(x); x = NULL;
3725 free(y); y = NULL;
3726 ope_opts_free(aopts); aopts = NULL;
3727 regress_1d_opts_destroy(regopts); regopts = NULL;
3728 c3opt_free(optimizer); optimizer = NULL;
3729 generic_function_free(gf); gf = NULL;;
3730 }
3731
3732
3733
PolyRegressionSuite()3734 CuSuite * PolyRegressionSuite(){
3735
3736 CuSuite * suite = CuSuiteNew();
3737 SUITE_ADD_TEST(suite, Test_LS_cheb_regress);
3738 SUITE_ADD_TEST(suite, Test_LS_leg_regress);
3739 SUITE_ADD_TEST(suite, Test_LS_herm_regress);
3740
3741 SUITE_ADD_TEST(suite, Test_RLS2_cheb_regress);
3742 SUITE_ADD_TEST(suite, Test_RLS2_leg_regress);
3743 SUITE_ADD_TEST(suite, Test_RLS2_herm_regress);
3744
3745 /* SUITE_ADD_TEST(suite, Test_RLSD2_cheb_regress); */
3746 /* SUITE_ADD_TEST(suite, Test_RLSD2_leg_regress); */
3747 /* SUITE_ADD_TEST(suite, Test_RLSD2_herm_regress); */
3748
3749 SUITE_ADD_TEST(suite, Test_RLSRKHS_alg_cheb_regress);
3750 SUITE_ADD_TEST(suite, Test_RLSRKHS_exp_cheb_regress);
3751 SUITE_ADD_TEST(suite, Test_RLSRKHS_alg_leg_regress);
3752 SUITE_ADD_TEST(suite, Test_RLSRKHS_exp_leg_regress);
3753 SUITE_ADD_TEST(suite, Test_RLSRKHS_alg_herm_regress);
3754 SUITE_ADD_TEST(suite, Test_RLSRKHS_exp_herm_regress);
3755
3756 return suite;
3757 }
3758
3759
3760