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
38
39
40 #include <stdlib.h>
41 #include <stdio.h>
42 #include <math.h>
43 #include <string.h>
44 #include <assert.h>
45 #include <float.h>
46
47 #include "CuTest.h"
48 #include "testfunctions.h"
49
50 #include "array.h"
51 #include "lib_linalg.h"
52
53 #include "lib_funcs.h"
54
55
Test_gauss_eval(CuTest * tc)56 void Test_gauss_eval(CuTest * tc){
57
58 printf("Testing functions: gauss_kernel_eval(deriv) \n");
59
60 double scale = 1.2;
61 double width = 0.2;
62 double center = -0.3;
63
64 double h = 1e-8;
65 double x = 0.2;
66 double xh = x+h;
67 double valh = gauss_kernel_eval(scale,width*width,center,xh);
68 double x2h = x-h;
69 double val2h = gauss_kernel_eval(scale,width*width,center,x2h);
70
71 double numerical_deriv = (valh-val2h)/(2.0*h);
72 double analytical_deriv = gauss_kernel_deriv(scale,width*width,center,x);
73 CuAssertDblEquals(tc,numerical_deriv,analytical_deriv,1e-5);
74 }
75
Test_gauss_integrate(CuTest * tc)76 void Test_gauss_integrate(CuTest * tc){
77
78 printf("Testing functions: gauss_kernel_integrate\n");
79
80 double scale = 1.2;
81 double width = 0.2;
82 double center = -0.3;
83
84 double lb = -1;
85 double ub = 1;
86 size_t N = 100;
87 double * x = linspace(-1,1,N);
88 double numerical_int = 0.0;
89 for (size_t ii = 0; ii < N-1; ii++){
90 numerical_int += (x[ii+1]-x[ii])*gauss_kernel_eval(scale,width*width,center,x[ii]);
91 }
92 free(x); x = NULL;
93
94 double analytical_int = gauss_kernel_integrate(scale,width*width,center,lb,ub);
95 /* printf("%3.15G,%3.15G\n",numerical_int,analytical_int); */
96 CuAssertDblEquals(tc,numerical_int,analytical_int,1e-5);
97 }
98
Test_gauss_inner(CuTest * tc)99 void Test_gauss_inner(CuTest * tc){
100
101 printf("Testing functions: gauss_kernel_inner\n");
102
103 double s1 = 1.2;
104 double w1 = 0.2;
105 double c1 = -0.3;
106 double s2 = 0.3;
107 double w2 = 0.5;
108 double c2 = -0.1;
109
110 double lb = -10;
111 double ub = 10;
112 size_t N = 10000;
113 double * x = linspace(lb,ub,N);
114 double numerical_int = 0.0;
115 for (size_t ii = 0; ii < N-1; ii++){
116 double v1 = gauss_kernel_eval(s1,w1*w1,c1,x[ii]);
117 double v2 = gauss_kernel_eval(s2,w2*w2,c2,x[ii]);
118 numerical_int += (x[ii+1]-x[ii])*v1*v2;
119 }
120 free(x); x= NULL;
121 double analytical_int = gauss_kernel_inner(s1,w1*w1,c1,s2,w2*w2,c2,-DBL_MAX,DBL_MAX);
122 /* printf("%3.15G,%3.15G\n",numerical_int,analytical_int); */
123 CuAssertDblEquals(tc,numerical_int,analytical_int,1e-5);
124 }
125
126
Test_kernel_expansion_mem(CuTest * tc)127 void Test_kernel_expansion_mem(CuTest * tc)
128 {
129
130 printf("Testing functions: kernel allocation and deallocation memory \n");
131
132 double scale = 1.2;
133 double width = 0.2;
134 size_t N = 20;
135 double * centers = linspace(-1,1,N);
136
137 struct KernelExpansion * ke = kernel_expansion_alloc(1);
138 CuAssertIntEquals(tc,1,ke!=NULL);
139 for (size_t ii = 0; ii < N; ii++){
140 struct Kernel * kern = kernel_gaussian(scale,width*width,centers[ii]);
141 kernel_expansion_add_kernel(ke,0.0,kern);
142 kernel_free(kern); kern = NULL;
143 }
144
145 kernel_expansion_free(ke); ke = NULL;
146
147 free(centers); centers = NULL;
148 }
149
Test_kernel_expansion_copy(CuTest * tc)150 void Test_kernel_expansion_copy(CuTest * tc)
151 {
152
153 printf("Testing functions: kernel_expansion_copy \n");
154
155 double scale = 1.2;
156 double width = 0.2;
157 size_t N = 20;
158 double * centers = linspace(-1,1,N);
159
160 struct KernelExpansion * ke = kernel_expansion_alloc(1);
161 for (size_t ii = 0; ii < N; ii++){
162 struct Kernel * kern = kernel_gaussian(scale,width*width,centers[ii]);
163 kernel_expansion_add_kernel(ke,randu(),kern);
164 kernel_free(kern); kern = NULL;
165 }
166
167 struct KernelExpansion * ke2 = kernel_expansion_copy(ke);
168
169 double * x = linspace(-1,1,200);
170 for (size_t ii = 0; ii < 200; ii++){
171 double val1 = kernel_expansion_eval(ke,x[ii]);
172 double val2 = kernel_expansion_eval(ke2,x[ii]);
173 CuAssertDblEquals(tc,val1,val2,1e-15);
174 }
175
176 free(x); x = NULL;
177
178 kernel_expansion_free(ke2); ke2 = NULL;
179 kernel_expansion_free(ke); ke = NULL;
180
181 free(centers); centers = NULL;
182 }
183
Test_serialize_kernel_expansion(CuTest * tc)184 void Test_serialize_kernel_expansion(CuTest * tc){
185
186 printf("Testing functions: serialize_kernel_expansion \n");
187 double scale = 1.2;
188 double width = 0.2;
189 size_t N = 18;
190 double * centers = linspace(-1,1,N);
191
192 struct KernelExpansion * ke = kernel_expansion_alloc(1);
193 for (size_t ii = 0; ii < N; ii++){
194 struct Kernel * kern = kernel_gaussian(scale,width*width,centers[ii]);
195 kernel_expansion_add_kernel(ke,randu(),kern);
196 kernel_free(kern); kern = NULL;
197 }
198
199 unsigned char * text = NULL;
200 size_t size_to_be;
201 serialize_kernel_expansion(text, ke, &size_to_be);
202 text = malloc(size_to_be * sizeof(char));
203
204 serialize_kernel_expansion(text, ke, NULL);
205
206 struct KernelExpansion * k2 = NULL;
207 deserialize_kernel_expansion(text, &k2);
208 free(text); text = NULL;
209
210 CuAssertIntEquals(tc,N,kernel_expansion_get_nkernels(k2));
211
212 double * x = linspace(-1,1,200);
213 for (size_t ii = 0; ii < 200; ii++){
214 double val1 = kernel_expansion_eval(ke,x[ii]);
215 double val2 = kernel_expansion_eval(k2,x[ii]);
216 /* printf("%G, %G\n",val1,val2); */
217 CuAssertDblEquals(tc,val1,val2,1e-15);
218 }
219 free(x); x = NULL;
220
221 kernel_expansion_free(k2); k2 = NULL;
222 kernel_expansion_free(ke); ke = NULL;
223
224 free(centers); centers = NULL;
225 }
226
Test_serialize_generic_function_kernel(CuTest * tc)227 void Test_serialize_generic_function_kernel (CuTest * tc){
228
229 printf("Testing functions: (de)serializing generic_function with kernels\n");
230 double scale = 1.2;
231 double width = 0.2;
232 size_t N = 18;
233 double lb = -1.0, ub = 2.0;
234 double * centers = linspace(lb,ub,N);
235
236 struct KernelExpansion * ke = kernel_expansion_alloc(1);
237 for (size_t ii = 0; ii < N; ii++){
238 struct Kernel * kern = kernel_gaussian(scale,width*width,centers[ii]);
239 kernel_expansion_add_kernel(ke,randu(),kern);
240 kernel_free(kern); kern = NULL;
241 }
242 kernel_expansion_set_bounds(ke,lb,ub);
243
244 struct GenericFunction * pl = generic_function_alloc(1,KERNEL);
245 pl->f = ke;
246
247 unsigned char * text = NULL;
248 size_t size_to_be;
249 serialize_generic_function(text, pl, &size_to_be);
250 text = malloc(size_to_be * sizeof(char));
251
252 serialize_generic_function(text, pl, NULL);
253
254 struct GenericFunction * pt = NULL;
255 deserialize_generic_function(text, &pt);
256
257 double * xtest = linspace(lb,ub,1000);
258 size_t ii;
259 double err = 0.0;
260 for (ii = 0; ii < 1000; ii++){
261 err += pow(generic_function_1d_eval(pl,xtest[ii]) -
262 generic_function_1d_eval(pt,xtest[ii]),2);
263 }
264 err = sqrt(err);
265 CuAssertDblEquals(tc, 0.0, err, 1e-15);
266 free(xtest);
267 free(text);
268
269 generic_function_free(pl);
270 generic_function_free(pt);
271 free(centers); centers = NULL;
272 }
273
Test_kernel_expansion_create_with_params_and_grad(CuTest * tc)274 void Test_kernel_expansion_create_with_params_and_grad(CuTest * tc){
275
276 printf("Testing functions: kernel_expansion_create_with_params and grad\n");
277
278 double scale = 1.2;
279 double width = 0.2;
280 size_t N = 18;
281 double * centers = linspace(-1,1,N);
282 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(N,centers,scale,width);
283
284 double * params = calloc_double(N);
285 for (size_t ii = 0; ii < N; ii++) params[ii]=0.1*randn();
286
287 struct KernelExpansion * ke = kernel_expansion_create_with_params(opts,N,params);
288
289 double grad[36];
290 double xloc[2] = {0.4, 0.8};
291 int res = kernel_expansion_param_grad_eval(ke,2,xloc,grad);
292 CuAssertIntEquals(tc,0,res);
293
294
295 // numerical derivative
296 struct KernelExpansion * ke1 = NULL;
297 struct KernelExpansion * ke2 = NULL;
298
299 size_t dim = N;
300 double * x1 = calloc_double(dim);
301 double * x2 = calloc_double(dim);
302 for (size_t ii = 0; ii < dim; ii++){
303 x1[ii] = params[ii];
304 x2[ii] = params[ii];
305 }
306
307 double diff = 0.0;
308 double v1,v2;
309 double norm = 0.0;
310 double eps = 1e-8;
311 for (size_t zz = 0; zz < 2; zz++){
312 /* size_t zz = 0; */
313 for (size_t ii = 0; ii < dim; ii++){
314 x1[ii] += eps;
315 x2[ii] -= eps;
316 ke1 = kernel_expansion_create_with_params(opts,dim,x1);
317 v1 = kernel_expansion_eval(ke1,xloc[zz]);
318
319 ke2 = kernel_expansion_create_with_params(opts,dim,x2);
320 v2 = kernel_expansion_eval(ke2,xloc[zz]);
321
322 double diff_iter = pow( (v1-v2)/2.0/eps - grad[zz*dim+ii], 2 );
323 /* printf("current diff = %G\n",diff_iter); */
324 /* printf("\t norm = %G\n",grad[zz*dim+ii]); */
325 diff += diff_iter;
326 norm += pow( (v1-v2)/2.0/eps,2);
327
328 x1[ii] -= eps;
329 x2[ii] += eps;
330
331 kernel_expansion_free(ke1); ke1 = NULL;
332 kernel_expansion_free(ke2); ke2 = NULL;
333 }
334 if (norm > 1){
335 diff /= norm;
336 }
337 CuAssertDblEquals(tc,0.0,diff,1e-7);
338 }
339 free(x1); x1 = NULL;
340 free(x2); x2 = NULL;
341
342 kernel_approx_opts_free(opts); opts = NULL;
343 free(centers); centers = NULL;
344 free(params); params = NULL;
345 kernel_expansion_free(ke); ke = NULL;
346 }
347
Test_kernel_expansion_integrate(CuTest * tc)348 void Test_kernel_expansion_integrate(CuTest * tc)
349 {
350
351 printf("Testing functions: kernel_expansion_integrate \n");
352
353 double scale = 1.2;
354 double width = 0.2;
355 size_t N = 20;
356 double * centers = linspace(-1,1,N);
357
358 struct KernelExpansion * ke = kernel_expansion_alloc(1);
359 for (size_t ii = 0; ii < N; ii++){
360 struct Kernel * kern = kernel_gaussian(scale,width*width,centers[ii]);
361 kernel_expansion_add_kernel(ke,randu(),kern);
362 kernel_free(kern); kern = NULL;
363 }
364
365
366 size_t nint = 10000;
367 double * x = linspace(-20,20,nint);
368 double num_int = 0.0;
369 for (size_t ii = 0; ii < nint-1; ii++){
370 double val1 = kernel_expansion_eval(ke,x[ii]);
371 num_int += val1*(x[ii+1]-x[ii]);
372 }
373 free(x); x = NULL;
374
375 double ana_int = kernel_expansion_integrate(ke);
376
377 CuAssertDblEquals(tc,num_int,ana_int,1e-5);
378
379 kernel_expansion_free(ke); ke = NULL;
380 free(centers); centers = NULL;
381 }
382
Test_kernel_expansion_inner(CuTest * tc)383 void Test_kernel_expansion_inner(CuTest * tc)
384 {
385
386 printf("Testing functions: kernel_expansion_inner\n");
387
388 double s1 = 1.2;
389 double w1 = 0.2;
390 size_t N1 = 20;
391 double * c1 = linspace(-1,1,N1);
392 struct KernelExpansion * ke = kernel_expansion_alloc(1);
393 for (size_t ii = 0; ii < N1; ii++){
394 struct Kernel * kern = kernel_gaussian(s1,w1*w1,c1[ii]);
395 kernel_expansion_add_kernel(ke,randu(),kern);
396 kernel_free(kern); kern = NULL;
397 }
398
399 double s2 = 3.2;
400 double w2 = 0.8;
401 size_t N2 = 14;
402 double * c2 = linspace(-3,0.2,N2);
403 struct KernelExpansion * ke2 = kernel_expansion_alloc(1);
404 for (size_t ii = 0; ii < N2; ii++){
405 struct Kernel * kern = kernel_gaussian(s2,w2*w2,c2[ii]);
406 kernel_expansion_add_kernel(ke2,randu(),kern);
407 kernel_free(kern); kern = NULL;
408 }
409
410 double ana_int = kernel_expansion_inner(ke,ke2);
411
412 size_t nint = 500000;
413 double * x = linspace(-20,20,nint);
414 double num_int = 0.0;
415 for (size_t ii = 0; ii < nint-1; ii++){
416 double val1 = kernel_expansion_eval(ke,x[ii]);
417 double val2 = kernel_expansion_eval(ke2,x[ii]);
418 num_int += val1*val2*(x[ii+1]-x[ii]);
419 }
420 free(x); x = NULL;
421
422 /* printf("%3.15G,%3.15G,%3.15G\n",num_int,ana_int,num_int-ana_int); */
423 CuAssertDblEquals(tc,num_int,ana_int,1e-13);
424
425 kernel_expansion_free(ke); ke = NULL;
426 kernel_expansion_free(ke2); ke2 = NULL;
427
428 free(c1); c1 = NULL;
429 free(c2); c2 = NULL;
430 }
431
Test_kernel_expansion_orth_basis(CuTest * tc)432 void Test_kernel_expansion_orth_basis(CuTest * tc)
433 {
434
435 printf("Testing functions: kernel_expansion_orth_basis\n");
436
437 double scale = 1.2;
438
439 size_t N = 30;
440 double width = pow((double) N, -0.2) * 1.0/12.0 * 2;
441 /* width *= 0.1; */
442 /* printf("width = %G\n",width); */
443 double * centers = linspace(-1,1,N);
444 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(N,centers,scale,width);
445 kernel_approx_opts_set_lb(opts,-10.0);
446 kernel_approx_opts_set_ub(opts,10.0);
447
448 size_t north = 30;
449 struct KernelExpansion * ke[200];
450 for (size_t ii = 0; ii < north; ii++){
451 ke[ii] = NULL;
452 }
453
454 kernel_expansion_orth_basis(north,ke,opts);
455 for (size_t ii = 0; ii < north; ii++){
456 /* printf("\n"); */
457 /* printf("ii = %zu\n",ii); */
458 for (size_t jj = 0; jj < north; jj++){
459 double val = kernel_expansion_inner(ke[ii],ke[jj]);
460 if (ii == jj){
461 /* printf("should be 1: %3.15G\n",val); */
462 CuAssertDblEquals(tc,1.0,val,1e-12);
463 }
464 else{
465 /* printf("should be 0: %3.15G\n",val); */
466 CuAssertDblEquals(tc,0.0,val,1e-13);
467 }
468 }
469 }
470
471
472 for (size_t ii = 0; ii < north; ii++){
473 kernel_expansion_free(ke[ii]); ke[ii] = NULL;
474 }
475
476 kernel_approx_opts_free(opts); opts = NULL;
477 free(centers); centers = NULL;
478 }
479
regress_func(size_t N,const double * x,double * out)480 static void regress_func(size_t N, const double * x, double * out)
481 {
482 for (size_t ii = 0; ii < N; ii++){
483 out[ii] = -1.0 * sin(2.0*x[ii]);
484 }
485 }
486
Test_kernel_LSregress(CuTest * tc)487 void Test_kernel_LSregress(CuTest * tc){
488
489 printf("Testing functions: least squares regression with kernel expansion\n");
490
491 // create data
492 size_t ndata = 200;
493 double * x = linspace(-1,1,ndata);
494 double * y = calloc_double(ndata);
495 regress_func(ndata,x,y);
496 // // add noise
497 double * params = calloc_double(ndata);
498 for (size_t ii =0 ; ii < ndata; ii++){
499 y[ii] += randn()*0.01;
500 params[ii] = y[ii];
501 }
502
503
504
505 size_t nparams = ndata;
506 double scale = 0.1;
507 double width = pow(ndata,-0.2)/12.0;
508 width *= 5.0;
509 /* printf("width = %G\n",width); */
510 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(nparams,x,scale,width);
511
512 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
513 c3opt_set_verbose(optimizer,0);
514 c3opt_set_maxiter(optimizer,2000);
515 c3opt_set_gtol(optimizer,1e-6);
516
517 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
518 regress_1d_opts_set_parametric_form(regopts,KERNEL,opts);
519 regress_1d_opts_set_initial_parameters(regopts,params);
520
521 // check derivative
522 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
523 double * deriv_diff = calloc_double(nparams);
524 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
525 for (size_t ii = 0; ii < nparams; ii++){
526 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
527 /* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); */
528 }
529 CuAssertDblEquals(tc,0.0,gerr,1e-3);
530 free(deriv_diff); deriv_diff = NULL;
531
532 int info;
533 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
534 /* CuAssertIntEquals(tc,1,info>-1); */
535
536 /* print_generic_function(gf,5,NULL); */
537
538 double * xtest = linspace(-1.0,1.0,1000);
539 double * vals = calloc_double(1000);
540 regress_func(1000,xtest,vals);
541 size_t ii;
542 double err = 0.0;
543 double norm = 0.0;
544 for (ii = 0; ii < 1000; ii++){
545 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
546 norm += vals[ii]*vals[ii];
547 }
548 err = sqrt(err);
549 norm = sqrt(norm);
550 double rat = err/norm;
551 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
552 CuAssertDblEquals(tc, 0.0, rat, 1e-2);
553 free(xtest); xtest = NULL;
554 free(vals); vals = NULL;
555
556 free(x); x = NULL;
557 free(y); y = NULL;
558 kernel_approx_opts_free(opts); opts = NULL;
559 regress_1d_opts_destroy(regopts); regopts = NULL;
560 c3opt_free(optimizer); optimizer = NULL;
561 generic_function_free(gf); gf = NULL;;
562 free(params); params = NULL;
563 }
564
Test_kernel_LSregress2(CuTest * tc)565 void Test_kernel_LSregress2(CuTest * tc){
566
567 printf("Testing functions: least squares regression with kernel expansion (2)\n");
568
569 // create data
570 size_t ndata = 200;
571 double * x = linspace(-1,1,ndata);
572 double * y = calloc_double(ndata);
573 regress_func(ndata,x,y);
574 // // add noise
575
576 for (size_t ii =0 ; ii < ndata; ii++){
577 y[ii] += randn()*0.01;
578 }
579
580 size_t nparams = 40;
581 double * params = linspace(-1.0,1.0,nparams);
582
583
584 double scale = 0.1;
585 double width = pow(nparams,-0.2)/12.0;
586 width *= 5.0;
587 /* printf("width = %G\n",width); */
588 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(nparams,params,scale,width);
589
590 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
591 c3opt_set_verbose(optimizer,0);
592 c3opt_set_maxiter(optimizer,2000);
593 c3opt_set_gtol(optimizer,1e-6);
594
595 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
596 regress_1d_opts_set_parametric_form(regopts,KERNEL,opts);
597 regress_1d_opts_set_initial_parameters(regopts,params);
598
599 // check derivative
600 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
601 double * deriv_diff = calloc_double(nparams);
602 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
603 for (size_t ii = 0; ii < nparams; ii++){
604 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
605 /* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); */
606 }
607 CuAssertDblEquals(tc,0.0,gerr,1e-3);
608 free(deriv_diff); deriv_diff = NULL;
609
610 int info;
611 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
612 /* CuAssertIntEquals(tc,1,info>-1); */
613
614 /* print_generic_function(gf,5,NULL); */
615
616 double * xtest = linspace(-1.0,1.0,1000);
617 double * vals = calloc_double(1000);
618 regress_func(1000,xtest,vals);
619 size_t ii;
620 double err = 0.0;
621 double norm = 0.0;
622 for (ii = 0; ii < 1000; ii++){
623 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
624 norm += vals[ii]*vals[ii];
625 }
626 err = sqrt(err);
627 norm = sqrt(norm);
628 double rat = err/norm;
629 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
630 CuAssertDblEquals(tc, 0.0, rat, 1e-2);
631 free(xtest); xtest = NULL;
632 free(vals); vals = NULL;
633
634 free(x); x = NULL;
635 free(y); y = NULL;
636 kernel_approx_opts_free(opts); opts = NULL;
637 regress_1d_opts_destroy(regopts); regopts = NULL;
638 c3opt_free(optimizer); optimizer = NULL;
639 generic_function_free(gf); gf = NULL;;
640 free(params); params = NULL;
641 }
642
Test_kernel_LSregress_with_centers(CuTest * tc)643 void Test_kernel_LSregress_with_centers(CuTest * tc){
644
645 printf("Testing functions: least squares regression with kernel expansion and moving centers\n");
646
647 // create data
648 size_t ndata = 200;
649 double * x = linspace(-1,1,ndata);
650 double * y = calloc_double(ndata);
651 regress_func(ndata,x,y);
652 // // add noise
653
654 for (size_t ii =0 ; ii < ndata; ii++){
655 y[ii] += randn()*0.01;
656 }
657
658 size_t nparams = 2*ndata;
659 double * params = calloc_double(2*ndata);
660 memmove(params,y,ndata*sizeof(double));
661 memmove(params+ndata,x,ndata*sizeof(double));
662
663 double scale = 0.1;
664 double width = pow(ndata,-0.2)/12.0;
665 width *= 5.0;
666 /* printf("width = %G\n",width); */
667 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(nparams/2,x,scale,width);
668 kernel_approx_opts_set_center_adapt(opts,1);
669
670 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
671 c3opt_set_verbose(optimizer,0);
672 c3opt_set_maxiter(optimizer,2000);
673 c3opt_set_gtol(optimizer,1e-6);
674
675 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
676 regress_1d_opts_set_parametric_form(regopts,KERNEL,opts);
677 regress_1d_opts_set_initial_parameters(regopts,params);
678
679 // check derivative
680 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
681 double * deriv_diff = calloc_double(nparams);
682 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
683 for (size_t ii = 0; ii < nparams; ii++){
684 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
685 /* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); */
686 }
687 CuAssertDblEquals(tc,0.0,gerr,1e-3);
688 free(deriv_diff); deriv_diff = NULL;
689
690 int info;
691 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
692 /* CuAssertIntEquals(tc,1,info>-1); */
693
694 /* print_generic_function(gf,5,NULL); */
695
696 double * xtest = linspace(-1.0,1.0,1000);
697 double * vals = calloc_double(1000);
698 regress_func(1000,xtest,vals);
699 size_t ii;
700 double err = 0.0;
701 double norm = 0.0;
702 for (ii = 0; ii < 1000; ii++){
703 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
704 norm += vals[ii]*vals[ii];
705 }
706 err = sqrt(err);
707 norm = sqrt(norm);
708 double rat = err/norm;
709 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
710 CuAssertDblEquals(tc, 0.0, rat, 1e-2);
711 free(xtest); xtest = NULL;
712 free(vals); vals = NULL;
713
714 free(x); x = NULL;
715 free(y); y = NULL;
716 kernel_approx_opts_free(opts); opts = NULL;
717 regress_1d_opts_destroy(regopts); regopts = NULL;
718 c3opt_free(optimizer); optimizer = NULL;
719 generic_function_free(gf); gf = NULL;;
720 free(params); params = NULL;
721 }
722
723
Test_kernel_LSregress_with_centers2(CuTest * tc)724 void Test_kernel_LSregress_with_centers2(CuTest * tc){
725
726 printf("Testing functions: least squares regression with kernel expansion and moving centers (2)\n");
727
728 // create data
729 size_t ndata = 200;
730 double * x = linspace(-1,1,ndata);
731 double * y = calloc_double(ndata);
732 regress_func(ndata,x,y);
733 // // add noise
734
735 for (size_t ii =0 ; ii < ndata; ii++){
736 y[ii] += randn()*0.01;
737 }
738
739 size_t nparams = 40;
740 double * params = calloc_double(nparams);
741 for (size_t ii = 0; ii < nparams; ii++){
742 params[ii] = randu()*2.0-1.0;
743 }
744
745 double scale = 0.1;
746 double width = pow(nparams,-0.2)/12.0;
747 width *= 5.0;
748 /* printf("width = %G\n",width); */
749 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(nparams/2,params+nparams/2,scale,width);
750 kernel_approx_opts_set_center_adapt(opts,1);
751
752 struct c3Opt * optimizer = c3opt_alloc(BFGS,nparams);
753 c3opt_set_verbose(optimizer,0);
754 c3opt_set_maxiter(optimizer,2000);
755 c3opt_set_gtol(optimizer,1e-6);
756
757 struct Regress1DOpts * regopts = regress_1d_opts_create(PARAMETRIC,LS,ndata,x,y);
758 regress_1d_opts_set_parametric_form(regopts,KERNEL,opts);
759 regress_1d_opts_set_initial_parameters(regopts,params);
760
761 // check derivative
762 c3opt_add_objective(optimizer,param_LSregress_cost,regopts);
763 double * deriv_diff = calloc_double(nparams);
764 double gerr = c3opt_check_deriv_each(optimizer,params,1e-8,deriv_diff);
765 for (size_t ii = 0; ii < nparams; ii++){
766 /* printf("ii = %zu, diff=%G\n",ii,deriv_diff[ii]); */
767 /* CuAssertDblEquals(tc,0.0,deriv_diff[ii],1e-3); */
768 }
769 CuAssertDblEquals(tc,0.0,gerr,1e-3);
770 free(deriv_diff); deriv_diff = NULL;
771
772 int info;
773 struct GenericFunction * gf = generic_function_regress1d(regopts,optimizer,&info);
774 /* CuAssertIntEquals(tc,1,info>-1); */
775
776 /* print_generic_function(gf,5,NULL); */
777
778 double * xtest = linspace(-1.0,1.0,1000);
779 double * vals = calloc_double(1000);
780 regress_func(1000,xtest,vals);
781 size_t ii;
782 double err = 0.0;
783 double norm = 0.0;
784 for (ii = 0; ii < 1000; ii++){
785 err += pow(generic_function_1d_eval(gf,xtest[ii]) - vals[ii],2);
786 norm += vals[ii]*vals[ii];
787 }
788 err = sqrt(err);
789 norm = sqrt(norm);
790 double rat = err/norm;
791 printf("\t error = %G, norm=%G, rat=%G\n",err,norm,rat);
792 CuAssertDblEquals(tc, 0.0, rat, 1e-2);
793 free(xtest); xtest = NULL;
794 free(vals); vals = NULL;
795
796 free(x); x = NULL;
797 free(y); y = NULL;
798 kernel_approx_opts_free(opts); opts = NULL;
799 regress_1d_opts_destroy(regopts); regopts = NULL;
800 c3opt_free(optimizer); optimizer = NULL;
801 generic_function_free(gf); gf = NULL;;
802 free(params); params = NULL;
803 }
804
Test_kernel_linear(CuTest * tc)805 void Test_kernel_linear(CuTest * tc){
806
807 printf("Testing function: kernel_expansion_linear\n");
808
809
810 double lb = -1.0;
811 double ub = 1.0;
812 size_t nparams = 10;
813 double * x = linspace(lb,ub,nparams);
814 double scale = 1.0;
815 double width = pow(nparams,-0.2)/12.0 * (ub-lb);
816 width *= 20.0;
817 /* printf("width = %G\n",width); */
818 struct KernelApproxOpts * opts = kernel_approx_opts_gauss(nparams,x,scale,width);
819
820 double a = 2.0, offset=3.0;
821
822 struct KernelExpansion * ke = kernel_expansion_linear(a,offset,opts);
823
824 size_t N = 100;
825 double * pts = linspace(lb,ub,N);
826 size_t ii;
827 for (ii = 0; ii < N; ii++){
828 double eval1 = kernel_expansion_eval(ke,pts[ii]);
829 double eval2 = a*pts[ii] + offset;
830 double diff= fabs(eval1-eval2);
831 /* printf("eval1 = %G, eval2=%G, diff=%G\n",eval1,eval2,diff); */
832 CuAssertDblEquals(tc, 0.0, diff, 1e-5);
833 }
834 free(pts); pts = NULL;
835
836 /* print_kernel_expansion(ke,5,NULL); */
837
838 kernel_expansion_free(ke); ke = NULL;
839 free(x); x = NULL;
840 kernel_approx_opts_free(opts); opts = NULL;
841 }
842
843
KernGetSuite()844 CuSuite * KernGetSuite(){
845
846 CuSuite * suite = CuSuiteNew();
847 SUITE_ADD_TEST(suite, Test_gauss_eval);
848 SUITE_ADD_TEST(suite, Test_gauss_integrate);
849 SUITE_ADD_TEST(suite, Test_gauss_inner);
850 SUITE_ADD_TEST(suite, Test_kernel_expansion_mem);
851 SUITE_ADD_TEST(suite, Test_kernel_expansion_copy);
852 SUITE_ADD_TEST(suite, Test_serialize_kernel_expansion);
853 SUITE_ADD_TEST(suite, Test_serialize_generic_function_kernel);
854 SUITE_ADD_TEST(suite, Test_kernel_expansion_create_with_params_and_grad);
855 SUITE_ADD_TEST(suite, Test_kernel_expansion_integrate);
856 SUITE_ADD_TEST(suite, Test_kernel_expansion_inner);
857 SUITE_ADD_TEST(suite, Test_kernel_expansion_orth_basis);
858 SUITE_ADD_TEST(suite, Test_kernel_LSregress);
859 SUITE_ADD_TEST(suite, Test_kernel_LSregress2);
860 SUITE_ADD_TEST(suite, Test_kernel_LSregress_with_centers);
861 SUITE_ADD_TEST(suite, Test_kernel_LSregress_with_centers2);
862 SUITE_ADD_TEST(suite, Test_kernel_linear);
863 return suite;
864 }
865