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 polynomials.c
40 * Provides routines for manipulating orthogonal polynomials
41 */
42
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <math.h>
46 #include <string.h>
47 #include <float.h>
48 #include <assert.h>
49 #include "futil.h"
50
51 //#define ZEROTHRESH 1e-20
52 /* #define ZEROTHRESH 1e0 * DBL_EPSILON */
53 /* #define ZEROTHRESH 1e-3 * DBL_EPSILON */
54 #define ZEROTHRESH 1e0 * DBL_EPSILON
55 /* #define ZEROTHRESH 1e0 * DBL_MIN */
56 /* #define ZEROTHRESH 1e-200 */
57 //#define ZEROTHRESH 0.0
58 /* #define ZEROTHRESH 1e0 * DBL_EPSILON */
59 //#define ZEROTHRESH 1e-12
60
61 #include "stringmanip.h"
62 #include "array.h"
63 #include "polynomials.h"
64 #include "hpoly.h"
65 #include "lib_quadrature.h"
66 #include "linalg.h"
67 #include "legtens.h"
68 #include "fourier.h"
69
70 enum SPACE_MAP {SM_LINEAR,SM_ROSENBLATT};
71
72 struct SpaceMapping
73 {
74 enum SPACE_MAP map;
75 int set;
76 double lin_slope;
77 double lin_offset;
78 double inv_lin_slope;
79 double inv_lin_offset;
80 };
81
82
space_mapping_create(enum SPACE_MAP map_type)83 static struct SpaceMapping * space_mapping_create(enum SPACE_MAP map_type)
84 {
85 struct SpaceMapping * map = malloc(sizeof(struct SpaceMapping));
86 if (map == NULL){
87 fprintf(stderr,"Failure to allocate space mapping for polynomials\n");
88 exit(1);
89 }
90
91 map->map = map_type;
92 map->set = 0;
93 if (map_type == SM_LINEAR){
94 map->lin_slope = 1;
95 map->lin_offset = 0;
96 map->inv_lin_slope = 1;
97 map->inv_lin_offset = 0;
98 }
99 else if (map_type == SM_ROSENBLATT){
100 fprintf(stderr,"Rosenblatt mapping not yet implemented\n");
101 exit(1);
102 }
103 else{
104 fprintf(stderr,"Mapping type %d is not defined\n",map_type);
105 exit(1);
106 }
107 return map;
108 }
109
space_mapping_copy(const struct SpaceMapping * map)110 static struct SpaceMapping * space_mapping_copy(const struct SpaceMapping * map)
111 {
112 struct SpaceMapping * map_out = space_mapping_create(map->map);
113 map_out->set = map->set;
114 map_out->lin_slope = map->lin_slope;
115 map_out->lin_offset = map->lin_offset;
116 map_out->inv_lin_slope = map->inv_lin_slope;
117 map_out->inv_lin_offset = map->inv_lin_offset;
118 return map_out;
119 }
120
space_mapping_free(struct SpaceMapping * map)121 static void space_mapping_free(struct SpaceMapping * map)
122 {
123 if (map != NULL){
124 free(map); map = NULL;
125 }
126 }
127
128 /********************************************************//**
129 * Serialize the mapping
130 *************************************************************/
131 static unsigned char *
serialize_space_mapping(unsigned char * ser,struct SpaceMapping * p,size_t * totSizeIn)132 serialize_space_mapping(unsigned char * ser,
133 struct SpaceMapping * p,
134 size_t * totSizeIn)
135 {
136 // order is ptype->lower_bound->upper_bound->orth_poly->coeff
137
138 size_t totsize =
139 sizeof(enum SPACE_MAP) +
140 sizeof(int) +
141 4*sizeof(double);
142
143 if (totSizeIn != NULL){
144 *totSizeIn = totsize;
145 return ser;
146 }
147
148 unsigned char * ptr = serialize_int(ser, p->map);
149 ptr = serialize_int(ptr,p->set);
150 ptr = serialize_double(ptr, p->lin_slope);
151 ptr = serialize_double(ptr, p->lin_offset);
152 ptr = serialize_double(ptr, p->inv_lin_slope);
153 ptr = serialize_double(ptr, p->inv_lin_offset);
154 return ptr;
155 }
156
157 /********************************************************//**
158 * Deserialize space mapping
159 *************************************************************/
160 static unsigned char *
deserialize_space_mapping(unsigned char * ser,struct SpaceMapping ** smap)161 deserialize_space_mapping(
162 unsigned char * ser,
163 struct SpaceMapping ** smap)
164 {
165 enum SPACE_MAP map;
166 int map_int;
167
168 unsigned char * ptr = deserialize_int(ser,&map_int);
169 map = (enum SPACE_MAP) map_int;
170
171 *smap = space_mapping_create(map);
172 ptr = deserialize_int(ptr, &((*smap)->set));
173 ptr = deserialize_double(ptr, &((*smap)->lin_slope));
174 ptr = deserialize_double(ptr, &((*smap)->lin_offset));
175 ptr = deserialize_double(ptr, &((*smap)->inv_lin_slope));
176 ptr = deserialize_double(ptr, &((*smap)->inv_lin_offset));
177
178 return ptr;
179 }
180
181
182 /********************************************************//**
183 Save a mapping in text format
184 ************************************************************/
space_mapping_savetxt(const struct SpaceMapping * f,FILE * stream,size_t prec)185 static void space_mapping_savetxt(const struct SpaceMapping * f,
186 FILE * stream, size_t prec)
187 {
188 assert (f != NULL);
189 fprintf(stream,"%d ",f->map);
190 fprintf(stream,"%d ",f->set);
191 fprintf(stream,"%3.*G ",(int)prec,f->lin_slope);
192 fprintf(stream,"%3.*G ",(int)prec,f->lin_offset);
193 fprintf(stream,"%3.*G ",(int)prec,f->inv_lin_slope);
194 fprintf(stream,"%3.*G ",(int)prec,f->inv_lin_offset);
195 }
196
197 /********************************************************//**
198 Load a mapping in text format
199 ************************************************************/
200 static struct SpaceMapping *
space_mapping_loadtxt(FILE * stream)201 space_mapping_loadtxt(FILE * stream)//l, size_t prec)
202 {
203
204 enum SPACE_MAP maptype;
205 int maptypeint;
206 int num = fscanf(stream,"%d ",&maptypeint);
207 maptype = (enum SPACE_MAP)maptypeint;
208 assert (num == 1);
209
210 struct SpaceMapping * map = space_mapping_create(maptype);
211 num = fscanf(stream,"%d ",&(map->set));
212 assert (num == 1);
213
214 num = fscanf(stream,"%lG ",&(map->lin_slope));
215 assert (num == 1);
216 num = fscanf(stream,"%lG ",&(map->lin_offset));
217 assert (num == 1);
218 num = fscanf(stream,"%lG ",&(map->inv_lin_slope));
219 assert (num == 1);
220 num = fscanf(stream,"%lG ",&(map->inv_lin_offset));
221 assert (num == 1);
222
223 return map;
224 }
225
226
227 // original space to normalized space
space_mapping_map(struct SpaceMapping * map,double x)228 double space_mapping_map(struct SpaceMapping * map, double x)
229 {
230 if (map->map == SM_LINEAR){
231 return map->lin_slope * x + map->lin_offset;
232 }
233 else{
234 fprintf(stderr, "NONLIENAR MAPPINGS NOT YET IMPLEMENTED\n");
235 exit(1);
236 }
237 }
238
space_mapping_map_deriv(struct SpaceMapping * map,double x)239 double space_mapping_map_deriv(struct SpaceMapping * map, double x)
240 {
241 (void)(x);
242 if (map->map == SM_LINEAR){
243 return map->lin_slope;
244 }
245 else{
246 fprintf(stderr, "NONLIENAR MAPPINGS NOT YET IMPLEMENTED\n");
247 exit(1);
248 }
249 }
250
251 // normalized space to original space
space_mapping_map_inverse(struct SpaceMapping * map,double x)252 double space_mapping_map_inverse(struct SpaceMapping * map, double x)
253 {
254 if (map->map == SM_LINEAR){
255 return map->inv_lin_slope * x + map->inv_lin_offset;
256 }
257 else{
258 fprintf(stderr, "NONLINEAR MAPPINGS NOT YET IMPLEMENTED\n");
259 exit(1);
260 }
261 }
262
space_mapping_map_inverse_deriv(struct SpaceMapping * map,double x)263 double space_mapping_map_inverse_deriv(struct SpaceMapping * map, double x)
264 {
265 (void)(x);
266 if (map->map == SM_LINEAR){
267 return map->inv_lin_slope;
268 }
269 else{
270 fprintf(stderr, "NONLIENAR MAPPINGS NOT YET IMPLEMENTED\n");
271 exit(1);
272 }
273 }
274
275
276 // Recurrence relationship sequences
zero_seq(size_t n)277 inline static double zero_seq(size_t n){ (void)n; return (0.0); }
278 /* inline static double one_seq(size_t n) { (void)n; return (1.0); } */
none_seq(size_t n)279 inline static double none_seq(size_t n){ (void)n; return (-1.0); }
two_seq(size_t n)280 inline static double two_seq(size_t n) { (void)n; return (2.0); }
281 /* inline static double n_seq(size_t n) { return ((double) n); } */
282 /* inline static double nn_seq (size_t n) { return -n_seq(n); } */
283 /* inline static double lega_seq (size_t n) { return ( (double)(2.0 * n -1.0) / (double) n);} */
284 /* inline static double legc_seq (size_t n) { return ( -((double)n - 1.0)/ (double) n );} */
285
286
287 static const double legaseqnorm[201] = {
288 0.0000000000000000000000000,
289 1.7320508075688772935737253,
290 1.9364916731037084426068906,
291 1.9720265943665386808867149,
292 1.9843134832984429429111536,
293 1.9899748742132399093648226,
294 1.9930434571835663367907546,
295 1.9948914348241344529156019,
296 1.9960899278339139998987573,
297 1.9969111950679364953821493,
298 1.9974984355438178915695749,
299 1.9979328159850827787737126,
300 1.9982631347136331423841246,
301 1.9985201625794738002940554,
302 1.9987240828047460811708186,
303 1.9988885800753266486651932,
304 1.9990231989649344737647318,
305 1.9991347609372268759927310,
306 1.9992282461607312885921994,
307 1.9993073592865872621432075,
308 1.9993749023132204957258970,
309 1.9994330262111444065009289,
310 1.9994834043566154743211405,
311 1.9995273543594641291838362,
312 1.9995659251169712031020662,
313 1.9995999599919979994094507,
314 1.9996301433163012802865163,
315 1.9996570350656427370343235,
316 1.9996810970242025987542861,
317 1.9997027127445487499745463,
318 1.9997222029294191164312117,
319 1.9997398373972733498548163,
320 1.9997558444720195391930417,
321 1.9997704184116869836722319,
322 1.9997837253305382904924534,
323 1.9997959079539561241378037,
324 1.9998070894618131628688376,
325 1.9998173766146948726477386,
326 1.9998268623119294251254049,
327 1.9998356276964555233591550,
328 1.9998437438960074911519643,
329 1.9998512734706997544817733,
330 1.9998582716222576653730680,
331 1.9998647872087154521263119,
332 1.9998708635995425604743592,
333 1.9998765393992465586080545,
334 1.9998818490620738438532863,
335 1.9998868234161445738743604,
336 1.9998914901119566046557366,
337 1.9998958740074785727074830,
338 1.9998999974998749922178512,
339 1.9999038808121516333851653,
340 1.9999075422415889522788357,
341 1.9999109983756762928906042,
342 1.9999142642803164142926986,
343 1.9999173536642966251576647,
344 1.9999202790233863698352718,
345 1.9999230517668953255793970,
346 1.9999256823290901509881587,
347 1.9999281802675053497939453,
348 1.9999305543498809755115567,
349 1.9999328126312063734840163,
350 1.9999349625221361967906952,
351 1.9999370108498654878071213,
352 1.9999389639123990025908889,
353 1.9999408275270214743832980,
354 1.9999426070736663534153227,
355 1.9999443075337875691619566,
356 1.9999459335252594405609242,
357 1.9999474893337618725376939,
358 1.9999489789410496206473189,
359 1.9999504060504542197682087,
360 1.9999517741099239061773607,
361 1.9999530863328694717815162,
362 1.9999543457170516339753766,
363 1.9999555550617174208343982,
364 1.9999567169831686592367317,
365 1.9999578339289244014115657,
366 1.9999589081906205670141272,
367 1.9999599419157738604399185,
368 1.9999609371185228226443925,
369 1.9999618956894464111197912,
370 1.9999628194045495514557728,
371 1.9999637099334954661120772,
372 1.9999645688471560844769875,
373 1.9999653976245443358231504,
374 1.9999661976591854849498106,
375 1.9999669702649787970689249,
376 1.9999677166815955984928432,
377 1.9999684380794551698127051,
378 1.9999691355643157882884461,
379 1.9999698101815145638701296,
380 1.9999704629198864438586677,
381 1.9999710947153898319165674,
382 1.9999717064544636596331098,
383 1.9999722989771384023096126,
384 1.9999728730799214406138944,
385 1.9999734295184752761939673,
386 1.9999739690101054293697169,
387 1.9999744922360733178625866,
388 1.9999749998437480468264568,
389 1.9999754924486098081387250,
390 1.9999759706361164681223674,
391 1.9999764349634439192238133,
392 1.9999768859611098577311328,
393 1.9999773241344898295265148,
394 1.9999777499652336322380833,
395 1.9999781639125894913602138,
396 1.9999785664146428101331163,
397 1.9999789578894757322230849,
398 1.9999793387362532517655470,
399 1.9999797093362411396682068,
400 1.9999800700537605352685028,
401 1.9999804212370836657046771,
402 1.9999807632192748093908491,
403 1.9999810963189802913647988,
404 1.9999814208411710151919866,
405 1.9999817370778407593125137,
406 1.9999820453086632256753732,
407 1.9999823458016106016892424,
408 1.9999826388135361899785550,
409 1.9999829245907234727581969,
410 1.9999832033694038013488942,
411 1.9999834753762447451218287,
412 1.9999837408288109826979759,
413 1.9999839999359994880149746,
414 1.9999842528984506348542841,
415 1.9999844999089367322906569,
416 1.9999847411527293949973227,
417 1.9999849768079470566052233,
418 1.9999852070458838431332368,
419 1.9999854320313209394806633,
420 1.9999856519228215081390720,
421 1.9999858668730101443622438,
422 1.9999860770288377853545064,
423 1.9999862825318329381286964,
424 1.9999864835183400217297200,
425 1.9999866801197455758263427,
426 1.9999868724626930385594725,
427 1.9999870606692867456861268,
428 1.9999872448572857672795949,
429 1.9999874251402881543361253,
430 1.9999876016279061367387007,
431 1.9999877744259327734793058,
432 1.9999879436365005309960211,
433 1.9999881093582322320878486,
434 1.9999882716863847930419476,
435 1.9999884307129861387439598,
436 1.9999885865269656637496429,
437 1.9999887392142785871298680,
438 1.9999888888580245199528412,
439 1.9999890355385605579790348,
440 1.9999891793336091811361319,
441 1.9999893203183612335350344,
442 1.9999894585655742371881408,
443 1.9999895941456662814238177,
444 1.9999897271268057108006139,
445 1.9999898575749968286869099,
446 1.9999899855541618157823633,
447 1.9999901111262190523407481,
448 1.9999902343511580255896315,
449 1.9999903552871109908319058,
450 1.9999904739904215434384907,
451 1.9999905905157102581825787,
452 1.9999907049159375333922603,
453 1.9999908172424637773983633,
454 1.9999909275451070651049429,
455 1.9999910358721983849204437,
456 1.9999911422706345906497027,
457 1.9999912467859291663333302,
458 1.9999913494622609072505146,
459 1.9999914503425206137960865,
460 1.9999915494683558901721845,
461 1.9999916468802141375580433,
462 1.9999917426173838194951998,
463 1.9999918367180340843811473,
464 1.9999919292192528146772171,
465 1.9999920201570831760143343,
466 1.9999921095665587343929645,
467 1.9999921974817372009999503,
468 1.9999922839357328720796131,
469 1.9999923689607478141661012,
470 1.9999924525881018556091456,
471 1.9999925348482614335075824,
472 1.9999926157708673449481598,
473 1.9999926953847614525313497,
474 1.9999927737180123875522508,
475 1.9999928507979402920362658,
476 1.9999929266511406438650011,
477 1.9999930013035072003373810,
478 1.9999930747802541001730339,
479 1.9999931471059371559419182,
480 1.9999932183044743736746382,
481 1.9999932883991657291437513,
482 1.9999933574127122330168702,
483 1.9999934253672343136129172,
484 1.9999934922842895436076449,
485 1.9999935581848897396366194,
486 1.9999936230895174589733757,
487 1.9999936870181419158341504,
488 1.9999937499902343445226660,
489 };
490
491 static const double legcseqnorm[201] = {
492 0.0000000000000000000000000,
493 0.0000000000000000000000000,
494 -1.1180339887498948482072100,
495 -1.0183501544346311125984611,
496 -1.0062305898749053633539630,
497 -1.0028530728448139498158037,
498 -1.0015420209622192481389857,
499 -1.0009272139219581055366928,
500 -1.0006007810695147948465422,
501 -1.0004114379931337589987872,
502 -1.0002940744071803442008517,
503 -1.0002174622185106379413169,
504 -1.0001653302482984141553307,
505 -1.0001286256356210525864103,
506 -1.0001020356106936058139534,
507 -1.0000823011400101476917751,
508 -1.0000673468701305763445264,
509 -1.0000558082429209263942635,
510 -1.0000467628422711159142180,
511 -1.0000395718327849260524676,
512 -1.0000337832131310390878717,
513 -1.0000290710350803458344310,
514 -1.0000251962155324265491343,
515 -1.0000219806789858287190270,
516 -1.0000192899374059475950394,
517 -1.0000170211317362819157528,
518 -1.0000150946813898723700489,
519 -1.0000134483616539496234413,
520 -1.0000120330427357868583946,
521 -1.0000108095837772925361039,
522 -1.0000097465411964244106496,
523 -1.0000088184587981272660392,
524 -1.0000080045786190982050334,
525 -1.0000072878595193167500232,
526 -1.0000066542232692270346719,
527 -1.0000060919704783676505502,
528 -1.0000055913245009481336162,
529 -1.0000051440726137630963036,
530 -1.0000047432817343602117627,
531 -1.0000043830717004995434563,
532 -1.0000040584333230009967883,
533 -1.0000037650815046152063745,
534 -1.0000034993360010045953049,
535 -1.0000032580241061291089627,
536 -1.0000030384008288828010663,
537 -1.0000028380831019170736640,
538 -1.0000026549953073004908810,
539 -1.0000024873239751952141291,
540 -1.0000023334799536602361322,
541 -1.0000021920666914396246330,
542 -1.0000020618535444830643141,
543 -1.0000019417532284092108713,
544 -1.0000018308027062940946880,
545 -1.0000017281469339661820195,
546 -1.0000016330249909989784721,
547 -1.0000015447582105839689076,
548 -1.0000014627399899225625515,
549 -1.0000013864270181358339523,
550 -1.0000013153317036379693836,
551 -1.0000012490156195652532892,
552 -1.0000011870838158338720483,
553 -1.0000011291798710141289983,
554 -1.0000010749815775035198584,
555 -1.0000010241971702458106327,
556 -1.0000009765620231633317591,
557 -1.0000009318357490399825985,
558 -1.0000008897996482624395340,
559 -1.0000008502544599036081924,
560 -1.0000008130183754348615718,
561 -1.0000007779252810548829833,
562 -1.0000007448231994552217961,
563 -1.0000007135729059221974924,
564 -1.0000006840466971318745829,
565 -1.0000006561272939441861046,
566 -1.0000006297068620102598893,
567 -1.0000006046861361493852810,
568 -1.0000005809736362878291024,
569 -1.0000005584849643282480472,
570 -1.0000005371421726703363639,
571 -1.0000005168731962740692049,
572 -1.0000004976113411596806008,
573 -1.0000004792948231127072342,
574 -1.0000004618663511196359830,
575 -1.0000004452727507103229274,
576 -1.0000004294646229590628680,
577 -1.0000004143960353902593347,
578 -1.0000004000242414690848708,
579 -1.0000003863094257375342469,
580 -1.0000003732144719883643755,
581 -1.0000003607047521588966840,
582 -1.0000003487479338864325415,
583 -1.0000003373138048878843165,
584 -1.0000003263741125275609889,
585 -1.0000003159024171116037866,
586 -1.0000003058739576006566652,
587 -1.0000002962655285725437904,
588 -1.0000002870553673827358140,
589 -1.0000002782230505849869057,
590 -1.0000002697493987622364559,
591 -1.0000002616163890110023357,
592 -1.0000002538070743900383910,
593 -1.0000002463055097187303796,
594 -1.0000002390966831671914941,
595 -1.0000002321664531315182148,
596 -1.0000002255014899446962720,
597 -1.0000002190892220043294189,
598 -1.0000002129177859488875363,
599 -1.0000002069759805409503856,
600 -1.0000002012532239508346352,
601 -1.0000001957395141607715799,
602 -1.0000001904253922358238246,
603 -1.0000001853019082297385073,
604 -1.0000001803605895158355218,
605 -1.0000001755934113488585510,
606 -1.0000001709927694851839236,
607 -1.0000001665514546978896088,
608 -1.0000001622626290424854581,
609 -1.0000001581198037385383656,
610 -1.0000001541168185448943406,
611 -1.0000001502478225146562651,
612 -1.0000001465072560282191724,
613 -1.0000001428898340098206168,
614 -1.0000001393905302394604986,
615 -1.0000001360045626810435851,
616 -1.0000001327273797537779215,
617 -1.0000001295546474781991350,
618 -1.0000001264822374353463674,
619 -1.0000001235062154811934418,
620 -1.0000001206228311637514566,
621 -1.0000001178285077943789713,
622 -1.0000001151198331268959291,
623 -1.0000001124935506041689973,
624 -1.0000001099465511333538870,
625 -1.0000001074758653540159800,
626 -1.0000001050786563654105762,
627 -1.0000001027522128845166638,
628 -1.0000001004939428029486687,
629 -1.0000000983013671212789805,
630 -1.0000000961721142318230568,
631 -1.0000000941039145299377849,
632 -1.0000000920945953319322191,
633 -1.0000000901420760810508342,
634 -1.0000000882443638214715567,
635 -1.0000000863995489252481624,
636 -1.0000000846058010559340087,
637 -1.0000000828613653539251108,
638 -1.0000000811645588297531945,
639 -1.0000000795137669539446035,
640 -1.0000000779074404293504316,
641 -1.0000000763440921371658435,
642 -1.0000000748222942452544601,
643 -1.0000000733406754682610490,
644 -1.0000000718979184725736242,
645 -1.0000000704927574156181957,
646 -1.0000000691239756129809549,
647 -1.0000000677904033262021616,
648 -1.0000000664909156620260150,
649 -1.0000000652244305800707413,
650 -1.0000000639899069999200210,
651 -1.0000000627863430035157863,
652 -1.0000000616127741272145396,
653 -1.0000000604682717384114402,
654 -1.0000000593519414925037722,
655 -1.0000000582629218640138419,
656 -1.0000000572003827506786799,
657 -1.0000000561635241443265987,
658 -1.0000000551515748656132565,
659 -1.0000000541637913591477838,
660 -1.0000000531994565462984639,
661 -1.0000000522578787309074810,
662 -1.0000000513383905583484138,
663 -1.0000000504403480213128433,
664 -1.0000000495631295122176546,
665 -1.0000000487061349200646268,
666 -1.0000000478687847678491846,
667 -1.0000000470505193895425292,
668 -1.0000000462507981442619037,
669 -1.0000000454690986663279506,
670 -1.0000000447049161481733945,
671 -1.0000000439577626562114715,
672 -1.0000000432271664758693969,
673 -1.0000000425126714863289720,
674 -1.0000000418138365622638258,
675 -1.0000000411302350023564495,
676 -1.0000000404614539811255780,
677 -1.0000000398070940262323231,
678 -1.0000000391667685166029872,
679 -1.0000000385401032031032834,
680 -1.0000000379267357488366139,
681 -1.0000000373263152893916691,
682 -1.0000000367385020107625221,
683 -1.0000000361629667461338400,
684 -1.0000000355993905882786080,
685 -1.0000000350474645183273051,
686 -1.0000000345068890502580120,
687 -1.0000000339773738891558860,
688 -1.0000000334586376040009459,
689 -1.0000000329504073138999626,
690 -1.0000000324524183859193141,
691 -1.0000000319644141456030054,
692 -1.0000000314861455999590131,
693 };
694
lega_seq_norm(size_t n)695 inline static double lega_seq_norm (size_t n) {
696 if (n <= 200){
697 return (double) legaseqnorm[n];
698 }
699 return ( sqrt( 4 * (double) (n * n) - 1) / (double) n);
700 }
701
legc_seq_norm(size_t n)702 inline static double legc_seq_norm (size_t n) {
703 if (n <= 200){
704 return (double) legcseqnorm[n];
705 }
706 return ( -(sqrt(2 * (double) n + 1) * ((double) n -1)) / (double)n / sqrt( 2 * (double) n - 3));}
707
708
709 // Orthonormality functions
chebortho(size_t n)710 double chebortho(size_t n) {
711 if (n == 0){
712 return M_PI;
713 }
714 else{
715 return M_PI/2.0;
716 }
717 }
718
719 /* static const double legorthoarr[201] = */
720 /* {1.000000000000000e+00, 3.333333333333333e-01,2.000000000000000e-01, */
721 /* 1.428571428571428e-01,1.111111111111111e-01,9.090909090909091e-02,7.692307692307693e-02, */
722 /* 6.666666666666667e-02,5.882352941176471e-02,5.263157894736842e-02,4.761904761904762e-02, */
723 /* 4.347826086956522e-02,4.000000000000000e-02,3.703703703703703e-02,3.448275862068965e-02, */
724 /* 3.225806451612903e-02,3.030303030303030e-02,2.857142857142857e-02,2.702702702702703e-02, */
725 /* 2.564102564102564e-02,2.439024390243903e-02,2.325581395348837e-02,2.222222222222222e-02, */
726 /* 2.127659574468085e-02,2.040816326530612e-02,1.960784313725490e-02,1.886792452830189e-02, */
727 /* 1.818181818181818e-02,1.754385964912281e-02,1.694915254237288e-02,1.639344262295082e-02 */
728 /* ,1.587301587301587e-02,1.538461538461539e-02,1.492537313432836e-02,1.449275362318841e-02 */
729 /* ,1.408450704225352e-02,1.369863013698630e-02,1.333333333333333e-02,1.298701298701299e-02 */
730 /* ,1.265822784810127e-02,1.234567901234568e-02,1.204819277108434e-02,1.176470588235294e-02 */
731 /* ,1.149425287356322e-02,1.123595505617977e-02,1.098901098901099e-02,1.075268817204301e-02 */
732 /* ,1.052631578947368e-02,1.030927835051546e-02,1.010101010101010e-02,9.900990099009901e-03 */
733 /* ,9.708737864077669e-03,9.523809523809525e-03,9.345794392523364e-03,9.174311926605505e-03 */
734 /* ,9.009009009009009e-03,8.849557522123894e-03,8.695652173913044e-03,8.547008547008548e-03 */
735 /* ,8.403361344537815e-03,8.264462809917356e-03,8.130081300813009e-03,8.000000000000000e-03 */
736 /* ,7.874015748031496e-03,7.751937984496124e-03,7.633587786259542e-03,7.518796992481203e-03 */
737 /* ,7.407407407407408e-03,7.299270072992700e-03,7.194244604316547e-03,7.092198581560284e-03 */
738 /* ,6.993006993006993e-03,6.896551724137931e-03,6.802721088435374e-03,6.711409395973154e-03 */
739 /* ,6.622516556291391e-03,6.535947712418301e-03,6.451612903225806e-03,6.369426751592357e-03 */
740 /* ,6.289308176100629e-03,6.211180124223602e-03,6.134969325153374e-03,6.060606060606061e-03 */
741 /* ,5.988023952095809e-03,5.917159763313609e-03,5.847953216374269e-03,5.780346820809248e-03 */
742 /* ,5.714285714285714e-03,5.649717514124294e-03,5.586592178770950e-03,5.524861878453038e-03 */
743 /* ,5.464480874316940e-03,5.405405405405406e-03,5.347593582887700e-03,5.291005291005291e-03 */
744 /* ,5.235602094240838e-03,5.181347150259068e-03,5.128205128205128e-03,5.076142131979695e-03 */
745 /* ,5.025125628140704e-03,4.975124378109453e-03,4.926108374384237e-03,4.878048780487805e-03 */
746 /* ,4.830917874396135e-03,4.784688995215311e-03,4.739336492890996e-03,4.694835680751174e-03 */
747 /* ,4.651162790697674e-03,4.608294930875576e-03,4.566210045662100e-03,4.524886877828055e-03 */
748 /* ,4.484304932735426e-03,4.444444444444444e-03,4.405286343612335e-03,4.366812227074236e-03 */
749 /* ,4.329004329004329e-03,4.291845493562232e-03,4.255319148936170e-03,4.219409282700422e-03 */
750 /* ,4.184100418410041e-03,4.149377593360996e-03,4.115226337448560e-03,4.081632653061225e-03 */
751 /* ,4.048582995951417e-03,4.016064257028112e-03,3.984063745019920e-03,3.952569169960474e-03 */
752 /* ,3.921568627450980e-03,3.891050583657588e-03,3.861003861003861e-03,3.831417624521073e-03 */
753 /* ,3.802281368821293e-03,3.773584905660377e-03,3.745318352059925e-03,3.717472118959108e-03 */
754 /* ,3.690036900369004e-03,3.663003663003663e-03,3.636363636363636e-03,3.610108303249098e-03 */
755 /* ,3.584229390681004e-03,3.558718861209964e-03,3.533568904593640e-03,3.508771929824561e-03 */
756 /* ,3.484320557491289e-03,3.460207612456748e-03,3.436426116838488e-03,3.412969283276451e-03 */
757 /* ,3.389830508474576e-03,3.367003367003367e-03,3.344481605351171e-03,3.322259136212625e-03 */
758 /* ,3.300330033003300e-03,3.278688524590164e-03,3.257328990228013e-03,3.236245954692557e-03 */
759 /* ,3.215434083601286e-03,3.194888178913738e-03,3.174603174603175e-03,3.154574132492113e-03 */
760 /* ,3.134796238244514e-03,3.115264797507788e-03,3.095975232198143e-03,3.076923076923077e-03 */
761 /* ,3.058103975535168e-03,3.039513677811550e-03,3.021148036253776e-03,3.003003003003003e-03 */
762 /* ,2.985074626865672e-03,2.967359050445104e-03,2.949852507374631e-03,2.932551319648094e-03 */
763 /* ,2.915451895043732e-03,2.898550724637681e-03,2.881844380403458e-03,2.865329512893983e-03 */
764 /* ,2.849002849002849e-03,2.832861189801700e-03,2.816901408450704e-03,2.801120448179272e-03 */
765 /* ,2.785515320334262e-03,2.770083102493075e-03,2.754820936639119e-03,2.739726027397260e-03 */
766 /* ,2.724795640326975e-03,2.710027100271003e-03,2.695417789757413e-03,2.680965147453083e-03 */
767 /* ,2.666666666666667e-03,2.652519893899204e-03,2.638522427440633e-03,2.624671916010499e-03 */
768 /* ,2.610966057441253e-03,2.597402597402597e-03,2.583979328165375e-03,2.570694087403599e-03 */
769 /* ,2.557544757033248e-03,2.544529262086514e-03,2.531645569620253e-03,2.518891687657431e-03 */
770 /* ,2.506265664160401e-03,2.493765586034913e-03}; */
771
772 // Helper functions
orth_poly_expansion_eval2(double x,void * p)773 double orth_poly_expansion_eval2(double x, void * p){
774
775 struct OrthPolyExpansion * temp = p;
776 return orth_poly_expansion_eval(temp,x);
777 }
778
orth_poly_expansion_eval3(double x,void * p)779 double orth_poly_expansion_eval3(double x, void * p)
780 {
781 struct OrthPolyExpansion ** temp = p;
782 double out = orth_poly_expansion_eval(temp[0],x);
783 out *= orth_poly_expansion_eval(temp[1],x);
784
785 return out;
786 }
787
788 struct lin_func
789 {
790 double slope;
791 double offset;
792 };
793
eval_lin_func(double x,void * args)794 double eval_lin_func(double x, void * args)
795 {
796 struct lin_func * lf = args;
797 double m = lf->slope;
798 double b = lf->offset;
799 return m*x + b;
800 }
801
802 struct quad_func
803 {
804 double scale;
805 double offset;
806 };
807
eval_quad_func(double x,void * args)808 double eval_quad_func(double x, void * args)
809 {
810 struct quad_func * qf = args;
811 double m = qf->scale;
812 double b = qf->offset;
813 return m*(x-b)*(x-b);
814 }
815
816 struct OpeOpts{
817 enum poly_type ptype;
818
819 size_t start_num;
820 size_t max_num;
821 size_t coeffs_check;
822 double tol;
823
824 double lb;
825 double ub;
826
827 // for hermite
828 double mean;
829 double std;
830
831
832 // kristoffel weighting for least squares
833 int kristoffel_eval;
834 enum quad_rule qrule;
835
836 };
837
ope_opts_alloc(enum poly_type ptype)838 struct OpeOpts * ope_opts_alloc(enum poly_type ptype)
839 {
840 struct OpeOpts * ao;
841 if ( NULL == (ao = malloc(sizeof(struct OpeOpts)))){
842 fprintf(stderr, "failed to allocate memory for struct OpeOpts in ope_opts_alloc.\n");
843 exit(1);
844 }
845
846 ao->start_num = 5;
847 ao->coeffs_check = 2;
848 ao->tol = 1e-10;
849 ao->max_num = 100;
850
851 ao->ptype = ptype;
852 if (ptype == HERMITE){
853 ao->lb = -DBL_MAX;
854 ao->ub = DBL_MAX;
855 }
856 else if (ptype == FOURIER){
857 ao->lb = 0.0;
858 ao->ub = 2.0*M_PI;
859 ao->start_num = 8;
860 }
861 else{
862 ao->lb = -1.0;
863 ao->ub = 1.0;
864 }
865
866
867 // for hermite
868 ao->mean = 0.0;
869 ao->std = 1.0;
870
871
872 // for least squares applications
873 ao->kristoffel_eval = 0;
874 ao->qrule = C3_GAUSS_QUAD;
875
876 return ao;
877 }
878
ope_opts_free(struct OpeOpts * ope)879 void ope_opts_free(struct OpeOpts * ope)
880 {
881 if (ope != NULL){
882 free(ope); ope = NULL;
883 }
884 }
885
ope_opts_free_deep(struct OpeOpts ** ope)886 void ope_opts_free_deep(struct OpeOpts ** ope)
887 {
888 if (*ope != NULL){
889 free(*ope); *ope = NULL;
890 }
891 }
892
ope_opts_set_start(struct OpeOpts * ope,size_t start)893 void ope_opts_set_start(struct OpeOpts * ope, size_t start)
894 {
895 assert (ope != NULL);
896 ope->start_num = start;
897 }
898
ope_opts_get_start(const struct OpeOpts * ope)899 size_t ope_opts_get_start(const struct OpeOpts * ope)
900 {
901 assert (ope != NULL);
902 return ope->start_num;
903 }
904
ope_opts_set_maxnum(struct OpeOpts * ope,size_t maxnum)905 void ope_opts_set_maxnum(struct OpeOpts * ope, size_t maxnum)
906 {
907 assert (ope != NULL);
908 ope->max_num = maxnum;
909 }
910
ope_opts_get_maxnum(const struct OpeOpts * ope)911 size_t ope_opts_get_maxnum(const struct OpeOpts * ope)
912 {
913 assert (ope != NULL);
914 assert (ope->max_num < 1000); // really a check if it exists
915
916 return ope->max_num;
917 }
918
ope_opts_set_coeffs_check(struct OpeOpts * ope,size_t num)919 void ope_opts_set_coeffs_check(struct OpeOpts * ope, size_t num)
920 {
921 assert (ope != NULL);
922 ope->coeffs_check = num;
923 }
924
ope_opts_set_tol(struct OpeOpts * ope,double tol)925 void ope_opts_set_tol(struct OpeOpts * ope, double tol)
926 {
927 assert (ope != NULL);
928 ope->tol = tol;
929 }
930
ope_opts_set_lb(struct OpeOpts * ope,double lb)931 void ope_opts_set_lb(struct OpeOpts * ope, double lb)
932 {
933 assert (ope != NULL);
934 ope->lb = lb;
935 }
936
ope_opts_get_lb(const struct OpeOpts * ope)937 double ope_opts_get_lb(const struct OpeOpts * ope)
938 {
939 assert (ope != NULL);
940 return ope->lb;
941 }
942
ope_opts_set_ub(struct OpeOpts * ope,double ub)943 void ope_opts_set_ub(struct OpeOpts * ope, double ub)
944 {
945 assert (ope != NULL);
946 ope->ub = ub;
947 }
948
ope_opts_get_ub(const struct OpeOpts * ope)949 double ope_opts_get_ub(const struct OpeOpts * ope)
950 {
951 assert (ope != NULL);
952 return ope->ub;
953 }
954
ope_opts_set_mean_and_std(struct OpeOpts * ope,double mean,double std)955 void ope_opts_set_mean_and_std(struct OpeOpts * ope, double mean, double std)
956 {
957 if (ope->ptype == HERMITE){
958 ope->mean = mean;
959 ope->std = std;
960 }
961 else{
962 fprintf(stderr,"Warning: setting mean and variance for non hermite polynomials\n");
963 fprintf(stderr," does not do anything\n");
964 }
965
966 }
967
ope_opts_set_ptype(struct OpeOpts * ope,enum poly_type ptype)968 void ope_opts_set_ptype(struct OpeOpts * ope, enum poly_type ptype)
969 {
970 assert (ope != NULL);
971 ope->ptype = ptype;
972 }
973
ope_opts_get_ptype(const struct OpeOpts * ope)974 enum poly_type ope_opts_get_ptype(const struct OpeOpts * ope)
975 {
976 assert (ope != NULL);
977 return ope->ptype;
978 }
979
ope_opts_set_qrule(struct OpeOpts * ope,enum quad_rule qrule)980 void ope_opts_set_qrule(struct OpeOpts * ope, enum quad_rule qrule)
981 {
982 assert (ope != NULL);
983 if ((qrule != C3_GAUSS_QUAD) && (qrule != C3_CC_QUAD)){
984 fprintf(stderr,"Specified qrule %d is not known\n",qrule);
985 exit(1);
986 }
987 ope->qrule = qrule;
988 }
989
990 /********************************************************//**
991 * Get number of free parameters
992 *************************************************************/
ope_opts_get_nparams(const struct OpeOpts * opts)993 size_t ope_opts_get_nparams(const struct OpeOpts * opts)
994 {
995 assert (opts != NULL);
996 return opts->start_num;
997 }
998
999 /********************************************************//**
1000 * Set number of free parameters
1001 *************************************************************/
ope_opts_set_nparams(struct OpeOpts * opts,size_t num)1002 void ope_opts_set_nparams(struct OpeOpts * opts, size_t num)
1003 {
1004 assert (opts != NULL);
1005 opts->start_num = num;
1006 }
1007
1008
1009 /********************************************************//**
1010 * Set kristoffel weighting for evaluation of polynomials
1011 * This is typically only used in the context of least squares
1012 *
1013 * \param[in] opts - options to modify
1014 * \param[in] kristoffel_weight - 1 to use, 0 to not use
1015 *************************************************************/
ope_opts_set_kristoffel_weight(struct OpeOpts * opts,int kristoffel_weight)1016 void ope_opts_set_kristoffel_weight(struct OpeOpts * opts, int kristoffel_weight)
1017 {
1018 assert (opts != NULL);
1019 opts->kristoffel_eval = kristoffel_weight;
1020 }
1021
1022 /********************************************************//**
1023 * Initialize a standard basis polynomial
1024 *
1025 * \param[in] num_poly - number of basis
1026 * \param[in] lb - lower bound
1027 * \param[in] ub - upper bound
1028 *
1029 * \return standard polynomial
1030 *************************************************************/
1031 struct StandardPoly *
standard_poly_init(size_t num_poly,double lb,double ub)1032 standard_poly_init(size_t num_poly, double lb, double ub){
1033
1034 struct StandardPoly * p;
1035 if ( NULL == (p = malloc(sizeof(struct StandardPoly)))){
1036 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1037 exit(1);
1038 }
1039 p->ptype = STANDARD;
1040 p->num_poly = num_poly;
1041 p->lower_bound = lb;
1042 p->upper_bound = ub;
1043 p->coeff = calloc_double(num_poly);
1044 return p;
1045 }
1046
1047 /********************************************************//**
1048 * Create the polynomial representing the derivative
1049 * of the standard polynomial
1050 *
1051 * \param[in] p - polynomial
1052 *
1053 * \return derivative polynomial
1054 *************************************************************/
1055 struct StandardPoly *
standard_poly_deriv(struct StandardPoly * p)1056 standard_poly_deriv(struct StandardPoly * p){
1057
1058 struct StandardPoly * dp;
1059 if ( NULL == (dp = malloc(sizeof(struct StandardPoly)))){
1060 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1061 exit(1);
1062 }
1063 dp->ptype = STANDARD;
1064 dp->lower_bound = p->lower_bound;
1065 dp->upper_bound = p->upper_bound;
1066 if (p->num_poly > 1){
1067 dp->num_poly = p->num_poly-1;
1068 dp->coeff = calloc_double(dp->num_poly);
1069 size_t ii;
1070 for (ii = 1; ii < p->num_poly; ii++){
1071 dp->coeff[ii-1] = p->coeff[ii] * (double) ii;
1072 }
1073 }
1074 else{
1075 dp->num_poly = 1;
1076 dp->coeff = calloc_double(dp->num_poly);
1077 }
1078
1079 return dp;
1080 }
1081
1082 /********************************************************//**
1083 * free memory allocated to a standard polynomial
1084 *
1085 * \param[in,out] p - polynomial structure
1086 *
1087 *************************************************************/
1088 void
standard_poly_free(struct StandardPoly * p)1089 standard_poly_free(struct StandardPoly * p)
1090 {
1091 free(p->coeff);
1092 free(p);
1093 }
1094
1095 /********************************************************//**
1096 * Initialize a Chebyshev polynomial
1097 *
1098 * \return p - polynomial
1099 *************************************************************/
init_cheb_poly()1100 struct OrthPoly * init_cheb_poly(){
1101
1102 struct OrthPoly * p;
1103 if ( NULL == (p = malloc(sizeof(struct OrthPoly)))){
1104 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1105 exit(1);
1106 }
1107 p->ptype = CHEBYSHEV;
1108 p->an = &two_seq;
1109 p->bn = &zero_seq;
1110 p->cn = &none_seq;
1111
1112 p->lower = -1.0;
1113 p->upper = 1.0;
1114
1115 p->const_term = 1.0;
1116 p->lin_coeff = 1.0;
1117 p->lin_const = 0.0;
1118
1119 p->norm = &chebortho;
1120
1121 return p;
1122 }
1123
legortho(size_t n)1124 inline static double legortho(size_t n){
1125 (void) n;
1126 return 1;
1127 /* if (n < 201){ */
1128 /* return legorthoarr[n]; */
1129 /* } */
1130 /* else{ */
1131 /* // printf("here?! n=%zu\n",n); */
1132 /* return (1.0 / (2.0 * (double) n + 1.0)); */
1133 /* } */
1134 }
1135
1136 /********************************************************//**
1137 * Initialize a Legendre polynomial
1138 *
1139 * \return p - polynomial
1140 *************************************************************/
init_leg_poly()1141 struct OrthPoly * init_leg_poly(){
1142
1143 struct OrthPoly * p;
1144 if ( NULL == (p = malloc(sizeof(struct OrthPoly)))){
1145 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1146 exit(1);
1147 }
1148 p->ptype = LEGENDRE;
1149 p->bn = &zero_seq;
1150
1151
1152 /* orthogonal */
1153 /* p->an = &lega_seq; */
1154 /* p->cn = &legc_seq; */
1155 /* p->lin_coeff = 1.0; */
1156
1157 /* orthonormal */
1158 p->an = &lega_seq_norm;
1159 p->cn = &legc_seq_norm;
1160 p->lin_coeff = 1.0 * sqrt(3.0);
1161
1162 p->lower = -1.0;
1163 p->upper = 1.0;
1164
1165 p->const_term = 1.0;
1166
1167
1168 p->lin_const = 0.0;
1169
1170 p->norm = &legortho;
1171
1172 return p;
1173 }
1174
1175 /********************************************************//**
1176 * free memory allocated to a polynomial
1177 *
1178 * \param[in,out] p - polynomial structure to free
1179 *************************************************************/
free_orth_poly(struct OrthPoly * p)1180 void free_orth_poly(struct OrthPoly * p)
1181 {
1182 if (p != NULL){
1183 free(p); p = NULL;
1184 }
1185 }
1186
1187 /********************************************************//**
1188 * serialize an orthonormal polynomial
1189 *
1190 * \param p - polynomial structure to serialize
1191 *
1192 * \return ser - serialize polynomial
1193 *
1194 * \note
1195 * This is actually pretty stupid because only need to
1196 * serialize the type. But hey. Its good practice.
1197 *************************************************************/
1198 unsigned char *
serialize_orth_poly(struct OrthPoly * p)1199 serialize_orth_poly(struct OrthPoly * p)
1200 {
1201
1202 /*
1203 char start[]= "type=";
1204 char * ser = NULL;
1205 concat_string_ow(&ser,start);
1206
1207 char temp[3];
1208 snprintf(temp,2,"%d",p->ptype);
1209 concat_string_ow(&ser,temp);
1210 */
1211
1212 unsigned char * ser = malloc(sizeof(int) * sizeof(unsigned char)) ;
1213 serialize_int(ser, p->ptype);
1214 return ser;
1215 }
1216
1217 /********************************************************//**
1218 * deserialize an orthonormal polynomial
1219 *
1220 * \param[in] ser - string to deserialize
1221 *
1222 * \return poly - orthonormal polynomial
1223 *************************************************************/
1224 struct OrthPoly *
deserialize_orth_poly(unsigned char * ser)1225 deserialize_orth_poly(unsigned char * ser)
1226 {
1227
1228 struct OrthPoly * poly = NULL;
1229
1230 int type;
1231 deserialize_int(ser, &type);
1232 if (type == LEGENDRE) {
1233 poly = init_leg_poly();
1234 }
1235 else if (type == HERMITE){
1236 poly = init_hermite_poly();
1237 }
1238 else if (type == CHEBYSHEV){
1239 poly = init_cheb_poly();
1240 }
1241 else if (type == FOURIER){
1242 poly = init_fourier_poly();
1243 }
1244 else{
1245 fprintf(stderr,"Cannot desrialize polynomial of type %d\n",type);
1246 }
1247 return poly;
1248 }
1249
1250 /********************************************************//**
1251 * Convert an orthogonal family polynomial of order *n*
1252 * to a standard_polynomial
1253 *
1254 * \param[in] p - polynomial
1255 * \param[in] n - polynomial order
1256 *
1257 * \return sp - standard polynomial
1258 *************************************************************/
orth_to_standard_poly(struct OrthPoly * p,size_t n)1259 struct StandardPoly * orth_to_standard_poly(struct OrthPoly * p, size_t n)
1260 {
1261 struct StandardPoly * sp = standard_poly_init(n+1,p->lower,p->upper);
1262 size_t ii, jj;
1263 if (n == 0){
1264 sp->coeff[n] = p->const_term;
1265 }
1266 else if (n == 1){
1267 sp->coeff[0] = p->lin_const;
1268 sp->coeff[1] = p->lin_coeff;
1269 }
1270 else{
1271
1272 double * a = calloc_double(n+1); //n-2 poly
1273 a[0] = p->const_term;
1274 double * b = calloc_double(n+1); // n- 1poly
1275 b[0] = p->lin_const;
1276 b[1] = p->lin_coeff;
1277 for (ii = 2; ii < n+1; ii++){
1278 sp->coeff[0] = p->bn(ii) * b[0] + p->cn(ii) * a[0];
1279 for (jj = 1; jj < ii-1; jj++){
1280 sp->coeff[jj] = p->an(ii)*b[jj-1] + p->bn(ii) * b[jj] +
1281 p->cn(ii) * a[jj];
1282 }
1283 sp->coeff[ii-1] = p->an(ii)*b[ii-2] + p->bn(ii) * b[ii-1];
1284 sp->coeff[ii] = p->an(ii) * b[ii-1];
1285
1286 memcpy(a, b, ii * sizeof(double));
1287 memcpy(b, sp->coeff, (ii+1) * sizeof(double));
1288 }
1289
1290 free(a);
1291 free(b);
1292 }
1293
1294 return sp;
1295 }
1296
1297 /********************************************************//**
1298 * Evaluate an orthogonal polynomial with previous two
1299 * polynomial orders specified
1300 *
1301 * \param[in] rec - orthogonal polynomial
1302 * \param[in] p2 - if evaluating polynomial P_n(x), then p2 is P_{n-2}(x)
1303 * \param[in] p1 - if evaluating polynomial P_n(x), then p1 is P_{n-1}(x)
1304 * \param[in] n - order
1305 * \param[in] x - location at which to evaluate
1306 *
1307 * \return out - polynomial value
1308 *************************************************************/
1309 inline static double
eval_orth_poly_wp(const struct OrthPoly * rec,double p2,double p1,size_t n,double x)1310 eval_orth_poly_wp(const struct OrthPoly * rec, double p2, double p1,
1311 size_t n, double x)
1312 {
1313 return (rec->an(n) * x + rec->bn(n)) * p1 + rec->cn(n) * p2;
1314 }
1315
1316 /********************************************************//**
1317 * Evaluate the derivative of a legendre polynomial up to a certain
1318 * order
1319 *
1320 * \param[in] x - location at which to evaluate
1321 * \param[in] order - maximum order;
1322 *
1323 * \return out - derivatives
1324 *************************************************************/
deriv_legen_upto(double x,size_t order)1325 double * deriv_legen_upto(double x, size_t order){
1326
1327 double * out = calloc_double(order+1);
1328 if (order == 0){
1329 return out;
1330 }
1331 else if( fabs(x-1.0) <= DBL_EPSILON) {
1332 size_t ii;
1333 for (ii = 1; ii < order+1; ii++){
1334 out[ii] = (double) (order * (order+1)/2.0);
1335
1336 out[ii] *= sqrt(2*ii+1); // for orthonormal
1337 }
1338 }
1339 else if (fabs(x-(-1.0)) <= DBL_EPSILON){
1340 size_t ii;
1341 for (ii = 1; ii < order+1; ii+=2){
1342 out[ii] = (double) (order * (order+1)/2.0);
1343 out[ii] *= sqrt(2*ii+1); // for orthonormal
1344 }
1345 for (ii = 2; ii < order+1; ii+=2){
1346 out[ii] = -(double) (order * (order+1)/2.0);
1347 out[ii] *= sqrt(2*ii+1); // for orthonormal
1348 }
1349 }
1350 else if (order == 1){
1351 struct OrthPoly * p = init_leg_poly();
1352 out[1] = x * orth_poly_eval(p,order,x) - orth_poly_eval(p,order-1,x);
1353 out[1] = (double) order * out[1] / ( x * x - 1.0);
1354 free_orth_poly(p);
1355 }
1356 else{
1357 struct OrthPoly * p = init_leg_poly();
1358 double eval0 = orth_poly_eval(p,0,x);
1359 double eval1 = orth_poly_eval(p,1,x);
1360 double evaltemp;
1361
1362 out[1] = x * eval1 - eval0;
1363 //printf("out[1]=%G\n",out[1]);
1364 out[1] = 1.0 * out[1] / ( x * x - 1.0);
1365
1366 size_t ii;
1367 for (ii = 2; ii < order+1; ii++){
1368 evaltemp = eval_orth_poly_wp(p, eval0, eval1, ii, x);
1369 eval0 = eval1;
1370 eval1 = evaltemp;
1371 out[ii] = x * eval1 - eval0;
1372 out[ii] = (double) ii * out[ii] / ( x * x - 1.0);
1373 }
1374 free_orth_poly(p);
1375 }
1376
1377 return out;
1378 }
1379
1380 /********************************************************//**
1381 * Evaluate the derivative of a legendre polynomial of a certain order
1382 *
1383 * \param[in] x - location at which to evaluate
1384 * \param[in] order - order of the polynomial
1385 *
1386 * \return out - derivative
1387 *************************************************************/
deriv_legen(double x,size_t order)1388 double deriv_legen(double x, size_t order){
1389
1390 if (order == 0){
1391 return 0.0;
1392 }
1393 double out;
1394
1395 if ( fabs(x-1.0) <= DBL_EPSILON) {
1396 out = (double) (order * (order+1)/2.0);
1397 }
1398 else if (fabs(x-(-1.0)) <= DBL_EPSILON){
1399 if (order % 2){ // odd
1400 out = (double) (order * (order+1)/2.0);
1401 }
1402 else{
1403 out = -(double) (order * (order+1)/2.0);
1404 }
1405 }
1406 else{
1407 struct OrthPoly * p = init_leg_poly();
1408 out = x * orth_poly_eval(p,order,x) - orth_poly_eval(p,order-1,x);
1409 //printf("out in plain = %G\n",out);
1410 out = (double) order * out / ( x * x - 1.0);
1411 free_orth_poly(p);
1412 }
1413 return out;
1414 }
1415
1416 /********************************************************//**
1417 * Evaluate the derivative of an orthogonal polynomial
1418 *
1419 * \param[in] ptype - polynomial type
1420 * \param[in] order - order of the polynomial
1421 * \param[in] x - location at which to evaluate
1422 *
1423 * \return out - orthonormal polynomial expansion
1424 *************************************************************/
orth_poly_deriv(enum poly_type ptype,size_t order,double x)1425 double orth_poly_deriv(enum poly_type ptype, size_t order, double x)
1426 {
1427 assert (1 == 0);
1428 double out = 0.0;
1429 if (ptype == LEGENDRE){
1430 out = deriv_legen(x,order);
1431 }
1432 else {
1433 fprintf(stderr,"Have not implemented orth_poly_deriv for %d\n",ptype);
1434 exit(1);
1435 }
1436 return out;
1437 }
1438
1439 /********************************************************//**
1440 * Evaluate an orthogonal polynomial of a given order
1441 *
1442 * \param[in] rec - orthogonal polynomial
1443 * \param[in] n - order
1444 * \param[in] x - location at which to evaluate
1445 *
1446 * \return out - polynomial value
1447 *************************************************************/
1448 double
orth_poly_eval(const struct OrthPoly * rec,size_t n,double x)1449 orth_poly_eval(const struct OrthPoly * rec, size_t n, double x)
1450 {
1451 if (n == 0){
1452 return rec->const_term;
1453 }
1454 else if (n == 1){
1455 return rec->lin_coeff * x + rec->lin_const;
1456 }
1457 else {
1458 double out = (rec->an(n)*x + rec->bn(n)) * orth_poly_eval(rec,n-1,x) +
1459 rec->cn(n) * orth_poly_eval(rec,n-2,x);
1460 return out;
1461 }
1462 }
1463
1464 /********************************************************//**
1465 * Get number of polynomials
1466 *************************************************************/
orth_poly_expansion_get_num_poly(const struct OrthPolyExpansion * ope)1467 size_t orth_poly_expansion_get_num_poly(const struct OrthPolyExpansion * ope)
1468 {
1469 assert (ope != NULL);
1470 return ope->num_poly;
1471 }
1472
1473 /********************************************************//**
1474 * Get number of parameters
1475 *************************************************************/
orth_poly_expansion_get_num_params(const struct OrthPolyExpansion * ope)1476 size_t orth_poly_expansion_get_num_params(const struct OrthPolyExpansion * ope)
1477 {
1478 assert (ope != NULL);
1479 return ope->num_poly;
1480 }
1481
1482 /********************************************************//**
1483 * Get lower bounds
1484 *************************************************************/
orth_poly_expansion_get_lb(const struct OrthPolyExpansion * ope)1485 double orth_poly_expansion_get_lb(const struct OrthPolyExpansion * ope)
1486 {
1487 assert (ope != NULL);
1488 return ope->lower_bound;
1489 }
1490
1491 /********************************************************//**
1492 * Get upper bounds
1493 *************************************************************/
orth_poly_expansion_get_ub(const struct OrthPolyExpansion * ope)1494 double orth_poly_expansion_get_ub(const struct OrthPolyExpansion * ope)
1495 {
1496 assert (ope != NULL);
1497 return ope->upper_bound;
1498 }
1499
1500
1501 /********************************************************//**
1502 * Initialize an expansion of a certain orthogonal polynomial family
1503 *
1504 * \param[in] ptype - type of polynomial
1505 * \param[in] num_poly - number of polynomials
1506 * \param[in] lb - lower bound
1507 * \param[in] ub - upper bound
1508 *
1509 * \return p - orthogonal polynomial expansion
1510 *************************************************************/
1511 struct OrthPolyExpansion *
orth_poly_expansion_init(enum poly_type ptype,size_t num_poly,double lb,double ub)1512 orth_poly_expansion_init(enum poly_type ptype, size_t num_poly,
1513 double lb, double ub)
1514 {
1515
1516 struct OrthPolyExpansion * p;
1517 if ( NULL == (p = malloc(sizeof(struct OrthPolyExpansion)))){
1518 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1519 exit(1);
1520 }
1521
1522 p->num_poly = num_poly;
1523
1524 //p->coeff = calloc_double(num_poly);
1525 p->lower_bound = lb;
1526 p->upper_bound = ub;
1527 p->kristoffel_eval = 0;
1528
1529 double m, off;
1530 switch (ptype) {
1531 case LEGENDRE:
1532 p->p = init_leg_poly();
1533 p->space_transform = space_mapping_create(SM_LINEAR);
1534 m = (p->p->upper - p->p->lower) / (ub - lb);
1535 off = p->p->upper - m * ub;
1536 p->space_transform->set = 1;
1537 p->space_transform->lin_slope = m;
1538 p->space_transform->lin_offset = off;
1539 p->space_transform->inv_lin_slope = 1.0/m;
1540 p->space_transform->inv_lin_offset = -off/m;
1541 break;
1542 case CHEBYSHEV:
1543 p->p = init_cheb_poly();
1544 p->space_transform = space_mapping_create(SM_LINEAR);
1545 m = (p->p->upper - p->p->lower) / (ub - lb);
1546 off = p->p->upper - m * ub;
1547 p->space_transform->set = 1;
1548 p->space_transform->lin_slope = m;
1549 p->space_transform->lin_offset = off;
1550 p->space_transform->inv_lin_slope = 1.0/m;
1551 p->space_transform->inv_lin_offset = -off/m;
1552 break;
1553 case FOURIER:
1554 p->p = init_fourier_poly();
1555 p->space_transform = space_mapping_create(SM_LINEAR);
1556 m = (p->p->upper - p->p->lower) / (ub - lb);
1557 off = p->p->upper - m * ub;
1558 p->space_transform->set = 1;
1559 p->space_transform->lin_slope = m;
1560 p->space_transform->lin_offset = off;
1561 p->space_transform->inv_lin_slope = 1.0/m;
1562 p->space_transform->inv_lin_offset = -off/m;
1563 break;
1564 case HERMITE:
1565 p->p = init_hermite_poly();
1566 p->space_transform = space_mapping_create(SM_LINEAR);;
1567 break;
1568 case STANDARD:
1569 p->space_transform = space_mapping_create(SM_LINEAR);
1570 break;
1571 //default:
1572 // fprintf(stderr, "Polynomial type does not exist: %d\n ", ptype);
1573 }
1574
1575 p->nalloc = num_poly+OPECALLOC;
1576 p->coeff = calloc_double(p->nalloc);
1577 p->ccoeff = NULL;
1578 if (ptype == FOURIER){
1579 p->ccoeff = malloc(p->nalloc * sizeof(double complex));
1580 for (size_t ii = 0; ii < p->nalloc; ii++){
1581 p->ccoeff[ii] = 0.0;
1582 }
1583 }
1584
1585 return p;
1586 }
1587
1588 /********************************************************//**
1589 * Initialize an expansion of a certain orthogonal polynomial family
1590 *
1591 * \param[in] opts - approximation options
1592 * \param[in] num_poly - number of polynomials
1593 *
1594 * \return p - orthogonal polynomial expansion
1595 *************************************************************/
1596 struct OrthPolyExpansion *
orth_poly_expansion_init_from_opts(const struct OpeOpts * opts,size_t num_poly)1597 orth_poly_expansion_init_from_opts(const struct OpeOpts * opts, size_t num_poly)
1598 {
1599
1600 struct OrthPolyExpansion * p =
1601 orth_poly_expansion_init(opts->ptype,num_poly,opts->lb,opts->ub);
1602 if (opts->ptype == HERMITE)
1603 {
1604 p->space_transform->set = 1;
1605 p->space_transform->lin_slope = 1.0/opts->std;
1606 p->space_transform->lin_offset = -opts->mean/opts->std;
1607 p->space_transform->inv_lin_slope = opts->std;
1608 p->space_transform->inv_lin_offset = opts->mean;
1609 }
1610
1611 p->kristoffel_eval = opts->kristoffel_eval;
1612
1613 return p;
1614 }
1615
1616 /********************************************************//**
1617 * Initialize an expanion of a certain orthogonal polynomial family
1618 *
1619 * \param[in] opts - approximation options
1620 * \param[in] nparams - number of polynomials
1621 * \param[in] param - parameters
1622 *
1623 * \return p - orthogonal polynomial expansion
1624 *************************************************************/
1625 struct OrthPolyExpansion *
orth_poly_expansion_create_with_params(struct OpeOpts * opts,size_t nparams,const double * param)1626 orth_poly_expansion_create_with_params(struct OpeOpts * opts,
1627 size_t nparams, const double * param)
1628 {
1629
1630 struct OrthPolyExpansion * poly = NULL;
1631 poly = orth_poly_expansion_init_from_opts(opts,nparams);
1632 for (size_t ii = 0; ii < nparams; ii++){
1633 poly->coeff[ii] = param[ii];
1634 }
1635 return poly;
1636 }
1637
1638 /********************************************************//**
1639 * Get parameters defining polynomial (for now just coefficients)
1640 *************************************************************/
orth_poly_expansion_get_params(const struct OrthPolyExpansion * ope,double * param)1641 size_t orth_poly_expansion_get_params(const struct OrthPolyExpansion * ope, double * param)
1642 {
1643 assert (ope != NULL);
1644 memmove(param,ope->coeff,ope->num_poly * sizeof(double));
1645 return ope->num_poly;
1646 }
1647
1648 /********************************************************//**
1649 * Get parameters defining polynomial (for now just coefficients)
1650 *************************************************************/
orth_poly_expansion_get_params_ref(const struct OrthPolyExpansion * ope,size_t * nparam)1651 double * orth_poly_expansion_get_params_ref(
1652 const struct OrthPolyExpansion * ope, size_t *nparam)
1653 {
1654 assert (ope != NULL);
1655 *nparam = ope->num_poly;
1656 return ope->coeff;
1657 }
1658
1659 /********************************************************//**
1660 * Update an expansion's parameters
1661 *
1662 * \param[in] ope - expansion to update
1663 * \param[in] nparams - number of polynomials
1664 * \param[in] param - parameters
1665
1666 * \returns 0 if successful
1667 *************************************************************/
1668 int
orth_poly_expansion_update_params(struct OrthPolyExpansion * ope,size_t nparams,const double * param)1669 orth_poly_expansion_update_params(struct OrthPolyExpansion * ope,
1670 size_t nparams, const double * param)
1671 {
1672
1673 size_t nold = ope->num_poly;
1674 ope->num_poly = nparams;
1675 if (nold >= nparams){
1676 for (size_t ii = 0; ii < nparams; ii++){
1677 ope->coeff[ii] = param[ii];
1678 }
1679 for (size_t ii = nparams; ii < nold; ii++){
1680 ope->coeff[ii] = 0.0;
1681 }
1682 }
1683 else{
1684 if (nparams <= ope->nalloc){
1685
1686 for (size_t ii = 0; ii < nparams; ii++){
1687 ope->coeff[ii] = param[ii];
1688 }
1689 }
1690 else{
1691 free(ope->coeff); ope->coeff = NULL;
1692 ope->nalloc = nparams;
1693 ope->coeff = calloc_double(ope->nalloc);
1694 for (size_t ii = 0; ii < nparams; ii++){
1695 ope->coeff[ii] = param[ii];
1696 }
1697 }
1698 }
1699 return 0;
1700 }
1701
1702 /********************************************************//**
1703 * Copy an orthogonal polynomial expansion
1704 *
1705 * \param[in] pin - polynomial to copy
1706 *
1707 * \return p - orthogonal polynomial expansion
1708 *************************************************************/
orth_poly_expansion_copy(const struct OrthPolyExpansion * pin)1709 struct OrthPolyExpansion * orth_poly_expansion_copy(const struct OrthPolyExpansion * pin)
1710 {
1711 struct OrthPolyExpansion * p = NULL;
1712 if (pin != NULL){
1713 assert (pin->nalloc >= pin->num_poly);
1714
1715 if ( NULL == (p = malloc(sizeof(struct OrthPolyExpansion)))){
1716 fprintf(stderr, "failed to allocate memory for poly exp.\n");
1717 exit(1);
1718 }
1719 // printf("copying polynomial\n");
1720 //printf("pin->num_poly = %zu, pin->nalloc = %zu\n",pin->num_poly,pin->nalloc);
1721 p->num_poly = pin->num_poly;
1722 p->nalloc = pin->nalloc;
1723 p->coeff = NULL;
1724 p->ccoeff = NULL;
1725 p->lower_bound = pin->lower_bound;
1726 p->upper_bound = pin->upper_bound;
1727 p->space_transform = space_mapping_copy(pin->space_transform);
1728 p->kristoffel_eval = pin->kristoffel_eval;
1729
1730 switch (pin->p->ptype) {
1731 case LEGENDRE:
1732 p->coeff = calloc_double(pin->nalloc);
1733 memmove(p->coeff,pin->coeff, p->num_poly * sizeof(double));
1734 p->p = init_leg_poly();
1735 break;
1736 case CHEBYSHEV:
1737 p->coeff = calloc_double(pin->nalloc);
1738 memmove(p->coeff,pin->coeff, p->num_poly * sizeof(double));
1739 p->p = init_cheb_poly();
1740 break;
1741 case HERMITE:
1742 p->coeff = calloc_double(pin->nalloc);
1743 memmove(p->coeff,pin->coeff, p->num_poly * sizeof(double));
1744 p->p = init_hermite_poly();
1745 break;
1746 case FOURIER:
1747 p->p = init_fourier_poly();
1748 p->ccoeff = malloc(p->nalloc * sizeof(double complex));
1749 for (size_t ii = 0; ii < p->nalloc; ii++){
1750 p->ccoeff[ii] = 0.0;
1751 }
1752 for (size_t ii = 0; ii < p->num_poly; ii++){
1753 p->ccoeff[ii] = pin->ccoeff[ii];
1754 }
1755 break;
1756 case STANDARD:
1757 break;
1758 //default:
1759 // fprintf(stderr, "Polynomial type does not exist: %d\n ", ptype);
1760 }
1761
1762
1763 }
1764 return p;
1765 }
1766
1767 enum poly_type
orth_poly_expansion_get_ptype(const struct OrthPolyExpansion * ope)1768 orth_poly_expansion_get_ptype(const struct OrthPolyExpansion * ope)
1769 {
1770 assert (ope != NULL);
1771 return ope->p->ptype;
1772 }
1773
1774 /********************************************************//**
1775 Return a zero function
1776
1777 \param[in] opts - extra arguments depending on function_class, sub_type, etc.
1778 \param[in] force_nparam - if == 1 then approximation will have the number of parameters
1779 defined by *get_nparams, for each approximation type
1780 if == 0 then it may be more compressed
1781
1782 \return p - zero function
1783 ************************************************************/
1784 struct OrthPolyExpansion *
orth_poly_expansion_zero(struct OpeOpts * opts,int force_nparam)1785 orth_poly_expansion_zero(struct OpeOpts * opts, int force_nparam)
1786 {
1787 assert (opts != NULL);
1788 struct OrthPolyExpansion * p = NULL;
1789
1790 if (opts->ptype == FOURIER){
1791 p = orth_poly_expansion_init_from_opts(opts, ope_opts_get_start(opts));
1792 return p;
1793 }
1794
1795
1796 if (force_nparam == 0){
1797 p = orth_poly_expansion_init_from_opts(opts,1);
1798 p->coeff[0] = 0.0;
1799 }
1800 else{
1801 size_t nparams = ope_opts_get_nparams(opts);
1802 p = orth_poly_expansion_init_from_opts(opts,nparams);
1803 for (size_t ii = 0; ii < nparams; ii++){
1804 p->coeff[ii] = 0.0;
1805 }
1806 }
1807
1808 return p;
1809 }
1810
1811
1812 /********************************************************//**
1813 * Generate a constant orthonormal polynomial expansion
1814 *
1815 * \param[in] a - value of the function
1816 * \param[in] opts - opts
1817 *
1818 * \return p - orthogonal polynomial expansion
1819 *************************************************************/
1820 struct OrthPolyExpansion *
orth_poly_expansion_constant(double a,struct OpeOpts * opts)1821 orth_poly_expansion_constant(double a, struct OpeOpts * opts)
1822 {
1823 assert (isnan(a) == 0);
1824 assert (isinf(a) == 0);
1825 assert (fabs(a) < 1e100);
1826 if (fabs(a) < ZEROTHRESH){
1827 return orth_poly_expansion_zero(opts, 0);
1828 }
1829
1830 assert (opts != NULL);
1831
1832
1833 struct OrthPolyExpansion * p = NULL;
1834 if (opts->ptype != FOURIER){
1835 p = orth_poly_expansion_init_from_opts(opts,1);
1836 p->coeff[0] = a / p->p->const_term;
1837 }
1838 else{
1839 p = orth_poly_expansion_init_from_opts(opts,ope_opts_get_start(opts));
1840 p->ccoeff[0] = a;
1841 }
1842
1843
1844
1845 return p;
1846 }
1847
1848 /********************************************************//**
1849 * Generate a linear orthonormal polynomial expansion
1850 *
1851 * \param[in] a - value of the slope function
1852 * \param[in] offset - offset
1853 * \param[in] opts - options
1854 *
1855 * \return p - orthogonal polynomial expansion
1856 *************************************************************/
1857 struct OrthPolyExpansion *
orth_poly_expansion_linear(double a,double offset,struct OpeOpts * opts)1858 orth_poly_expansion_linear(double a, double offset, struct OpeOpts * opts)
1859 {
1860 assert (isnan(a) == 0);
1861 assert (isinf(a) == 0);
1862 assert (isnan(offset) == 0);
1863 assert (isinf(offset) == 0);
1864
1865 if (opts->ptype == FOURIER){
1866 fprintf(stderr,
1867 "Cannot initialize linear function for fourier basis\n");
1868 return NULL;
1869 }
1870 struct OrthPolyExpansion * p =
1871 orth_poly_expansion_init_from_opts(opts,2);
1872 p->coeff[1] = a / (p->p->lin_coeff * p->space_transform->lin_slope);
1873 p->coeff[0] = (offset - p->p->lin_const -
1874 p->coeff[1] * p->p->lin_coeff *
1875 p->space_transform->lin_offset) /
1876 p->p->const_term;
1877
1878 return p;
1879 }
1880
1881 /********************************************************//**
1882 * Update a linear orthonormal polynomial expansion
1883 *
1884 * \param[in] p - existing linear polynomial
1885 * \param[in] a - value of the slope function
1886 * \param[in] offset - offset
1887 *
1888 * \return 0 if succesfull, 1 otherwise
1889 *
1890 *************************************************************/
1891 int
orth_poly_expansion_linear_update(struct OrthPolyExpansion * p,double a,double offset)1892 orth_poly_expansion_linear_update(struct OrthPolyExpansion * p, double a, double offset)
1893 {
1894 assert (isnan(a) == 0);
1895 assert (isinf(a) == 0);
1896 assert (isnan(offset) == 0);
1897 assert (isinf(offset) == 0);
1898
1899 if (p->p->ptype == FOURIER){
1900 fprintf(stderr,
1901 "Cannot update linear function for fourier basis\n");
1902 return 1;
1903 }
1904
1905 p->coeff[1] = a / (p->p->lin_coeff * p->space_transform->lin_slope);
1906 p->coeff[0] = (offset - p->p->lin_const -
1907 p->coeff[1] * p->p->lin_coeff * p->space_transform->lin_offset) /
1908 p->p->const_term;
1909
1910 return 0;
1911 }
1912
1913 /********************************************************//**
1914 * Generate a quadratic orthonormal polynomial expansion
1915 a * (x-offset)^2
1916 *
1917 * \param[in] a - value of the slope function
1918 * \param[in] offset - offset
1919 * \param[in] opts - options
1920 *
1921 * \return quadratic polynomial
1922 *************************************************************/
1923 struct OrthPolyExpansion *
orth_poly_expansion_quadratic(double a,double offset,struct OpeOpts * opts)1924 orth_poly_expansion_quadratic(double a, double offset, struct OpeOpts * opts)
1925 {
1926
1927 assert (isnan(a) == 0);
1928 assert (isinf(a) == 0);
1929 assert (isnan(offset) == 0);
1930 assert (isinf(offset) == 0);
1931
1932 if (opts->ptype == FOURIER){
1933 fprintf(stderr, "Cannot initialize quadratic function for fourier basis\n");
1934 return NULL;
1935 }
1936
1937 struct OrthPolyExpansion * p = orth_poly_expansion_init_from_opts(opts, 3);
1938
1939 struct quad_func qf;
1940 qf.scale = a;
1941 qf.offset = offset;
1942 orth_poly_expansion_approx(eval_quad_func, &qf, p);
1943
1944 return p;
1945 }
1946
1947 /********************************************************//**
1948 * Generate a polynomial expansion with only the
1949 * *order* coefficient being nonzero
1950 *
1951 * \param[in] order - order of the polynomial
1952 * \param[in] opts - options for building polynomial
1953 *
1954 * \return p - orthogonal polynomial expansion
1955 *
1956 * \note Not sure about hermite!
1957 *************************************************************/
1958 struct OrthPolyExpansion *
orth_poly_expansion_genorder(size_t order,struct OpeOpts * opts)1959 orth_poly_expansion_genorder(size_t order, struct OpeOpts * opts)
1960 {
1961
1962 struct OrthPolyExpansion * p =
1963 orth_poly_expansion_init_from_opts(opts, order+1);
1964
1965 double m = space_mapping_map_inverse_deriv(p->space_transform,0);
1966 switch (opts->ptype){
1967 case LEGENDRE:
1968 p->coeff[order] = 1.0 / sqrt(p->p->norm(order)) / sqrt(2.0) / sqrt(m);
1969 break;
1970 case HERMITE:
1971 p->coeff[order] = 1.0 / sqrt(p->p->norm(order)) / sqrt(m);
1972 break;
1973 case CHEBYSHEV:
1974 p->coeff[order] = 1.0 / sqrt(p->p->norm(order)) / sqrt(m);
1975 /* double norm = orth_poly_expansion_inner(p, p); */
1976 /* p->coeff[order] /= (norm * sqrt(m)); */
1977 /* break; */
1978 /* fprintf(stderr,"cannot generate orthonormal polynomial of\n"); */
1979 /* fprintf(stderr,"order for CHEBYSHEV type\n;"); */
1980 /* exit(1); */
1981 break;
1982 case FOURIER:
1983 /* p->coeff[order] = 1.0 / sqrt(p->p->norm(order)) / sqrt(m); */
1984 /* break; */
1985 fprintf(stderr,"cannot generate orthonormal polynomial of certin\n");
1986 fprintf(stderr,"order for FOURIER type\n;");
1987 exit(1);
1988 case STANDARD:
1989 fprintf(stderr,"cannot generate orthonormal polynomial of certin\n");
1990 fprintf(stderr,"order for STANDARD type\n;");
1991 exit(1);
1992 }
1993 //printf("wherhhh\n");
1994
1995
1996 //printf("there \n");
1997 return p;
1998 }
1999
2000 /********************************************************//**
2001 Generate an orthonormal basis
2002
2003 \param[in] n - number of basis function
2004 \param[in,out] f - space to write polynomials
2005 \param[in] opts - approximation options
2006
2007 \note
2008 Uses modified gram schmidt to determine function coefficients
2009 Each function f[ii] must have the same nodes
2010 *************************************************************/
2011 void
orth_poly_expansion_orth_basis(size_t n,struct OrthPolyExpansion ** f,struct OpeOpts * opts)2012 orth_poly_expansion_orth_basis(size_t n, struct OrthPolyExpansion ** f,
2013 struct OpeOpts * opts)
2014 {
2015
2016 if (opts->ptype == CHEBYSHEV){
2017 for (size_t ii = 0; ii < n; ii++){
2018 f[ii] = orth_poly_expansion_init_from_opts(opts, ii+1);
2019 f[ii]->coeff[ii] = 1.0;
2020 }
2021 // now gram schmidt
2022 double norm, proj;
2023 for (size_t ii = 0; ii < n; ii++){
2024 norm = sqrt(orth_poly_expansion_inner(f[ii], f[ii]));
2025 if (norm > 1e-200){
2026 orth_poly_expansion_scale(1.0/norm, f[ii]);
2027 for (size_t jj = ii+1; jj < n; jj++){
2028 proj = orth_poly_expansion_inner(f[ii],f[jj]);
2029 orth_poly_expansion_axpy(-proj,f[ii],f[jj]);
2030 }
2031 }
2032 }
2033 }
2034 else if ((opts->ptype == HERMITE) || (opts->ptype == LEGENDRE)){
2035 for (size_t ii = 0; ii < n; ii++){
2036 f[ii] = orth_poly_expansion_init_from_opts(opts, ii+1);
2037 double m = space_mapping_map_inverse_deriv(f[ii]->space_transform,0);
2038 f[ii]->coeff[ii] = 1.0 / sqrt(f[ii]->p->norm(ii)) / sqrt(m);
2039 if (opts->ptype == LEGENDRE){
2040 f[ii]->coeff[ii] /= sqrt(2.0);
2041 }
2042 }
2043 }
2044 else if (opts->ptype == FOURIER) {
2045 size_t N = ope_opts_get_start(opts);
2046
2047 if (n > N){
2048 fprintf(stderr, "Cannot look for a rank so large for fourier\n");
2049
2050 /* printf("n = %zu\n", n); */
2051 /* N = N-1; */
2052 /* while (n > N){ */
2053 /* N = 2 * N; */
2054 /* } */
2055 /* N = N+1; */
2056 /* printf("N = %zu\n", N); */
2057 }
2058 for (size_t ii = 0; ii < n; ii++){
2059 f[ii] = orth_poly_expansion_init_from_opts(opts, N);
2060 if (n < N){
2061 f[ii]->ccoeff[ii] = 1.0;
2062 }
2063 }
2064 // now gram schmidt
2065 double norm, proj;
2066 for (size_t ii = 0; ii < n; ii++){
2067 /* printf("ii = %zu %zu\n", ii, f[ii]->num_poly); */
2068 norm = sqrt(orth_poly_expansion_inner(f[ii], f[ii]));
2069 if (norm > 1e-200){
2070 orth_poly_expansion_scale(1.0/norm, f[ii]);
2071 for (size_t jj = ii+1; jj < n; jj++){
2072 /* printf("jj = %zu before %zu\n", jj, f[jj]->num_poly); */
2073 proj = orth_poly_expansion_inner(f[ii],f[jj]);
2074 /* printf("jj = %zu after 1 %zu\n", jj, f[jj]->num_poly); */
2075 orth_poly_expansion_axpy(-proj,f[ii],f[jj]);
2076 /* printf("jj = %zu after 2%zu\n", jj, f[jj]->num_poly); */
2077 /* if (f[jj]->num_poly != 33){ */
2078 /* printf("WARNING\n"); */
2079 /* exit(1); */
2080 /* } */
2081 }
2082 }
2083 }
2084 }
2085 else{
2086 fprintf(stderr, "Cannot generate an orthonormal basis for polytype %d\n", opts->ptype);
2087 exit(1);
2088 }
2089 }
2090
2091
2092 /********************************************************//**
2093 * Evaluate the derivative of an orthogonal polynomial expansion
2094 *
2095 * \param[in] poly - pointer to orth poly expansion
2096 * \param[in] x - location at which to evaluate
2097 *
2098 *
2099 * \return out - value of derivative
2100 *************************************************************/
orth_poly_expansion_deriv_eval(const struct OrthPolyExpansion * poly,double x)2101 double orth_poly_expansion_deriv_eval(const struct OrthPolyExpansion * poly, double x)
2102 {
2103 assert (poly != NULL);
2104 assert (poly->kristoffel_eval == 0);
2105
2106 if (poly->p->ptype == FOURIER){
2107 return fourier_expansion_deriv_eval(poly, x);
2108 }
2109
2110 double x_normalized = space_mapping_map(poly->space_transform,x);
2111
2112 //values
2113 double p[2];
2114 double pnew;
2115
2116 //gradients
2117 double pg[2];
2118 double pgnew;
2119
2120 size_t iter = 0;
2121 double out = 0.0;
2122 p[0] = poly->p->const_term;
2123 pg[0] = 0.0;
2124 out += pg[0] * poly->coeff[iter];
2125 iter++;
2126 if (poly->num_poly > 1){
2127 p[1] = poly->p->lin_const + poly->p->lin_coeff * x_normalized;
2128 pg[1] = poly->p->lin_coeff;
2129 out += pg[1] * poly->coeff[iter];
2130 iter++;
2131 }
2132
2133 double a,b,c;
2134 for (iter = 2; iter < poly->num_poly; iter++){
2135 a = poly->p->an(iter);
2136 b = poly->p->bn(iter);
2137 c = poly->p->cn(iter);
2138 pnew = eval_orth_poly_wp(poly->p, p[0], p[1], iter, x_normalized);
2139
2140 pgnew = (a*x_normalized + b) * pg[1] + a * p[1] + c*pg[0];
2141 out += poly->coeff[iter] * pgnew;
2142
2143 p[0] = p[1];
2144 p[1] = pnew;
2145
2146 pg[0] = pg[1];
2147 pg[1] = pgnew;
2148 }
2149
2150 out *= space_mapping_map_deriv(poly->space_transform,x);
2151 return out;
2152 }
2153
2154
orth_poly_expansion_deriv_eval_for_approx(double x,void * poly)2155 static inline double orth_poly_expansion_deriv_eval_for_approx(double x, void* poly){
2156 return orth_poly_expansion_deriv_eval(poly, x);
2157 }
2158
2159 /********************************************************//**
2160 * Evaluate the derivative of an orthogonal polynomial expansion
2161 *
2162 * \param[in] poly - pointer to orth poly expansion
2163 * \param[in] x - location at which to evaluate
2164 *
2165 *
2166 * \return out - value of derivative
2167 *************************************************************/
cheb_expansion_deriv_eval(const struct OrthPolyExpansion * poly,double x)2168 double cheb_expansion_deriv_eval(const struct OrthPolyExpansion * poly, double x)
2169 {
2170 assert (poly != NULL);
2171 assert (poly->kristoffel_eval == 0);
2172
2173 double dmult = space_mapping_map_deriv(poly->space_transform,x);
2174 if (poly->num_poly == 1){
2175 return 0.0;
2176 }
2177 else if (poly->num_poly == 2){
2178 return poly->coeff[1] * dmult;
2179 }
2180
2181 double x_norm = space_mapping_map(poly->space_transform,x);
2182
2183 if (poly->num_poly == 3){
2184 return (poly->coeff[1] + poly->coeff[2] * 4 * x_norm) * dmult;
2185 }
2186
2187 double * cheb_eval = calloc_double(poly->num_poly);
2188 double * cheb_evald = calloc_double(poly->num_poly);
2189 cheb_eval[0] = 1.0;
2190 cheb_eval[1] = x_norm;
2191 cheb_eval[2] = 2.0*x_norm*cheb_eval[1] - cheb_eval[0];
2192
2193 cheb_evald[0] = 0.0;
2194 cheb_evald[1] = 1.0;
2195 cheb_evald[2] = 4.0*x_norm;
2196
2197 double out = poly->coeff[1]*cheb_evald[1] + poly->coeff[2]*cheb_evald[2];
2198 for (size_t ii = 3; ii < poly->num_poly; ii++){
2199 cheb_eval[ii] = 2.0 * x_norm * cheb_eval[ii-1] - cheb_eval[ii-2];
2200 cheb_evald[ii] = 2.0 * cheb_eval[ii-1] + 2.0 * x_norm * cheb_evald[ii-1] -
2201 cheb_evald[ii-2];
2202 out += poly->coeff[ii]*cheb_evald[ii];
2203 }
2204
2205 out *= dmult;
2206 free(cheb_eval); cheb_eval = NULL;
2207 free(cheb_evald); cheb_evald = NULL;
2208 return out;
2209 }
2210
cheb_expansion_deriv_eval_for_approx(double x,void * poly)2211 static inline double cheb_expansion_deriv_eval_for_approx(double x, void* poly){
2212 return cheb_expansion_deriv_eval(poly, x);
2213 }
2214
2215 /********************************************************//**
2216 * Evaluate the derivative of an orth poly expansion
2217 *
2218 * \param[in] pin - orthogonal polynomial expansion
2219 *
2220 * \return derivative
2221 *
2222 * \note
2223 * Could speed this up slightly by using partial sum
2224 * to keep track of sum of coefficients
2225 *************************************************************/
2226 struct OrthPolyExpansion *
orth_poly_expansion_deriv(const struct OrthPolyExpansion * pin)2227 orth_poly_expansion_deriv(const struct OrthPolyExpansion * pin)
2228 {
2229 if (pin == NULL) return NULL;
2230 assert (pin->kristoffel_eval == 0);
2231
2232 if (pin->p->ptype == FOURIER){
2233 return fourier_expansion_deriv(pin);
2234 }
2235
2236 struct OrthPolyExpansion * p = orth_poly_expansion_copy(pin);
2237 orth_poly_expansion_round(&p);
2238
2239 struct OrthPolyExpansion * out = NULL;
2240
2241 out = orth_poly_expansion_copy(p);
2242 for (size_t ii = 0; ii < out->nalloc; ii++){
2243 out->coeff[ii] = 0.0;
2244 }
2245 if (p->num_poly == 1){
2246 orth_poly_expansion_round(&out);
2247 orth_poly_expansion_free(p); p = NULL;
2248 return out;
2249 }
2250
2251 out->num_poly -= 1;
2252 if (p->p->ptype == LEGENDRE){
2253 /* if (1 == 0){ */
2254
2255 double dtransform_dx = space_mapping_map_deriv(p->space_transform,0.0);
2256 for (size_t ii = 0; ii < p->num_poly-1; ii++){ // loop over coefficients
2257 for (size_t jj = ii+1; jj < p->num_poly; jj+=2){
2258 /* out->coeff[ii] += p->coeff[jj]; */
2259 out->coeff[ii] += p->coeff[jj]*sqrt(2*jj+1);
2260 }
2261 /* out->coeff[ii] *= (double) ( 2 * (ii) + 1) * dtransform_dx; */
2262 out->coeff[ii] *= sqrt((double) ( 2 * (ii) + 1))* dtransform_dx;
2263 }
2264 }
2265 else if (p->p->ptype == CHEBYSHEV){
2266 orth_poly_expansion_approx(cheb_expansion_deriv_eval_for_approx, p, out);
2267 }
2268 else{
2269 orth_poly_expansion_approx(orth_poly_expansion_deriv_eval_for_approx, p, out);
2270 }
2271
2272 orth_poly_expansion_round(&out);
2273 orth_poly_expansion_free(p); p = NULL;
2274 return out;
2275 }
2276
2277 /********************************************************//**
2278 * Evaluate the second derivative of a chebyshev expansion
2279 *
2280 * \param[in] poly - pointer to orth poly expansion
2281 * \param[in] x - location at which to evaluate
2282 *
2283 *
2284 * \return out - value of derivative
2285 *************************************************************/
cheb_expansion_dderiv_eval(const struct OrthPolyExpansion * poly,double x)2286 double cheb_expansion_dderiv_eval(const struct OrthPolyExpansion * poly, double x)
2287 {
2288 assert (poly != NULL);
2289 assert (poly->kristoffel_eval == 0);
2290
2291
2292 double dmult = space_mapping_map_deriv(poly->space_transform,x);
2293 if (poly->num_poly <= 2){
2294 return 0.0;
2295 }
2296 else if (poly->num_poly == 3){
2297 return poly->coeff[2] * 4.0;
2298 }
2299
2300 double x_norm = space_mapping_map(poly->space_transform,x);
2301 /* printf("x_norm = %3.15G\n", x_norm); */
2302
2303 if (poly->num_poly == 4){
2304 /* printf("here yo!\n"); */
2305 return (poly->coeff[2] * 4 + poly->coeff[3]*24.0*x_norm) * dmult * dmult;
2306 }
2307
2308 double * cheb_eval = calloc_double(poly->num_poly);
2309 double * cheb_evald = calloc_double(poly->num_poly);
2310 double * cheb_evaldd = calloc_double(poly->num_poly);
2311 cheb_eval[0] = 1.0;
2312 cheb_eval[1] = x_norm;
2313 cheb_eval[2] = 2.0*x_norm*cheb_eval[1] - cheb_eval[0];
2314 cheb_eval[3] = 2.0*x_norm*cheb_eval[2] - cheb_eval[1];
2315
2316 cheb_evald[0] = 0.0;
2317 cheb_evald[1] = 1.0;
2318 cheb_evald[2] = 4.0*x_norm;
2319 cheb_evald[3] = 2.0 * cheb_eval[2] + 2.0 * x_norm * cheb_evald[2] - cheb_evald[1];
2320
2321 cheb_evaldd[0] = 0.0;
2322 cheb_evaldd[1] = 0.0;
2323 cheb_evaldd[2] = 4.0;
2324 cheb_evaldd[3] = 24.0 * x_norm;
2325
2326 double out = poly->coeff[2]*cheb_evaldd[2] + poly->coeff[3]*cheb_evaldd[3];
2327 for (size_t ii = 4; ii < poly->num_poly; ii++){
2328 cheb_eval[ii] = 2.0 * x_norm * cheb_eval[ii-1] - cheb_eval[ii-2];
2329 cheb_evald[ii] = 2.0 * cheb_eval[ii-1] + 2.0 * x_norm * cheb_evald[ii-1] -
2330 cheb_evald[ii-2];
2331 cheb_evaldd[ii] = 4.0 * cheb_evald[ii-1] + 2.0 * x_norm * cheb_evaldd[ii-1] -
2332 cheb_evaldd[ii-2];
2333
2334 out += poly->coeff[ii]*cheb_evaldd[ii];
2335 }
2336
2337 out *= dmult*dmult;
2338 /* if (fabs(x_norm) > 0.999){ */
2339 /* printf("out = %3.15G\n", out); */
2340 /* } */
2341 free(cheb_eval); cheb_eval = NULL;
2342 free(cheb_evald); cheb_evald = NULL;
2343 free(cheb_evaldd); cheb_evaldd = NULL;
2344 return out;
2345 }
2346
cheb_expansion_dderiv_eval_for_approx(double x,void * poly)2347 static inline double cheb_expansion_dderiv_eval_for_approx(double x, void* poly){
2348 return cheb_expansion_dderiv_eval(poly, x);
2349 }
2350 /********************************************************//**
2351 * Evaluate the second derivative of an orth poly expansion
2352 *
2353 * \param[in] pin - orthogonal polynomial expansion
2354 *
2355 * \return derivative
2356 *
2357 * \note
2358 * Could speed this up slightly by using partial sum
2359 * to keep track of sum of coefficients
2360 *************************************************************/
2361 struct OrthPolyExpansion *
orth_poly_expansion_dderiv(const struct OrthPolyExpansion * pin)2362 orth_poly_expansion_dderiv(const struct OrthPolyExpansion * pin)
2363 {
2364 if (pin == NULL) return NULL;
2365 assert (pin->kristoffel_eval == 0);
2366
2367 if (pin->p->ptype == FOURIER){
2368 return fourier_expansion_dderiv(pin);
2369 }
2370
2371 struct OrthPolyExpansion * p = orth_poly_expansion_copy(pin);
2372 orth_poly_expansion_round(&p);
2373
2374 struct OrthPolyExpansion * out = NULL;
2375
2376 out = orth_poly_expansion_copy(p);
2377 for (size_t ii = 0; ii < out->nalloc; ii++){
2378 out->coeff[ii] = 0.0;
2379 }
2380 if (p->num_poly < 2){
2381 orth_poly_expansion_free(p); p = NULL;
2382 return out;
2383 }
2384
2385 /* printf("lets go!\n"); */
2386 /* dprint(p->num_poly, p->coeff); */
2387 out->num_poly -= 2;
2388 if (p->p->ptype == CHEBYSHEV){
2389 orth_poly_expansion_approx(cheb_expansion_dderiv_eval_for_approx, p, out);
2390 }
2391 else{
2392
2393 struct OrthPolyExpansion * temp = orth_poly_expansion_deriv(p);
2394 orth_poly_expansion_free(out);
2395 out = orth_poly_expansion_deriv(temp);
2396 orth_poly_expansion_free(temp); temp = NULL;
2397 /* fprintf(stderr, "Cannot yet take second derivative for polynomial of type %d\n", */
2398 /* p->p->ptype); */
2399 /* exit(1); */
2400 }
2401
2402 orth_poly_expansion_round(&out);
2403 orth_poly_expansion_free(p); p = NULL;
2404 return out;
2405 }
2406
2407 /********************************************************//**
2408 Take a second derivative and enforce periodic bc
2409 **************************************************************/
orth_poly_expansion_dderiv_periodic(const struct OrthPolyExpansion * f)2410 struct OrthPolyExpansion * orth_poly_expansion_dderiv_periodic(const struct OrthPolyExpansion * f)
2411 {
2412 if (f->p->ptype == FOURIER){
2413 return orth_poly_expansion_dderiv(f);
2414 }
2415 else{
2416 NOT_IMPLEMENTED_MSG("orth_poly_expansion_dderiv_periodic");
2417 exit(1);
2418 }
2419 }
2420
2421 /********************************************************//**
2422 * free the memory of an orthonormal polynomial expansion
2423 *
2424 * \param[in,out] p - orthogonal polynomial expansion
2425 *************************************************************/
orth_poly_expansion_free(struct OrthPolyExpansion * p)2426 void orth_poly_expansion_free(struct OrthPolyExpansion * p){
2427 if (p != NULL){
2428 if (p->p->ptype == FOURIER){
2429 free(p->ccoeff); p->ccoeff = NULL;
2430 }
2431 free_orth_poly(p->p); p->p = NULL;
2432 space_mapping_free(p->space_transform); p->space_transform = NULL;
2433 free(p->coeff); p->coeff = NULL;
2434
2435 free(p); p = NULL;
2436 }
2437 }
2438
2439 /********************************************************//**
2440 * Serialize orth_poly_expansion
2441 *
2442 * \param[in] ser - location to which to serialize
2443 * \param[in] p - polynomial
2444 * \param[in] totSizeIn - if not null then only return total size of
2445 * array without serialization! if NULL then serialiaze
2446 *
2447 * \return ptr : pointer to end of serialization
2448 *************************************************************/
2449 unsigned char *
serialize_orth_poly_expansion(unsigned char * ser,struct OrthPolyExpansion * p,size_t * totSizeIn)2450 serialize_orth_poly_expansion(unsigned char * ser,
2451 struct OrthPolyExpansion * p,
2452 size_t * totSizeIn)
2453 {
2454 // order is ptype->lower_bound->upper_bound->orth_poly->coeff
2455
2456 size_t totsize;
2457 if (p->p->ptype != FOURIER){
2458 totsize = sizeof(int) + 2*sizeof(double) +
2459 p->num_poly * sizeof(double) + sizeof(size_t);
2460 }
2461 else{
2462 /* fprintf(stderr, "Cannot serialized fourier polynomials yet\n"); */
2463 /* exit(1); */
2464 totsize = sizeof(int) + 2*sizeof(double) +
2465 2*p->num_poly * sizeof(double) + sizeof(size_t);
2466 }
2467
2468 size_t size_mapping;
2469 serialize_space_mapping(NULL,p->space_transform,&size_mapping);
2470 totsize += size_mapping;
2471
2472 totsize += sizeof(int); // for kristoffel flag
2473 if (totSizeIn != NULL){
2474 *totSizeIn = totsize;
2475 return ser;
2476 }
2477 unsigned char * ptr = serialize_int(ser, p->p->ptype);
2478 ptr = serialize_double(ptr, p->lower_bound);
2479 ptr = serialize_double(ptr, p->upper_bound);
2480 if (p->p->ptype != FOURIER){
2481 ptr = serialize_doublep(ptr, p->coeff, p->num_poly);
2482 }
2483 else{
2484 ptr = serialize_size_t(ptr, p->num_poly);
2485 for (size_t ii = 0; ii < p->num_poly; ii++){
2486 ptr = serialize_double(ptr, creal(p->ccoeff[ii]));
2487 ptr = serialize_double(ptr, cimag(p->ccoeff[ii]));
2488 }
2489 }
2490
2491 ptr = serialize_space_mapping(ptr,p->space_transform,NULL);
2492 ptr = serialize_int(ptr,p->kristoffel_eval);
2493 return ptr;
2494 }
2495
2496 /********************************************************//**
2497 * Deserialize orth_poly_expansion
2498 *
2499 * \param[in] ser - input string
2500 * \param[in,out] poly - poly expansion
2501 *
2502 * \return ptr - ser + number of bytes of poly expansion
2503 *************************************************************/
2504 unsigned char *
deserialize_orth_poly_expansion(unsigned char * ser,struct OrthPolyExpansion ** poly)2505 deserialize_orth_poly_expansion(
2506 unsigned char * ser,
2507 struct OrthPolyExpansion ** poly)
2508 {
2509
2510 size_t num_poly = 0;
2511 //size_t npoly_check = 0;
2512 double lower_bound = 0;
2513 double upper_bound = 0;
2514 double * coeff = NULL;
2515 complex double * ccoeff = NULL;
2516 struct OrthPoly * p = NULL;
2517 struct SpaceMapping * map = NULL;
2518 // order is ptype->lower_bound->upper_bound->orth_poly->coeff
2519 p = deserialize_orth_poly(ser);
2520 unsigned char * ptr = ser + sizeof(int);
2521 ptr = deserialize_double(ptr,&lower_bound);
2522 ptr = deserialize_double(ptr,&upper_bound);
2523
2524 if (p->ptype != FOURIER){
2525 ptr = deserialize_doublep(ptr, &coeff, &num_poly);
2526 }
2527 else{
2528 ptr = deserialize_size_t(ptr, &num_poly);
2529 /* printf("num_poly = %zu\n", num_poly); */
2530 ccoeff = malloc( num_poly * sizeof(complex double));
2531 for (size_t ii = 0; ii < num_poly; ii++){
2532 double cr;
2533 double ci;
2534 ptr = deserialize_double(ptr, &cr);
2535 ptr = deserialize_double(ptr, &ci);
2536 ccoeff[ii] = cr + ci * (complex double) I;
2537 /* printf("%3.5f, %3.5f\n", cr, ci); */
2538 }
2539
2540 }
2541 ptr = deserialize_space_mapping(ptr, &map);
2542
2543 int kristoffel_eval;
2544 ptr = deserialize_int(ptr,&kristoffel_eval);
2545 if ( NULL == (*poly = malloc(sizeof(struct OrthPolyExpansion)))){
2546 fprintf(stderr, "failed to allocate memory for poly exp.\n");
2547 exit(1);
2548 }
2549 (*poly)->num_poly = num_poly;
2550 (*poly)->lower_bound = lower_bound;
2551 (*poly)->upper_bound = upper_bound;
2552 (*poly)->coeff = coeff;
2553 (*poly)->ccoeff = ccoeff;
2554 (*poly)->nalloc = num_poly;//+OPECALLOC;
2555 (*poly)->p = p;
2556 (*poly)->space_transform = map;
2557 (*poly)->kristoffel_eval = kristoffel_eval;
2558 return ptr;
2559 }
2560
2561 /********************************************************//**
2562 Save an orthonormal polynomial expansion in text format
2563
2564 \param[in] f - function to save
2565 \param[in] stream - stream to save it to
2566 \param[in] prec - precision with which to save it
2567 ************************************************************/
orth_poly_expansion_savetxt(const struct OrthPolyExpansion * f,FILE * stream,size_t prec)2568 void orth_poly_expansion_savetxt(const struct OrthPolyExpansion * f,
2569 FILE * stream, size_t prec)
2570 {
2571 assert (f != NULL);
2572 if (f->p->ptype == FOURIER){
2573 fprintf(stderr, "Cannot savetxt fourier polynomials yet\n");
2574 exit(1);
2575 }
2576
2577 fprintf(stream,"%d ",f->p->ptype);
2578 fprintf(stream,"%3.*G ",(int)prec,f->lower_bound);
2579 fprintf(stream,"%3.*G ",(int)prec,f->upper_bound);
2580 fprintf(stream,"%zu ",f->num_poly);
2581 for (size_t ii = 0; ii < f->num_poly; ii++){
2582 if (prec < 100){
2583 fprintf(stream, "%3.*G ",(int)prec,f->coeff[ii]);
2584 }
2585 }
2586 space_mapping_savetxt(f->space_transform,stream,prec);
2587 fprintf(stream,"%d ",f->kristoffel_eval);
2588 }
2589
2590 /********************************************************//**
2591 Load an orthonormal polynomial expansion in text format
2592
2593 \param[in] stream - stream to save it to
2594
2595 \return Orthonormal polynomial expansion
2596 ************************************************************/
2597 struct OrthPolyExpansion *
orth_poly_expansion_loadtxt(FILE * stream)2598 orth_poly_expansion_loadtxt(FILE * stream)//l, size_t prec)
2599 {
2600
2601 enum poly_type ptype;
2602 double lower_bound = 0;
2603 double upper_bound = 0;
2604 size_t num_poly;
2605
2606 int ptypeint;
2607 int num = fscanf(stream,"%d ",&ptypeint);
2608 ptype = (enum poly_type)ptypeint;
2609 assert (num == 1);
2610 num = fscanf(stream,"%lG ",&lower_bound);
2611 assert (num == 1);
2612 num = fscanf(stream,"%lG ",&upper_bound);
2613 assert (num == 1);
2614 num = fscanf(stream,"%zu ",&num_poly);
2615 assert (num == 1);
2616
2617 struct OrthPolyExpansion * ope =
2618 orth_poly_expansion_init(ptype,num_poly,lower_bound,upper_bound);
2619
2620 for (size_t ii = 0; ii < ope->num_poly; ii++){
2621 num = fscanf(stream, "%lG ",ope->coeff+ii);
2622 assert (num == 1);
2623 }
2624
2625 space_mapping_free(ope->space_transform); ope->space_transform = NULL;
2626 ope->space_transform = space_mapping_loadtxt(stream);
2627
2628 int kristoffel_eval;
2629 num = fscanf(stream,"%d ",&kristoffel_eval);
2630 assert (num == 1);
2631 ope->kristoffel_eval = kristoffel_eval;
2632
2633 return ope;
2634 }
2635
2636 /********************************************************//**
2637 * Convert an orthogonal polynomial expansion to a standard_polynomial
2638 *
2639 * \param[in] p - polynomial
2640 *
2641 * \return sp - standard polynomial
2642 *************************************************************/
2643 struct StandardPoly *
orth_poly_expansion_to_standard_poly(struct OrthPolyExpansion * p)2644 orth_poly_expansion_to_standard_poly(struct OrthPolyExpansion * p)
2645 {
2646 if (p->p->ptype == FOURIER){
2647 fprintf(stderr, "Cannot convert fourier poly to a standard poly");
2648 exit(1);
2649 }
2650 struct StandardPoly * sp =
2651 standard_poly_init(p->num_poly,p->lower_bound,p->upper_bound);
2652
2653 double m = (p->p->upper - p->p->lower) / (p->upper_bound - p->lower_bound);
2654 double off = p->p->upper - m * p->upper_bound;
2655
2656 size_t ii, jj;
2657 size_t n = p->num_poly-1;
2658
2659 sp->coeff[0] = p->coeff[0]*p->p->const_term;
2660
2661 if (n > 0){
2662 sp->coeff[0]+=p->coeff[1] * (p->p->lin_const + p->p->lin_coeff * off);
2663 sp->coeff[1]+=p->coeff[1] * p->p->lin_coeff * m;
2664 }
2665 if (n > 1){
2666
2667 double * a = calloc_double(n+1); //n-2 poly
2668 a[0] = p->p->const_term;
2669 double * b = calloc_double(n+1); // n- 1poly
2670 double * c = calloc_double(n+1); // n- 1poly
2671 b[0] = p->p->lin_const + p->p->lin_coeff * off;
2672 b[1] = p->p->lin_coeff * m;
2673 for (ii = 2; ii < n+1; ii++){ // starting at the order 2 polynomial
2674 c[0] = (p->p->bn(ii) + p->p->an(ii)*off) * b[0] +
2675 p->p->cn(ii) * a[0];
2676 sp->coeff[0] += p->coeff[ii] * c[0];
2677 for (jj = 1; jj < ii-1; jj++){
2678 c[jj] = (p->p->an(ii) * m) * b[jj-1] +
2679 (p->p->bn(ii) + p->p->an(ii) * off) * b[jj] +
2680 p->p->cn(ii) * a[jj];
2681 sp->coeff[jj] += p->coeff[ii] * c[jj];
2682 }
2683 c[ii-1] = (p->p->an(ii) * m) * b[ii-2] +
2684 (p->p->bn(ii) + p->p->an(ii) * off) * b[ii-1];
2685 c[ii] = (p->p->an(ii) * m) * b[ii-1];
2686
2687 sp->coeff[ii-1] += p->coeff[ii] * c[ii-1];
2688 sp->coeff[ii] += p->coeff[ii] * c[ii];
2689
2690 memcpy(a, b, ii * sizeof(double));
2691 memcpy(b, c, (ii+1) * sizeof(double));
2692 }
2693
2694 free(a);
2695 free(b);
2696 free(c);
2697 }
2698
2699 // Need to do something with lower and upper bounds!!
2700 return sp;
2701 }
2702
2703 /********************************************************//**
2704 * Evaluate each orthonormal polynomial expansion that is in an
2705 * array of generic functions
2706 *
2707 * \param[in] n - number of polynomials
2708 * \param[in] parr - polynomial expansions
2709 * \param[in] x - location at which to evaluate
2710 * \param[in,out] out - evaluations
2711 *
2712 * \return 0 - successful
2713 *
2714 * \note
2715 * Assumes all functions have the same bounds
2716 *************************************************************/
orth_poly_expansion_arr_eval(size_t n,struct OrthPolyExpansion ** parr,double x,double * out)2717 int orth_poly_expansion_arr_eval(size_t n,
2718 struct OrthPolyExpansion ** parr,
2719 double x, double * out)
2720 {
2721
2722
2723 if (parr[0]->kristoffel_eval == 1){
2724 for (size_t ii = 0; ii < n; ii++){
2725 out[ii] = orth_poly_expansion_eval(parr[ii],x);
2726 }
2727 return 0;
2728 }
2729
2730 int all_same = 1;
2731 enum poly_type ptype = parr[0]->p->ptype;
2732 for (size_t ii = 1; ii < n; ii++){
2733 if (parr[ii]->p->ptype != ptype){
2734 all_same = 0;
2735 break;
2736 }
2737 }
2738
2739 if ((all_same == 0) || (ptype == CHEBYSHEV) || (ptype == FOURIER)){
2740 for (size_t ii = 0; ii < n; ii++){
2741 out[ii] = orth_poly_expansion_eval(parr[ii],x);
2742 }
2743 return 0;
2744 }
2745
2746 // all the polynomials are of the same type
2747 size_t maxpoly = 0;
2748 for (size_t ii = 0; ii < n; ii++){
2749 if (parr[ii]->num_poly > maxpoly){
2750 maxpoly = parr[ii]->num_poly;
2751 }
2752 out[ii] = 0.0;
2753 }
2754
2755
2756 double x_norm = space_mapping_map(parr[0]->space_transform,x);
2757
2758 // double out = 0.0;
2759 double p[2];
2760 double pnew;
2761 p[0] = parr[0]->p->const_term;
2762 size_t iter = 0;
2763 for (size_t ii = 0; ii < n; ii++){
2764 out[ii] += p[0] * parr[ii]->coeff[iter];
2765 }
2766 iter++;
2767 p[1] = parr[0]->p->lin_const + parr[0]->p->lin_coeff * x_norm;
2768 for (size_t ii = 0; ii < n; ii++){
2769 if (parr[ii]->num_poly > iter){
2770 out[ii] += p[1] * parr[ii]->coeff[iter];
2771 }
2772 }
2773 iter++;
2774 for (iter = 2; iter < maxpoly; iter++){
2775 pnew = eval_orth_poly_wp(parr[0]->p, p[0], p[1], iter, x_norm);
2776 for (size_t ii = 0; ii < n; ii++){
2777 if (parr[ii]->num_poly > iter){
2778 out[ii] += parr[ii]->coeff[iter] * pnew;
2779 }
2780 }
2781 p[0] = p[1];
2782 p[1] = pnew;
2783 }
2784
2785 return 0;
2786 }
2787
2788 /********************************************************//**
2789 * Evaluate each orthonormal polynomial expansion that is in an
2790 * array of generic functions at an array of points
2791 *
2792 * \param[in] n - number of polynomials
2793 * \param[in] parr - polynomial expansions (all have the same bounds)
2794 * \param[in] N - number of evaluations
2795 * \param[in] x - locations at which to evaluate
2796 * \param[in] incx - increment between locations
2797 * \param[in,out] y - evaluations
2798 * \param[in] incy - increment between evaluations of array (at least n)
2799 *
2800 * \return 0 - successful
2801 *
2802 * \note
2803 * Assumes all functions have the same bounds
2804 *************************************************************/
orth_poly_expansion_arr_evalN(size_t n,struct OrthPolyExpansion ** parr,size_t N,const double * x,size_t incx,double * y,size_t incy)2805 int orth_poly_expansion_arr_evalN(size_t n,
2806 struct OrthPolyExpansion ** parr,
2807 size_t N,
2808 const double * x, size_t incx,
2809 double * y, size_t incy)
2810 {
2811 if (parr[0]->kristoffel_eval == 1){
2812 for (size_t jj = 0; jj < N; jj++){
2813 for (size_t ii = 0; ii < n; ii++){
2814 y[ii + jj * incy] = orth_poly_expansion_eval(parr[ii],x[jj*incx]);
2815 /* printf("y = %G\n",y[ii+jj*incy]); */
2816 }
2817 }
2818 return 0;
2819 }
2820
2821
2822 for (size_t jj = 0; jj < N; jj++){
2823 for (size_t ii = 0; ii < n; ii++){
2824 y[ii + jj * incy] = 0.0;
2825 }
2826 }
2827
2828 int res;
2829 for (size_t jj = 0; jj < N; jj++){
2830 res = orth_poly_expansion_arr_eval(n, parr,
2831 x[jj*incx], y + jj*incy);
2832 if (res != 0){
2833 return res;
2834 }
2835 }
2836
2837
2838 for (size_t jj = 0; jj < N; jj++){
2839 for (size_t ii = 0; ii < n; ii++){
2840 if (isnan(y[ii + jj* incy]) || y[ii+jj * incy] > 1e100){
2841 fprintf(stderr,"Warning, evaluation in orth_poly_expansion_array_eval is nan\n");
2842 fprintf(stderr,"Polynomial %zu, evaluation %zu\n",ii,jj);
2843 print_orth_poly_expansion(parr[ii],0,NULL,stderr);
2844 exit(1);
2845 }
2846 else if (isinf(y[ii + jj * incy])){
2847 fprintf(stderr,"Warning, evaluation in orth_poly_expansion_array_eval inf\n");
2848 exit(1);
2849 }
2850 }
2851 }
2852
2853
2854 return 0;
2855 }
2856
2857 /********************************************************//**
2858 * Evaluate a Chebyshev polynomial expansion using clenshaw algorithm
2859 *
2860 * \param[in] poly - polynomial expansion
2861 * \param[in] x - location at which to evaluate
2862 *
2863 * \return out - polynomial value
2864 *************************************************************/
chebyshev_poly_expansion_eval(const struct OrthPolyExpansion * poly,double x)2865 double chebyshev_poly_expansion_eval(const struct OrthPolyExpansion * poly, double x)
2866 {
2867
2868 assert (poly->kristoffel_eval == 0);
2869
2870 double p[2] = {0.0, 0.0};
2871 double pnew;
2872
2873 double x_norm = space_mapping_map(poly->space_transform,x);
2874 size_t n = poly->num_poly-1;
2875 while (n > 0){
2876 pnew = poly->coeff[n] + 2.0 * x_norm * p[0] - p[1];
2877 p[1] = p[0];
2878 p[0] = pnew;
2879 n--;
2880 }
2881 double out = poly->coeff[0] + x_norm * p[0] - p[1];
2882 return out;
2883 }
2884
2885 /********************************************************//**
2886 * Evaluate a polynomial expansion consisting of sequentially increasing
2887 * order polynomials from the same family.
2888 *
2889 * \param[in] poly - polynomial expansion
2890 * \param[in] x - location at which to evaluate
2891 *
2892 * \return out - polynomial value
2893 *************************************************************/
orth_poly_expansion_eval(const struct OrthPolyExpansion * poly,double x)2894 double orth_poly_expansion_eval(const struct OrthPolyExpansion * poly, double x)
2895 {
2896 double out = 0.0;
2897 if (poly->p->ptype == FOURIER){
2898 out = fourier_expansion_eval(poly,x);
2899 }
2900 else if (poly->p->ptype != CHEBYSHEV){
2901 size_t iter = 0;
2902 double p [2];
2903 double pnew;
2904
2905 double x_normalized = space_mapping_map(poly->space_transform,x);
2906
2907 double den = 0.0;
2908
2909 p[0] = poly->p->const_term;
2910 out += p[0] * poly->coeff[iter];
2911
2912 den += p[0]*p[0];
2913
2914 iter++;
2915 if (poly->num_poly > 1){
2916 p[1] = poly->p->lin_const + poly->p->lin_coeff * x_normalized;
2917 out += p[1] * poly->coeff[iter];
2918
2919 den += p[1]*p[1];
2920 iter++;
2921 }
2922 for (iter = 2; iter < poly->num_poly; iter++){
2923 pnew = eval_orth_poly_wp(poly->p, p[0], p[1], iter, x_normalized);
2924 out += poly->coeff[iter] * pnew;
2925 p[0] = p[1];
2926 p[1] = pnew;
2927
2928 den += pnew*pnew;
2929 }
2930
2931
2932 if (poly->kristoffel_eval == 1){
2933 /* printf("normalizing for kristoffel out = %G, %G\n",out,den); */
2934 /* dprint(poly->num_poly,poly->coeff); */
2935 /* out /= sqrt(den); */
2936 out /= sqrt( den / ( (double) poly->num_poly )) ;
2937 }
2938 }
2939 else{
2940 out = chebyshev_poly_expansion_eval(poly,x);
2941 }
2942 return out;
2943 }
2944
2945 /********************************************************//**
2946 * Get the kristoffel weight of an orthonormal polynomial expansion
2947 *
2948 * \param[in] poly - polynomial expansion
2949 * \param[in] x - location at which to evaluate
2950 *
2951 * \return out - weight
2952 *************************************************************/
orth_poly_expansion_get_kristoffel_weight(const struct OrthPolyExpansion * poly,double x)2953 double orth_poly_expansion_get_kristoffel_weight(const struct OrthPolyExpansion * poly, double x)
2954 {
2955 assert (poly != NULL);
2956 if (poly->p->ptype == FOURIER){
2957 fprintf(stderr,
2958 "Cannot get kristoffel_weight for fourier basis\n");
2959 exit(1);
2960 }
2961 size_t iter = 0;
2962 double p [2];
2963 double pnew;
2964
2965 double x_normalized = space_mapping_map(poly->space_transform,x);
2966 double den = 0.0;
2967
2968 p[0] = poly->p->const_term;
2969
2970
2971 den += p[0]*p[0];
2972
2973 iter++;
2974 if (poly->num_poly > 1){
2975 p[1] = poly->p->lin_const + poly->p->lin_coeff * x_normalized;
2976
2977 den += p[1]*p[1];
2978 iter++;
2979 }
2980 for (iter = 2; iter < poly->num_poly; iter++){
2981 pnew = eval_orth_poly_wp(poly->p, p[0], p[1], iter, x_normalized);
2982
2983 p[0] = p[1];
2984 p[1] = pnew;
2985
2986 den += pnew*pnew;
2987 }
2988
2989 // Normalize by number of functions
2990 return sqrt(den / ( (double) poly->num_poly ) );
2991 //return sqrt(den);
2992 }
2993
2994 /********************************************************//**
2995 * Evaluate a polynomial expansion consisting of sequentially increasing
2996 * order polynomials from the same family.
2997 *
2998 * \param[in] poly - function
2999 * \param[in] N - number of evaluations
3000 * \param[in] x - location at which to evaluate
3001 * \param[in] incx - increment of x
3002 * \param[in,out] y - allocated space for evaluations
3003 * \param[in] incy - increment of y*
3004 *
3005 * \note Currently just calls the single evaluation code
3006 * Note sure if this is optimal, cache-wise
3007 *************************************************************/
orth_poly_expansion_evalN(const struct OrthPolyExpansion * poly,size_t N,const double * x,size_t incx,double * y,size_t incy)3008 void orth_poly_expansion_evalN(const struct OrthPolyExpansion * poly, size_t N,
3009 const double * x, size_t incx, double * y, size_t incy)
3010 {
3011 for (size_t ii = 0; ii < N; ii++){
3012 y[ii*incy] = orth_poly_expansion_eval(poly,x[ii*incx]);
3013 }
3014 }
3015
3016 /********************************************************//*
3017 * Evaluate the gradient of an orthonormal polynomial expansion
3018 * with respect to the parameters
3019 *
3020 * \param[in] poly - polynomial expansion
3021 * \param[in] nx - number of x points
3022 * \param[in] x - location at which to evaluate
3023 * \param[in,out] grad - gradient values (N,nx)
3024 *
3025 * \return 0 success, 1 otherwise
3026 *************************************************************/
orth_poly_expansion_param_grad_eval(const struct OrthPolyExpansion * poly,size_t nx,const double * x,double * grad)3027 int orth_poly_expansion_param_grad_eval(
3028 const struct OrthPolyExpansion * poly, size_t nx, const double * x, double * grad)
3029 {
3030 if (poly->p->ptype == FOURIER){
3031 fprintf(stderr, "Cannot perform param_grad_eval with fourier basis\n");
3032 exit(1);
3033 }
3034 size_t nparam = orth_poly_expansion_get_num_params(poly);
3035 for (size_t ii = 0; ii < nx; ii++){
3036
3037 double p[2];
3038 double pnew;
3039
3040 double x_norm = space_mapping_map(poly->space_transform,x[ii]);
3041
3042 size_t iter = 0;
3043 p[0] = poly->p->const_term;
3044 double den = p[0]*p[0];
3045
3046 grad[ii*nparam] = p[0];
3047 iter++;
3048 if (poly->num_poly > 1){
3049 p[1] = poly->p->lin_const + poly->p->lin_coeff * x_norm;
3050 grad[ii*nparam + iter] = p[1];
3051 iter++;
3052 den += p[1]*p[1];
3053
3054 for (iter = 2; iter < poly->num_poly; iter++){
3055 pnew = (poly->p->an(iter)*x_norm + poly->p->bn(iter)) * p[1] + poly->p->cn(iter) * p[0];
3056 grad[ii*nparam + iter] = pnew;
3057 den += pnew * pnew;
3058 p[0] = p[1];
3059 p[1] = pnew;
3060 }
3061 }
3062
3063 if (poly->kristoffel_eval == 1){
3064 /* printf("gradient normalized by kristoffel %G\n",den); */
3065 for (size_t jj = 0; jj < poly->num_poly; jj++){
3066 // Normalize by number of functions
3067 grad[ii*nparam+jj] /= sqrt(den / ( (double) poly->num_poly ));
3068 //grad[ii*nparam+jj] /= sqrt(den);
3069 }
3070 }
3071 }
3072 return 0;
3073 }
3074
3075
3076 /********************************************************//*
3077 * Evaluate the gradient of an orthonormal polynomial expansion
3078 * with respect to the parameters
3079 *
3080 * \param[in] poly - polynomial expansion
3081 * \param[in] x - location at which to evaluate
3082 * \param[in,out] grad - gradient values (N,nx)
3083 *
3084 * \return evaluation
3085 *************************************************************/
orth_poly_expansion_param_grad_eval2(const struct OrthPolyExpansion * poly,double x,double * grad)3086 double orth_poly_expansion_param_grad_eval2(
3087 const struct OrthPolyExpansion * poly, double x, double * grad)
3088 {
3089 if (poly->p->ptype == FOURIER){
3090 fprintf(stderr, "Cannot perform param_grad_eval with fourier basis\n");
3091 exit(1);
3092 }
3093
3094 double out = 0.0;
3095
3096 double p[2];
3097 double pnew;
3098
3099 double x_norm = space_mapping_map(poly->space_transform,x);
3100
3101 double den = 0.0;
3102 size_t iter = 0;
3103 p[0] = poly->p->const_term;
3104 grad[0] = p[0];
3105 den += p[0]*p[0];
3106
3107 out += p[0]*poly->coeff[0];
3108 iter++;
3109 if (poly->num_poly > 1){
3110 p[1] = poly->p->lin_const + poly->p->lin_coeff * x_norm;
3111 grad[iter] = p[1];
3112 den += p[1]*p[1];
3113 out += p[1]*poly->coeff[1];
3114 iter++;
3115
3116 for (iter = 2; iter < poly->num_poly; iter++){
3117 pnew = (poly->p->an(iter)*x_norm + poly->p->bn(iter)) * p[1] + poly->p->cn(iter) * p[0];
3118 grad[iter] = pnew;
3119 out += pnew*poly->coeff[iter];
3120 p[0] = p[1];
3121 p[1] = pnew;
3122
3123 den += pnew * pnew;
3124 }
3125 }
3126 if (poly->kristoffel_eval == 1){
3127 /* printf("gradient normalized by kristoffel %G\n",den); */
3128 for (size_t jj = 0; jj < poly->num_poly; jj++){
3129 // Normalize by number of functions
3130 grad[jj] /= sqrt(den / ( (double) poly->num_poly ));
3131 /* grad[jj] /= sqrt(den); */
3132 }
3133 }
3134 return out;
3135 }
3136
3137 /********************************************************//**
3138 Take a gradient of the squared norm
3139 with respect to its parameters, and add a scaled version
3140 of this gradient to *grad*
3141
3142 \param[in] poly - polynomial
3143 \param[in] scale - scaling for additional gradient
3144 \param[in,out] grad - gradient, on output adds scale * new_grad
3145
3146 \return 0 - success, 1 -failure
3147
3148 ************************************************************/
3149 int
orth_poly_expansion_squared_norm_param_grad(const struct OrthPolyExpansion * poly,double scale,double * grad)3150 orth_poly_expansion_squared_norm_param_grad(const struct OrthPolyExpansion * poly,
3151 double scale, double * grad)
3152 {
3153 if (poly->p->ptype == FOURIER){
3154 fprintf(stderr, "Cannot perform squared_norm_param_grad with fourier basis\n");
3155 exit(1);
3156 }
3157
3158 assert (poly->kristoffel_eval == 0);
3159 int res = 1;
3160
3161
3162 // assuming linear transformation
3163 double dtransform_dx = space_mapping_map_deriv(poly->space_transform,0.0);
3164 if (poly->p->ptype == LEGENDRE){
3165 for (size_t ii = 0; ii < poly->num_poly; ii++){
3166 //the extra 2 is for the weight
3167 grad[ii] += 2.0*scale * poly->coeff[ii] * poly->p->norm(ii) * /* *dtransform_dx * */
3168 (poly->upper_bound-poly->lower_bound);
3169 }
3170 res = 0;
3171 }
3172 else if (poly->p->ptype == HERMITE){
3173 for (size_t ii = 0; ii < poly->num_poly; ii++){
3174 grad[ii] += 2.0*scale * poly->coeff[ii] * poly->p->norm(ii) * dtransform_dx;
3175 }
3176 res = 0;
3177 }
3178 else if (poly->p->ptype == CHEBYSHEV){
3179 double * temp = calloc_double(poly->num_poly);
3180 for (size_t ii = 0; ii < poly->num_poly; ii++){
3181 temp[ii] = 2.0/(1.0 - (double)2*ii*2*ii);
3182 }
3183 for (size_t ii = 0; ii < poly->num_poly; ii++){
3184 for (size_t jj = 0; jj < poly->num_poly; jj++){
3185 size_t n1 = ii+jj;
3186 size_t n2;
3187 if (ii > jj){
3188 n2 = ii-jj;
3189 }
3190 else{
3191 n2 = jj-ii;
3192 }
3193 if (n1 % 2 == 0){
3194 grad[ii] += 2.0*scale*temp[n1/2] * dtransform_dx;
3195 }
3196 if (n2 % 2 == 0){
3197 grad[ii] += 2.0*scale*temp[n2/2] * dtransform_dx;
3198 }
3199 }
3200 }
3201 free(temp); temp = NULL;
3202 res = 0;
3203 }
3204 else{
3205 fprintf(stderr,
3206 "Cannot evaluate derivative with respect to parameters for polynomial type %d\n",
3207 poly->p->ptype);
3208 exit(1);
3209 }
3210 return res;
3211 }
3212
3213 /********************************************************//**
3214 Squared norm of a function in RKHS
3215
3216 \param[in] poly - polynomial
3217 \param[in] decay_type - type of decay
3218 \param[in] decay_param - parameter of decay
3219
3220 \return 0 - success, 1 -failure
3221 ************************************************************/
3222 double
orth_poly_expansion_rkhs_squared_norm(const struct OrthPolyExpansion * poly,enum coeff_decay_type decay_type,double decay_param)3223 orth_poly_expansion_rkhs_squared_norm(const struct OrthPolyExpansion * poly,
3224 enum coeff_decay_type decay_type,
3225 double decay_param)
3226 {
3227
3228 assert (poly->kristoffel_eval == 0);
3229 if (poly->p->ptype == FOURIER){
3230 fprintf(stderr, "Cannot perform rkhs_squared_norm with fourier basis\n");
3231 exit(1);
3232 }
3233
3234 // assuming linear transformation
3235 double m = space_mapping_map_deriv(poly->space_transform,0.0);
3236 /* double m = (poly->upper_bound-poly->lower_bound) /(poly->p->upper - poly->p->lower); */
3237 double out = 0.0;
3238 if (poly->p->ptype == LEGENDRE){
3239 if (decay_type == ALGEBRAIC){
3240 for (size_t ii = 0; ii < poly->num_poly; ii++){
3241 out += poly->coeff[ii] * poly->coeff[ii]*pow(decay_param,ii) * poly->p->norm(ii)*2.0 * m;
3242 }
3243 }
3244 else if (decay_type == EXPONENTIAL){
3245 for (size_t ii = 0; ii < poly->num_poly; ii++){
3246 out += poly->coeff[ii] * poly->coeff[ii]*pow((double)ii+1.0,-decay_param)*m*poly->p->norm(ii)*2.0;
3247 }
3248 }
3249 else{
3250 for (size_t ii = 0; ii < poly->num_poly; ii++){
3251 out += poly->coeff[ii] * poly->coeff[ii] * poly->p->norm(ii)*2.0*m;
3252 }
3253 }
3254 }
3255 else if (poly->p->ptype == HERMITE){
3256 if (decay_type == ALGEBRAIC){
3257 for (size_t ii = 0; ii < poly->num_poly; ii++){
3258 out += poly->coeff[ii] * poly->coeff[ii]*pow(decay_param,ii) * poly->p->norm(ii);
3259 }
3260 }
3261 else if (decay_type == EXPONENTIAL){
3262 for (size_t ii = 0; ii < poly->num_poly; ii++){
3263 out += poly->coeff[ii] * poly->coeff[ii]*pow((double)ii+1.0,-decay_param) * poly->p->norm(ii);
3264 }
3265 }
3266 else{
3267 for (size_t ii = 0; ii < poly->num_poly; ii++){
3268 out += poly->coeff[ii] * poly->coeff[ii] * poly->p->norm(ii);
3269 }
3270 }
3271 }
3272 else if (poly->p->ptype == CHEBYSHEV){
3273 double * temp = calloc_double(poly->num_poly);
3274 for (size_t ii = 0; ii < poly->num_poly; ii++){
3275 temp[ii] = 2.0/(1.0 - (double)2*ii*2*ii);
3276 }
3277 for (size_t ii = 0; ii < poly->num_poly; ii++){
3278 double temp_sum = 0.0;
3279 for (size_t jj = 0; jj < poly->num_poly; jj++){
3280 size_t n1 = ii+jj;
3281 size_t n2;
3282 if (ii > jj){
3283 n2 = ii-jj;
3284 }
3285 else{
3286 n2 = jj-ii;
3287 }
3288 if (n1 % 2 == 0){
3289 temp_sum += poly->coeff[jj]*temp[n1/2];
3290 }
3291 if (n2 % 2 == 0){
3292 temp_sum += poly->coeff[jj]*temp[n2/2];
3293 }
3294 }
3295 if (decay_type == ALGEBRAIC){
3296 out += temp_sum*temp_sum * pow(decay_param,ii)*m;
3297 }
3298 else if (decay_type == EXPONENTIAL){
3299 out += temp_sum*temp_sum * pow((double)ii+1.0,-decay_param)*m;
3300 }
3301 else{
3302 out += temp_sum * temp_sum*m;
3303 }
3304 }
3305 free(temp); temp = NULL;
3306 }
3307 else{
3308 fprintf(stderr, "Cannot evaluate derivative with respect to parameters for polynomial type %d\n",poly->p->ptype);
3309 exit(1);
3310 }
3311 return out;
3312 }
3313
3314 /********************************************************//**
3315 Take a gradient of the squared norm
3316 with respect to its parameters, and add a scaled version
3317 of this gradient to *grad*
3318
3319 \param[in] poly - polynomial
3320 \param[in] scale - scaling for additional gradient
3321 \param[in] decay_type - type of decay
3322 \param[in] decay_param - parameter of decay
3323 \param[in,out] grad - gradient, on output adds scale * new_grad
3324
3325 \return 0 - success, 1 -failure
3326
3327 \note
3328 NEED TO DO SOME TESTS FOR CHEBYSHEV (dont use for now)
3329 ************************************************************/
3330 int
orth_poly_expansion_rkhs_squared_norm_param_grad(const struct OrthPolyExpansion * poly,double scale,enum coeff_decay_type decay_type,double decay_param,double * grad)3331 orth_poly_expansion_rkhs_squared_norm_param_grad(const struct OrthPolyExpansion * poly,
3332 double scale, enum coeff_decay_type decay_type,
3333 double decay_param, double * grad)
3334 {
3335 assert (poly->kristoffel_eval == 0);
3336 if (poly->p->ptype == FOURIER){
3337 fprintf(stderr, "Cannot perform rkhs_squared_norm with fourier basis\n");
3338 exit(1);
3339 }
3340 int res = 1;
3341 if ((poly->p->ptype == LEGENDRE) || (poly->p->ptype == HERMITE)){
3342 if (decay_type == ALGEBRAIC){
3343 for (size_t ii = 0; ii < poly->num_poly; ii++){
3344 grad[ii] += 2.0*scale * poly->coeff[ii] * pow(decay_param,ii);
3345 }
3346 }
3347 else if (decay_type == EXPONENTIAL){
3348 for (size_t ii = 0; ii < poly->num_poly; ii++){
3349 grad[ii] += 2.0*scale * poly->coeff[ii] * pow((double)ii+1.0,-decay_param);
3350 }
3351 }
3352 else{
3353 for (size_t ii = 0; ii < poly->num_poly; ii++){
3354 grad[ii] += 2.0*scale * poly->coeff[ii];
3355 }
3356 }
3357 res = 0;
3358 }
3359 else if (poly->p->ptype == CHEBYSHEV){
3360 // THIS COULD BE WRONG!!
3361 double * temp = calloc_double(poly->num_poly);
3362 for (size_t ii = 0; ii < poly->num_poly; ii++){
3363 temp[ii] = 2.0/(1.0 - (double)2*ii*2*ii);
3364 }
3365 for (size_t ii = 0; ii < poly->num_poly; ii++){
3366 for (size_t jj = 0; jj < poly->num_poly; jj++){
3367 size_t n1 = ii+jj;
3368 size_t n2;
3369 if (ii > jj){
3370 n2 = ii-jj;
3371 }
3372 else{
3373 n2 = jj-ii;
3374 }
3375 if (decay_type == ALGEBRAIC){
3376 if (n1 % 2 == 0){
3377 grad[ii] += 2.0*scale*temp[n1/2]* pow(decay_param,ii);
3378 }
3379 if (n2 % 2 == 0){
3380 grad[ii] += 2.0*scale*temp[n2/2]* pow(decay_param,ii);
3381 }
3382 }
3383 else if (decay_type == EXPONENTIAL){
3384 if (n1 % 2 == 0){
3385 grad[ii] += 2.0*scale*temp[n1/2]*pow((double)ii+1.0,-decay_param);
3386 }
3387 if (n2 % 2 == 0){
3388 grad[ii] += 2.0*scale*temp[n2/2]*pow((double)ii+1.0,-decay_param);
3389 }
3390 }
3391 else {
3392 if (n1 % 2 == 0){
3393 grad[ii] += 2.0*scale*temp[n1/2];
3394 }
3395 if (n2 % 2 == 0){
3396 grad[ii] += 2.0*scale*temp[n2/2];
3397 }
3398 }
3399 }
3400 }
3401 free(temp); temp = NULL;
3402 res = 0;
3403 }
3404 else{
3405 fprintf(stderr, "Cannot evaluate derivative with respect to parameters for polynomial type %d\n",poly->p->ptype);
3406 exit(1);
3407 }
3408 return res;
3409 }
3410
3411 /********************************************************//**
3412 * Round an orthogonal polynomial expansion
3413 *
3414 * \param[in,out] p - orthogonal polynomial expansion
3415 *
3416 * \note
3417 * (UNTESTED, use with care!!!!
3418 *************************************************************/
orth_poly_expansion_round(struct OrthPolyExpansion ** p)3419 void orth_poly_expansion_round(struct OrthPolyExpansion ** p)
3420 {
3421 if ((0 == 0) && ((*p)->p->ptype != FOURIER)){
3422 /* double thresh = 1e-3*ZEROTHRESH; */
3423 double thresh = ZEROTHRESH;
3424 /* double thresh = 1e-30; */
3425 /* double thresh = 10.0*DBL_EPSILON; */
3426 /* printf("thresh = %G\n",thresh); */
3427 size_t jj = 0;
3428 //
3429 int allzero = 1;
3430 double maxcoeff = fabs((*p)->coeff[0]);
3431 for (size_t ii = 1; ii < (*p)->num_poly; ii++){
3432 double val = fabs((*p)->coeff[ii]);
3433 if (val > maxcoeff){
3434 maxcoeff = val;
3435 }
3436 }
3437 maxcoeff = maxcoeff * (*p)->num_poly;
3438 /* printf("maxcoeff = %3.15G\n", maxcoeff); */
3439 for (jj = 0; jj < (*p)->num_poly;jj++){
3440 if (fabs((*p)->coeff[jj]) < thresh){
3441 (*p)->coeff[jj] = 0.0;
3442 }
3443 if (fabs((*p)->coeff[jj])/maxcoeff < thresh){
3444 (*p)->coeff[jj] = 0.0;
3445 }
3446 else{
3447 allzero = 0;
3448 }
3449
3450 }
3451 if (allzero == 1){
3452 (*p)->num_poly = 1;
3453 }
3454 else {
3455 jj = 0;
3456 size_t end = (*p)->num_poly;
3457 if ((*p)->num_poly > 2){
3458 while (fabs((*p)->coeff[end-1]) < thresh){
3459 end-=1;
3460 if (end == 0){
3461 break;
3462 }
3463 }
3464
3465 if (end > 0){
3466 //printf("SHOULD NOT BE HERE\n");
3467 size_t num_poly = end;
3468 //
3469 //double * new_coeff = calloc_double(num_poly);
3470 //for (jj = 0; jj < num_poly; jj++){
3471 // new_coeff[jj] = (*p)->coeff[jj];
3472 // }
3473 //free((*p)->coeff); (*p)->coeff=NULL;
3474 //(*p)->coeff = new_coeff;
3475 (*p)->num_poly = num_poly;
3476 }
3477 }
3478 }
3479
3480 /* printf("rounded coeffs = "); dprint((*p)->num_poly, (*p)->coeff); */
3481
3482 /* orth_poly_expansion_roundt(p,thresh); */
3483
3484 }
3485 }
3486
3487 /********************************************************//**
3488 * Round an orthogonal polynomial expansion to a threshold
3489 *
3490 * \param[in,out] p - orthogonal polynomial expansion
3491 * \param[in] thresh - threshold (relative) to round to
3492 *
3493 * \note
3494 * (UNTESTED, use with care!!!!
3495 *************************************************************/
orth_poly_expansion_roundt(struct OrthPolyExpansion ** p,double thresh)3496 void orth_poly_expansion_roundt(struct OrthPolyExpansion ** p, double thresh)
3497 {
3498
3499 size_t jj = 0;
3500 double sum = 0.0;
3501 /* double maxval = fabs((*p)->coeff[0]); */
3502 /* for (jj = 1; jj < (*p)->num_poly;jj++){ */
3503 /* sum += pow((*p)->coeff[jj],2); */
3504 /* if (fabs((*p)->coeff[jj]) > maxval){ */
3505 /* maxval = fabs((*p)->coeff[jj]); */
3506 /* } */
3507 /* } */
3508 size_t keep = (*p)->num_poly;
3509 if (sum <= thresh){
3510 keep = 1;
3511 }
3512 else{
3513 double sumrun = 0.0;
3514 for (jj = 0; jj < (*p)->num_poly; jj++){
3515 /* if ((fabs((*p)->coeff[jj]) / maxval) < thresh){ */
3516 /* (*p)->coeff[jj] = 0.0; */
3517 /* } */
3518 sumrun += pow((*p)->coeff[jj],2);
3519 if ( (sumrun / sum) > (1.0-thresh)){
3520 keep = jj+1;
3521 break;
3522 }
3523 }
3524 }
3525 /* dprint((*p)->num_poly, (*p)->coeff); */
3526 /* printf("number keep = %zu\n",keep); */
3527 //printf("tolerance = %G\n",thresh);
3528 double * new_coeff = calloc_double(keep + OPECALLOC);
3529 memmove(new_coeff,(*p)->coeff, keep * sizeof(double));
3530 free((*p)->coeff);
3531 (*p)->num_poly = keep;
3532 (*p)->nalloc = (*p)->num_poly + OPECALLOC;
3533 (*p)->coeff = new_coeff;
3534 }
3535
3536
3537
3538 /********************************************************//**
3539 * Approximate a function with an orthogonal polynomial
3540 * series with a fixed number of basis
3541 *
3542 * \param[in] A - function to approximate
3543 * \param[in] args - arguments to function
3544 * \param[in] poly - orthogonal polynomial expansion
3545 *
3546 * \note
3547 * Wont work for polynomial expansion with only the constant
3548 * term.
3549 *************************************************************/
3550 void
orth_poly_expansion_approx(double (* A)(double,void *),void * args,struct OrthPolyExpansion * poly)3551 orth_poly_expansion_approx(double (*A)(double,void *), void *args,
3552 struct OrthPolyExpansion * poly)
3553 {
3554
3555 size_t ii, jj;
3556 double p[2];
3557 double pnew;
3558
3559 /* double m = 1.0; */
3560 /* double off = 0.0; */
3561
3562 double * fvals = NULL;
3563 double * pt_un = NULL; // unormalized point
3564 double * pt = NULL;
3565 double * wt = NULL;
3566
3567 size_t nquad = poly->num_poly+1;
3568
3569 switch (poly->p->ptype) {
3570 case FOURIER:
3571 fprintf(stderr, "Cannot perform orth_poly_expansion_approx with Fourier basis\n");
3572 exit(1);
3573 case CHEBYSHEV:
3574 /* m = (poly->upper_bound - poly->lower_bound) / */
3575 /* (poly->p->upper - poly->p->lower); */
3576 /* off = poly->upper_bound - m * poly->p->upper; */
3577 pt = calloc_double(nquad);
3578 wt = calloc_double(nquad);
3579 cheb_gauss(poly->num_poly,pt,wt);
3580
3581 /* clenshaw_curtis(nquad,pt,wt); */
3582 /* for (ii = 0; ii < nquad; ii++){wt[ii] *= 0.5;} */
3583
3584 break;
3585 case LEGENDRE:
3586 /* m = (poly->upper_bound - poly->lower_bound) / */
3587 /* (poly->p->upper - poly->p->lower); */
3588 /* off = poly->upper_bound - m * poly->p->upper; */
3589 // nquad = poly->num_poly*2.0-1.0;//*2.0;
3590 pt = calloc_double(nquad);
3591 wt = calloc_double(nquad);
3592
3593 // uncomment next two for cc
3594 // clenshaw_curtis(nquad,pt,wt);
3595 // for (ii = 0; ii < nquad; ii++){wt[ii] *= 0.5;}
3596
3597 gauss_legendre(nquad,pt,wt);
3598 break;
3599 case HERMITE:
3600 pt = calloc_double(nquad);
3601 wt = calloc_double(nquad);
3602 gauss_hermite(nquad,pt,wt);
3603 // printf("point = ");
3604 // dprint(nquad,pt);
3605 break;
3606 case STANDARD:
3607 fprintf(stderr, "Cannot call orth_poly_expansion_approx for STANDARD type\n");
3608 break;
3609 //default:
3610 // fprintf(stderr, "Polynomial type does not exist: %d\n ",
3611 // poly->p->ptype);
3612 }
3613
3614 fvals = calloc_double(nquad);
3615 pt_un = calloc_double(nquad);
3616 for (ii = 0; ii < nquad; ii++){
3617 /* pt_un[ii] = m * pt[ii] + off; */
3618 pt_un[ii] = space_mapping_map_inverse(poly->space_transform,pt[ii]);
3619 fvals[ii] = A(pt_un[ii],args) * wt[ii];
3620 }
3621
3622 if (poly->num_poly > 1){
3623 for (ii = 0; ii < nquad; ii++){ // loop over all points
3624 p[0] = poly->p->const_term;
3625 poly->coeff[0] += fvals[ii] * poly->p->const_term;
3626
3627 p[1] = poly->p->lin_const + poly->p->lin_coeff * pt[ii];
3628 poly->coeff[1] += fvals[ii] * p[1] ;
3629 // loop over all coefficients
3630 for (jj = 2; jj < poly->num_poly; jj++){
3631 pnew = eval_orth_poly_wp(poly->p, p[0], p[1], jj, pt[ii]);
3632 poly->coeff[jj] += fvals[ii] * pnew;
3633 p[0] = p[1];
3634 p[1] = pnew;
3635 }
3636 }
3637
3638 for (ii = 0; ii < poly->num_poly; ii++){
3639 poly->coeff[ii] /= poly->p->norm(ii);
3640 }
3641
3642 }
3643 else{
3644 for (ii = 0; ii < nquad; ii++){
3645
3646 poly->coeff[0] += fvals[ii] *poly->p->const_term;
3647 }
3648 poly->coeff[0] /= poly->p->norm(0);
3649 }
3650 free(fvals); fvals = NULL;
3651 free(wt); wt = NULL;
3652 free(pt); pt = NULL;
3653 free(pt_un); pt_un = NULL;
3654
3655 }
3656
3657 /********************************************************//**
3658 * Construct an orthonormal polynomial expansion from (weighted) function
3659 * evaluations and quadrature nodes
3660 *
3661 * \param[in,out] poly - orthogonal polynomial expansion
3662 * \param[in] num_nodes - number of nodes
3663 * \param[in] fvals - function values (multiplied by a weight if necessary)
3664 * \param[in] nodes - locations of evaluations
3665 *************************************************************/
3666 static void
orth_poly_expansion_construct(struct OrthPolyExpansion * poly,size_t num_nodes,double * fvals,double * nodes)3667 orth_poly_expansion_construct(struct OrthPolyExpansion * poly,
3668 size_t num_nodes, double * fvals,
3669 double * nodes)
3670 {
3671 assert (poly->p->ptype != FOURIER);
3672
3673 double p[2];
3674 double pnew;
3675 size_t ii,jj;
3676 if (poly->num_poly > 1){
3677 for (ii = 0; ii < num_nodes; ii++){ // loop over all points
3678 p[0] = poly->p->const_term;
3679 poly->coeff[0] += fvals[ii] * poly->p->const_term;
3680 p[1] = poly->p->lin_const + poly->p->lin_coeff * nodes[ii];
3681 poly->coeff[1] += fvals[ii] * p[1] ;
3682 // loop over all coefficients
3683 for (jj = 2; jj < poly->num_poly; jj++){
3684 pnew = eval_orth_poly_wp(poly->p, p[0], p[1], jj, nodes[ii]);
3685 poly->coeff[jj] += fvals[ii] * pnew;
3686 p[0] = p[1];
3687 p[1] = pnew;
3688 }
3689 }
3690
3691 for (ii = 0; ii < poly->num_poly; ii++){
3692 poly->coeff[ii] /= poly->p->norm(ii);
3693 }
3694
3695 }
3696 else{
3697 for (ii = 0; ii < num_nodes; ii++){
3698 poly->coeff[0] += fvals[ii] * poly->p->const_term;
3699 }
3700 poly->coeff[0] /= poly->p->norm(0);
3701 }
3702 }
3703
3704 /********************************************************//**
3705 * Approximating a function that can take a vector of points as
3706 * input
3707 *
3708 * \param[in,out] poly - orthogonal polynomial expansion
3709 * \param[in] f - wrapped function
3710 * \param[in] opts - approximation options
3711 *
3712 * \return 0 - no problems, > 0 problem
3713 *
3714 * \note Maximum quadrature limited to 200 nodes
3715 *************************************************************/
3716 int
orth_poly_expansion_approx_vec(struct OrthPolyExpansion * poly,struct Fwrap * f,const struct OpeOpts * opts)3717 orth_poly_expansion_approx_vec(struct OrthPolyExpansion * poly,
3718 struct Fwrap * f,
3719 const struct OpeOpts * opts)
3720 {
3721 assert (poly != NULL);
3722 assert (f != NULL);
3723
3724 enum quad_rule qrule = C3_GAUSS_QUAD;
3725 size_t nquad = poly->num_poly+1;
3726 if (opts != NULL){
3727 qrule = opts->qrule;
3728 if (qrule == C3_CC_QUAD){
3729 nquad = 2*poly->num_poly-1;
3730 }
3731 }
3732
3733 if ((nquad < 1 || nquad > 200) && (poly->p->ptype != FOURIER)){
3734 return 1;
3735 }
3736
3737 double fvals[200];
3738 double pt_un[200];
3739 double qpt[200];
3740 double wt[200];
3741
3742 double * quadpt = NULL;
3743 double * quadwt = NULL;
3744
3745 int return_val = 0;
3746 switch (poly->p->ptype) {
3747 case FOURIER:
3748 return fourier_expansion_approx_vec(poly, f, opts);
3749 case CHEBYSHEV:
3750 assert (qrule == C3_GAUSS_QUAD);
3751 return_val = cheb_gauss(nquad,qpt,wt);
3752 break;
3753 case LEGENDRE:
3754 if (qrule == C3_GAUSS_QUAD){
3755 return_val = getLegPtsWts2(nquad,&quadpt,&quadwt);
3756 }
3757 else if (qrule == C3_CC_QUAD){
3758 clenshaw_curtis(nquad,qpt,wt);
3759 for (size_t ii = 0; ii < nquad; ii++){wt[ii] *= 0.5;}
3760 }
3761 else{
3762 fprintf(stderr,"Specified quadrature rule not valid for legendre poly\n");
3763 exit(1);
3764 }
3765 break;
3766 case HERMITE:
3767 assert (qrule == C3_GAUSS_QUAD);
3768 return_val = gauss_hermite(nquad,qpt,wt);
3769 break;
3770 case STANDARD:
3771 fprintf(stderr, "Cannot call orth_poly_expansion_approx_vec");
3772 fprintf(stderr," for STANDARD type\n");
3773 return 1;
3774 }
3775
3776 if (return_val != 0){
3777 return return_val;
3778 }
3779
3780 // something other than legendre
3781 if (quadpt == NULL){
3782 quadpt = qpt;
3783 quadwt = wt;
3784 }
3785
3786 for (size_t ii = 0; ii < nquad; ii++){
3787 pt_un[ii] = space_mapping_map_inverse(poly->space_transform,quadpt[ii]);
3788 }
3789
3790
3791 // Evaluate functions
3792 return_val = fwrap_eval(nquad,pt_un,fvals,f);
3793 if (return_val != 0){
3794 return return_val;
3795 }
3796 for (size_t ii = 0; ii < nquad; ii++){
3797 fvals[ii] *= quadwt[ii];
3798 }
3799
3800 orth_poly_expansion_construct(poly,nquad,fvals,quadpt);
3801
3802 /* printf("constructed\n"); */
3803 /* printf("pts = "); dprint(nquad, quadpt); */
3804 /* printf("vals = "); dprint(nquad, fvals); */
3805 /* for (size_t ii = 0; ii < nquad; ii++){ */
3806 /* double peval = orth_poly_expansion_eval(poly, pt_un[ii]); */
3807
3808 /* printf("peval = %3.5G, fval = %3.5G\n", peval*quadwt[ii], fvals[ii]); */
3809 /* } */
3810
3811 /* exit(1); */
3812 return return_val;
3813 }
3814
3815 /********************************************************//**
3816 * Create an approximation adaptively
3817 *
3818 * \param[in] aopts - approximation options
3819 * \param[in] fw - wrapped function
3820 *
3821 * \return poly
3822 *
3823 * \note
3824 * Follows general scheme that trefethan outlines about
3825 * Chebfun in his book Approximation Theory and practice
3826 *************************************************************/
3827 struct OrthPolyExpansion *
orth_poly_expansion_approx_adapt(const struct OpeOpts * aopts,struct Fwrap * fw)3828 orth_poly_expansion_approx_adapt(const struct OpeOpts * aopts,
3829 struct Fwrap * fw)
3830 {
3831 assert (aopts != NULL);
3832 assert (fw != NULL);
3833 if (aopts->ptype == FOURIER){
3834
3835 size_t N = ope_opts_get_start(aopts);
3836 struct OrthPolyExpansion * poly = orth_poly_expansion_init_from_opts(aopts, N);
3837 orth_poly_expansion_approx_vec(poly, fw, NULL);
3838
3839 return poly;
3840 /* fprintf(stderr, "Cannot perform approx_adapt with fourier basis\n"); */
3841 /* exit(1); */
3842 }
3843
3844 size_t ii;
3845 size_t N = aopts->start_num;
3846 struct OrthPolyExpansion * poly = NULL;
3847 poly = orth_poly_expansion_init_from_opts(aopts,N);
3848 orth_poly_expansion_approx_vec(poly,fw,aopts);
3849
3850 size_t coeffs_too_big = 0;
3851 for (ii = 0; ii < aopts->coeffs_check; ii++){
3852 if (fabs(poly->coeff[N-1-ii]) > aopts->tol){
3853 coeffs_too_big = 1;
3854 break;
3855 }
3856 }
3857
3858
3859
3860 size_t maxnum = ope_opts_get_maxnum(aopts);
3861 /* printf("TOL SPECIFIED IS %G\n",aopts->tol); */
3862 /* printf("Ncoeffs check=%zu \n",aopts->coeffs_check); */
3863 /* printf("maxnum = %zu\n", maxnum); */
3864 while ((coeffs_too_big == 1) && (N < maxnum)){
3865 /* printf("N = %zu\n",N); */
3866 coeffs_too_big = 0;
3867
3868 free(poly->coeff); poly->coeff = NULL;
3869 if (aopts->qrule == C3_CC_QUAD){
3870 N = N * 2 - 1; // for nested cc
3871 }
3872 else{
3873 N = N + 7;
3874 }
3875 poly->num_poly = N;
3876 poly->nalloc = N + OPECALLOC;
3877 poly->coeff = calloc_double(poly->nalloc);
3878 // printf("Number of coefficients to check = %zu\n",aopts->coeffs_check);
3879 orth_poly_expansion_approx_vec(poly,fw,aopts);
3880 double sum_coeff_squared = 0.0;
3881 for (ii = 0; ii < N; ii++){
3882 sum_coeff_squared += pow(poly->coeff[ii],2);
3883 }
3884 sum_coeff_squared = fmax(sum_coeff_squared,ZEROTHRESH);
3885 /* sum_coeff_squared = 1.0; */
3886 for (ii = 0; ii < aopts->coeffs_check; ii++){
3887 /* printf("aopts->tol=%3.15G last coefficients %3.15G\n", */
3888 /* aopts->tol * sum_coeff_squared, */
3889 /* poly->coeff[N-1-ii]); */
3890 if (fabs(poly->coeff[N-1-ii]) > (aopts->tol * sum_coeff_squared) ){
3891 coeffs_too_big = 1;
3892 break;
3893 }
3894 }
3895 if (N > 100){
3896 //printf("Warning: num of poly is %zu: last coeff = %G \n",N,poly->coeff[N-1]);
3897 //printf("tolerance is %3.15G\n", aopts->tol * sum_coeff_squared);
3898 //printf("Considering using piecewise polynomials\n");
3899 /*
3900 printf("first 5 coeffs\n");
3901
3902 size_t ll;
3903 for (ll = 0; ll<5;ll++){
3904 printf("%3.10G ",poly->coeff[ll]);
3905 }
3906 printf("\n");
3907
3908 printf("Last 10 coeffs\n");
3909 for (ll = 0; ll<10;ll++){
3910 printf("%3.10G ",poly->coeff[N-10+ll]);
3911 }
3912 printf("\n");
3913 */
3914 coeffs_too_big = 0;
3915 }
3916
3917 }
3918
3919 orth_poly_expansion_round(&poly);
3920
3921 // verify
3922 /* double pt = (upper - lower)*randu() + lower; */
3923 /* double val_true = A(pt,args); */
3924 /* double val_test = orth_poly_expansion_eval(poly,pt); */
3925 /* double diff = val_true-val_test; */
3926 /* double err = fabs(diff); */
3927 /* if (fabs(val_true) > 1.0){ */
3928 /* //if (fabs(val_true) > ZEROTHRESH){ */
3929 /* err /= fabs(val_true); */
3930 /* } */
3931 /* if (err > 100.0*aopts->tol){ */
3932 /* //fprintf(stderr, "Approximating at point %G in (%3.15G,%3.15G)\n",pt,lower,upper); */
3933 /* //fprintf(stderr, "leads to error %G, while tol is %G \n",err,aopts->tol); */
3934 /* //fprintf(stderr, "actual value is %G \n",val_true); */
3935 /* //fprintf(stderr, "predicted value is %3.15G \n",val_test); */
3936 /* //fprintf(stderr, "%zu N coeffs, last coeffs are %3.15G,%3.15G \n",N,poly->coeff[N-2],poly->coeff[N-1]); */
3937 /* //exit(1); */
3938 /* } */
3939
3940 /* if (default_opts == 1){ */
3941 /* ope_opts_free(aopts); */
3942 /* } */
3943 return poly;
3944 }
3945
3946 /********************************************************//**
3947 * Generate an orthonormal polynomial with pseudorandom coefficients
3948 * between [-1,1]
3949 *
3950 * \param[in] ptype - polynomial type
3951 * \param[in] maxorder - maximum order of the polynomial
3952 * \param[in] lower - lower bound of input
3953 * \param[in] upper - upper bound of input
3954 *
3955 * \return poly
3956 *************************************************************/
3957 struct OrthPolyExpansion *
orth_poly_expansion_randu(enum poly_type ptype,size_t maxorder,double lower,double upper)3958 orth_poly_expansion_randu(enum poly_type ptype, size_t maxorder,
3959 double lower, double upper)
3960 {
3961 struct OrthPolyExpansion * poly =
3962 orth_poly_expansion_init(ptype,maxorder+1, lower, upper);
3963
3964 size_t ii;
3965 for (ii = 0; ii < poly->num_poly; ii++){
3966 poly->coeff[ii] = randu()*2.0-1.0;
3967 }
3968 return poly;
3969 }
3970
3971 /********************************************************//**
3972 * Integrate a Chebyshev approximation
3973 *
3974 * \param[in] poly - polynomial to integrate
3975 *
3976 * \return out - integral of approximation
3977 *************************************************************/
3978 double
cheb_integrate2(const struct OrthPolyExpansion * poly)3979 cheb_integrate2(const struct OrthPolyExpansion * poly)
3980 {
3981 size_t ii;
3982 double out = 0.0;
3983
3984 double m = space_mapping_map_inverse_deriv(poly->space_transform,0.0);
3985 /* double m = */
3986 for (ii = 0; ii < poly->num_poly; ii+=2){
3987 out += poly->coeff[ii] * 2.0 / (1.0 - (double) (ii*ii));
3988 }
3989 out = out * m;
3990 return out;
3991 }
3992
3993 /********************************************************//**
3994 * Integrate a Legendre approximation
3995 *
3996 * \param[in] poly - polynomial to integrate
3997 *
3998 * \return out - integral of approximation
3999 *************************************************************/
4000 double
legendre_integrate(const struct OrthPolyExpansion * poly)4001 legendre_integrate(const struct OrthPolyExpansion * poly)
4002 {
4003 double out = 0.0;
4004
4005 double m = space_mapping_map_inverse_deriv(poly->space_transform,0.0);
4006 out = poly->coeff[0] * 2.0;
4007 out = out * m;
4008 return out;
4009 }
4010
4011 /********************************************************//**
4012 * Compute the product of two polynomial expansion
4013 *
4014 * \param[in] a - first polynomial
4015 * \param[in] b - second polynomial
4016 *
4017 * \return c - polynomial expansion
4018 *
4019 * \note
4020 * Computes c(x) = a(x)b(x) where c is same form as a
4021 * Lower and upper bounds of both polynomials must be the same
4022 *************************************************************/
4023 struct OrthPolyExpansion *
orth_poly_expansion_prod(const struct OrthPolyExpansion * a,const struct OrthPolyExpansion * b)4024 orth_poly_expansion_prod(const struct OrthPolyExpansion * a,
4025 const struct OrthPolyExpansion * b)
4026 {
4027
4028 struct OrthPolyExpansion * c = NULL;
4029 double lb = a->lower_bound;
4030 double ub = a->upper_bound;
4031 assert ( fabs(a->lower_bound - b->lower_bound) < DBL_EPSILON );
4032 assert ( fabs(a->upper_bound - b->upper_bound) < DBL_EPSILON );
4033
4034 enum poly_type p = a->p->ptype;
4035 if (p == FOURIER){
4036 return fourier_expansion_prod(a, b);
4037 }
4038 if ( (p == LEGENDRE) && (a->num_poly < 100) && (b->num_poly < 100)){
4039 //printf("in special prod\n");
4040 //double lb = a->lower_bound;
4041 //double ub = b->upper_bound;
4042
4043 size_t ii,jj;
4044 c = orth_poly_expansion_init(p, a->num_poly + b->num_poly+1, lb, ub);
4045 double * allprods = calloc_double(a->num_poly * b->num_poly);
4046 for (ii = 0; ii < a->num_poly; ii++){
4047 for (jj = 0; jj < b->num_poly; jj++){
4048 allprods[jj + ii * b->num_poly] = a->coeff[ii] * b->coeff[jj];
4049 }
4050 }
4051
4052 //printf("A = \n");
4053 //print_orth_poly_expansion(a,1,NULL);
4054
4055 //printf("B = \n");
4056 //print_orth_poly_expansion(b,1,NULL);
4057
4058 //dprint2d_col(b->num_poly, a->num_poly, allprods);
4059
4060 size_t kk;
4061 for (kk = 0; kk < c->num_poly; kk++){
4062 for (ii = 0; ii < a->num_poly; ii++){
4063 for (jj = 0; jj < b->num_poly; jj++){
4064 c->coeff[kk] += lpolycoeffs[ii+jj*100+kk*10000] *
4065 allprods[jj+ii*b->num_poly];
4066 }
4067 }
4068 //printf("c coeff[%zu]=%G\n",kk,c->coeff[kk]);
4069 }
4070 orth_poly_expansion_round(&c);
4071 free(allprods); allprods=NULL;
4072 }
4073 /* else if (p == CHEBYSHEV){ */
4074 /* c = orth_poly_expansion_init(p,a->num_poly+b->num_poly+1,lb,ub); */
4075 /* for (size_t ii = 0; ii) */
4076 /* } */
4077 else{
4078 // printf("OrthPolyExpansion product greater than order 100 is slow\n");
4079 const struct OrthPolyExpansion * comb[2];
4080 comb[0] = a;
4081 comb[1] = b;
4082
4083 double norma = 0.0, normb = 0.0;
4084 size_t ii;
4085 for (ii = 0; ii < a->num_poly; ii++){
4086 norma += pow(a->coeff[ii],2);
4087 }
4088 for (ii = 0; ii < b->num_poly; ii++){
4089 normb += pow(b->coeff[ii],2);
4090 }
4091
4092 if ( (norma < ZEROTHRESH) || (normb < ZEROTHRESH) ){
4093 //printf("in here \n");
4094 // c = orth_poly_expansion_constant(0.0,a->p->ptype,lb,ub);
4095 c = orth_poly_expansion_init(p,1, lb, ub);
4096 space_mapping_free(c->space_transform); c->space_transform = NULL;
4097 c->space_transform = space_mapping_copy(a->space_transform);
4098 c->coeff[0] = 0.0;
4099 }
4100 else{
4101 //printf(" total order of product = %zu\n",a->num_poly+b->num_poly);
4102 c = orth_poly_expansion_init(p, a->num_poly + b->num_poly+5, lb, ub);
4103 space_mapping_free(c->space_transform); c->space_transform = NULL;
4104 c->space_transform = space_mapping_copy(a->space_transform);
4105 orth_poly_expansion_approx(&orth_poly_expansion_eval3,comb,c);
4106 /* printf("num_coeff pre_round = %zu\n", c->num_poly); */
4107 /* orth_poly_expansion_approx(&orth_poly_expansion_eval3,comb,c); */
4108 orth_poly_expansion_round(&c);
4109 /* printf("num_coeff post_round = %zu\n", c->num_poly); */
4110 }
4111 }
4112
4113 //*
4114 //printf("compute product\n");
4115 //struct OpeOpts ao;
4116 //ao.start_num = 3;
4117 //ao.coeffs_check = 2;
4118 //ao.tol = 1e-13;
4119 //c = orth_poly_expansion_approx_adapt(&orth_poly_expansion_eval3,comb,
4120 // p, lb, ub, &ao);
4121
4122 //orth_poly_expansion_round(&c);
4123 //printf("done\n");
4124 //*/
4125 return c;
4126 }
4127
4128 /********************************************************//**
4129 * Compute the sum of the product between the functions in two arraysarrays
4130 *
4131 * \param[in] n - number of functions
4132 * \param[in] lda - stride of first array
4133 * \param[in] a - array of orthonormal polynomial expansions
4134 * \param[in] ldb - stride of second array
4135 * \param[in] b - array of orthonormal polynomial expansions
4136 *
4137 * \return c - polynomial expansion
4138 *
4139 * \note
4140 * All the functions need to have the same lower
4141 * and upper bounds and be of the same type
4142 *
4143 * If the maximum order of the polynomials is greater than 25 then this is
4144 * inefficient because I haven't precomputed triple product rules
4145 *************************************************************/
4146 struct OrthPolyExpansion *
orth_poly_expansion_sum_prod(size_t n,size_t lda,struct OrthPolyExpansion ** a,size_t ldb,struct OrthPolyExpansion ** b)4147 orth_poly_expansion_sum_prod(size_t n, size_t lda,
4148 struct OrthPolyExpansion ** a, size_t ldb,
4149 struct OrthPolyExpansion ** b)
4150 {
4151
4152 enum poly_type ptype;
4153 ptype = a[0]->p->ptype;
4154 struct OrthPolyExpansion * c = NULL;
4155 double lb = a[0]->lower_bound;
4156 double ub = a[0]->upper_bound;
4157
4158 size_t ii;
4159 size_t maxordera = 0;
4160 size_t maxorderb = 0;
4161 size_t maxorder = 0;
4162 //int legen = 1;
4163 for (ii = 0; ii < n; ii++){
4164
4165 if (a[ii*lda]->p->ptype != b[ii*ldb]->p->ptype){
4166 return c; // cant do it
4167 }
4168 else if (a[ii*lda]->p->ptype != ptype){
4169 return c;
4170 }
4171 size_t neworder = a[ii*lda]->num_poly + b[ii*ldb]->num_poly;
4172 if (neworder > maxorder){
4173 maxorder = neworder;
4174 }
4175 if (a[ii*lda]->num_poly > maxordera){
4176 maxordera = a[ii*lda]->num_poly;
4177 }
4178 if (b[ii*ldb]->num_poly > maxorderb){
4179 maxorderb = b[ii*ldb]->num_poly;
4180 }
4181 }
4182 if ( (maxordera > 99) || (maxorderb > 99) || (ptype != LEGENDRE)){
4183 printf("OrthPolyExpansion sum_product greater than order 100 is slow\n");
4184 c = orth_poly_expansion_prod(a[0],b[0]);
4185 for (ii = 1; ii< n; ii++){
4186 struct OrthPolyExpansion * temp =
4187 orth_poly_expansion_prod(a[ii*lda],b[ii*ldb]);
4188 orth_poly_expansion_axpy(1.0,temp,c);
4189 orth_poly_expansion_free(temp);
4190 temp = NULL;
4191 }
4192 }
4193 else{
4194 enum poly_type p = LEGENDRE;
4195 c = orth_poly_expansion_init(p, maxorder, lb, ub);
4196 size_t kk,jj,ll;
4197 double * allprods = calloc_double( maxorderb * maxordera);
4198 for (kk = 0; kk < n; kk++){
4199 for (ii = 0; ii < a[kk*lda]->num_poly; ii++){
4200 for (jj = 0; jj < b[kk*ldb]->num_poly; jj++){
4201 allprods[jj + ii * maxorderb] +=
4202 a[kk*lda]->coeff[ii] * b[kk*ldb]->coeff[jj];
4203 }
4204 }
4205 }
4206
4207 for (ll = 0; ll < c->num_poly; ll++){
4208 for (ii = 0; ii < maxordera; ii++){
4209 for (jj = 0; jj < maxorderb; jj++){
4210 c->coeff[ll] += lpolycoeffs[ii+jj*100+ll*10000] *
4211 allprods[jj+ii*maxorderb];
4212 }
4213 }
4214 }
4215 free(allprods); allprods=NULL;
4216 orth_poly_expansion_round(&c);
4217 }
4218 return c;
4219 }
4220
4221 /********************************************************//**
4222 * Compute a linear combination of generic functions
4223 *
4224 * \param[in] n - number of functions
4225 * \param[in] ldx - stride of first array
4226 * \param[in] x - functions
4227 * \param[in] ldc - stride of coefficients
4228 * \param[in] c - scaling coefficients
4229 *
4230 * \return out = \f$ \sum_{i=1}^n coeff[ldc[i]] * gfa[ldgf[i]] \f$
4231 *
4232 * \note
4233 * If both arrays do not consist of only LEGENDRE polynomials
4234 * return NULL. All the functions need to have the same lower
4235 * and upper bounds
4236 *************************************************************/
4237 struct OrthPolyExpansion *
orth_poly_expansion_lin_comb(size_t n,size_t ldx,struct OrthPolyExpansion ** x,size_t ldc,const double * c)4238 orth_poly_expansion_lin_comb(size_t n, size_t ldx,
4239 struct OrthPolyExpansion ** x, size_t ldc,
4240 const double * c )
4241 {
4242
4243 struct OrthPolyExpansion * out = NULL;
4244 double lb = x[0]->lower_bound;
4245 double ub = x[0]->upper_bound;
4246 enum poly_type ptype = x[0]->p->ptype;
4247 size_t ii;
4248 size_t maxorder = 0;
4249 //int legen = 1;
4250 for (ii = 0; ii < n; ii++){
4251 if (x[ii*ldx]->p->ptype != ptype){
4252 //legen = 0;
4253 return out; // cant do it
4254 }
4255 size_t neworder = x[ii*ldx]->num_poly;
4256 if (neworder > maxorder){
4257 maxorder = neworder;
4258 }
4259 }
4260
4261
4262 out = orth_poly_expansion_init(ptype, maxorder, lb, ub);
4263 space_mapping_free(out->space_transform);
4264 out->space_transform = space_mapping_copy(x[0]->space_transform);
4265 size_t kk;
4266 if (ptype != FOURIER){
4267 for (kk = 0; kk < n; kk++){
4268 for (ii = 0; ii < x[kk*ldx]->num_poly; ii++){
4269 out->coeff[ii] += c[kk*ldc]*x[kk*ldx]->coeff[ii];
4270 }
4271 }
4272 }
4273 else{
4274 for (kk = 0; kk < n; kk++){
4275 for (ii = 0; ii < x[kk*ldx]->num_poly; ii++){
4276 out->ccoeff[ii] += c[kk*ldc]*x[kk*ldx]->ccoeff[ii];
4277 }
4278 }
4279 }
4280 orth_poly_expansion_round(&out);
4281 return out;
4282 }
4283
4284 /********************************************************//**
4285 * Integrate an orthogonal polynomial expansion
4286 *
4287 * \param[in] poly - polynomial to integrate
4288 *
4289 * \return out - Integral of approximation
4290 *
4291 * \note
4292 * Need to an 'else' or default behavior to switch case
4293 * int_{lb}^ub f(x) dx
4294 * For Hermite polynomials this integrates with respec to
4295 * the Gaussian weight
4296 *************************************************************/
4297 double
orth_poly_expansion_integrate(const struct OrthPolyExpansion * poly)4298 orth_poly_expansion_integrate(const struct OrthPolyExpansion * poly)
4299 {
4300 double out = 0.0;
4301 switch (poly->p->ptype){
4302 case LEGENDRE: out = legendre_integrate(poly); break;
4303 case HERMITE: out = hermite_integrate(poly); break;
4304 case CHEBYSHEV: out = cheb_integrate2(poly); break;
4305 case FOURIER: out = fourier_integrate(poly); break;
4306 case STANDARD: fprintf(stderr, "Cannot integrate STANDARD type\n"); break;
4307 }
4308 return out;
4309 }
4310
4311 /********************************************************//**
4312 * Integrate an orthogonal polynomial expansion
4313 *
4314 * \param[in] poly - polynomial to integrate
4315 *
4316 * \return out - Integral of approximation
4317 *
4318 \note Computes \f$ \int f(x) w(x) dx \f$ for every univariate function
4319 in the qmarray
4320
4321 w(x) depends on underlying parameterization
4322 for example, it is 1/2 for legendre (and default for others),
4323 gauss for hermite,etc
4324 *************************************************************/
4325 double
orth_poly_expansion_integrate_weighted(const struct OrthPolyExpansion * poly)4326 orth_poly_expansion_integrate_weighted(const struct OrthPolyExpansion * poly)
4327 {
4328 double out = 0.0;
4329 switch (poly->p->ptype){
4330 case LEGENDRE: out = poly->coeff[0]; break;
4331 case HERMITE: out = poly->coeff[0]; break;
4332 case CHEBYSHEV: out = poly->coeff[0]; break;
4333 case FOURIER: out = creal(poly->ccoeff[0]); break;
4334 case STANDARD: fprintf(stderr, "Cannot integrate STANDARD type\n"); break;
4335 }
4336
4337 return out;
4338 }
4339
4340
4341 /********************************************************//**
4342 * Weighted inner product between two polynomial
4343 * expansions of the same type
4344 *
4345 * \param[in] a - first polynomial
4346 * \param[in] b - second polynomai
4347 *
4348 * \return inner product
4349 *
4350 * \note
4351 * Computes \f[ \int_{lb}^ub a(x)b(x) w(x) dx \f]
4352 *
4353 *************************************************************/
4354 double
orth_poly_expansion_inner_w(const struct OrthPolyExpansion * a,const struct OrthPolyExpansion * b)4355 orth_poly_expansion_inner_w(const struct OrthPolyExpansion * a,
4356 const struct OrthPolyExpansion * b)
4357 {
4358 assert(a->p->ptype == b->p->ptype);
4359 assert (a->p->ptype != FOURIER);
4360 double out = 0.0;
4361 size_t N = a->num_poly < b->num_poly ? a->num_poly : b->num_poly;
4362 size_t ii;
4363 for (ii = 0; ii < N; ii++){
4364 out += a->coeff[ii] * b->coeff[ii] * a->p->norm(ii);
4365 }
4366
4367 return out;
4368 }
4369
4370 /********************************************************//**
4371 * Inner product between two polynomial expansions of the same type
4372 *
4373 * \param[in] a - first polynomial
4374 * \param[in] b - second polynomai
4375 *
4376 * \return inner product
4377 *
4378 * \note
4379 * If the polynomial is NOT HERMITE then
4380 * Computes \f$ \int_{lb}^ub a(x)b(x) dx \f$ by first
4381 * converting each polynomial to a Legendre polynomial
4382 * Otherwise it computes the error with respect to gaussia weight
4383 *************************************************************/
4384 double
orth_poly_expansion_inner(const struct OrthPolyExpansion * a,const struct OrthPolyExpansion * b)4385 orth_poly_expansion_inner(const struct OrthPolyExpansion * a,
4386 const struct OrthPolyExpansion * b)
4387 {
4388 struct OrthPolyExpansion * t1 = NULL;
4389 struct OrthPolyExpansion * t2 = NULL;
4390
4391 // assert (a->ptype == b->ptype);
4392 // enum poly_type ptype = a->ptype;
4393 // switch (ptype)
4394 if ((a->p->ptype == HERMITE) && (b->p->ptype == HERMITE)){
4395 return orth_poly_expansion_inner_w(a,b);
4396 }
4397 else if ((a->p->ptype == CHEBYSHEV) && (b->p->ptype == CHEBYSHEV)){
4398
4399 struct OrthPolyExpansion * prod = orth_poly_expansion_prod(a, b);
4400 double int_val = orth_poly_expansion_integrate(prod);
4401 orth_poly_expansion_free(prod); prod = NULL;
4402 return int_val;
4403 /* // can possibly make this more efficient */
4404 /* double out = 0.0; */
4405 /* size_t N = a->num_poly < b->num_poly ? a->num_poly : b->num_poly; */
4406 /* for (size_t ii = 0; ii < N; ii++){ */
4407 /* for (size_t jj = 0; jj < ii; jj++){ */
4408 /* if ( ((ii+jj) % 2) == 0){ */
4409 /* out += (a->coeff[ii]*b->coeff[jj] * */
4410 /* (1.0 / (1.0 - (double) (ii-jj)*(ii-jj)) */
4411 /* + 1.0 / (1.0 - (double) (ii+jj)*(ii+jj)))); */
4412 /* } */
4413 /* } */
4414 /* for (size_t jj = ii; jj < N; jj++){ */
4415 /* if ( ((ii+jj) % 2) == 0){ */
4416 /* out += (a->coeff[ii]*b->coeff[jj] * */
4417 /* (1.0 / (1.0 - (double) (jj-ii)*(jj-ii)) */
4418 /* + 1.0 / (1.0 - (double) (ii+jj)*(ii+jj)))); */
4419 /* } */
4420 /* } */
4421 /* } */
4422 /* double m = (a->upper_bound - a->lower_bound) / (a->p->upper - a->p->lower); */
4423 /* out *= m; */
4424 /* return out */;
4425 }
4426 else if ((a->p->ptype == FOURIER) && (b->p->ptype == FOURIER)){
4427 struct OrthPolyExpansion * prod = orth_poly_expansion_prod(a, b);
4428 double int_val = orth_poly_expansion_integrate(prod);
4429 orth_poly_expansion_free(prod); prod = NULL;
4430 return int_val;
4431 }
4432 else{
4433 if (a->p->ptype == CHEBYSHEV){
4434 t1 = orth_poly_expansion_init(LEGENDRE, a->num_poly,
4435 a->lower_bound, a->upper_bound);
4436 orth_poly_expansion_approx(&orth_poly_expansion_eval2, (void*)a, t1);
4437 orth_poly_expansion_round(&t1);
4438 }
4439 else if (a->p->ptype != LEGENDRE){
4440 fprintf(stderr, "Don't know how to take inner product using polynomial type. \n");
4441 fprintf(stderr, "type1 = %d, and type2= %d\n",a->p->ptype,b->p->ptype);
4442 exit(1);
4443 }
4444
4445 if (b->p->ptype == CHEBYSHEV){
4446 t2 = orth_poly_expansion_init(LEGENDRE, b->num_poly,
4447 b->lower_bound, b->upper_bound);
4448 orth_poly_expansion_approx(&orth_poly_expansion_eval2, (void*)b, t2);
4449 orth_poly_expansion_round(&t2);
4450 }
4451 else if (b->p->ptype != LEGENDRE){
4452 fprintf(stderr, "Don't know how to take inner product using polynomial type. \n");
4453 fprintf(stderr, "type1 = %d, and type2= %d\n",a->p->ptype,b->p->ptype);
4454 exit(1);
4455 }
4456
4457 double out;
4458 if ((t1 == NULL) && (t2 == NULL)){
4459 out = orth_poly_expansion_inner_w(a,b) * (a->upper_bound - a->lower_bound);
4460 }
4461 else if ((t1 == NULL) && (t2 != NULL)){
4462 out = orth_poly_expansion_inner_w(a,t2) * (a->upper_bound - a->lower_bound);
4463 orth_poly_expansion_free(t2); t2 = NULL;
4464 }
4465 else if ((t2 == NULL) && (t1 != NULL)){
4466 out = orth_poly_expansion_inner_w(t1,b) * (b->upper_bound - b->lower_bound);
4467 orth_poly_expansion_free(t1); t1 = NULL;
4468 }
4469 else{
4470 out = orth_poly_expansion_inner_w(t1,t2) * (t1->upper_bound - t1->lower_bound);
4471 orth_poly_expansion_free(t1); t1 = NULL;
4472 orth_poly_expansion_free(t2); t2 = NULL;
4473 }
4474 return out;
4475 }
4476 }
4477
4478 /********************************************************//**
4479 * Compute the norm of an orthogonal polynomial
4480 * expansion with respect to family weighting
4481 * function
4482 *
4483 * \param[in] p - polynomial to integrate
4484 *
4485 * \return out - norm of function
4486 *
4487 * \note
4488 * Computes \f$ \sqrt(\int f(x)^2 w(x) dx) \f$
4489 *************************************************************/
orth_poly_expansion_norm_w(const struct OrthPolyExpansion * p)4490 double orth_poly_expansion_norm_w(const struct OrthPolyExpansion * p){
4491
4492 double out = sqrt(orth_poly_expansion_inner_w(p,p));
4493 return sqrt(out);
4494 }
4495
4496 /********************************************************//**
4497 * Compute the norm of an orthogonal polynomial
4498 * expansion with respect to uniform weighting
4499 * (except in case of HERMITE, then do gaussian weighting)
4500 *
4501 * \param[in] p - polynomial of which to obtain norm
4502 *
4503 * \return out - norm of function
4504 *
4505 * \note
4506 * Computes \f$ \sqrt(\int_a^b f(x)^2 dx) \f$
4507 *************************************************************/
orth_poly_expansion_norm(const struct OrthPolyExpansion * p)4508 double orth_poly_expansion_norm(const struct OrthPolyExpansion * p){
4509
4510 double out = 0.0;
4511 out = sqrt(orth_poly_expansion_inner(p,p));
4512 return out;
4513 }
4514
4515 /********************************************************//**
4516 * Multiply polynomial expansion by -1
4517 *
4518 * \param[in,out] p - polynomial multiply by -1
4519 *************************************************************/
4520 void
orth_poly_expansion_flip_sign(struct OrthPolyExpansion * p)4521 orth_poly_expansion_flip_sign(struct OrthPolyExpansion * p)
4522 {
4523
4524 if (p->p->ptype != FOURIER){
4525 for (size_t ii = 0; ii < p->num_poly; ii++){
4526 p->coeff[ii] *= -1.0;
4527 }
4528 }
4529 else{
4530 for (size_t ii = 0; ii < p->num_poly; ii++){
4531 p->ccoeff[ii] *= -1.0;
4532 }
4533 }
4534 }
4535
4536 /********************************************************//**
4537 * Multiply by scalar and overwrite expansion
4538 *
4539 * \param[in] a - scaling factor for first polynomial
4540 * \param[in] x - polynomial to scale
4541 *************************************************************/
orth_poly_expansion_scale(double a,struct OrthPolyExpansion * x)4542 void orth_poly_expansion_scale(double a, struct OrthPolyExpansion * x)
4543 {
4544
4545 size_t ii;
4546 if (x->p->ptype != FOURIER){
4547 for (ii = 0; ii < x->num_poly; ii++){
4548 x->coeff[ii] *= a;
4549 }
4550 }
4551 else{
4552 for (ii = 0; ii < x->num_poly; ii++){
4553 x->ccoeff[ii] *= a;
4554 }
4555 }
4556 orth_poly_expansion_round(&x);
4557 }
4558
4559 /********************************************************//**
4560 * Multiply and add 3 expansions \f$ z \leftarrow ax + by + cz \f$
4561 *
4562 * \param[in] a - scaling factor for first polynomial
4563 * \param[in] x - first polynomial
4564 * \param[in] b - scaling factor for second polynomial
4565 * \param[in] y - second polynomial
4566 * \param[in] c - scaling factor for third polynomial
4567 * \param[in] z - third polynomial
4568 *
4569 *************************************************************/
4570 void
orth_poly_expansion_sum3_up(double a,struct OrthPolyExpansion * x,double b,struct OrthPolyExpansion * y,double c,struct OrthPolyExpansion * z)4571 orth_poly_expansion_sum3_up(double a, struct OrthPolyExpansion * x,
4572 double b, struct OrthPolyExpansion * y,
4573 double c, struct OrthPolyExpansion * z)
4574 {
4575 assert (x->p->ptype == y->p->ptype);
4576 assert (y->p->ptype == z->p->ptype);
4577 assert ( x != NULL );
4578 assert ( y != NULL );
4579 assert ( z != NULL );
4580
4581 assert (x->p->ptype != FOURIER);
4582
4583 size_t ii;
4584 if ( (z->num_poly >= x->num_poly) && (z->num_poly >= y->num_poly) ){
4585
4586 if (x->num_poly > y->num_poly){
4587 for (ii = 0; ii < y->num_poly; ii++){
4588 z->coeff[ii] = c*z->coeff[ii] + a*x->coeff[ii] + b*y->coeff[ii];
4589 }
4590 for (ii = y->num_poly; ii < x->num_poly; ii++){
4591 z->coeff[ii] = c*z->coeff[ii] + a*x->coeff[ii];
4592 }
4593 for (ii = x->num_poly; ii < z->num_poly; ii++){
4594 z->coeff[ii] = c*z->coeff[ii];
4595 }
4596 }
4597 else{
4598 for (ii = 0; ii < x->num_poly; ii++){
4599 z->coeff[ii] = c*z->coeff[ii] + a*x->coeff[ii] + b*y->coeff[ii];
4600 }
4601 for (ii = x->num_poly; ii < y->num_poly; ii++){
4602 z->coeff[ii] = c*z->coeff[ii] + b*y->coeff[ii];
4603 }
4604 for (ii = x->num_poly; ii < z->num_poly; ii++){
4605 z->coeff[ii] = c*z->coeff[ii];
4606 }
4607 }
4608 }
4609 else if ((z->num_poly >= x->num_poly) && ( z->num_poly < y->num_poly)) {
4610 double * temp = realloc(z->coeff, (y->num_poly)*sizeof(double));
4611 if (temp == NULL){
4612 fprintf(stderr,"cannot allocate new size fo z-coeff in sum3_up\n");
4613 exit(1);
4614 }
4615 else{
4616 z->coeff = temp;
4617 }
4618 for (ii = 0; ii < x->num_poly; ii++){
4619 z->coeff[ii] = c*z->coeff[ii]+a*x->coeff[ii]+b*y->coeff[ii];
4620 }
4621 for (ii = x->num_poly; ii < z->num_poly; ii++){
4622 z->coeff[ii] = c*z->coeff[ii] + b*y->coeff[ii];
4623 }
4624 for (ii = z->num_poly; ii < y->num_poly; ii++){
4625 z->coeff[ii] = b*y->coeff[ii];
4626 }
4627 z->num_poly = y->num_poly;
4628 }
4629 else if ( (z->num_poly < x->num_poly) && ( z->num_poly >= y->num_poly) ){
4630 double * temp = realloc(z->coeff, (x->num_poly)*sizeof(double));
4631 if (temp == NULL){
4632 fprintf(stderr,"cannot allocate new size fo z-coeff in sum3_up\n");
4633 exit(1);
4634 }
4635 else{
4636 z->coeff = temp;
4637 }
4638 for (ii = 0; ii < y->num_poly; ii++){
4639 z->coeff[ii] = c*z->coeff[ii]+a*x->coeff[ii]+b*y->coeff[ii];
4640 }
4641 for (ii = y->num_poly; ii < z->num_poly; ii++){
4642 z->coeff[ii] = c*z->coeff[ii] + a*x->coeff[ii];
4643 }
4644 for (ii = z->num_poly; ii < x->num_poly; ii++){
4645 z->coeff[ii] = a*x->coeff[ii];
4646 }
4647 z->num_poly = x->num_poly;
4648 }
4649 else if ( x->num_poly <= y->num_poly){
4650 double * temp = realloc(z->coeff, (y->num_poly)*sizeof(double));
4651 if (temp == NULL){
4652 fprintf(stderr,"cannot allocate new size fo z-coeff in sum3_up\n");
4653 exit(1);
4654 }
4655 for (ii = 0; ii < z->num_poly; ii++){
4656 z->coeff[ii] = c*z->coeff[ii]+a*x->coeff[ii]+b*y->coeff[ii];
4657 }
4658 for (ii = z->num_poly; ii < x->num_poly; ii++){
4659 z->coeff[ii] = a*x->coeff[ii] + b*y->coeff[ii];
4660 }
4661 for (ii = x->num_poly; ii < y->num_poly; ii++){
4662 z->coeff[ii] = b*y->coeff[ii];
4663 }
4664 z->num_poly = y->num_poly;
4665 }
4666 else if (y->num_poly <= x->num_poly) {
4667 double * temp = realloc(z->coeff, (x->num_poly)*sizeof(double));
4668 if (temp == NULL){
4669 fprintf(stderr,"cannot allocate new size fo z-coeff in sum3_up\n");
4670 exit(1);
4671 }
4672 for (ii = 0; ii < z->num_poly; ii++){
4673 z->coeff[ii] = c*z->coeff[ii]+a*x->coeff[ii]+b*y->coeff[ii];
4674 }
4675 for (ii = z->num_poly; ii < y->num_poly; ii++){
4676 z->coeff[ii] = a*x->coeff[ii] + b*y->coeff[ii];
4677 }
4678 for (ii = y->num_poly; ii < x->num_poly; ii++){
4679 z->coeff[ii] = a*x->coeff[ii];
4680 }
4681 z->num_poly = x->num_poly;
4682 }
4683 else{
4684 fprintf(stderr,"Haven't accounted for anything else?! %zu %zu %zu\n",
4685 x->num_poly, y->num_poly, z->num_poly);
4686 exit(1);
4687 }
4688 // z->nalloc = z->num_poly + OPECALLOC;
4689 orth_poly_expansion_round(&z);
4690 }
4691
4692 /********************************************************//**
4693 * Multiply by scalar and add two orthgonal
4694 * expansions of the same family together \f[ y \leftarrow ax + y \f]
4695 *
4696 * \param[in] a - scaling factor for first polynomial
4697 * \param[in] x - first polynomial
4698 * \param[in] y - second polynomial
4699 *
4700 * \return 0 if successfull 1 if error with allocating more space for y
4701 *
4702 * \note
4703 * Computes z=ax+by, where x and y are polynomial expansionx
4704 * Requires both polynomials to have the same upper
4705 * and lower bounds
4706 *
4707 **************************************************************/
orth_poly_expansion_axpy(double a,struct OrthPolyExpansion * x,struct OrthPolyExpansion * y)4708 int orth_poly_expansion_axpy(double a, struct OrthPolyExpansion * x,
4709 struct OrthPolyExpansion * y)
4710 {
4711 assert (y != NULL);
4712 assert (x != NULL);
4713 assert (x->p->ptype == y->p->ptype);
4714
4715 assert ( fabs(x->lower_bound - y->lower_bound) < DBL_EPSILON );
4716 assert ( fabs(x->upper_bound - y->upper_bound) < DBL_EPSILON );
4717
4718 if (x->p->ptype == FOURIER){
4719 return fourier_expansion_axpy(a, x, y);
4720 }
4721
4722 if (x->num_poly < y->num_poly){
4723 // shouldnt need rounding here
4724 size_t ii;
4725 for (ii = 0; ii < x->num_poly; ii++){
4726 y->coeff[ii] += a * x->coeff[ii];
4727 if (fabs(y->coeff[ii]) < ZEROTHRESH){
4728 y->coeff[ii] = 0.0;
4729 }
4730 }
4731 }
4732 else{
4733 size_t ii;
4734 if (x->num_poly > y->nalloc){
4735 //printf("hereee\n");
4736 y->nalloc = x->num_poly+10;
4737 double * temp = realloc(y->coeff, (y->nalloc)*sizeof(double));
4738 if (temp == NULL){
4739 return 0;
4740 }
4741 else{
4742 y->coeff = temp;
4743 for (ii = y->num_poly; ii < y->nalloc; ii++){
4744 y->coeff[ii] = 0.0;
4745 }
4746 }
4747 //printf("finished\n");
4748 }
4749 for (ii = y->num_poly; ii < x->num_poly; ii++){
4750 y->coeff[ii] = a * x->coeff[ii];
4751 if (fabs(y->coeff[ii]) < ZEROTHRESH){
4752 y->coeff[ii] = 0.0;
4753 }
4754 }
4755 for (ii = 0; ii < y->num_poly; ii++){
4756 y->coeff[ii] += a * x->coeff[ii];
4757 if (fabs(y->coeff[ii]) < ZEROTHRESH){
4758 y->coeff[ii] = 0.0;
4759 }
4760 }
4761 y->num_poly = x->num_poly;
4762 size_t nround = y->num_poly;
4763 for (ii = 0; ii < y->num_poly-1;ii++){
4764 if (fabs(y->coeff[y->num_poly-1-ii]) > ZEROTHRESH){
4765 break;
4766 }
4767 else{
4768 nround = nround-1;
4769 }
4770 }
4771 y->num_poly = nround;
4772 }
4773
4774 return 0;
4775 }
4776
4777
4778 /********************************************************//**
4779 * Multiply by scalar and add two orthgonal
4780 * expansions of the same family together
4781 *
4782 * \param[in] a - scaling factor for first polynomial
4783 * \param[in] x - first polynomial
4784 * \param[in] b - scaling factor for second polynomial
4785 * \param[in] y second polynomial
4786 *
4787 * \return p - orthogonal polynomial expansion
4788 *
4789 * \note
4790 * Computes z=ax+by, where x and y are polynomial expansionx
4791 * Requires both polynomials to have the same upper
4792 * and lower bounds
4793 *
4794 *************************************************************/
4795 struct OrthPolyExpansion *
orth_poly_expansion_daxpby(double a,struct OrthPolyExpansion * x,double b,struct OrthPolyExpansion * y)4796 orth_poly_expansion_daxpby(double a, struct OrthPolyExpansion * x,
4797 double b, struct OrthPolyExpansion * y)
4798 {
4799 /* assert (x->p->ptype != FOURIER); */
4800 struct OrthPolyExpansion * p ;
4801 if (x != NULL && y != NULL){
4802 p = orth_poly_expansion_copy(y);
4803 orth_poly_expansion_scale(b, p);
4804 orth_poly_expansion_axpy(a, x, p);
4805 }
4806 else if (x != NULL){
4807 p = orth_poly_expansion_copy(x);
4808 orth_poly_expansion_scale(a, p);
4809 }
4810 else if (y != NULL){
4811 p = orth_poly_expansion_copy(y);
4812 orth_poly_expansion_scale(b, p);
4813 }
4814 else{
4815 fprintf(stderr, "Error in orth_poly_expansion_daxpby\n");
4816 fprintf(stderr, "trying to run with all NULL arguments\n");
4817 exit(1);
4818 }
4819 return p;
4820 }
4821
4822 ////////////////////////////////////////////////////////////////////////////
4823 // Algorithms
4824
4825 /********************************************************//**
4826 * Obtain the real roots of a standard polynomial
4827 *
4828 * \param[in] p - standard polynomial
4829 * \param[in,out] nkeep - returns how many real roots tehre are
4830 *
4831 * \return real_roots - real roots of a standard polynomial
4832 *
4833 * \note
4834 * Only roots within the bounds are returned
4835 *************************************************************/
4836 double *
standard_poly_real_roots(struct StandardPoly * p,size_t * nkeep)4837 standard_poly_real_roots(struct StandardPoly * p, size_t * nkeep)
4838 {
4839 if (p->num_poly == 1) // constant function
4840 {
4841 double * real_roots = NULL;
4842 *nkeep = 0;
4843 return real_roots;
4844 }
4845 else if (p->num_poly == 2){ // linear
4846 double root = -p->coeff[0] / p->coeff[1];
4847
4848 if ((root > p->lower_bound) && (root < p->upper_bound)){
4849 *nkeep = 1;
4850 }
4851 else{
4852 *nkeep = 0;
4853 }
4854 double * real_roots = NULL;
4855 if (*nkeep == 1){
4856 real_roots = calloc_double(1);
4857 real_roots[0] = root;
4858 }
4859 return real_roots;
4860 }
4861
4862 size_t nrows = p->num_poly-1;
4863 //printf("coeffs = \n");
4864 //dprint(p->num_poly, p->coeff);
4865 while (fabs(p->coeff[nrows]) < ZEROTHRESH ){
4866 //while (fabs(p->coeff[nrows]) < DBL_MIN){
4867 nrows--;
4868 if (nrows == 1){
4869 break;
4870 }
4871 }
4872
4873 //printf("nrows left = %zu \n", nrows);
4874 if (nrows == 1) // linear
4875 {
4876 double root = -p->coeff[0] / p->coeff[1];
4877 if ((root > p->lower_bound) && (root < p->upper_bound)){
4878 *nkeep = 1;
4879 }
4880 else{
4881 *nkeep = 0;
4882 }
4883 double * real_roots = NULL;
4884 if (*nkeep == 1){
4885 real_roots = calloc_double(1);
4886 real_roots[0] = root;
4887 }
4888 return real_roots;
4889 }
4890 else if (nrows == 0)
4891 {
4892 double * real_roots = NULL;
4893 *nkeep = 0;
4894 return real_roots;
4895 }
4896
4897 // transpose of the companion matrix
4898 double * t_companion = calloc_double((p->num_poly-1)*(p->num_poly-1));
4899 size_t ii;
4900
4901
4902 // size_t m = nrows;
4903 t_companion[nrows-1] = -p->coeff[0]/p->coeff[nrows];
4904 for (ii = 1; ii < nrows; ii++){
4905 t_companion[ii * nrows + ii-1] = 1.0;
4906 t_companion[ii * nrows + nrows-1] = -p->coeff[ii]/p->coeff[nrows];
4907 }
4908 double * real = calloc_double(nrows);
4909 double * img = calloc_double(nrows);
4910 int info;
4911 int lwork = 8 * nrows;
4912 double * iwork = calloc_double(8 * nrows);
4913 //double * vl;
4914 //double * vr;
4915 int n = nrows;
4916
4917 //printf("hello! n=%d \n",n);
4918 dgeev_("N","N", &n, t_companion, &n, real, img, NULL, &n,
4919 NULL, &n, iwork, &lwork, &info);
4920
4921 //printf("info = %d", info);
4922
4923 free (iwork);
4924
4925 int * keep = calloc_int(nrows);
4926 *nkeep = 0;
4927 // the 1e-10 is kinda hacky
4928 for (ii = 0; ii < nrows; ii++){
4929 //printf("real[ii] - p->lower_bound = %G\n",real[ii]-p->lower_bound);
4930 //printf("real root = %3.15G, imag = %G \n",real[ii],img[ii]);
4931 //printf("lower thresh = %3.20G\n",p->lower_bound-1e-8);
4932 //printf("zero thresh = %3.20G\n",1e-8);
4933 //printf("upper thresh = %G\n",p->upper_bound+ZEROTHRESH);
4934 //printf("too low? %d \n", real[ii] < (p->lower_bound-1e-8));
4935 if ((fabs(img[ii]) < 1e-7) &&
4936 (real[ii] > (p->lower_bound-1e-8)) &&
4937 //(real[ii] >= (p->lower_bound-1e-7)) &&
4938 (real[ii] < (p->upper_bound+1e-8))) {
4939 //(real[ii] <= (p->upper_bound+1e-7))) {
4940
4941 //*
4942 if (real[ii] < p->lower_bound){
4943 real[ii] = p->lower_bound;
4944 }
4945 if (real[ii] > p->upper_bound){
4946 real[ii] = p->upper_bound;
4947 }
4948 //*/
4949
4950 keep[ii] = 1;
4951 *nkeep = *nkeep + 1;
4952 //printf("keep\n");
4953 }
4954 else{
4955 keep[ii] = 0;
4956 }
4957 }
4958
4959 /*
4960 printf("real portions roots = ");
4961 dprint(nrows, real);
4962 printf("imag portions roots = ");
4963 for (ii = 0; ii < nrows; ii++) printf("%E ",img[ii]);
4964 printf("\n");
4965 //dprint(nrows, img);
4966 */
4967
4968 double * real_roots = calloc_double(*nkeep);
4969 size_t counter = 0;
4970 for (ii = 0; ii < nrows; ii++){
4971 if (keep[ii] == 1){
4972 real_roots[counter] = real[ii];
4973 counter++;
4974 }
4975
4976 }
4977
4978 free(t_companion);
4979 free(real);
4980 free(img);
4981 free(keep);
4982
4983 return real_roots;
4984 }
4985
dblcompare(const void * a,const void * b)4986 static int dblcompare(const void * a, const void * b)
4987 {
4988 const double * aa = a;
4989 const double * bb = b;
4990 if ( *aa < *bb){
4991 return -1;
4992 }
4993 return 1;
4994 }
4995
4996 /********************************************************//**
4997 * Obtain the real roots of a legendre polynomial expansion
4998 *
4999 * \param[in] p - orthogonal polynomial expansion
5000 * \param[in,out] nkeep - returns how many real roots tehre are
5001 *
5002 * \return real_roots - real roots of an orthonormal polynomial expansion
5003 *
5004 * \note
5005 * Only roots within the bounds are returned
5006 * Algorithm is based on eigenvalues of non-standard companion matrix from
5007 * Roots of Polynomials Expressed in terms of orthogonal polynomials
5008 * David Day and Louis Romero 2005
5009 *
5010 * Multiplying by a factor of sqrt(2*N+1) because using orthonormal,
5011 * rather than orthogonal polynomials
5012 *************************************************************/
5013 double *
legendre_expansion_real_roots(struct OrthPolyExpansion * p,size_t * nkeep)5014 legendre_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep)
5015 {
5016
5017 double * real_roots = NULL; // output
5018 *nkeep = 0;
5019
5020 double m = (p->upper_bound - p->lower_bound) /
5021 (p->p->upper - p->p->lower);
5022 double off = p->upper_bound - m * p->p->upper;
5023
5024 orth_poly_expansion_round(&p);
5025 // print_orth_poly_expansion(p,3,NULL);
5026 //printf("last 2 = %G\n",p->coeff[p->num_poly-1]);
5027 size_t N = p->num_poly-1;
5028 //printf("N = %zu\n",N);
5029 if (N == 0){
5030 return real_roots;
5031 }
5032 else if (N == 1){
5033 if (fabs(p->coeff[N]) <= ZEROTHRESH){
5034 return real_roots;
5035 }
5036 else{
5037 double root = -p->coeff[0] / p->coeff[1];
5038 if ( (root >= -1.0-ZEROTHRESH) && (root <= 1.0 - ZEROTHRESH)){
5039 if (root <-1.0){
5040 root = -1.0;
5041 }
5042 else if (root > 1.0){
5043 root = 1.0;
5044 }
5045 *nkeep = 1;
5046 real_roots = calloc_double(1);
5047 real_roots[0] = m*root+off;
5048 }
5049 }
5050 }
5051 else{
5052 /* printf("I am here\n"); */
5053 double * nscompanion = calloc_double(N*N); // nonstandard companion
5054 size_t ii;
5055 double hnn1 = - (double) (N) / (2.0 * (double) (N) - 1.0);
5056 /* double hnn1 = - 1.0 / p->p->an(N); */
5057 nscompanion[1] = 1.0;
5058 /* nscompanion[(N-1)*N] += hnn1 * p->coeff[0] / p->coeff[N]; */
5059 nscompanion[(N-1)*N] += hnn1 * p->coeff[0] / (p->coeff[N] * sqrt(2*N+1));
5060 for (ii = 1; ii < N-1; ii++){
5061 assert (fabs(p->p->bn(ii)) < 1e-14);
5062 double in = (double) ii;
5063 nscompanion[ii*N+ii-1] = in / ( 2.0 * in + 1.0);
5064 nscompanion[ii*N+ii+1] = (in + 1.0) / ( 2.0 * in + 1.0);
5065
5066 /* nscompanion[(N-1)*N + ii] += hnn1 * p->coeff[ii] / p->coeff[N]; */
5067 nscompanion[(N-1)*N + ii] += hnn1 * p->coeff[ii] * sqrt(2*ii+1)/ p->coeff[N] / sqrt(2*N+1);
5068 }
5069 nscompanion[N*N-2] += (double) (N-1) / (2.0 * (double) (N-1) + 1.0);
5070
5071
5072 /* nscompanion[N*N-1] += hnn1 * p->coeff[N-1] / p->coeff[N]; */
5073 nscompanion[N*N-1] += hnn1 * p->coeff[N-1] * sqrt(2*(N-1)+1)/ p->coeff[N] / sqrt(2*N+1);
5074
5075 //printf("good up to here!\n");
5076 //dprint2d_col(N,N,nscompanion);
5077
5078 int info;
5079 double * scale = calloc_double(N);
5080 //*
5081 //Balance
5082 int ILO, IHI;
5083 //printf("am I here? N=%zu \n",N);
5084 //dprint(N*N,nscompanion);
5085 dgebal_("S", (int*)&N, nscompanion, (int *)&N,&ILO,&IHI,scale,&info);
5086 //printf("yep\n");
5087 if (info < 0){
5088 fprintf(stderr, "Calling dgebl had error in %d-th input in the legendre_expansion_real_roots function\n",info);
5089 exit(1);
5090 }
5091
5092 //printf("balanced!\n");
5093 //dprint2d_col(N,N,nscompanion);
5094
5095 //IHI = M1;
5096 //printf("M1=%zu\n",M1);
5097 //printf("ilo=%zu\n",ILO);
5098 //printf("IHI=%zu\n",IHI);
5099 //*/
5100
5101 double * real = calloc_double(N);
5102 double * img = calloc_double(N);
5103 //printf("allocated eigs N = %zu\n",N);
5104 int lwork = 8 * (int)N;
5105 //printf("got lwork\n");
5106 double * iwork = calloc_double(8*N);
5107 //printf("go here");
5108
5109 //dgeev_("N","N", &N, nscompanion, &N, real, img, NULL, &N,
5110 // NULL, &N, iwork, &lwork, &info);
5111 dhseqr_("E","N",(int*)&N,&ILO,&IHI,nscompanion,(int*)&N,real,img,NULL,(int*)&N,iwork,&lwork,&info);
5112 //printf("done here");
5113
5114 if (info < 0){
5115 fprintf(stderr, "Calling dhesqr had error in %d-th input in the legendre_expansion_real_roots function\n",info);
5116 exit(1);
5117 }
5118 else if(info > 0){
5119 //fprintf(stderr, "Eigenvalues are still uncovered in legendre_expansion_real_roots function\n");
5120 // printf("coeffs are \n");
5121 // dprint(p->num_poly, p->coeff);
5122 // printf("last 2 = %G\n",p->coeff[p->num_poly-1]);
5123 // exit(1);
5124 }
5125
5126 // printf("eigenvalues \n");
5127 size_t * keep = calloc_size_t(N);
5128 for (ii = 0; ii < N; ii++){
5129 /* printf("(%3.15G, %3.15G)\n",real[ii],img[ii]); */
5130 if ((fabs(img[ii]) < 1e-6) && (real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){
5131 if (real[ii] < -1.0){
5132 real[ii] = -1.0;
5133 }
5134 else if (real[ii] > 1.0){
5135 real[ii] = 1.0;
5136 }
5137 keep[ii] = 1;
5138 *nkeep = *nkeep + 1;
5139 }
5140 }
5141
5142
5143 if (*nkeep > 0){
5144 real_roots = calloc_double(*nkeep);
5145 size_t counter = 0;
5146 for (ii = 0; ii < N; ii++){
5147 if (keep[ii] == 1){
5148 real_roots[counter] = real[ii]*m+off;
5149 counter++;
5150 }
5151 }
5152 }
5153
5154
5155 free(keep); keep = NULL;
5156 free(iwork); iwork = NULL;
5157 free(real); real = NULL;
5158 free(img); img = NULL;
5159 free(nscompanion); nscompanion = NULL;
5160 free(scale); scale = NULL;
5161 }
5162
5163 if (*nkeep > 1){
5164 qsort(real_roots, *nkeep, sizeof(double), dblcompare);
5165 }
5166 return real_roots;
5167 }
5168
5169 /********************************************************//**
5170 * Obtain the real roots of a chebyshev polynomial expansion
5171 *
5172 * \param[in] p - orthogonal polynomial expansion
5173 * \param[in,out] nkeep - returns how many real roots tehre are
5174 *
5175 * \return real_roots - real roots of an orthonormal polynomial expansion
5176 *
5177 * \note
5178 * Only roots within the bounds are returned
5179 * Algorithm is based on eigenvalues of non-standard companion matrix from
5180 * Roots of Polynomials Expressed in terms of orthogonal polynomials
5181 * David Day and Louis Romero 2005
5182 *
5183 * Multiplying by a factor of sqrt(2*N+1) because using orthonormal,
5184 * rather than orthogonal polynomials
5185 *************************************************************/
5186 double *
chebyshev_expansion_real_roots(struct OrthPolyExpansion * p,size_t * nkeep)5187 chebyshev_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep)
5188 {
5189 /* fprintf(stderr, "Chebyshev real_roots not finished yet\n"); */
5190 /* exit(1); */
5191 double * real_roots = NULL; // output
5192 *nkeep = 0;
5193
5194 double m = (p->upper_bound - p->lower_bound) / (p->p->upper - p->p->lower);
5195 double off = p->upper_bound - m * p->p->upper;
5196
5197
5198 /* printf("coeff pre truncate = "); dprint(p->num_, p->coeff); */
5199 /* for (size_t ii = 0; ii < p->num_poly; ii++){ */
5200 /* if (fabs(p->coeff[ii]) < 1e-13){ */
5201 /* p->coeff[ii] = 0.0; */
5202 /* } */
5203 /* } */
5204 orth_poly_expansion_round(&p);
5205
5206 size_t N = p->num_poly-1;
5207 if (N == 0){
5208 return real_roots;
5209 }
5210 else if (N == 1){
5211 if (fabs(p->coeff[N]) <= ZEROTHRESH){
5212 return real_roots;
5213 }
5214 else{
5215 double root = -p->coeff[0] / p->coeff[1];
5216 if ( (root >= -1.0-ZEROTHRESH) && (root <= 1.0 - ZEROTHRESH)){
5217 if (root <-1.0){
5218 root = -1.0;
5219 }
5220 else if (root > 1.0){
5221 root = 1.0;
5222 }
5223 *nkeep = 1;
5224 real_roots = calloc_double(1);
5225 real_roots[0] = m*root+off;
5226 }
5227 }
5228 }
5229 else{
5230 /* printf("I am heare\n"); */
5231 /* dprint(N+1, p->coeff); */
5232 double * nscompanion = calloc_double(N*N); // nonstandard companion
5233 size_t ii;
5234
5235 double hnn1 = 0.5;
5236 double gamma = p->coeff[N];
5237
5238 nscompanion[1] = 1.0;
5239 nscompanion[(N-1)*N] -= hnn1*p->coeff[0] / gamma;
5240 for (ii = 1; ii < N-1; ii++){
5241 assert (fabs(p->p->bn(ii)) < 1e-14);
5242
5243 nscompanion[ii*N+ii-1] = 0.5; // ii-th column
5244 nscompanion[ii*N+ii+1] = 0.5;
5245
5246 // update last column
5247 nscompanion[(N-1)*N + ii] -= hnn1 * p->coeff[ii] / gamma;
5248 }
5249 nscompanion[N*N-2] += 0.5;
5250 nscompanion[N*N-1] -= hnn1 * p->coeff[N-1] / gamma;
5251
5252 //printf("good up to here!\n");
5253 /* dprint2d_col(N,N,nscompanion); */
5254
5255 int info;
5256 double * scale = calloc_double(N);
5257 //*
5258 //Balance
5259 int ILO, IHI;
5260 //printf("am I here? N=%zu \n",N);
5261 //dprint(N*N,nscompanion);
5262 dgebal_("S", (int*)&N, nscompanion, (int *)&N,&ILO,&IHI,scale,&info);
5263 //printf("yep\n");
5264 if (info < 0){
5265 fprintf(stderr, "Calling dgebl had error in %d-th input in the chebyshev_expansion_real_roots function\n",info);
5266 exit(1);
5267 }
5268
5269 //printf("balanced!\n");
5270 //dprint2d_col(N,N,nscompanion);
5271
5272 //IHI = M1;
5273 //printf("M1=%zu\n",M1);
5274 //printf("ilo=%zu\n",ILO);
5275 //printf("IHI=%zu\n",IHI);
5276 //*/
5277
5278 double * real = calloc_double(N);
5279 double * img = calloc_double(N);
5280 //printf("allocated eigs N = %zu\n",N);
5281 int lwork = 8 * (int)N;
5282 //printf("got lwork\n");
5283 double * iwork = calloc_double(8*N);
5284 //printf("go here");
5285
5286 //dgeev_("N","N", &N, nscompanion, &N, real, img, NULL, &N,
5287 // NULL, &N, iwork, &lwork, &info);
5288 dhseqr_("E","N",(int*)&N,&ILO,&IHI,nscompanion,(int*)&N,real,img,NULL,(int*)&N,iwork,&lwork,&info);
5289 //printf("done here");
5290
5291 if (info < 0){
5292 fprintf(stderr, "Calling dhesqr had error in %d-th input in the legendre_expansion_real_roots function\n",info);
5293 exit(1);
5294 }
5295 else if(info > 0){
5296 //fprintf(stderr, "Eigenvalues are still uncovered in legendre_expansion_real_roots function\n");
5297 // printf("coeffs are \n");
5298 // dprint(p->num_poly, p->coeff);
5299 // printf("last 2 = %G\n",p->coeff[p->num_poly-1]);
5300 // exit(1);
5301 }
5302
5303 /* printf("eigenvalues \n"); */
5304 size_t * keep = calloc_size_t(N);
5305 for (ii = 0; ii < N; ii++){
5306 /* printf("(%3.15G, %3.15G)\n",real[ii],img[ii]); */
5307 if ((fabs(img[ii]) < 1e-6) && (real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){
5308 /* if ((real[ii] > -1.0-1e-12) && (real[ii] < 1.0+1e-12)){ */
5309 if (real[ii] < -1.0){
5310 real[ii] = -1.0;
5311 }
5312 else if (real[ii] > 1.0){
5313 real[ii] = 1.0;
5314 }
5315 keep[ii] = 1;
5316 *nkeep = *nkeep + 1;
5317 }
5318 }
5319
5320 /* printf("nkeep = %zu\n", *nkeep); */
5321
5322 if (*nkeep > 0){
5323 real_roots = calloc_double(*nkeep);
5324 size_t counter = 0;
5325 for (ii = 0; ii < N; ii++){
5326 if (keep[ii] == 1){
5327 real_roots[counter] = real[ii]*m+off;
5328 counter++;
5329 }
5330 }
5331 }
5332
5333
5334 free(keep); keep = NULL;
5335 free(iwork); iwork = NULL;
5336 free(real); real = NULL;
5337 free(img); img = NULL;
5338 free(nscompanion); nscompanion = NULL;
5339 free(scale); scale = NULL;
5340 }
5341
5342 if (*nkeep > 1){
5343 qsort(real_roots, *nkeep, sizeof(double), dblcompare);
5344 }
5345 return real_roots;
5346 }
5347
5348 /********************************************************//**
5349 * Obtain the real roots of a orthogonal polynomial expansion
5350 *
5351 * \param[in] p - orthogonal polynomial expansion
5352 * \param[in] nkeep - returns how many real roots tehre are
5353 *
5354 * \return real_roots - real roots of an orthonormal polynomial expansion
5355 *
5356 * \note
5357 * Only roots within the bounds are returned
5358 *************************************************************/
5359 double *
orth_poly_expansion_real_roots(struct OrthPolyExpansion * p,size_t * nkeep)5360 orth_poly_expansion_real_roots(struct OrthPolyExpansion * p, size_t * nkeep)
5361 {
5362 double * real_roots = NULL;
5363 enum poly_type ptype = p->p->ptype;
5364 switch(ptype){
5365 case LEGENDRE:
5366 real_roots = legendre_expansion_real_roots(p,nkeep);
5367 break;
5368 case STANDARD:
5369 assert (1 == 0);
5370 //x need to convert polynomial to standard polynomial first
5371 //real_roots = standard_poly_real_roots(sp,nkeep);
5372 break;
5373 case CHEBYSHEV:
5374 real_roots = chebyshev_expansion_real_roots(p,nkeep);
5375 break;
5376 case HERMITE:
5377 assert (1 == 0);
5378 break;
5379 case FOURIER:
5380 assert (1 == 0);
5381 break;
5382 }
5383 return real_roots;
5384 }
5385
5386 /********************************************************//**
5387 * Obtain the maximum of an orthogonal polynomial expansion
5388 *
5389 * \param[in] p - orthogonal polynomial expansion
5390 * \param[in] x - location of maximum value
5391 *
5392 * \return maxval - maximum value
5393 *
5394 * \note
5395 * if constant function then just returns the left most point
5396 *************************************************************/
orth_poly_expansion_max(struct OrthPolyExpansion * p,double * x)5397 double orth_poly_expansion_max(struct OrthPolyExpansion * p, double * x)
5398 {
5399
5400 double maxval;
5401 double tempval;
5402
5403 assert(p->p->ptype != FOURIER);
5404 maxval = orth_poly_expansion_eval(p,p->lower_bound);
5405 *x = p->lower_bound;
5406
5407 tempval = orth_poly_expansion_eval(p,p->upper_bound);
5408 if (tempval > maxval){
5409 maxval = tempval;
5410 *x = p->upper_bound;
5411 }
5412
5413 if (p->num_poly > 2){
5414 size_t nroots;
5415 struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p);
5416 double * roots = orth_poly_expansion_real_roots(deriv,&nroots);
5417 if (nroots > 0){
5418 size_t ii;
5419 for (ii = 0; ii < nroots; ii++){
5420 tempval = orth_poly_expansion_eval(p, roots[ii]);
5421 if (tempval > maxval){
5422 *x = roots[ii];
5423 maxval = tempval;
5424 }
5425 }
5426 }
5427
5428 free(roots); roots = NULL;
5429 orth_poly_expansion_free(deriv); deriv = NULL;
5430 }
5431 return maxval;
5432 }
5433
5434 /********************************************************//**
5435 * Obtain the minimum of an orthogonal polynomial expansion
5436 *
5437 * \param[in] p - orthogonal polynomial expansion
5438 * \param[in,out] x - location of minimum value
5439 *
5440 * \return minval - minimum value
5441 *************************************************************/
orth_poly_expansion_min(struct OrthPolyExpansion * p,double * x)5442 double orth_poly_expansion_min(struct OrthPolyExpansion * p, double * x)
5443 {
5444 assert(p->p->ptype != FOURIER);
5445
5446 double minval;
5447 double tempval;
5448
5449 minval = orth_poly_expansion_eval(p,p->lower_bound);
5450 *x = p->lower_bound;
5451
5452 tempval = orth_poly_expansion_eval(p,p->upper_bound);
5453 if (tempval < minval){
5454 minval = tempval;
5455 *x = p->upper_bound;
5456 }
5457
5458 if (p->num_poly > 2){
5459 size_t nroots;
5460 struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p);
5461 double * roots = orth_poly_expansion_real_roots(deriv,&nroots);
5462 if (nroots > 0){
5463 size_t ii;
5464 for (ii = 0; ii < nroots; ii++){
5465 tempval = orth_poly_expansion_eval(p, roots[ii]);
5466 if (tempval < minval){
5467 *x = roots[ii];
5468 minval = tempval;
5469 }
5470 }
5471 }
5472 free(roots); roots = NULL;
5473 orth_poly_expansion_free(deriv); deriv = NULL;
5474 }
5475 return minval;
5476 }
5477
5478 /********************************************************//**
5479 * Obtain the maximum in absolute value of an orthogonal polynomial expansion
5480 *
5481 * \param[in] p - orthogonal polynomial expansion
5482 * \param[in,out] x - location of maximum
5483 * \param[in] oargs - optimization arguments
5484 * required for HERMITE otherwise can set NULL
5485 *
5486 * \return maxval : maximum value (absolute value)
5487 *
5488 * \note
5489 * if no roots then either lower or upper bound
5490 *************************************************************/
orth_poly_expansion_absmax(struct OrthPolyExpansion * p,double * x,void * oargs)5491 double orth_poly_expansion_absmax(
5492 struct OrthPolyExpansion * p, double * x, void * oargs)
5493 {
5494
5495 //printf("in absmax\n");
5496 // print_orth_poly_expansion(p,3,NULL);
5497 //printf("%G\n", orth_poly_expansion_norm(p));
5498
5499 enum poly_type ptype = p->p->ptype;
5500 if (oargs != NULL){
5501
5502 struct c3Vector * optnodes = oargs;
5503 double mval = fabs(orth_poly_expansion_eval(p,optnodes->elem[0]));
5504 *x = optnodes->elem[0];
5505 double cval = mval;
5506 if (ptype == HERMITE){
5507 mval *= exp(-pow(optnodes->elem[0],2)/2.0);
5508 }
5509 *x = optnodes->elem[0];
5510 for (size_t ii = 0; ii < optnodes->size; ii++){
5511 double val = fabs(orth_poly_expansion_eval(p,optnodes->elem[ii]));
5512 double tval = val;
5513 if (ptype == HERMITE){
5514 val *= exp(-pow(optnodes->elem[ii],2)/2.0);
5515 //printf("ii=%zu, x = %G. val=%G, tval=%G\n",ii,optnodes->elem[ii],val,tval);
5516 }
5517 if (val > mval){
5518 // printf("min achieved\n");
5519 mval = val;
5520 cval = tval;
5521 *x = optnodes->elem[ii];
5522 }
5523 }
5524 // printf("optloc=%G .... cval = %G\n",*x,cval);
5525 return cval;
5526 }
5527 else if ((ptype == HERMITE) || (ptype == FOURIER)){
5528 fprintf(stderr,"Must specify optimizatino arguments\n");
5529 fprintf(stderr,"In the form of candidate points for \n");
5530 fprintf(stderr,"finding the absmax of hermite or fourier expansion\n");
5531 exit(1);
5532
5533 }
5534 double maxval;
5535 double norm = orth_poly_expansion_norm(p);
5536
5537 if (norm < ZEROTHRESH) {
5538 *x = p->lower_bound;
5539 maxval = 0.0;
5540 }
5541 else{
5542 //printf("nroots=%zu\n", nroots);
5543 double tempval;
5544
5545 maxval = fabs(orth_poly_expansion_eval(p,p->lower_bound));
5546 *x = p->lower_bound;
5547
5548 tempval = fabs(orth_poly_expansion_eval(p,p->upper_bound));
5549 if (tempval > maxval){
5550 maxval = tempval;
5551 *x = p->upper_bound;
5552 }
5553 if (p->num_poly > 2){
5554 size_t nroots;
5555 struct OrthPolyExpansion * deriv = orth_poly_expansion_deriv(p);
5556 double * roots = orth_poly_expansion_real_roots(deriv,&nroots);
5557 if (nroots > 0){
5558 size_t ii;
5559 for (ii = 0; ii < nroots; ii++){
5560 tempval = fabs(orth_poly_expansion_eval(p, roots[ii]));
5561 if (tempval > maxval){
5562 *x = roots[ii];
5563 maxval = tempval;
5564 }
5565 }
5566 }
5567
5568 free(roots); roots = NULL;
5569 orth_poly_expansion_free(deriv); deriv = NULL;
5570 }
5571 }
5572 //printf("done\n");
5573 return maxval;
5574 }
5575
5576
5577 /////////////////////////////////////////////////////////
5578 // Utilities
convert_ptype_to_char(enum poly_type ptype)5579 char * convert_ptype_to_char(enum poly_type ptype)
5580 {
5581 char * out = NULL;
5582 switch (ptype) {
5583 case LEGENDRE:
5584 out = "Legendre";
5585 break;
5586 case HERMITE:
5587 out = "Hermite";
5588 break;
5589 case CHEBYSHEV:
5590 out = "Chebyshev";
5591 break;
5592 case STANDARD:
5593 out = "Standard";
5594 break;
5595 case FOURIER:
5596 out = "Fourier";
5597 break;
5598 //default:
5599 // fprintf(stderr, "Polynomial type does not exist: %d\n ", ptype);
5600 }
5601 return out;
5602
5603 }
print_orth_poly_expansion(struct OrthPolyExpansion * p,size_t prec,void * args,FILE * fp)5604 void print_orth_poly_expansion(struct OrthPolyExpansion * p, size_t prec,
5605 void * args, FILE *fp)
5606 {
5607
5608 if (args == NULL){
5609 fprintf(fp, "Orthogonal Polynomial Expansion:\n");
5610 fprintf(fp, "--------------------------------\n");
5611 fprintf(fp, "Polynomial basis is %s\n",convert_ptype_to_char(p->p->ptype));
5612 fprintf(fp, "Coefficients = ");
5613 size_t ii;
5614 if (p->p->ptype == FOURIER){
5615 for (ii = 0; ii < p->num_poly; ii++){
5616 fprintf(fp, "%3.3f %3.3f\n",
5617 creal(p->ccoeff[ii]), cimag(p->ccoeff[ii]));
5618 }
5619 }
5620 else{
5621 for (ii = 0; ii < p->num_poly; ii++){
5622 if (prec == 0){
5623 fprintf(fp, "%3.1f ", p->coeff[ii]);
5624 }
5625 else if (prec == 1){
5626 fprintf(fp, "%3.3f ", p->coeff[ii]);
5627 }
5628 else if (prec == 2){
5629 fprintf(fp, "%3.15f ", p->coeff[ii]);
5630 }
5631 else{
5632 fprintf(fp, "%3.15E ", p->coeff[ii]);
5633 }
5634 }
5635 }
5636 printf("\n");
5637 }
5638 }
5639
5640