1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3 // Copyright (c) 2017 NTESS, LLC.
4
5 // This file is part of the Compressed Continuous Computation (C3) Library
6 // Author: Alex A. Gorodetsky
7 // Contact: alex@alexgorodetsky.com
8
9 // All rights reserved.
10
11 // Redistribution and use in source and binary forms, with or without modification,
12 // are permitted provided that the following conditions are met:
13
14 // 1. Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16
17 // 2. Redistributions in binary form must reproduce the above copyright notice,
18 // this list of conditions and the following disclaimer in the documentation
19 // and/or other materials provided with the distribution.
20
21 // 3. Neither the name of the copyright holder nor the names of its contributors
22 // may be used to endorse or promote products derived from this software
23 // without specific prior written permission.
24
25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 //Code
37
38
39 /** \file functions.c
40 * Provides basic routines for interfacing specific functions to the outside world through
41 * generic functions
42 */
43
44 #include <stdlib.h>
45 #include <stdio.h>
46 #include <string.h>
47 #include <math.h>
48 #include <assert.h>
49
50 #include "stringmanip.h"
51 #include "functions.h"
52 #include "futil.h"
53
54 /** \struct Regress1DOpts
55 * \brief One dimensional regression options
56 * \var Regress1DOpts:atype
57 * approximation type
58 * \var Regress1DOpts:rtype
59 * regression problem
60 * \var Regress1DOpts:fc
61 * function class of the approximation
62 * \var Regress1DOpts:reg_param_set
63 * indicator of whethe the regularization parameter is set
64 * \var Regress1DOpts:lambda
65 * regularization parameter
66 * \var Regress1DOpts:decay_type
67 * decay type (used for regularized RKHS regression)
68 * \var Regress1DOpts:coeff_decay_opt
69 * parameter specifying decay rate
70 * \var Regress1DOpts:N
71 * number of training samples
72 * \var Regress1DOpts:x
73 * location of training samples
74 * \var Regress1DOpts:y
75 * value of training samples
76 * \var Regress1DOpts:aopts
77 * approximation options for function class
78 * \var Regress1DOpts:nparam
79 * number of parameters for parametric regression
80 * \var Regress1DOpts:init_param
81 * initial parameters
82 * \var Regress1DOpts:gf
83 * Generic function currently being worked with
84 * \var Regress1DOpts:eval
85 * Storage locations for evaluation of current guess
86 * \var Regress1DOpts:grad
87 * Storage location for gradient
88 * \var Regress1DOpts:resid
89 * Storage location for residual
90 */
91 struct Regress1DOpts
92 {
93 enum approx_type atype;
94 enum regress_type rtype;
95 enum function_class fc;
96
97 // Regularization options
98 int reg_param_set;
99 double lambda;
100 enum coeff_decay_type decay_type;
101 double coeff_decay_param;
102
103 size_t N;
104 const double * x;
105 const double * y;
106
107 void * aopts; // approximation options
108
109 // parameteric stuff
110 size_t nparam; // for parametric
111 const double * init_param;
112
113 // store current generic funciton
114 struct GenericFunction * gf;
115
116 // stuff to speed up storage
117 double * eval;
118 double * grad;
119 double * resid;
120
121 };
122
123
124 /********************************************************//**
125 Allocate memory for a generic function without specifying class or sub_type
126
127 \param[in] dim - dimension of functions
128
129 \return out - generic function
130 ************************************************************/
generic_function_alloc_base(size_t dim)131 struct GenericFunction * generic_function_alloc_base(size_t dim)
132 {
133 struct GenericFunction * out;
134 if (NULL == ( out = malloc(sizeof(struct GenericFunction)))){
135 fprintf(stderr, "failed to allocate for a generic function.\n");
136 exit(1);
137 }
138 out->dim = dim;
139 out->fargs = NULL;
140 out->f = NULL;
141 return out;
142 }
143
144 /********************************************************//**
145 Allocate memory for a generic function array
146
147 \param[in] size - size of array
148
149 \return out - generic function
150 ************************************************************/
generic_function_array_alloc(size_t size)151 struct GenericFunction ** generic_function_array_alloc(size_t size)
152 {
153 struct GenericFunction ** out;
154 if (NULL == ( out = malloc( size * sizeof(struct GenericFunction *)))){
155 fprintf(stderr, "failed to allocate a generic function array.\n");
156 exit(1);
157 }
158 size_t ii;
159 for (ii = 0; ii < size; ii++){
160 out[ii] = NULL;
161 }
162 return out;
163 }
164
165 /********************************************************//**
166 Allocate memory for a generic function of a particular
167 function class
168
169 \param[in] dim - dimension of functions
170 \param[in] fc - function class
171
172 \return generic function
173 ************************************************************/
174 struct GenericFunction *
generic_function_alloc(size_t dim,enum function_class fc)175 generic_function_alloc(size_t dim, enum function_class fc)
176 {
177
178 struct GenericFunction * out;
179 if (NULL == ( out = malloc(sizeof(struct GenericFunction)))){
180 fprintf(stderr, "failed to allocate for a generic function.\n");
181 exit(1);
182 }
183 out->dim = dim;
184 out->fc = fc;
185 out->f = NULL;
186 out->fargs = NULL;
187
188 return out;
189 }
190
191
192 /********************************************************//**
193 Free memory for generic function
194
195 \param[in,out] gf - generic function
196 ************************************************************/
generic_function_free(struct GenericFunction * gf)197 void generic_function_free(struct GenericFunction * gf){
198 if (gf != NULL){
199 if (gf->f != NULL) {
200 GF_SWITCH_NO_OUT(free)
201 gf->f = NULL;
202 }
203 free(gf); gf = NULL;
204 }
205 }
206
207 /********************************************************//**
208 Free memory for generic function array
209
210 \param[in,out] gf - generic function array
211 \param[in] n - number of generic functions in the array
212 ************************************************************/
generic_function_array_free(struct GenericFunction ** gf,size_t n)213 void generic_function_array_free(struct GenericFunction ** gf, size_t n){
214 if (gf != NULL){
215 size_t ii;
216 for (ii = 0; ii < n; ii++){
217 generic_function_free(gf[ii]);
218 gf[ii] = NULL;
219 }
220 free(gf); gf = NULL;
221 }
222 }
223
224 /* (Deep) Copy a generic function */
GF_IN_OUT(copy)225 GF_IN_OUT(copy)
226
227 /********************************************************//**
228 Copy a generic function to a preallocated generic function
229
230 \param[in] gf - generic function
231 \param[in,out] gfpa - preallocated function
232
233 ************************************************************/
234 void generic_function_copy_pa(const struct GenericFunction * gf,
235 struct GenericFunction * gfpa)
236 {
237 gfpa->fc = gf->fc;
238 void * temp = NULL;
239 GF_SWITCH_TEMPOUT(copy)
240 gfpa->f = temp;
241 }
242
243
244
245 /********************************************************//**
246 * Serialize a generic function
247 *
248 * \param[in,out] ser - location to serialize to
249 * \param[in] gf - generic function
250 * \param[in] totSizeIn - if not null then only total size
251 in bytes of generic function si returned
252 * if NULL then serialization occurs
253 *
254 * \return ptr - ser + num_bytes
255 ************************************************************/
256 unsigned char *
serialize_generic_function(unsigned char * ser,const struct GenericFunction * gf,size_t * totSizeIn)257 serialize_generic_function(unsigned char * ser,
258 const struct GenericFunction * gf,
259 size_t * totSizeIn)
260 {
261 // order =
262 // function_class -> sub_type -> function
263 unsigned char * ptr = ser;
264 size_t totSize = sizeof(int) + sizeof(size_t); // for function class and dim
265 if (totSizeIn != NULL){
266 size_t sizef = 0;
267 GF_SWITCH_THREE_FRONT(serialize,gf->fc,NULL,gf->f,&sizef)
268 totSize += sizef;
269 *totSizeIn = totSize;
270 return ptr;
271 }
272 else{
273 ptr = serialize_size_t(ptr, gf->dim);
274 ptr = serialize_int(ptr, gf->fc);
275 GF_SWITCH_THREEOUT_FRONT(serialize,gf->fc,ptr,ptr,gf->f, NULL)
276 }
277 return ptr;
278
279 }
280
281
282
283 /********************************************************//**
284 Save a generic function in text format
285
286 \param[in] gf - generic function to save
287 \param[in] stream - stream to save it to
288 \param[in] prec - precision with which to save it
289
290 ************************************************************/
generic_function_savetxt(const struct GenericFunction * gf,FILE * stream,size_t prec)291 void generic_function_savetxt(const struct GenericFunction * gf,
292 FILE * stream, size_t prec)
293 {
294 assert (gf != NULL);
295 fprintf(stream,"%zu ",gf->dim);
296 fprintf(stream,"%d ",(int)(gf->fc));
297 GF_SWITCH_NO_THREEOUT(savetxt,gf->fc,gf->f,stream,prec)
298 }
299
300 /********************************************************//**
301 Load a generic function in text format
302
303 \param[in] stream - stream to save it to
304
305 \return Generic function
306 ************************************************************/
307 struct GenericFunction *
generic_function_loadtxt(FILE * stream)308 generic_function_loadtxt(FILE * stream)
309 {
310 size_t dim;
311 int num = fscanf(stream,"%zu ",&dim);
312 assert (num == 1);
313 struct GenericFunction * gf = generic_function_alloc_base(dim);
314 int fcint;
315 num = fscanf(stream,"%d ",&fcint);
316 gf->fc = (enum function_class)fcint;
317 assert (num = 1);
318
319 GF_SWITCH_ONEOUT(loadtxt, gf->fc, gf->f, stream)
320 return gf;
321 }
322
print_generic_function(const struct GenericFunction * gf,size_t prec,void * args,FILE * fp)323 void print_generic_function(const struct GenericFunction * gf, size_t prec,void * args, FILE * fp)
324 {
325 GF_SWITCH_NO_FOUROUT_FRONT(print,gf->fc,gf->f,prec,args,fp)
326 }
327
328 /*******************************************************//**
329 Update a linear function
330
331 \param[in] gf - existing linear function
332 \param[in] a - slope of the function
333 \param[in] offset - offset of the function
334
335 \returns 0 if successfull, 1 otherwise
336 \note
337 Existing function must be linear
338 ***********************************************************/
339 int
generic_function_linear_update(struct GenericFunction * gf,double a,double offset)340 generic_function_linear_update(struct GenericFunction * gf,
341 double a, double offset)
342 {
343 int temp = 0;
344 GF_SWITCH_THREEOUT(linear_update,gf->fc,temp,gf->f,a,offset)
345 return temp;
346 }
347
348
349 /********************************************************//**
350 Create a generic function through regression of data
351
352 \return gf - generic function
353 ************************************************************/
354 struct GenericFunction *
generic_function_regress1d(struct Regress1DOpts * opts,struct c3Opt * optimizer,int * info)355 generic_function_regress1d(struct Regress1DOpts * opts, struct c3Opt * optimizer, int *info)
356 {
357
358 struct GenericFunction * func = NULL;
359 // perform linear regression to generate the starting point
360
361 // Initialize generic function to this linear function
362
363 double val;
364 if (opts->atype == PARAMETRIC){
365 double * start = calloc_double(opts->nparam);
366 memmove(start,opts->init_param,opts->nparam*sizeof(double));
367 if (opts->rtype == LS) {
368 c3opt_add_objective(optimizer,param_LSregress_cost,opts);
369 }
370 else if (opts->rtype == RLS2){
371 if (opts->reg_param_set == 0){
372 printf("Must set regularization parameter for RLS2 regression\n");
373 free(start); start = NULL;
374 return NULL;
375 }
376 c3opt_add_objective(optimizer,param_RLS2regress_cost,opts);
377 }
378 else if (opts->rtype == RLSD2){
379 if (opts->reg_param_set == 0){
380 printf("Must set regularization parameter for RLSD2 regression\n");
381 free(start); start = NULL;
382 return NULL;
383 }
384 c3opt_add_objective(optimizer,param_RLSD2regress_cost,opts);
385 }
386 else if (opts->rtype == RLSRKHS){
387 if (opts->reg_param_set == 0){
388 printf("Must set regularization parameter for RLSRKHS regression\n");
389 free(start); start = NULL;
390 return NULL;
391 }
392 else if (opts->decay_type == NONE){
393 printf("Must set decay type for parameter for RLSRKHS regression\n");
394 free(start); start = NULL;
395 return NULL;
396 }
397 c3opt_add_objective(optimizer,param_RLSRKHSregress_cost,opts);
398 }
399 else if (opts->rtype == RLS1){
400 printf("L1 regularization not yet implemented\n");
401 free(start); start = NULL;
402 return NULL;
403 /* c3opt_add_objective(optimizer,param_RLS1regress_cost,opts); */
404 }
405 else{
406 printf("Parameteric regression type %d is not recognized\n",opts->rtype);
407 free(start); start = NULL;
408 return NULL;
409 }
410
411
412 *info = c3opt_minimize(optimizer,start,&val);
413 /* if (*info > -1){ */
414 func = generic_function_create_with_params(opts->fc,opts->aopts,opts->nparam,start);
415 /* } */
416 free(start); start = NULL;
417 }
418 else if (opts->atype == NONPARAMETRIC){
419 printf("Non-parametric regression is not yet implemented\n");
420 return NULL;
421 }
422 else{
423 printf("Regression of type %d is not recognized\n",opts->atype);
424 return NULL;
425 }
426
427 return func;
428 }
429
430 /********************************************************//**
431 Create a generic function with particular parameters
432
433 \param[in] fc - function class
434 \param[in] aopts - approximation options
435 \param[in] dim - number of parameters
436 \param[in] param - parameter values to set
437
438 \return generic function
439 ************************************************************/
440 struct GenericFunction *
generic_function_create_with_params(enum function_class fc,void * aopts,size_t dim,const double * param)441 generic_function_create_with_params(enum function_class fc, void * aopts, size_t dim,
442 const double * param)
443 {
444
445 struct GenericFunction * gf = generic_function_alloc(1,fc);
446 GF_SWITCH_THREEOUT(create_with_params, fc, gf->f, aopts, dim, param);
447 return gf;
448 }
449
450 /********************************************************//**
451 Return a zero function
452
453 \param[in] fc - function class
454 \param[in] aopts - extra arguments depending on function_class, sub_type, etc.
455 \param[in] force_nparam - if == 1 then approximation will have the number of parameters
456 defined by *get_nparams, for each approximation type
457 if == 0 then it may be more compressed
458
459 \return gf - zero function
460 ************************************************************/
461 struct GenericFunction *
generic_function_zero(enum function_class fc,void * aopts,int force_nparam)462 generic_function_zero(enum function_class fc, void * aopts, int force_nparam)
463 {
464 struct GenericFunction * gf = generic_function_alloc(1,fc);
465 if (force_nparam == 0){ GF_SWITCH_TWOOUT(constant, fc, gf->f, 0, aopts) }
466 else { GF_SWITCH_TWOOUT(zero, fc, gf->f, aopts, 1) }
467 return gf;
468 }
469
470
471 /********************************************************//**
472 Return a constant function
473
474 \param[in] a - value of the function
475 \param[in] fc - function class
476 \param[in] aopts - extra arguments depending on function_class, sub_type, etc.
477
478 \return gf - constant function
479 ************************************************************/
480 struct GenericFunction *
generic_function_constant(double a,enum function_class fc,void * aopts)481 generic_function_constant(double a, enum function_class fc, void * aopts)
482 {
483 struct GenericFunction * gf = generic_function_alloc(1,fc);
484 GF_SWITCH_TWOOUT(constant, fc, gf->f, a, aopts)
485 return gf;
486 }
487
488 /*******************************************************//**
489 Return a linear function
490
491 \param[in] a - slope of the function
492 \param[in] offset - offset of the function
493 \param[in] fc - function class
494 \param[in] aopts - extra arguments depending on function_class,
495 sub_type, etc.
496
497 \return gf - linear function
498
499 \note
500 For kernel, this is only approximate
501 ***********************************************************/
502 struct GenericFunction *
generic_function_linear(double a,double offset,enum function_class fc,void * aopts)503 generic_function_linear(double a, double offset,
504 enum function_class fc, void * aopts)
505 {
506 struct GenericFunction * gf = generic_function_alloc(1,fc);
507 GF_SWITCH_THREEOUT(linear, fc, gf->f, a, offset, aopts)
508 return gf;
509 }
510
511 /*******************************************************//**
512 Return a quadratic function a * (x - offset)^2 = a (x^2 - 2offset x + offset^2)
513
514 \param[in] a - quadratic coefficients
515 \param[in] offset - shift of the function
516 \param[in] fc - function class
517 \param[in] aopts - extra arguments depending on function_class, sub_type, etc.
518
519 \return gf - quadratic
520 ************************************************************/
521 struct GenericFunction *
generic_function_quadratic(double a,double offset,enum function_class fc,void * aopts)522 generic_function_quadratic(double a, double offset,
523 enum function_class fc, void * aopts)
524 {
525 struct GenericFunction * gf = generic_function_alloc(1,fc);
526 GF_SWITCH_THREEOUT(quadratic, fc, gf->f, a, offset, aopts);
527 return gf;
528 }
529
530 /********************************************************//**
531 Create a pseudo-random polynomial generic function
532
533 * \param[in] ptype - polynomial type
534 * \param[in] maxorder - maximum order of the polynomial
535 * \param[in] lower - lower bound of input
536 * \param[in] upper - upper bound of input
537
538 \return gf - generic function
539 ************************************************************/
540 struct GenericFunction *
generic_function_poly_randu(enum poly_type ptype,size_t maxorder,double lower,double upper)541 generic_function_poly_randu(enum poly_type ptype,
542 size_t maxorder, double lower,
543 double upper)
544 {
545 enum function_class fc = POLYNOMIAL;
546 struct GenericFunction * gf = generic_function_alloc(1,fc);
547 gf->f = orth_poly_expansion_randu(ptype,maxorder,lower,upper);
548 gf->fargs = NULL;
549 return gf;
550 }
551
552 struct GenericFunction *
generic_function_onezero(enum function_class fc,double one,size_t nz,double * zeros,double lb,double ub)553 generic_function_onezero(enum function_class fc, double one, size_t nz,
554 double * zeros, double lb, double ub)
555 {
556 assert ((fc == LINELM) || (fc == CONSTELM));
557
558 struct GenericFunction * f =
559 generic_function_alloc(1, fc);
560
561
562 if (fc == LINELM){
563 struct LinElemExp * lexp = lin_elem_exp_alloc();
564 lexp->num_nodes = nz+3;
565 lexp->nodes = calloc_double(nz+3);
566 lexp->coeff = calloc_double(nz+3);
567
568 lexp->nodes[0] = lb;
569 size_t ind = 1;
570 int alloc = 0;
571 for (size_t ii = 0; ii < nz; ii++){
572 if (zeros[ii] < one){
573 lexp->nodes[ind] = zeros[ii];
574 ind++;
575 }
576 else if (alloc == 0){
577 // printf("lets go\n");
578 lexp->nodes[ind] = one;
579 lexp->coeff[ind] = 1.0;
580 ind++;
581 lexp->nodes[ind] = zeros[ii];
582 ind++;
583 alloc = 1;
584 }
585 else{
586 lexp->nodes[ind] = zeros[ii];
587 ind++;
588 }
589 }
590 if (alloc == 0){
591 lexp->nodes[ind] = one;
592 lexp->coeff[ind] = 1.0;
593 ind++;
594 }
595 assert (ind == nz+2);
596 lexp->nodes[nz+2] = ub;
597 f->f = lexp;
598 }
599 else if (fc == CONSTELM){
600 struct ConstElemExp * lexp = const_elem_exp_alloc();
601 lexp->num_nodes = nz+3;
602 lexp->nodes = calloc_double(nz+3);
603 lexp->coeff = calloc_double(nz+3);
604
605 lexp->nodes[0] = lb;
606 size_t ind = 1;
607 int alloc = 0;
608 for (size_t ii = 0; ii < nz; ii++){
609 if (zeros[ii] < one){
610 lexp->nodes[ind] = zeros[ii];
611 ind++;
612 }
613 else if (alloc == 0){
614 // printf("lets go\n");
615 lexp->nodes[ind] = one;
616 lexp->coeff[ind] = 1.0;
617 ind++;
618 lexp->nodes[ind] = zeros[ii];
619 ind++;
620 alloc = 1;
621 }
622 else{
623 lexp->nodes[ind] = zeros[ii];
624 ind++;
625 }
626 }
627 if (alloc == 0){
628 lexp->nodes[ind] = one;
629 lexp->coeff[ind] = 1.0;
630 ind++;
631 }
632 assert (ind == nz+2);
633 lexp->nodes[nz+2] = ub;
634 f->f = lexp;
635 }
636 return f;
637 }
638
639 /********************************************************//**
640 * Compute a linear combination of generic functions
641 *
642 * \param[in] n - number of functions
643 * \param[in] gfarray - array of functions
644 * \param[in] coeffs - scaling coefficients
645 *
646 * \return out = sum_i=1^n coeff[i] * gfarray[i]
647 ************************************************************/
648 struct GenericFunction *
generic_function_lin_comb(size_t n,struct GenericFunction ** gfarray,const double * coeffs)649 generic_function_lin_comb(size_t n,struct GenericFunction ** gfarray,
650 const double * coeffs)
651 {
652 // this function is not optimal
653 struct GenericFunction * out = NULL;
654 struct GenericFunction * temp1 = NULL;
655 struct GenericFunction * temp2 = NULL;
656 size_t ii;
657 if (n == 1){
658 out = generic_function_daxpby(coeffs[0],gfarray[0], 0.0, NULL);
659 }
660 else{
661 temp1 = generic_function_daxpby(coeffs[0],gfarray[0],coeffs[1],gfarray[1]);
662 for (ii = 2; ii < n; ii++){
663 if (ii % 2 == 0){
664 temp2 = generic_function_daxpby(coeffs[ii],gfarray[ii],1.0,temp1);
665 generic_function_free(temp1);
666 temp1 = NULL;
667 }
668 else{
669 temp1 = generic_function_daxpby(coeffs[ii],gfarray[ii],1.0,temp2);
670 generic_function_free(temp2);
671 temp2 = NULL;
672 }
673 }
674 }
675 if (temp1 != NULL){
676 return temp1;
677 }
678 else if (temp2 != NULL){
679 return temp2;
680 }
681 else{
682 assert (out != NULL);
683 return out;
684 }
685 }
686
687 /********************************************************//**
688 * Compute a linear combination of generic functions
689 *
690 * \param[in] n - number of functions
691 * \param[in] ldgf - stride of array to use
692 * \param[in] gfa - array of functions
693 * \param[in] ldc - stride of coefficents
694 * \param[in] c - scaling coefficients
695 *
696 * \return function representing
697 * \f$ \sum_{i=1}^n coeff[ldc[i]] * gfa[ldgf[i]] \f$
698 ************************************************************/
699 struct GenericFunction *
generic_function_lin_comb2(size_t n,size_t ldgf,struct GenericFunction ** gfa,size_t ldc,const double * c)700 generic_function_lin_comb2(size_t n, size_t ldgf,
701 struct GenericFunction ** gfa,
702 size_t ldc, const double * c)
703 {
704 // this function is not optimal
705 struct GenericFunction * out = NULL;
706 struct GenericFunction * temp1 = NULL;
707 struct GenericFunction * temp2 = NULL;
708 size_t ii;
709 if (n == 1){
710 out = generic_function_daxpby(c[0],gfa[0], 0.0, NULL);
711 }
712 else{
713
714 int allpoly = 1;
715 for (ii = 0; ii < n; ii++){
716 if (gfa[ii*ldgf]->fc != POLYNOMIAL){
717 allpoly = 0;
718 break;
719 }
720 }
721
722 if (allpoly == 1){
723 struct OrthPolyExpansion ** xx = NULL;
724 if (NULL == (xx = malloc(n * sizeof(struct OrthPolyExpansion *)))){
725 fprintf(stderr, "failed to allocate memmory in generic_function_lin_comb2\n");
726 exit(1);
727 }
728 for (ii = 0; ii < n; ii++){
729 xx[ii] = gfa[ii*ldgf]->f;
730 }
731
732 struct GenericFunction * gf = generic_function_alloc(1,gfa[0]->fc);
733 gf->f = orth_poly_expansion_lin_comb(n,1,xx,ldc,c);
734 gf->fargs = NULL;
735 free(xx); xx = NULL;
736
737 assert (gf->f != NULL);
738 return gf;
739 }
740
741 temp1 = generic_function_daxpby(c[0],gfa[0],c[ldc],gfa[ldgf]);
742 for (ii = 2; ii < n; ii++){
743 if (ii % 2 == 0){
744 temp2 = generic_function_daxpby(c[ii*ldc],gfa[ii*ldgf],
745 1.0,temp1);
746 generic_function_free(temp1);
747 temp1 = NULL;
748 }
749 else{
750 temp1 = generic_function_daxpby(c[ii*ldc],gfa[ii*ldgf],
751 1.0,temp2);
752 generic_function_free(temp2);
753 temp2 = NULL;
754 }
755 }
756 }
757 if (temp1 != NULL){
758 return temp1;
759 }
760 else if (temp2 != NULL){
761 return temp2;
762 }
763 else{
764 assert (out != NULL);
765 return out;
766 }
767 }
768
769
770 /* Take the derivative of a generic function */
771 GF_IN_OUT(deriv)
GF_IN_OUT(dderiv)772 GF_IN_OUT(dderiv)
773 GF_IN_OUT(dderiv_periodic)
774
775
776 /********************************************************//**
777 * Create a nodal basis at particular points
778 *
779 * \param[in] f - function to interpolate
780 * \param[in] N - number of nodes
781 * \param[in] x - locations of nodes
782 *
783 * \return nodal basis (LINELM) function
784 ************************************************************/
785 struct GenericFunction *
786 generic_function_create_nodal(struct GenericFunction * f,size_t N, double * x)
787 {
788 struct GenericFunction * out = NULL;
789 /* out = generic_function_alloc(f->dim,LINELM); */
790 out = generic_function_alloc(f->dim,f->fc);
791 out->fargs = NULL;
792 double * fvals = calloc_double(N);
793 for (size_t ii = 0; ii < N; ii++){
794 fvals[ii] = generic_function_1d_eval(f,x[ii]);
795 }
796 if (f->fc == LINELM){
797 out->f = lin_elem_exp_init(N,x,fvals);
798 }
799 else if (f->fc == CONSTELM){
800 out->f = const_elem_exp_init(N,x,fvals);
801 }
802 else{
803 fprintf(stderr,"Cannot create nodal function of this type\n");
804 exit(1);
805 }
806 free(fvals); fvals = NULL;
807
808 return out;
809 }
810
811
812 /********************************************************//**
813 * Compute axpby for a an array of generic functions
814 *
815 * \param[in] n - number of functions
816 * \param[in] a - scaling for x
817 * \param[in] ldx - stride of functions to use in a
818 * \param[in] x - functions
819 * \param[in] b - scaling for y
820 * \param[in] ldy - stride of functions to use in a
821 * \param[in] y - functions
822 *
823 * \return fout - array of generic functions
824 *************************************************************/
825 struct GenericFunction **
generic_function_array_daxpby(size_t n,double a,size_t ldx,struct GenericFunction ** x,double b,size_t ldy,struct GenericFunction ** y)826 generic_function_array_daxpby(size_t n, double a, size_t ldx,
827 struct GenericFunction ** x, double b, size_t ldy,
828 struct GenericFunction ** y)
829 {
830 struct GenericFunction ** fout = NULL;
831 if (NULL == ( fout = malloc(n*sizeof(struct GenericFunction *)))){
832 fprintf(stderr, "failed to allocate in generic_function_array_daxpby.\n");
833 exit(1);
834 }
835 //printf("in daxpby here!\n");
836 size_t ii;
837 if ( y == NULL){
838 for (ii = 0; ii < n ;ii++){
839 fout[ii] = generic_function_daxpby(a,x[ii*ldx],0.0, NULL);
840 }
841 }
842 else if (x == NULL){
843 for (ii = 0; ii < n ;ii++){
844 fout[ii] = generic_function_daxpby(b,y[ii*ldy],0.0, NULL);
845 }
846 }
847 else{
848 for (ii = 0; ii < n ;ii++){
849 //printf("array daxpby ii=(%zu/%zu)!\n",ii,n);
850 fout[ii] = generic_function_daxpby(a,x[ii*ldx],b, y[ii*ldy]);
851 //printf("outhere ii=(%zu/%zu)!\n",ii,n);
852 }
853 }
854 //printf("return \n");
855 return fout;
856 }
857
858 GF_IN_GENOUT(get_num_params, size_t, 0) // Get the number of parameters describing the generic function
859 GF_IN_GENOUT(get_lb, double, -0.123456789) // Get the lower bound of the input space
860 GF_IN_GENOUT(get_ub, double, 0.123456789) // Get the lower bound of the input space
861
862 /********************************************************//**
863 * Get the function class
864 ************************************************************/
generic_function_get_fc(const struct GenericFunction * f)865 enum function_class generic_function_get_fc(const struct GenericFunction * f)
866 {
867 assert (f != NULL);
868 return f->fc;
869 }
870
871 /***********************************************************//**
872 Determine whether kristoffel weighting is active
873
874 \param[in] gf - generic function
875
876 \return 1 if active, 0 otherwise
877 ***************************************************************/
generic_function_is_kristoffel_active(const struct GenericFunction * gf)878 int generic_function_is_kristoffel_active(const struct GenericFunction * gf)
879 {
880 if (gf->fc != POLYNOMIAL){
881 return 0;
882 }
883 else{
884 struct OrthPolyExpansion * ope = gf->f;
885 return ope->kristoffel_eval;
886 }
887 }
888
generic_function_activate_kristoffel(struct GenericFunction * gf)889 void generic_function_activate_kristoffel(struct GenericFunction * gf)
890 {
891 if (gf->fc != POLYNOMIAL){
892 fprintf(stderr,"Cannot activate kristoffel for non polynomial basis\n");
893 exit(1);
894 }
895 else{
896 struct OrthPolyExpansion * ope = gf->f;
897 ope->kristoffel_eval = 1;
898 }
899 }
900
generic_function_deactivate_kristoffel(struct GenericFunction * gf)901 void generic_function_deactivate_kristoffel(struct GenericFunction * gf)
902 {
903 if (gf->fc == POLYNOMIAL){
904 struct OrthPolyExpansion * ope = gf->f;
905 ope->kristoffel_eval = 0;
906 }
907 }
908
909
910 /***********************************************************//**
911 Get the kristoffel normalization factor
912
913 \param[in] gf - generic function
914 \param[in] x - location at which to obtain normalization
915
916 \return normalization factor
917 ***************************************************************/
generic_function_get_kristoffel_weight(const struct GenericFunction * gf,double x)918 double generic_function_get_kristoffel_weight(const struct GenericFunction * gf,
919 double x)
920 {
921 if (gf->fc != POLYNOMIAL){
922 fprintf(stderr, "Cannot get the kristoffel weight of a function that is not a polynomial\n");
923 exit(1);
924 }
925 else{
926 double weight = orth_poly_expansion_get_kristoffel_weight(gf->f,x);
927 return weight;
928 }
929 }
930
931 /********************************************************//**
932 Get the parameters of generic function
933
934 \param[in] gf - generic function
935 \param[in,out] params - location to write parameters
936
937 \returns number of parameters
938 ************************************************************/
generic_function_get_params(const struct GenericFunction * gf,double * params)939 size_t generic_function_get_params(const struct GenericFunction * gf, double * params)
940 {
941
942 assert (gf != NULL);
943 size_t nparam = 0;
944 GF_SWITCH_TWOOUT(get_params, gf->fc, nparam, gf->f, params)
945 return nparam;
946 }
947
948 /********************************************************//**
949 Get the parameters of generic function
950
951 \param[in] gf - generic function
952 \param[in,out] nparam - location to write parameters
953
954 \returns reference to parameters
955 ************************************************************/
generic_function_get_params_ref(const struct GenericFunction * gf,size_t * nparam)956 double * generic_function_get_params_ref(const struct GenericFunction * gf, size_t * nparam)
957 {
958
959 assert (gf != NULL);
960 double * params = NULL;
961 GF_SWITCH_TWOOUT(get_params_ref, gf->fc, params, gf->f, nparam)
962
963 return params;
964 }
965
966
967 /********************************************************//**
968 Update a generic function with particular parameters
969
970 \param[in] f - function to update
971 \param[in] dim - number of parameters
972 \param[in] param - parameter values to set
973
974 \returns 0 if successfull, 1 otherwise
975 ************************************************************/
976 int
generic_function_update_params(struct GenericFunction * f,size_t dim,const double * param)977 generic_function_update_params(struct GenericFunction * f, size_t dim,
978 const double * param)
979 {
980
981 for (size_t ii = 0; ii < dim; ii++){
982 if (isnan(param[ii])){
983 fprintf(stderr,"Updating generic functions with params that are NaN\n");
984 exit(1);
985 }
986 else if (isinf(param[ii])){
987 fprintf(stderr,"Updating generic functions with params that are inf\n");
988 exit(1);
989 }
990 }
991
992 int res = 0;
993 GF_SWITCH_THREEOUT(update_params, f->fc, res, f->f, dim, param)
994 return res;
995 }
996
997
998 /********************************************************//**
999 * Round an generic function to some tolerance
1000 *
1001 * \param[in,out] gf - generic function
1002 * \param[in] thresh - threshold (relative) to round to
1003 *
1004 * \note
1005 * (UNTESTED, use with care!!!!
1006 *************************************************************/
generic_function_roundt(struct GenericFunction ** gf,double thresh)1007 void generic_function_roundt(struct GenericFunction ** gf, double thresh)
1008 {
1009 struct OrthPolyExpansion * ope = NULL;
1010 assert ( (*gf)->fc == POLYNOMIAL);
1011 ope = (*gf)->f;
1012 orth_poly_expansion_roundt(&ope,thresh);
1013 }
1014
1015 static struct GenericFunction *
generic_function_onezero2(enum function_class fc,size_t nzeros,double * zero_locations,void * opts)1016 generic_function_onezero2(
1017 enum function_class fc,
1018 size_t nzeros,
1019 double * zero_locations,
1020 void * opts
1021 )
1022 {
1023 struct GenericFunction * gf = generic_function_alloc(1,fc);
1024 if (fc == LINELM){
1025 gf->f = lin_elem_exp_onezero(nzeros, zero_locations, opts);
1026 }
1027 else if (fc == CONSTELM){
1028 gf->f = const_elem_exp_onezero(nzeros, zero_locations, opts);
1029 }
1030 else{
1031 fprintf(stderr,"Cannot create a onezero generic function for non-nodal basis\n");
1032 exit(1);
1033 }
1034 return gf;
1035 }
1036
generic_function_array_onezero(struct GenericFunction ** L,size_t n,enum function_class fc,size_t upto,size_t * piv,double * px,void * opts)1037 void generic_function_array_onezero(
1038 struct GenericFunction ** L,
1039 size_t n,
1040 enum function_class fc,
1041 size_t upto,
1042 size_t * piv,
1043 double * px,
1044 void * opts)
1045 // void * opt_args)
1046 {
1047 //create an arbitrary array that has zeros at piv[:upto-1],px[:upto-1]
1048 // and one at piv[upto],piv[upto] less than one every else
1049
1050 // note that need to set piv[upto] and px[upto] in this function
1051 // number of pivots per array
1052
1053 size_t * npiv = calloc_size_t(n);
1054 double ** x = malloc_dd(n); // pivots per function
1055 for (size_t ii = 0; ii < upto; ii++){
1056 npiv[piv[ii]]++;
1057 }
1058
1059 for (size_t ii = 0; ii < n; ii++){
1060 x[ii] = calloc_double(npiv[ii]);
1061 size_t on = 0;
1062 for (size_t jj = 0; jj < upto; jj++){
1063 if (piv[jj] == ii){
1064 x[ii][on] = px[jj];
1065 on++;
1066 }
1067 }
1068 }
1069
1070 for (size_t ii = 0; ii < n; ii++){
1071 L[ii] = generic_function_onezero2(fc,npiv[ii],x[ii],opts);
1072 }
1073
1074 double xval;
1075 size_t amind;
1076 generic_function_array_absmax(n, 1, L,&amind, &xval,NULL);//optargs);
1077 px[upto] = xval;
1078 piv[upto] = amind;
1079 double val = generic_function_1d_eval(L[piv[upto]],px[upto]);
1080 generic_function_array_scale(1.0/val,L,n);
1081
1082 free_dd(n,x);
1083 free(npiv); npiv = NULL;
1084 }
1085
1086 /********************************************************//**
1087 * Flip the sign of a generic function f(x) to -f(x)
1088 *
1089 * \param[in,out] gf - number of functions
1090 ************************************************************/
generic_function_flip_sign(struct GenericFunction * gf)1091 void generic_function_flip_sign(struct GenericFunction * gf)
1092 {
1093 GF_SWITCH_NO_OUT(flip_sign)
1094 }
1095
1096 /********************************************************//**
1097 * Flip the sign of each generic function in an array
1098 *
1099 * \param[in] n - number of functions
1100 * \param[in] lda - stride of array
1101 * \param[in,out] a - array of functions
1102 ************************************************************/
1103 void
generic_function_array_flip_sign(size_t n,size_t lda,struct GenericFunction ** a)1104 generic_function_array_flip_sign(size_t n, size_t lda,
1105 struct GenericFunction ** a){
1106 size_t ii;
1107 for (ii = 0; ii < n; ii++){
1108 generic_function_flip_sign(a[ii*lda]);
1109 }
1110 }
1111
1112
1113 /********************************************************//**
1114 * Multiply and add 3 functions \f$ z \leftarrow ax + by + cz \f$
1115 *
1116 * \param a [in] - first scaling factor
1117 * \param x [in] - first function
1118 * \param b [in] - second scaling factor
1119 * \param y [in] - second function
1120 * \param c [in] - third scaling factor
1121 * \param z [in] - third function
1122 *
1123 *************************************************************/
1124 void
generic_function_sum3_up(double a,struct GenericFunction * x,double b,struct GenericFunction * y,double c,struct GenericFunction * z)1125 generic_function_sum3_up(double a, struct GenericFunction * x,
1126 double b, struct GenericFunction * y,
1127 double c, struct GenericFunction * z)
1128 {
1129 if ( x->fc != POLYNOMIAL){
1130 fprintf(stderr, "Have not yet implemented generic_function_sum3_up \n");
1131 fprintf(stderr, "for functions other than polynomials\n");
1132 exit(1);
1133 }
1134 assert (x->fc == y->fc);
1135 assert (y->fc == z->fc);
1136
1137 orth_poly_expansion_sum3_up(a,x->f,b,y->f,c,z->f);
1138 }
1139
1140 /********************************************************//**
1141 * Add two generic functions \f$ y \leftarrow ax + y \f$
1142 *
1143 * \param[in] a - scaling of first function
1144 * \param[in] x - first function
1145 * \param[in,out] y - second function
1146 *
1147 * \return 0 if successfull, 1 if error
1148 *
1149 * \note
1150 * Handling the function class of the output is not very smart
1151 ************************************************************/
generic_function_axpy(double a,const struct GenericFunction * x,struct GenericFunction * y)1152 int generic_function_axpy(double a, const struct GenericFunction * x,
1153 struct GenericFunction * y)
1154 {
1155 //printf("in here! a =%G b = %G\n",a,b);
1156
1157 assert (y != NULL);
1158 assert (x != NULL);
1159 assert (x->fc == y->fc);
1160
1161 int out = 1;
1162 GF_SWITCH_THREEOUT(axpy, x->fc, out, a, x->f, y->f);
1163 return out;
1164 }
1165
1166 /********************************************************//**
1167 * Add generic functions \f$ y[i] \leftarrow a x[i] + y[i] \f$
1168 *
1169 * \param[in] n - number of functions
1170 * \param[in] a - scaling of the first functions
1171 * \param[in] x - first function array
1172 * \param[in,out] y - second function array
1173 *
1174 * \return 0 if successfull, 1 if error
1175 *
1176 * \note
1177 * Handling the function class of the output is not very smart
1178 ************************************************************/
generic_function_array_axpy(size_t n,double a,struct GenericFunction ** x,struct GenericFunction ** y)1179 int generic_function_array_axpy(size_t n, double a,
1180 struct GenericFunction ** x,
1181 struct GenericFunction ** y)
1182 {
1183 //printf("in here! a =%G b = %G\n",a,b);
1184
1185 assert (y != NULL);
1186 assert (x != NULL);
1187
1188 int success = 1;
1189 size_t ii;
1190 for (ii = 0; ii < n; ii++){
1191 success = generic_function_axpy(a,x[ii],y[ii]);
1192 if (success == 1){
1193 break;
1194 }
1195 }
1196
1197 return success;
1198 }
1199
1200 /********************************************************//**
1201 Scale a generic function
1202
1203 \param[in] a - value with which to scale the functions
1204 \param[in,out] gf - function to scale
1205 ************************************************************/
generic_function_scale(double a,struct GenericFunction * gf)1206 void generic_function_scale(double a, struct GenericFunction * gf)
1207 {
1208 GF_SWITCH_NO_ONEOUT(scale, gf->fc, a, gf->f)
1209 }
1210
1211 /********************************************************//**
1212 Scale a generic function array
1213
1214 \param[in] a - value with which to scale the functions
1215 \param[in,out] gf - functions to scale
1216 \param[in] N - number of functions
1217 ************************************************************/
generic_function_array_scale(double a,struct GenericFunction ** gf,size_t N)1218 void generic_function_array_scale(double a, struct GenericFunction ** gf,
1219 size_t N)
1220 {
1221 size_t ii;
1222 for (ii = 0; ii < N; ii++){
1223 generic_function_scale(a,gf[ii]);
1224 }
1225 }
1226
1227
1228
1229 /********************************************************//**
1230 Helper function for computing the first part of
1231 \f[
1232 a kron(\cdot,c)
1233 \f]
1234 if left = 1,
1235 otherwise
1236 \f[
1237 kron(\cdot,c)a
1238 \f]
1239
1240 \param left [in]
1241 \param r [in]
1242 \param m [in]
1243 \param n [in]
1244 \param l [in]
1245 \param a [in] - if left = 1 (r, m * n) otherwise (l * m,r)
1246 \param c [in] - (n,l)
1247 \param d [inout] - (rm,l)
1248
1249 ************************************************************/
generic_function_kronh(int left,size_t r,size_t m,size_t n,size_t l,const double * a,struct GenericFunction ** c,struct GenericFunction ** d)1250 void generic_function_kronh(int left,
1251 size_t r, size_t m, size_t n, size_t l,
1252 const double * a,
1253 struct GenericFunction ** c,
1254 struct GenericFunction ** d)
1255 {
1256 size_t ii,jj,kk;
1257 if (left == 1){
1258 for (kk = 0; kk < l; kk++){
1259 for (jj = 0; jj < m; jj++){
1260 for (ii = 0; ii < r; ii++){
1261 //printf("(%zu/%zu,%zu/%zu,%zu/%zu)\n",jj,m-1,kk,l-1,ii,r-1);
1262 d[kk*r*m + jj*r + ii] =
1263 generic_function_lin_comb2(
1264 n, 1, c + n*kk, r, a + ii + jj*n*r);
1265 }
1266 }
1267 }
1268 }
1269 else{
1270 for (ii = 0; ii < r; ii++){
1271 for (jj = 0; jj < n; jj++){
1272 for (kk = 0; kk < m; kk++){
1273 d[jj + kk*n + ii*n*m] =
1274 generic_function_lin_comb2(
1275 l, n, c + jj, 1, a + kk*l + ii*l*m);
1276 }
1277 }
1278 }
1279 }
1280 }
1281
generic_function_kronh2(int left,size_t r,size_t m,size_t n,size_t l,struct GenericFunction ** b,struct GenericFunction ** t,struct GenericFunction ** out)1282 void generic_function_kronh2(int left, size_t r, size_t m, size_t n, size_t l,
1283 struct GenericFunction ** b, struct GenericFunction ** t,
1284 struct GenericFunction ** out)
1285 {
1286 if (left == 1){
1287 size_t ii,jj, kk;
1288 for (jj = 0; jj < l; jj++){
1289 for (kk = 0; kk < m; kk++){
1290 for (ii = 0; ii < r; ii++){
1291 out[ii + kk*r + jj*r*m] =
1292 generic_function_sum_prod(
1293 n, 1, b + jj*n, r, t + ii + kk*r*n);
1294 }
1295 }
1296 }
1297 }
1298 else {
1299 size_t ii,jj, kk;
1300 for (ii = 0; ii < r; ii++){
1301 for (jj = 0; jj < n; jj++){
1302 for (kk = 0; kk < m; kk++){
1303 out[kk + jj*m + ii*n*m] =
1304 generic_function_sum_prod(
1305 l, n, b + jj, m, t + kk + ii * l * m);
1306 }
1307 }
1308 }
1309 }
1310 }
1311
1312
1313
1314 /********************************************************//**
1315 * Compute axpby for a an array of generic functions and overwrite into z
1316 *
1317 * \param[in] n - number of functions
1318 * \param[in] a - scaling for x
1319 * \param[in] ldx - stride of functions to use in a
1320 * \param[in] x - functions
1321 * \param[in] b - scaling for y
1322 * \param[in] ldy - stride of functions to use in a
1323 * \param[in] y - functions
1324 * \param[in] ldz - stride for z
1325 * \param[in,out] z - locations for resulting functions
1326 *************************************************************/
1327 void
generic_function_array_daxpby2(size_t n,double a,size_t ldx,struct GenericFunction ** x,double b,size_t ldy,struct GenericFunction ** y,size_t ldz,struct GenericFunction ** z)1328 generic_function_array_daxpby2(size_t n, double a, size_t ldx,
1329 struct GenericFunction ** x, double b, size_t ldy,
1330 struct GenericFunction ** y, size_t ldz,
1331 struct GenericFunction ** z)
1332 {
1333 size_t ii;
1334 if ( y == NULL){
1335 for (ii = 0; ii < n ;ii++){
1336 z[ii*ldz] = generic_function_daxpby(a,x[ii*ldx],0.0, NULL);
1337 }
1338 }
1339 else if (x == NULL){
1340 for (ii = 0; ii < n ;ii++){
1341 z[ii*ldz] = generic_function_daxpby(b,y[ii*ldy],0.0, NULL);
1342 }
1343 }
1344 else{
1345 for (ii = 0; ii < n ;ii++){
1346 z[ii*ldz] = generic_function_daxpby(a,x[ii*ldx],b, y[ii*ldy]);
1347 }
1348 }
1349 }
1350
1351
1352
1353 /********************************************************//**
1354 * Evaluate a generic function
1355 *
1356 * \param[in] f - function
1357 * \param[in] x - location at which to evaluate
1358 *
1359 * \return evaluation
1360 ************************************************************/
generic_function_1d_eval(const struct GenericFunction * f,double x)1361 double generic_function_1d_eval(const struct GenericFunction * f, double x){
1362 assert (f != NULL);
1363 double out = 0.1234567890;
1364
1365 GF_SWITCH_TWOOUT(eval, f->fc, out, f->f, x)
1366
1367 if (isnan(out)){
1368 fprintf(stderr,"Warning, evaluation of generic_function is nan\n");
1369 exit(1);
1370 }
1371 else if (isinf(out)){
1372 fprintf(stderr,"Warning, evaluation of generic_function is inf\n");
1373 exit(1);
1374 }
1375 return out;
1376 }
1377
1378 /********************************************************//**
1379 Evaluate the derivative of a generic function
1380
1381 \param[in] gf - generic function
1382 \param[in] x - location at which to evaluate
1383
1384 \return value of the derivative
1385 ************************************************************/
generic_function_deriv_eval(const struct GenericFunction * gf,double x)1386 double generic_function_deriv_eval(const struct GenericFunction * gf, double x)
1387 {
1388 double out = 0.1234567890;
1389 GF_SWITCH_TWOOUT(deriv_eval, gf->fc, out, gf->f, x);
1390 return out;
1391 }
1392
1393
1394 /********************************************************//**
1395 * Evaluate a generic function at multiple locations
1396 *
1397 * \param[in] f - function
1398 * \param[in] N - number of evaluations
1399 * \param[in] x - location at which to evaluate
1400 * \param[in] incx - increment of x
1401 * \param[in,out] y - allocated space for evaluations
1402 * \param[in] incy - increment of y
1403 ************************************************************/
generic_function_1d_evalN(const struct GenericFunction * f,size_t N,const double * x,size_t incx,double * y,size_t incy)1404 void generic_function_1d_evalN(const struct GenericFunction * f, size_t N,
1405 const double * x, size_t incx, double * y, size_t incy)
1406 {
1407 assert (f != NULL);
1408 assert (f->f != NULL);
1409 GF_SWITCH_SIX(evalN,f->fc,f->f,N,x,incx,y,incy)
1410 }
1411
1412 /********************************************************//**
1413 * Evaluate a generic function consisting of nodal
1414 * basis functions at some node
1415 *
1416 * \param[in] f - function
1417 * \param[in] ind - location at which to evaluate
1418 *
1419 * \return evaluation
1420 ************************************************************/
generic_function_1d_eval_ind(const struct GenericFunction * f,size_t ind)1421 double generic_function_1d_eval_ind(const struct GenericFunction * f, size_t ind)
1422 {
1423 assert (f != NULL);
1424 double out = 0.1234567890;
1425 if (f->fc == LINELM){
1426 out = lin_elem_exp_get_nodal_val(f->f,ind);
1427 }
1428 else if (f->fc == CONSTELM){
1429 out = const_elem_exp_get_nodal_val(f->f,ind);
1430 }
1431 else{
1432 assert (1 == 0);
1433 }
1434
1435 return out;
1436 }
1437
1438
1439 /********************************************************//**
1440 * Evaluate an array of generic functions
1441 *
1442 * \param[in] n - number of functions
1443 * \param[in] f - array of functions
1444 * \param[in] x - location at which to evaluate
1445 *
1446 * \return array of values
1447 ************************************************************/
1448 double *
generic_function_1darray_eval(size_t n,struct GenericFunction ** f,double x)1449 generic_function_1darray_eval(size_t n, struct GenericFunction ** f, double x)
1450 {
1451 double * out = calloc_double(n);
1452 size_t ii;
1453 for (ii = 0; ii < n; ii++){
1454 out[ii] = generic_function_1d_eval(f[ii],x);
1455 }
1456 return out;
1457 }
1458
1459 /********************************************************//**
1460 * Evaluate a generic function consisting of nodal
1461 * basis functions at some node
1462 *
1463 * \param[in] f - function
1464 * \param[in] x - location at which to Evaluate
1465 * \param[in] size - byte size of location (sizeof(double) or (sizeof(size_t)))
1466 *
1467 * \return evaluation
1468 ************************************************************/
generic_function_1d_eval_gen(const struct GenericFunction * f,void * x,size_t size)1469 static double generic_function_1d_eval_gen(const struct GenericFunction * f,
1470 void * x, size_t size)
1471 {
1472 assert (f != NULL);
1473
1474 size_t dsize = sizeof(double);
1475 size_t stsize = sizeof(size_t);
1476 double out;
1477 if (size == dsize){
1478 out = generic_function_1d_eval(f,*(double *)x);
1479 }
1480 else if (size == stsize){
1481 out = generic_function_1d_eval_ind(f,*(size_t *)x);
1482 }
1483 else{
1484 fprintf(stderr, "Cannot evaluate generic function at \n");
1485 fprintf(stderr, "input of byte size %zu\n ", size);
1486 exit(1);
1487 }
1488
1489 return out;
1490 }
1491
1492 /********************************************************//**
1493 Evaluate a generic function array at a given pivot
1494 ************************************************************/
generic_function_1darray_eval_piv(struct GenericFunction ** f,struct Pivot * piv)1495 double generic_function_1darray_eval_piv(struct GenericFunction ** f,
1496 struct Pivot * piv)
1497 {
1498 size_t size = pivot_get_size(piv);
1499 size_t ind = pivot_get_ind(piv);
1500 void * loc = pivot_get_loc(piv);
1501 double out = generic_function_1d_eval_gen(f[ind],loc,size);
1502 return out;
1503 }
1504
1505 /********************************************************//**
1506 * Evaluate an array of generic functions
1507 *
1508 * \param[in] n - number of functions
1509 * \param[in] f - array of functions
1510 * \param[in] x - location at which to evaluate
1511 * \param[in,out] out - array of values
1512 ************************************************************/
1513 void
generic_function_1darray_eval2(size_t n,struct GenericFunction ** f,double x,double * out)1514 generic_function_1darray_eval2(size_t n,
1515 struct GenericFunction ** f,
1516 double x, double * out)
1517 {
1518 int allpoly = 1;
1519 struct OrthPolyExpansion * parr[1000];
1520 for (size_t ii = 0; ii < n; ii++){
1521 if (f[ii]->fc != POLYNOMIAL){
1522 allpoly = 0;
1523 break;
1524 }
1525 parr[ii] = f[ii]->f;
1526 }
1527 if ((allpoly == 1) && (n <= 1000)){
1528 int res = orth_poly_expansion_arr_eval(n,parr,x,out);
1529 if (res == 1){ //something when wrong
1530 size_t ii;
1531 for (ii = 0; ii < n; ii++){
1532 out[ii] = generic_function_1d_eval(f[ii],x);
1533 }
1534 }
1535 }
1536 else{
1537 size_t ii;
1538 for (ii = 0; ii < n; ii++){
1539 out[ii] = generic_function_1d_eval(f[ii],x);
1540 }
1541 }
1542 }
1543
1544 /********************************************************//**
1545 * Evaluate an array of functions at an array of points
1546 *
1547 * \param[in] n - number of functions
1548 * \param[in] f - function
1549 * \param[in] N - number of evaluations
1550 * \param[in] x - location at which to evaluate
1551 * \param[in] incx - increment of x
1552 * \param[in,out] y - allocated space for evaluations
1553 * \param[in] incy - increment of y*
1554 *
1555 * \note Currently just calls the single evaluation code
1556 * Note sure if this is optimal, cache-wise
1557 *************************************************************/
1558 void
generic_function_1darray_eval2N(size_t n,struct GenericFunction ** f,size_t N,const double * x,size_t incx,double * y,size_t incy)1559 generic_function_1darray_eval2N(size_t n,
1560 struct GenericFunction ** f,
1561 size_t N, const double * x, size_t incx,
1562 double * y, size_t incy)
1563 {
1564
1565 int allpoly = 1;
1566 struct OrthPolyExpansion * parr[1000];
1567 for (size_t ii = 0; ii < n; ii++){
1568 if (f[ii]->fc != POLYNOMIAL){
1569 allpoly = 0;
1570 break;
1571 }
1572 parr[ii] = f[ii]->f;
1573 }
1574 if ((allpoly == 1) && (n <= 1000)){
1575 /* printf("generic function, kristoffel_active = %d\n",generic_function_is_kristoffel_active(f[0])); */
1576 int res = orth_poly_expansion_arr_evalN(n,parr,N,x,incx,y,incy);
1577 /* for (size_t ii = 0; ii < N*n; ii++){ */
1578 /* printf("y[%zu] = %G\n",ii,y[ii]); */
1579 /* } */
1580 if (res == 1){ //something when wrong
1581 size_t ii;
1582 for (ii = 0; ii < n; ii++){
1583 generic_function_1darray_eval2(n,f,x[ii*incx],y+ii*incy);
1584 }
1585 }
1586 }
1587 else{
1588 for (size_t ii = 0; ii < N; ii++){
1589 generic_function_1darray_eval2(n,f,x[ii*incx],y+ii*incy);
1590 }
1591 }
1592 }
1593
1594 /********************************************************//**
1595 * Evaluate an array of generic functions which should be
1596 * of nodal basis class at particular nodal locations
1597 *
1598 * \param[in] n - number of functions
1599 * \param[in] f - array of functions
1600 * \param[in] ind - location at which to evaluate
1601 * \param[in,out] out - array of values
1602 ************************************************************/
1603 void
generic_function_1darray_eval2_ind(size_t n,struct GenericFunction ** f,size_t ind,double * out)1604 generic_function_1darray_eval2_ind(size_t n,
1605 struct GenericFunction ** f,
1606 size_t ind, double * out)
1607 {
1608
1609 size_t ii;
1610 for (ii = 0; ii < n; ii++){
1611 out[ii] = generic_function_1d_eval_ind(f[ii],ind);
1612 }
1613 }
1614
1615 /********************************************************//**
1616 Take a gradient with respect to function parameters
1617
1618 \param[in] gf - generic function
1619 \param[in] nx - number of x values
1620 \param[in] x - x values
1621 \param[in,out] grad - gradient (N,nx)
1622
1623 \return 0 - success, 1 -failure
1624 ************************************************************/
generic_function_param_grad_eval(const struct GenericFunction * gf,size_t nx,const double * x,double * grad)1625 int generic_function_param_grad_eval(const struct GenericFunction * gf,
1626 size_t nx, const double * x,
1627 double * grad)
1628 {
1629
1630 enum function_class fc = generic_function_get_fc(gf);
1631 int res = 1;
1632 GF_SWITCH_FOUROUT(param_grad_eval, fc, res, gf->f, nx, x, grad)
1633 assert (res == 0);
1634 return res;
1635 }
1636
1637
1638 /********************************************************//**
1639 Take a gradient with respect to function parameters
1640
1641 \param[in] gf - generic function
1642 \param[in] x - x values
1643 \param[in,out] grad - gradient (N)
1644
1645 \return evaluation
1646 ************************************************************/
generic_function_param_grad_eval2(const struct GenericFunction * gf,double x,double * grad)1647 double generic_function_param_grad_eval2(const struct GenericFunction * gf,
1648 double x,double * grad)
1649
1650 {
1651
1652 enum function_class fc = generic_function_get_fc(gf);
1653 double ret = 0.1234;
1654 GF_SWITCH_THREEOUT(param_grad_eval2, fc, ret, gf->f, x, grad)
1655 return ret;
1656 }
1657
1658
1659 GF_IN_GENOUT(integrate, double, 0.0) // Compute an integral
1660 GF_IN_GENOUT(integrate_weighted, double, 0.0) // Take the derivative of a generic function
1661
1662 /********************************************************//**
1663 * Compute norm of a generic function
1664 *
1665 * \param[in] f - generic function
1666 *
1667 * \return out - norm
1668 ************************************************************/
generic_function_norm(const struct GenericFunction * f)1669 double generic_function_norm(const struct GenericFunction * f)
1670 {
1671 double out = generic_function_inner(f,f);
1672
1673 if (out < 0.0){
1674 fprintf(stderr, "Norm of a function cannot be negative %G\n",out);
1675 exit(1);
1676 /* generic_function_scale(0.0, f); */
1677 /* out = 0.0; */
1678 }
1679 //assert (out > -1e-15);
1680 return sqrt(out);
1681 }
1682
1683
1684 /********************************************************//**
1685 * Compute the weighted inner product between two generic functions
1686 *
1687 * \param[in] a - generic function
1688 * \param[in] b - generic function
1689 *
1690 * \return out - int a(x) b(x) w(x) dx
1691 ************************************************************/
generic_function_inner_weighted(const struct GenericFunction * a,const struct GenericFunction * b)1692 double generic_function_inner_weighted(const struct GenericFunction * a,
1693 const struct GenericFunction * b)
1694 {
1695 assert(a->fc == POLYNOMIAL);
1696 assert(b->fc == POLYNOMIAL);
1697 double out = orth_poly_expansion_inner_w(a->f,b->f);
1698
1699 return out;
1700 }
1701
1702 /********************************************************//**
1703 * Compute the sum of the inner products between
1704 * two arrays of generic functions
1705 *
1706 * \param[in] n - number of inner products
1707 * \param[in] lda - stride of functions to use in a
1708 * \param[in] a - first array of generic functions
1709 * \param[in] ldb - stride of functions to use in b
1710 * \param[in] b - second array of generic functions
1711 *
1712 * \return val - sum_{i=1^N} int a[ii*lda](x) b[ii*ldb](x) dx
1713 ************************************************************/
generic_function_inner_sum(size_t n,size_t lda,struct GenericFunction ** a,size_t ldb,struct GenericFunction ** b)1714 double generic_function_inner_sum(size_t n, size_t lda,
1715 struct GenericFunction ** a,
1716 size_t ldb,
1717 struct GenericFunction ** b)
1718 {
1719 double val = 0.0;
1720 for (size_t ii = 0; ii < n; ii++){
1721 val += generic_function_inner(a[ii*lda], b[ii*ldb]);
1722 }
1723 return val;
1724 }
1725
1726 /********************************************************//**
1727 * Compute the sum of the (weighted) inner products between
1728 * two arrays of generic functions
1729 *
1730 * \param[in] n - number of inner products
1731 * \param[in] lda - stride of functions to use in a
1732 * \param[in] a - first array of generic functions
1733 * \param[in] ldb - stride of functions to use in b
1734 * \param[in] b - second array of generic functions
1735 *
1736 * \return val - sum_{i=1^N} int a[ii*lda](x) b[ii*ldb](x) w(x) dx
1737 ************************************************************/
generic_function_inner_weighted_sum(size_t n,size_t lda,struct GenericFunction ** a,size_t ldb,struct GenericFunction ** b)1738 double generic_function_inner_weighted_sum(size_t n, size_t lda,
1739 struct GenericFunction ** a,
1740 size_t ldb,
1741 struct GenericFunction ** b)
1742 {
1743 double val = 0.0;
1744 size_t ii;
1745 for (ii = 0; ii < n; ii++){
1746 val += generic_function_inner_weighted(a[ii*lda], b[ii*ldb]);
1747 }
1748 return val;
1749 }
1750
1751
1752 /********************************************************//**
1753 * Compute the norm of the difference between two generic function
1754 *
1755 * \param[in] f1 - generic function
1756 * \param[in] f2 - generic function
1757 *
1758 * \return out - norm of difference
1759 ************************************************************/
generic_function_norm2diff(const struct GenericFunction * f1,const struct GenericFunction * f2)1760 double generic_function_norm2diff(const struct GenericFunction * f1,
1761 const struct GenericFunction * f2)
1762 {
1763 struct GenericFunction * f3 = generic_function_daxpby(1.0,f1,-1.0,f2);
1764 double out = generic_function_norm(f3);
1765 generic_function_free(f3); f3 = NULL;
1766 return out;
1767 }
1768
1769 /********************************************************//**
1770 * Compute the norm of the difference between two generic function arrays
1771 *
1772 * \param[in] n - number of elements
1773 * \param[in] f1 - generic function array
1774 * \param[in] inca - incremenent of first array
1775 * \param[in] f2 - generic function array
1776 * \param[in] incb - incremenent of second array
1777 *
1778 * \return out - norm of difference
1779 ************************************************************/
generic_function_array_norm2diff(size_t n,struct GenericFunction ** f1,size_t inca,struct GenericFunction ** f2,size_t incb)1780 double generic_function_array_norm2diff(
1781 size_t n, struct GenericFunction ** f1, size_t inca,
1782 struct GenericFunction ** f2, size_t incb)
1783 {
1784 double out = 0.0;
1785 size_t ii;
1786 for (ii = 0; ii < n; ii++){
1787 out += pow(generic_function_norm2diff(f1[ii*inca],f2[ii*incb]),2);
1788 }
1789 assert (out >= 0.0);
1790 return sqrt(out);
1791 }
1792
1793 /********************************************************//**
1794 Compute the integral of a generic function
1795
1796 \param[in] f - generic function
1797
1798 \return out - integral
1799
1800 \note Computes \f$ \int f(x) w(x) dx\f$ for every univariate function
1801 in the qmarray
1802
1803 w(x) depends on underlying parameterization
1804 for example, it is 1/2 for legendre (and default for others),
1805 gauss for hermite,etc
1806 ************************************************************/
generic_function_integral_weighted(const struct GenericFunction * f)1807 double generic_function_integral_weighted(
1808 const struct GenericFunction * f){
1809
1810 assert (f != NULL);
1811 assert (f->fc == POLYNOMIAL);
1812
1813 double out = orth_poly_expansion_integrate_weighted(f->f);
1814 return out;
1815 }
1816
1817 /********************************************************//**
1818 * Compute the integral of all the functions in a generic function array
1819 *
1820 * \param[in] n - number of functions
1821 * \param[in] lda - stride
1822 * \param[in] a - array of generic functions
1823 *
1824 * \return out - array of integrals
1825 ************************************************************/
1826 double *
generic_function_integral_array(size_t n,size_t lda,struct GenericFunction ** a)1827 generic_function_integral_array(size_t n,size_t lda,struct GenericFunction ** a)
1828 {
1829 double * out = calloc_double(n);
1830 size_t ii;
1831 for (ii = 0; ii < n; ii++){
1832 out[ii] = generic_function_integrate(a[ii*lda]);
1833 }
1834 return out;
1835 }
1836
1837 /********************************************************//**
1838 * Compute the norm of an array of generic functions
1839 *
1840 * \param[in] n - number of functions
1841 * \param[in] lda - stride of functions to use in a
1842 * \param[in] a - functions
1843 *
1844 * \return val -sqrt(sum_{i=1^N} int a[ii*lda](x)^2 ) dx)
1845 ***********************************************************/
generic_function_array_norm(size_t n,size_t lda,struct GenericFunction ** a)1846 double generic_function_array_norm(size_t n, size_t lda,
1847 struct GenericFunction ** a)
1848 {
1849
1850 double val = 0.0;
1851 size_t ii;
1852 for (ii = 0; ii < n; ii++){
1853 val += pow(generic_function_norm(a[lda*ii]),2.0);
1854 }
1855 //val = generic_function_inner_sum(n,lda,a,lda,a);
1856
1857 return sqrt(val);
1858 }
1859
1860
1861
1862 /********************************************************//**
1863 Compute the index, location and value of the maximum, in absolute value,
1864 element of a generic function array
1865
1866 \param[in] n - number of functions
1867 \param[in] lda - stride
1868 \param[in] a - array of functions
1869 \param[in,out] ind - index of maximum
1870 \param[in,out] x - location of maximum
1871 \param[in] optargs - optimization arguments
1872
1873 \return maxval - absolute value of the maximum
1874 ************************************************************/
1875 double
generic_function_array_absmax(size_t n,size_t lda,struct GenericFunction ** a,size_t * ind,double * x,void * optargs)1876 generic_function_array_absmax(size_t n, size_t lda,
1877 struct GenericFunction ** a,
1878 size_t * ind, double * x,
1879 void * optargs)
1880 {
1881 size_t ii = 0;
1882 *ind = ii;
1883 //printf("do absmax\n");
1884 //print_generic_function(a[ii],0,NULL);
1885 double maxval = generic_function_absmax(a[ii],x,optargs);
1886 //printf("maxval=%G\n",maxval);
1887 double tempval, tempx;
1888 for (ii = 1; ii < n; ii++){
1889 tempval = generic_function_absmax(a[ii*lda],&tempx,optargs);
1890 if (tempval > maxval){
1891 maxval = tempval;
1892 *x = tempx;
1893 *ind = ii;
1894 }
1895 }
1896 return maxval;
1897 }
1898
1899 /********************************************************//**
1900 Compute the index, location and value of the maximum, in absolute value,
1901 element of a generic function array (Pivot Based)
1902
1903 \param[in] n - number of functions
1904 \param[in] lda - stride
1905 \param[in] a - array of functions
1906 \param[in,out] piv - pivot
1907 \param[in] optargs - optimization arguments
1908
1909 \return maxval - absolute value of the maximum
1910 ************************************************************/
1911 double
generic_function_array_absmax_piv(size_t n,size_t lda,struct GenericFunction ** a,struct Pivot * piv,void * optargs)1912 generic_function_array_absmax_piv(size_t n, size_t lda,
1913 struct GenericFunction ** a,
1914 struct Pivot * piv,
1915 void * optargs)
1916 {
1917 size_t ii = 0;
1918 pivot_set_ind(piv,ii);
1919 //printf("do absmax\n");
1920 //print_generic_function(a[ii],0,NULL);
1921 size_t size = pivot_get_size(piv);
1922 void * x = pivot_get_loc(piv);
1923 double maxval = generic_function_absmax_gen(a[ii],x,size,optargs);
1924 //printf("maxval=%G\n",maxval);
1925 if (size == sizeof(double)){
1926 double tempval, tempx;
1927 for (ii = 1; ii < n; ii++){
1928 tempval = generic_function_absmax_gen(a[ii*lda],&tempx,size,optargs);
1929 if (tempval > maxval){
1930 maxval = tempval;
1931 *(double *)(x) = tempx;
1932 pivot_set_ind(piv,ii);
1933 }
1934 }
1935 }
1936 else if (size == sizeof(size_t)){
1937 double tempval;
1938 size_t tempx;
1939 for (ii = 1; ii < n; ii++){
1940 tempval = generic_function_absmax_gen(a[ii*lda],&tempx,size,optargs);
1941 if (tempval > maxval){
1942 maxval = tempval;
1943 *(size_t *)(x) = tempx;
1944 pivot_set_ind(piv,ii);
1945 }
1946 }
1947 }
1948 else{
1949 fprintf(stderr, "Cannot perform generic_function_array_absmax_piv\n");
1950 fprintf(stderr, "with the specified elements of size %zu\n",size);
1951 exit(1);
1952 }
1953 return maxval;
1954 }
1955
1956
1957
1958
1959
1960
1961
1962 ////////////////////////////////////////////////////////////////////
1963 // High dimensional helper functions
1964
1965 /********************************************************//**
1966 Allocate memory for a fiber cut
1967
1968 \param totdim [in] - total dimension of underlying function
1969 \param dim [in] - dimension along which fiber is obtained
1970 ************************************************************/
alloc_fiber_cut(size_t totdim,size_t dim)1971 struct FiberCut * alloc_fiber_cut(size_t totdim, size_t dim)
1972 {
1973 struct FiberCut * fcut;
1974 if (NULL == ( fcut = malloc(sizeof(struct FiberCut)))){
1975 fprintf(stderr, "failed to allocate fiber_cut.\n");
1976 exit(1);
1977 }
1978 fcut->totdim = totdim;
1979 fcut->dimcut = dim;
1980 fcut->vals = calloc_double(totdim);
1981
1982 return fcut;
1983 }
1984
1985 /********************************************************//**
1986 Free memory allocated to a fiber cut
1987
1988 \param fc [inout] - fiber cut
1989 ************************************************************/
fiber_cut_free(struct FiberCut * fc)1990 void fiber_cut_free(struct FiberCut * fc)
1991 {
1992 free(fc->vals); fc->vals = NULL;
1993 free(fc); fc = NULL;
1994 }
1995
1996 /********************************************************//**
1997 Free memory allocated to an array of fiber cuts
1998
1999 \param n [in] - number of fiber cuts
2000 \param fc [inout] - array of fiber cuts
2001 ************************************************************/
fiber_cut_array_free(size_t n,struct FiberCut ** fc)2002 void fiber_cut_array_free(size_t n, struct FiberCut ** fc)
2003 {
2004 size_t ii;
2005 for (ii = 0; ii < n; ii++){
2006 fiber_cut_free(fc[ii]);
2007 fc[ii] = NULL;
2008 }
2009 free(fc); fc = NULL;
2010 }
2011
2012 /********************************************************//**
2013 Generate a fibercut of a two dimensional function
2014
2015 \param f [in] - function to cut
2016 \param args [in] - function arguments
2017 \param dim [in] - dimension along which we obtain the cut
2018 \param val [in] - value of the input which is not *dim*
2019
2020 \return fcut - struct necessary for computing values in the cut
2021 ************************************************************/
2022 struct FiberCut *
fiber_cut_init2d(double (* f)(double,double,void *),void * args,size_t dim,double val)2023 fiber_cut_init2d( double (*f)(double, double, void *), void * args,
2024 size_t dim, double val)
2025 {
2026 struct FiberCut * fcut = alloc_fiber_cut(2,dim);
2027 fcut->fpoint.f2d = f;
2028 fcut->args = args;
2029 fcut->ftype_flag=0;
2030 if (dim == 0){
2031 fcut->vals[1] = val;
2032 }
2033 else{
2034 fcut->vals[0] = val;
2035 }
2036 return fcut;
2037 }
2038
2039 /********************************************************//**
2040 Generate an array fibercuts of a two dimensional function
2041
2042 \param[in] f - function to cut
2043 \param[in] args - function arguments
2044 \param[in] dim - dimension along which we obtain the cut
2045 \param[in] n - number of fibercuts
2046 \param[in] val - values of the input for each fibercut
2047
2048 \return fcut - array of struct necessary for computing values in the cut
2049 ***************************************************************/
2050 struct FiberCut **
fiber_cut_2darray(double (* f)(double,double,void *),void * args,size_t dim,size_t n,const double * val)2051 fiber_cut_2darray( double (*f)(double, double, void *), void * args,
2052 size_t dim, size_t n, const double * val)
2053 {
2054 struct FiberCut ** fcut;
2055 if (NULL == ( fcut = malloc(n *sizeof(struct FiberCut *)))){
2056 fprintf(stderr, "failed to allocate fiber_cut.\n");
2057 exit(1);
2058 }
2059 size_t ii;
2060 for (ii = 0; ii < n; ii++){
2061 fcut[ii] = alloc_fiber_cut(2,dim);
2062 fcut[ii]->fpoint.f2d = f;
2063 fcut[ii]->args = args;
2064 fcut[ii]->ftype_flag = 0;
2065 if (dim == 0){
2066 fcut[ii]->vals[1] = val[ii];
2067 }
2068 else{
2069 fcut[ii]->vals[0] = val[ii];
2070 }
2071
2072 }
2073 return fcut;
2074 }
2075
2076 /********************************************************//**
2077 Generate an array fibercuts of a n-dimensional functions
2078
2079 \param[in] f - function to cut
2080 \param[in] args - function arguments
2081 \param[in] totdim - total number of dimensions
2082 \param[in] dim - dimension along which we obtain the cut
2083 \param[in] n - number of fibercuts
2084 \param[in] val - array of values of the inputs for each fibercut
2085
2086 \return fcut -array of struct necessary for computing values in the cut
2087 ***************************************************************/
2088 struct FiberCut **
fiber_cut_ndarray(double (* f)(double *,void *),void * args,size_t totdim,size_t dim,size_t n,double ** val)2089 fiber_cut_ndarray( double (*f)(double *, void *), void * args,
2090 size_t totdim, size_t dim, size_t n, double ** val)
2091 {
2092 struct FiberCut ** fcut;
2093 if (NULL == ( fcut = malloc(n *sizeof(struct FiberCut *)))){
2094 fprintf(stderr, "failed to allocate fiber_cut.\n");
2095 exit(1);
2096 }
2097 size_t ii;
2098 // printf("vals are \n");
2099 for (ii = 0; ii < n; ii++){
2100 fcut[ii] = alloc_fiber_cut(totdim,dim);
2101 fcut[ii]->fpoint.fnd = f;
2102 fcut[ii]->args = args;
2103 fcut[ii]->ftype_flag = 1;
2104 memmove(fcut[ii]->vals, val[ii], totdim*sizeof(double));
2105 // dprint(totdim,val[ii]);
2106
2107 }
2108 return fcut;
2109 }
2110
2111 /********************************************************//**
2112 Evaluate a fiber of a two dimensional function
2113
2114 \param x [in] - value at which to evaluate
2115 \param vfcut [in] - void pointer to fiber_cut structure
2116
2117 \return val - value of the function
2118 ************************************************************/
fiber_cut_eval2d(double x,void * vfcut)2119 double fiber_cut_eval2d(double x, void * vfcut){
2120
2121 struct FiberCut * fcut = vfcut;
2122
2123 double val;
2124 if (fcut->dimcut == 0){
2125 val = fcut->fpoint.f2d(x, fcut->vals[1], fcut->args);
2126 }
2127 else{
2128 val = fcut->fpoint.f2d(fcut->vals[0], x, fcut->args);
2129 }
2130 return val;
2131 }
2132
2133 /********************************************************//**
2134 Evaluate a fiber of an n dimensional function
2135
2136 \param[in] x - value at which to evaluate
2137 \param[in] vfcut - void pointer to fiber_cut structure
2138
2139 \return val - value of the function
2140 ************************************************************/
fiber_cut_eval(double x,void * vfcut)2141 double fiber_cut_eval(double x, void * vfcut){
2142
2143 struct FiberCut * fcut = vfcut;
2144
2145 double val;
2146 fcut->vals[fcut->dimcut] = x;
2147 val = fcut->fpoint.fnd(fcut->vals, fcut->args);
2148 return val;
2149 }
2150
2151
2152 /////////////////////////////////////////////////////////
2153 // Utilities
2154
2155
2156
2157 /***********************************************************
2158 Generate a set of orthonormal arrays of functions for helping
2159 generate an orthonormal qmarray
2160
2161 \param[in] fc - function class
2162 \param[in] st - function class sub_type
2163 \param[in] nrows - number of rows
2164 \param[in] ncols - number of columns
2165 \param[in] lb - lower bound on 1d functions
2166 \param[in] ub - upper bound on 1d functions
2167
2168 \note
2169 - Not super efficient because of copies
2170 ***************************************************************/
2171 /* void */
2172 /* generic_function_array_orth1d_columns(struct GenericFunction ** f, */
2173 /* struct GenericFunction ** funcs, */
2174 /* enum function_class fc, */
2175 /* void * st, size_t nrows, */
2176 /* size_t ncols, double lb, */
2177 /* double ub) */
2178 /* { */
2179
2180 /* struct Interval ob; */
2181 /* ob.lb = lb; */
2182 /* ob.ub = ub; */
2183 /* size_t jj,kk; */
2184 /* generic_function_array_orth(ncols, fc, st, funcs, &ob); */
2185 /* struct GenericFunction * zero = */
2186 /* generic_function_constant(0.0,fc,st,lb,ub,NULL); */
2187 /* size_t onnon = 0; */
2188 /* size_t onorder = 0; */
2189 /* for (jj = 0; jj < ncols; jj++){ */
2190 /* f[jj*nrows+onnon] = generic_function_copy(funcs[onorder]); */
2191 /* for (kk = 0; kk < onnon; kk++){ */
2192 /* f[jj*nrows+kk] = generic_function_copy(zero); */
2193 /* } */
2194 /* for (kk = onnon+1; kk < nrows; kk++){ */
2195 /* f[jj*nrows+kk] = generic_function_copy(zero); */
2196 /* } */
2197 /* onnon = onnon+1; */
2198 /* if (onnon == nrows){ */
2199 /* //generic_function_free(funcs[onorder]); */
2200 /* //funcs[onorder] = NULL; */
2201 /* onorder = onorder+1; */
2202 /* onnon = 0; */
2203 /* } */
2204 /* } */
2205 /* generic_function_free(zero); zero = NULL; */
2206
2207 /* } */
2208
2209
2210
2211
2212 //////////////////////////////////////////////////
2213 //////////////////////////////////////////////////
2214 //////////////////////////////////////////////////
2215 ///// Regression //////////
2216 //////////////////////////////////////////////////
2217 //////////////////////////////////////////////////
2218 //////////////////////////////////////////////////
2219
2220
2221 /********************************************************//**
2222 Create a regression options
2223
2224 \param[in] atype - approximation type
2225 \param[in] rtype - regression problem type
2226 \param[in] N - number of training samples
2227 \param[in] x - location of training samples
2228 \param[in] y - values at training samples
2229
2230 \return opts - regression options
2231 ************************************************************/
2232 struct Regress1DOpts *
regress_1d_opts_create(enum approx_type atype,enum regress_type rtype,size_t N,const double * x,const double * y)2233 regress_1d_opts_create(enum approx_type atype, enum regress_type rtype,
2234 size_t N, const double * x, const double * y)
2235 {
2236 struct Regress1DOpts * opts = malloc(sizeof(struct Regress1DOpts));
2237 if (opts == NULL){
2238 fprintf(stderr, "Error allocating regression options\n");
2239 exit(1);
2240 }
2241 opts->atype = atype;
2242 opts->rtype = rtype;
2243 opts->N = N;
2244 opts->x = x;
2245 opts->y = y;
2246
2247 opts->aopts = NULL;
2248
2249 opts->nparam = 0;
2250 opts->init_param = NULL;
2251
2252 opts->gf = NULL;
2253
2254 opts->eval = calloc_double(N);
2255 opts->grad = NULL;
2256 opts->resid = calloc_double(N);
2257
2258 // regularization options
2259 opts->reg_param_set = 0;
2260 opts->lambda = 0.0;
2261 opts->decay_type = NONE;
2262 opts->coeff_decay_param = 1.0;
2263
2264 return opts;
2265 }
2266
2267 /********************************************************//**
2268 Destroy regression options
2269 ************************************************************/
regress_1d_opts_destroy(struct Regress1DOpts * opts)2270 void regress_1d_opts_destroy(struct Regress1DOpts * opts)
2271 {
2272 if (opts != NULL){
2273 generic_function_free(opts->gf); opts->gf = NULL;
2274 free(opts->eval); opts->eval = NULL;
2275 free(opts->grad); opts->grad = NULL;
2276 free(opts->resid); opts->resid = NULL;
2277 free(opts); opts = NULL;
2278 }
2279 }
2280
2281 /********************************************************//**
2282 Add a parametric form to learn
2283
2284 \param[in] opts - regression options structure
2285 \param[in] fc - regression problem type
2286 \param[in] aopts - parametric approximation options
2287 ************************************************************/
regress_1d_opts_set_parametric_form(struct Regress1DOpts * opts,enum function_class fc,void * aopts)2288 void regress_1d_opts_set_parametric_form(
2289 struct Regress1DOpts * opts, enum function_class fc, void * aopts)
2290 {
2291 assert(opts != NULL);
2292 assert(aopts != NULL);
2293
2294 opts->fc = fc;
2295 opts->aopts = aopts;
2296
2297 GF_OPTS_SWITCH_ONEOUT(get_nparams, fc, opts->nparam, aopts)
2298 opts->grad = calloc_double(opts->nparam);
2299 }
2300
2301 /********************************************************//**
2302 Add starting parameters for optimization
2303 ************************************************************/
regress_1d_opts_set_initial_parameters(struct Regress1DOpts * opts,const double * param)2304 void regress_1d_opts_set_initial_parameters(
2305 struct Regress1DOpts * opts, const double * param)
2306 {
2307 assert (opts != NULL);
2308 opts->init_param = param;
2309
2310 opts->gf = generic_function_create_with_params(opts->fc,opts->aopts,opts->nparam,param);
2311 }
2312
2313 /********************************************************//**
2314 Set regularization penalty
2315 ************************************************************/
regress_1d_opts_set_regularization_penalty(struct Regress1DOpts * opts,double lambda)2316 void regress_1d_opts_set_regularization_penalty(
2317 struct Regress1DOpts * opts, double lambda)
2318 {
2319 assert (opts != NULL);
2320 opts->reg_param_set = 1;
2321 opts->lambda = lambda;
2322 }
2323
2324 /********************************************************//**
2325 Set RKHS decay
2326 ************************************************************/
regress_1d_opts_set_RKHS_decay_rate(struct Regress1DOpts * opts,enum coeff_decay_type decay_type,double lambda)2327 void regress_1d_opts_set_RKHS_decay_rate(
2328 struct Regress1DOpts * opts, enum coeff_decay_type decay_type, double lambda)
2329 {
2330 assert (opts != NULL);
2331 opts->decay_type = decay_type;
2332
2333 if (decay_type == ALGEBRAIC){
2334 if ((lambda < 1e-15) || (lambda > 1)){
2335 fprintf(stderr,"For algebraic decay of RKHS must specify decay rate in (0,1)\n");
2336 fprintf(stderr,"\t Currently specified as %G\n",lambda);
2337 exit(1);
2338 }
2339 else{
2340 opts->coeff_decay_param = lambda;
2341 }
2342 }
2343 else if (decay_type == EXPONENTIAL){
2344 if (lambda < 0){
2345 fprintf(stderr,"For exponential decay of RKHS must specify decay rate > 0\n");
2346 fprintf(stderr,"\t Currently specified as %G\n",lambda);
2347 exit(1);
2348 }
2349 else{
2350 opts->coeff_decay_param = lambda;
2351 }
2352 }
2353 else{
2354 fprintf(stderr,"Do not recognized RKHS decay type %d\n",decay_type);
2355 exit(1);
2356 }
2357 }
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367 /********************************************************//**
2368 Take a gradient of the squared norm of a generic function
2369 with respect to its parameters, and add a scaled version
2370 of this gradient to *grad*
2371
2372 \param[in] gf - generic function
2373 \param[in] scale - scaling for additional gradient
2374 \param[in,out] grad - gradient, on output adds scale * new_grad
2375
2376 \return 0 - success, 1 -failure
2377 ************************************************************/
2378 int
generic_function_squared_norm_param_grad(const struct GenericFunction * gf,double scale,double * grad)2379 generic_function_squared_norm_param_grad(const struct GenericFunction * gf,
2380 double scale, double * grad)
2381 {
2382
2383 enum function_class fc = generic_function_get_fc(gf);
2384 int res = 1;
2385 GF_SWITCH_THREEOUT(squared_norm_param_grad, fc, res, gf->f, scale, grad)
2386 return res;
2387 }
2388
2389 /********************************************************//**
2390 Norm in the RKHS (instead of L2)
2391
2392 \param[in] gf - generic function
2393 \param[in] decay_type - type of decay
2394 \param[in] decay_param - parameter of decay
2395
2396 \return 0 - success, 1 -failure
2397 ************************************************************/
2398 double
generic_function_rkhs_squared_norm(const struct GenericFunction * gf,enum coeff_decay_type decay_type,double decay_param)2399 generic_function_rkhs_squared_norm(const struct GenericFunction * gf,
2400 enum coeff_decay_type decay_type,
2401 double decay_param)
2402 {
2403
2404 enum function_class fc = generic_function_get_fc(gf);
2405 assert (fc == POLYNOMIAL);
2406 double out = orth_poly_expansion_rkhs_squared_norm(gf->f,decay_type,decay_param);
2407
2408 return out;
2409 }
2410
2411
2412 /********************************************************//**
2413 Take a gradient of the norm in the RKHS (instead of L2)
2414
2415 \param[in] gf - generic function
2416 \param[in] scale - scaling for additional gradient
2417 \param[in] decay_type - type of decay
2418 \param[in] decay_param - parameter of decay
2419 \param[in,out] grad - gradient, on output adds scale * new_grad
2420
2421 \return 0 - success, 1 -failure
2422 ************************************************************/
2423 int
generic_function_rkhs_squared_norm_param_grad(const struct GenericFunction * gf,double scale,enum coeff_decay_type decay_type,double decay_param,double * grad)2424 generic_function_rkhs_squared_norm_param_grad(const struct GenericFunction * gf,
2425 double scale, enum coeff_decay_type decay_type,
2426 double decay_param, double * grad)
2427 {
2428
2429 enum function_class fc = generic_function_get_fc(gf);
2430 assert (fc == POLYNOMIAL);
2431 int res = orth_poly_expansion_rkhs_squared_norm_param_grad(
2432 gf->f,scale,decay_type,decay_param,grad);
2433
2434
2435 return res;
2436 }
2437
2438
2439
2440
2441 /********************************************************//**
2442 LS regression objective function
2443 ************************************************************/
param_LSregress_cost(size_t dim,const double * param,double * grad,void * arg)2444 double param_LSregress_cost(size_t dim, const double * param, double * grad, void * arg)
2445 {
2446
2447 struct Regress1DOpts * opts = arg;
2448
2449 assert (opts->nparam == dim);
2450 assert (opts->gf != NULL);
2451 // update function
2452 /* printf("update param\n"); */
2453 /* printf("\t old = "); */
2454 /* print_generic_function(opts->gf,0,NULL); */
2455 /* printf("\t param = "); dprint(dim,param); */
2456 generic_function_update_params(opts->gf,dim,param);
2457
2458 /* printf("evaluate\n"); */
2459 for (size_t ii = 0; ii < opts->N; ii++){
2460 opts->eval[ii] = generic_function_1d_eval(opts->gf,opts->x[ii]);
2461 }
2462
2463 /* printf("compute resid\n"); */
2464 double out = 0.0;
2465 for (size_t ii = 0; ii < opts->N; ii++){
2466 opts->resid[ii] = opts->y[ii]-opts->eval[ii];
2467 out += opts->resid[ii] * opts->resid[ii];
2468 }
2469 out *= 0.5;
2470
2471 if (grad != NULL){
2472 /* printf("grad is not null!\n"); */
2473 for (size_t ii = 0; ii < dim; ii++){
2474 grad[ii] = 0.0;
2475 }
2476 for (size_t jj = 0; jj < opts->N; jj++){
2477 int res = generic_function_param_grad_eval(opts->gf,1,
2478 opts->x+jj,
2479 opts->grad);
2480 assert (res == 0);
2481 for (size_t ii = 0; ii < dim; ii++){
2482 grad[ii] += opts->resid[jj] * (-1.0)*opts->grad[ii];
2483 }
2484 }
2485 /* printf("done\n"); */
2486 }
2487
2488 return out;
2489 }
2490
2491 /********************************************************//**
2492 Ridge regression
2493 ************************************************************/
param_RLS2regress_cost(size_t dim,const double * param,double * grad,void * arg)2494 double param_RLS2regress_cost(size_t dim, const double * param, double * grad, void * arg)
2495 {
2496
2497 struct Regress1DOpts * opts = arg;
2498
2499 // first part (recall this function updates parameters already!)
2500 double ls_portion = param_LSregress_cost(dim,param,grad,arg);
2501
2502 // second part
2503 double regularization = generic_function_inner(opts->gf,opts->gf);
2504
2505 double out = ls_portion + 0.5*opts->lambda * regularization;
2506
2507 if (grad != NULL){
2508 int res = generic_function_squared_norm_param_grad(opts->gf,0.5*opts->lambda,grad);
2509 assert (res == 0);
2510 }
2511
2512 return out;
2513 }
2514
2515 /********************************************************//**
2516 Ridge regression penalizing second derivative
2517 ************************************************************/
param_RLSD2regress_cost(size_t dim,const double * param,double * grad,void * arg)2518 double param_RLSD2regress_cost(size_t dim, const double * param, double * grad, void * arg)
2519 {
2520
2521 struct Regress1DOpts * opts = arg;
2522
2523 // first part (recall this function updates parameters already!)
2524 double ls_portion = param_LSregress_cost(dim,param,grad,arg);
2525
2526 // second part
2527 struct GenericFunction * gf1 = generic_function_deriv(opts->gf);
2528 struct GenericFunction * gf2 = generic_function_deriv(gf1);
2529 double regularization = generic_function_inner(gf2,gf2);
2530
2531 double out = ls_portion + 0.5*opts->lambda * regularization;
2532
2533 if (grad != NULL){
2534 int res = generic_function_squared_norm_param_grad(gf2,0.5*opts->lambda,grad);
2535 assert (res == 0);
2536 }
2537
2538 generic_function_free(gf1); gf1 = NULL;
2539 generic_function_free(gf2); gf2 = NULL;
2540 return out;
2541 }
2542
2543 /********************************************************//**
2544 Ridge regression with an RKHS penalty
2545 ************************************************************/
param_RLSRKHSregress_cost(size_t dim,const double * param,double * grad,void * arg)2546 double param_RLSRKHSregress_cost(size_t dim, const double * param, double * grad, void * arg)
2547 {
2548
2549 struct Regress1DOpts * opts = arg;
2550
2551 // first part (recall this function updates parameters already!)
2552 double ls_portion = param_LSregress_cost(dim,param,grad,arg);
2553
2554 // second part
2555 double regularization =
2556 generic_function_rkhs_squared_norm(opts->gf,
2557 opts->decay_type,
2558 opts->coeff_decay_param);
2559
2560 double out = ls_portion + 0.5*opts->lambda * regularization;
2561
2562 if (grad != NULL){
2563 int res = generic_function_rkhs_squared_norm_param_grad(opts->gf,opts->lambda,
2564 opts->decay_type,
2565 opts->coeff_decay_param,grad);
2566 assert (res == 0);
2567 }
2568
2569 return out;
2570 }
2571
2572
2573