Doxygen 1.15.0
Toolkit for Adaptive Stochastic Modeling and Non-Intrusive ApproximatioN: Tasmanian v8.2
Loading...
Searching...
No Matches
tsgRuleLocalPolynomial.hpp
1/*
2 * Copyright (c) 2017, Miroslav Stoyanov
3 *
4 * This file is part of
5 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
6 *
7 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
12 * and the following disclaimer in the documentation and/or other materials provided with the distribution.
13 *
14 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
15 * or promote products derived from this software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
18 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
19 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
20 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
21 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
22 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23 *
24 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
25 * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
26 * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
27 * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
28 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
29 */
30
31#ifndef __TSG_RULE_LOCAL_POLYNOMIAL_HPP
32#define __TSG_RULE_LOCAL_POLYNOMIAL_HPP
33
34#include "tsgCoreOneDimensional.hpp"
35
36namespace TasGrid{
37
38#ifndef __TASMANIAN_DOXYGEN_SKIP
39
40namespace RuleLocal {
41 // effective rule type
42 enum class erule {
43 pwc, localp, semilocalp, localp0, localpb
44 };
45
46 inline erule getEffectiveRule(int order, TypeOneDRule rule) {
47 if (order == 0) return erule::pwc;
48 switch(rule) {
49 case rule_semilocalp: return erule::semilocalp;
50 case rule_localp0: return erule::localp0;
51 case rule_localpb: return erule::localpb;
52 default:
53 return erule::localp;
54 }
55 }
56 inline TypeOneDRule getRule(erule effective_rule) {
57 switch(effective_rule) {
58 case erule::pwc:
59 case erule::localp: return rule_localp;
60 case erule::semilocalp: return rule_semilocalp;
61 case erule::localp0: return rule_localp0;
62 default: //case erule::localpb:
63 return rule_localpb;
64 };
65 }
66
67 template<erule effective_rule>
68 int getNumPoints(int level) {
69 switch(effective_rule) {
70 case erule::pwc: {
71 int n = 1;
72 while (level-- > 0) n *= 3;
73 return n;
74 }
75 case erule::localp:
76 case erule::semilocalp:
77 return (level == 0) ? 1 : ((1 << level) + 1);
78 case erule::localp0:
79 return (1 << (level+1)) -1;
80 default: // case erule::localpb:
81 return ((1 << level) + 1);
82 };
83 }
84
85 template<erule effective_rule>
86 int getMaxNumKids() { return (effective_rule == erule::pwc) ? 4 : 2; }
87 template<erule effrule>
88 int getMaxNumParents() {
89 return ((effrule == erule::pwc or effrule == erule::semilocalp or effrule == erule::localpb) ? 2 : 1);
90 }
91
92 template<erule effective_rule>
93 int getParent(int point) {
94 switch(effective_rule) {
95 case erule::pwc:
96 return (point == 0) ? -1 : point / 3;
97 case erule::localp:
98 case erule::semilocalp: {
99 int dad = (point + 1) / 2;
100 if (point < 4) dad--;
101 return dad;
102 }
103 case erule::localp0:
104 return (point == 0) ? -1 : (point - 1) / 2;
105 default: // case erule::localpb:
106 return (point < 2) ? -1 : ((point + 1) / 2);
107 };
108 }
109
110 template<erule effective_rule>
111 int getStepParent(int point) {
112 if (effective_rule == erule::pwc){
113 int i3l3 = Maths::int3log3(point);
114 if (point == i3l3/3) return -1;
115 if (point == i3l3-1) return -1;
116 int mod3 = point % 3;
117 int mod2 = point % 2;
118 if (mod3 == 2 and mod2 == 0) return point / 3 + 1;
119 if (mod3 == 0 and mod2 == 1) return point / 3 - 1;
120 return -1;
121 }else{
122 if (effective_rule == erule::semilocalp){
123 switch(point) {
124 case 3: return 2;
125 case 4: return 1;
126 default:
127 return -1;
128 };
129 }else if (effective_rule == erule::localpb){
130 return (point == 2) ? 0 : -1;
131 }
132 return -1;
133 }
134 }
135
136 template<erule effective_rule>
137 int getKid(int point, int kid_number) {
138 switch(effective_rule) {
139 case erule::pwc:
140 if (point == 0) return (kid_number == 0) ? 1 : (kid_number==1) ? 2 : -1;
141 if (kid_number == 3){
142 int i3l3 = Maths::int3log3(point);
143 if (point == i3l3/3) return -1;
144 if (point == i3l3-1) return -1;
145 return (point % 2 == 0) ? 3*point + 3 : 3*point - 1;
146 }
147 return 3*point + kid_number;
148 case erule::localp:
149 case erule::semilocalp:
150 if (kid_number == 0){
151 switch(point) {
152 case 0: return 1;
153 case 1: return 3;
154 case 2: return 4;
155 default:
156 return 2 * point - 1;
157 };
158 } else {
159 switch(point) {
160 case 0: return 2;
161 case 1:
162 case 2: return -1;
163 default:
164 return 2 * point;
165 };
166 }
167 case erule::localp0:
168 return 2 * point + ((kid_number == 0) ? 1 : 2);
169 default: // case erule::localpb:
170 switch(point) {
171 case 0:
172 case 1:
173 return (kid_number == 0) ? 2 : -1;
174 default:
175 return 2*point - ((kid_number == 0) ? 1 : 0);
176 };
177 };
178 }
179
180 template<erule effective_rule>
181 double getNode(int point) {
182 switch(effective_rule) {
183 case erule::pwc:
184 return -2.0 + (1.0 / ((double) Maths::int3log3(point))) * (3*point + 2 - point % 2);
185 case erule::localp:
186 case erule::semilocalp:
187 switch(point) {
188 case 0: return 0.0;
189 case 1: return -1.0;
190 case 2: return 1.0;
191 default:
192 return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0;
193 };
194 case erule::localp0:
195 return ((double)(2*point + 3) ) / ((double) Maths::int2log2(point + 1) ) - 3.0;
196 default: // case erule::localpb:
197 switch(point) {
198 case 0: return -1.0;
199 case 1: return 1.0;
200 case 2: return 0.0;
201 default:
202 return ((double)(2*point - 1)) / ((double) Maths::int2log2(point - 1)) - 3.0;
203 };
204 };
205 }
206
207 template<erule effective_rule>
208 int getLevel(int point) {
209 switch(effective_rule) {
210 case erule::pwc: {
211 int level = 0;
212 while(point >= 1){ point /= 3; level += 1; }
213 return level;
214 }
215 case erule::localp:
216 case erule::semilocalp:
217 return (point == 0) ? 0 : (point == 1) ? 1 : (Maths::intlog2(point - 1) + 1);
218 case erule::localp0:
219 return Maths::intlog2(point + 1);
220 default: // case erule::localpb:
221 return (point <= 1) ? 0 : (Maths::intlog2(point - 1) + 1);
222 };
223 }
224
225 template<erule effective_rule>
226 double getSupport(int point) {
227 switch(effective_rule) {
228 case erule::pwc:
229 return 1.0 / (double) Maths::int3log3(point);
230 case erule::localp:
231 case erule::semilocalp:
232 return (point == 0) ? 1.0 : 1.0 / ((double) Maths::int2log2(point - 1));
233 case erule::localp0:
234 return 1.0 / ((double) Maths::int2log2(point + 1));
235 default: // case erule::localpb:
236 return (point <= 1) ? 2.0 : 1.0 / ((double) Maths::int2log2(point - 1));
237 };
238 }
239
240 template<erule effective_rule>
241 double scaleDiffX(int point) {
242 switch(effective_rule) {
243 case erule::pwc:
244 return 0.0;
245 case erule::localp:
246 return (point <= 2) ? 1.0 : static_cast<double>(Maths::int2log2(point - 1));
247 case erule::semilocalp:
248 return static_cast<double>(Maths::int2log2(point - 1));
249 case erule::localp0:
250 return (point == 0) ? 1.0 : static_cast<double>(Maths::int2log2(point + 1));
251 default: // case erule::localpb:
252 switch(point) {
253 case 0:
254 case 1: return 0.5;
255 case 2: return 1.0;
256 default:
257 return static_cast<double>(Maths::int2log2(point - 1));
258 };
259 };
260 }
261
262 template<erule effective_rule>
263 double scaleX(int point, double x) {
264 switch(effective_rule) {
265 case erule::pwc:
266 return 0.0; // should never be called
267 case erule::localp:
268 switch(point) {
269 case 0: return x;
270 case 1: return (x + 1.0);
271 case 2: return (x - 1.0);
272 default:
273 return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point));
274 };
275 case erule::semilocalp:
276 return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point));
277 case erule::localp0:
278 return ((double) Maths::int2log2(point + 1) * (x + 3.0) - 3.0 - (double) (2*point));
279 default: // case erule::localpb:
280 switch(point) {
281 case 0: return (x + 1.0) / 2.0;
282 case 1: return (x - 1.0) / 2.0;
283 case 2: return x;
284 default:
285 return ((double) Maths::int2log2(point - 1) * (x + 3.0) + 1.0 - (double) (2*point));
286 };
287 };
288 }
289
290 template<erule effective_rule>
291 double evalPWQuadratic(int point, double x) {
292 if (effective_rule == erule::localp){
293 switch(point) {
294 case 1: return 1.0 - x;
295 case 2: return 1.0 + x;
296 default:
297 return (1.0 - x) * (1.0 + x);
298 };
299 }else if (effective_rule == erule::localpb){
300 switch(point) {
301 case 0: return 1.0 - x;
302 case 1: return 1.0 + x;
303 default:
304 return (1.0 - x) * (1.0 + x);
305 };
306 }
307 return (1.0 - x) * (1.0 + x);
308 }
309 template<erule effective_rule>
310 double evalPWCubic(int point, double x) {
311 if (effective_rule == erule::localp){
312 switch(point) {
313 case 0: return 1.0;
314 case 1: return 1.0 - x;
315 case 2: return 1.0 + x;
316 case 3:
317 case 4: return (1.0 - x) * (1.0 + x);
318 default:
319 return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0;
320 };
321 }else if (effective_rule == erule::localp0){
322 if (point == 0) return (1.0 - x) * (1.0 + x);
323 }else if (effective_rule == erule::localpb){
324 switch(point) {
325 case 0: return 1.0 - x;
326 case 1: return 1.0 + x;
327 case 2: return (1.0 - x) * (1.0 + x);
328 default:
329 return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0;
330 };
331 }
332 return (point % 2 == 0) ? (1.0 - x) * (1.0 + x) * (3.0 + x) / 3.0 : (1.0 - x) * (1.0 + x) * (3.0 - x) / 3.0;
333 }
334
335 template<erule effective_rule>
336 double evalPWPower(int max_order, int point, double x) {
337 // use the cubic implementation until we have enough points
338 if (effective_rule == erule::localp) if (point <= 8) return evalPWCubic<effective_rule>(point, x); // if order is cubic or less, use the hard-coded functions
339 if (effective_rule == erule::semilocalp) if (point <= 4) return evalPWCubic<effective_rule>(point, x);
340 if (effective_rule == erule::localpb) if (point <= 4) return evalPWCubic<effective_rule>(point, x);
341 if (effective_rule == erule::localp0) if (point <= 2) return evalPWCubic<effective_rule>(point, x);
342 int level = getLevel<effective_rule>(point);
343 int most_turns = 1;
344 double value = (1.0 - x)*(1.0 + x);
345 double phantom_distance = 1.0;
346 int max_ancestors = [&]()->int{ // maximum number of ancestors to consider, first set the possible max then constrain by max_order
347 switch(effective_rule) {
348 case erule::pwc: return 0; // should never happen
349 case erule::localp: return max_ancestors = level-2;
350 case erule::semilocalp: return max_ancestors = level-1;
351 case erule::localpb: return max_ancestors = level-1;
352 default: return max_ancestors = level;
353 };
354 }();
355 if (max_order > 0) max_ancestors = std::min(max_ancestors, max_order - 2); // use the minimum of the available ancestors or the order restriction
356
357 for(int j=0; j < max_ancestors; j++){
358 // Lagrange polynomial needs to be constructed using normalized x (basis support (-1, 1)).
359 // The support of the basis is equal to 2, thus we use units of "half-support"
360 // The first two nodes are used in the initialization of value, those are the nearest ancestors (in spacial distance, not hierarchy level).
361 // The other nodes are "phantoms" that lay strictly outside of the support [-1, 1] (i.e., not on the edge)
362 // The walking distance more than doubles and is an odd number (due to the half-support)
363 // The walk through the ancestors can take a left or right turn at each step,
364 // most_turns is the total number of turns possible for the current ancestor, turns (or most_turns - 1 - turns) is the actual number of turns
365 // Every time we turn, we backtrack and we lose 2 units, the phantom node is at maximum distance minus 2 time the number of turns (to the left or right)
366 most_turns *= 2;
367 phantom_distance = 2.0 * phantom_distance + 1.0;
368 int turns = (effective_rule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns);
369 double node = (turns < most_turns / 2) ? (phantom_distance - 2.0 * ((double) turns)) : (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns)));
370 value *= - ( x - node ) / node;
371 }
372 return value;
373 }
374
375 template<erule effective_rule>
376 double evalSupport(int max_order, int point, double x, bool &isSupported) {
377 switch(effective_rule) {
378 case erule::pwc: {
379 double distance = std::abs(x - getNode<effective_rule>(point));
380 double support = getSupport<effective_rule>(point);
381 isSupported = (distance <= 2.0 * support);
382 return (distance <= support) ? 1.0 : 0.0;
383 }
384 case erule::localp:
385 isSupported = true;
386 if (point == 0) {
387 return 1.0;
388 } else {
389 double xn = scaleX<effective_rule>(point, x);
390 if (std::abs(xn) <= 1.0) {
391 switch(max_order) {
392 case 1: return 1.0 - std::abs(xn);
393 case 2: return evalPWQuadratic<effective_rule>(point, xn);
394 case 3: return evalPWCubic<effective_rule>(point, xn);
395 default:
396 return evalPWPower<effective_rule>(max_order, point, xn);
397 };
398 } else {
399 isSupported = false;
400 return 0.0;
401 }
402 }
403 case erule::semilocalp:
404 isSupported = true;
405 switch(point) {
406 case 0: return 1.0;
407 case 1: return 0.5 * x * (x - 1.0);
408 case 2: return 0.5 * x * (x + 1.0);
409 default: {
410 double xn = scaleX<effective_rule>(point, x);
411 if (std::abs(xn) <= 1.0) {
412 switch(max_order) {
413 case 1: return 1.0 - std::abs(xn);
414 case 2: return evalPWQuadratic<effective_rule>(point, xn);
415 case 3: return evalPWCubic<effective_rule>(point, xn);
416 default:
417 return evalPWPower<effective_rule>(max_order, point, xn);
418 };
419 } else {
420 isSupported = false;
421 return 0.0;
422 }
423 }
424 };
425 case erule::localp0:
426 default: { // case erule::localpb:
427 double xn = scaleX<effective_rule>(point, x);
428 if (std::abs(xn) <= 1.0) {
429 isSupported = true;
430 switch(max_order) {
431 case 1: return 1.0 - std::abs(xn);
432 case 2: return evalPWQuadratic<effective_rule>(point, xn);
433 case 3: return evalPWCubic<effective_rule>(point, xn);
434 default:
435 return evalPWPower<effective_rule>(max_order, point, xn);
436 };
437 } else {
438 isSupported = false;
439 return 0.0;
440 }
441 }
442 };
443 }
444
445 template<erule effective_rule>
446 double evalRaw(int max_order, int point, double x) {
447 switch(effective_rule) {
448 case erule::pwc:
449 return (std::abs(x - getNode<effective_rule>(point)) <= getSupport<effective_rule>(point)) ? 1.0 : 0.0;
450 case erule::localp:
451 if (point == 0) {
452 return 1.0;
453 } else {
454 double xn = scaleX<effective_rule>(point, x);
455 if (std::abs(xn) <= 1.0) {
456 switch(max_order) {
457 case 1: return 1.0 - std::abs(xn);
458 case 2: return evalPWQuadratic<effective_rule>(point, xn);
459 case 3: return evalPWCubic<effective_rule>(point, xn);
460 default:
461 return evalPWPower<effective_rule>(max_order, point, xn);
462 };
463 } else {
464 return 0.0;
465 }
466 }
467 case erule::semilocalp:
468 switch(point) {
469 case 0: return 1.0;
470 case 1: return 0.5 * x * (x - 1.0);
471 case 2: return 0.5 * x * (x + 1.0);
472 default: {
473 double xn = scaleX<effective_rule>(point, x);
474 if (std::abs(xn) <= 1.0) {
475 switch(max_order) {
476 case 1: return 1.0 - std::abs(xn);
477 case 2: return evalPWQuadratic<effective_rule>(point, xn);
478 case 3: return evalPWCubic<effective_rule>(point, xn);
479 default:
480 return evalPWPower<effective_rule>(max_order, point, xn);
481 };
482 } else {
483 return 0.0;
484 }
485 }
486 };
487 case erule::localp0:
488 default: { // case erule::localpb:
489 double xn = scaleX<effective_rule>(point, x);
490 if (std::abs(xn) <= 1.0) {
491 switch(max_order) {
492 case 1: return 1.0 - std::abs(xn);
493 case 2: return evalPWQuadratic<effective_rule>(point, xn);
494 case 3: return evalPWCubic<effective_rule>(point, xn);
495 default:
496 return evalPWPower<effective_rule>(max_order, point, xn);
497 };
498 } else {
499 return 0.0;
500 }
501 }
502 };
503 }
504
505 template<erule effective_rule>
506 double diffPWQuadratic(int point, double x) {
507 if (effective_rule == erule::localp) {
508 switch(point) {
509 case 1: return -1.0;
510 case 2: return 1.0;
511 default:
512 return -2.0 * x;
513 };
514 } else if (effective_rule == erule::localpb) {
515 switch(point) {
516 case 0: return -1.0;
517 case 1: return 1.0;
518 default:
519 return -2.0 * x;
520 };
521 }
522 return -2.0 * x;
523 }
524 template<erule effective_rule>
525 double diffPWCubic(int point, double x) {
526 if (effective_rule == erule::localp) {
527 switch(point) {
528 case 0: return 0.0;
529 case 1: return -1.0;
530 case 2: return 1.0;
531 case 3:
532 case 4: return -2.0 * x;
533 default:
534 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
535 };
536 } else if (effective_rule == erule::localpb) {
537 switch(point) {
538 case 0: return -1.0;
539 case 1: return 1.0;
540 case 2: return -2.0 * x;
541 default:
542 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
543 };
544 } else if (effective_rule == erule::localp0) {
545 return (point == 0) ? -2.0 * x : ((point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0));
546 }
547 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
548 }
549
550 template<erule effrule>
551 double diffPWPower(int max_order, int point, double x) {
552 // use the cubic implementation until we have enough points
553 if (effrule == erule::localp and point <= 8) return diffPWCubic<effrule>(point, x);
554 if (effrule == erule::semilocalp and point <= 4) return diffPWCubic<effrule>(point, x);
555 if (effrule == erule::localpb and point <= 4) return diffPWCubic<effrule>(point, x);
556 if (effrule == erule::localp0 and point <= 2) return diffPWCubic<effrule>(point, x);
557 int level = getLevel<effrule>(point);
558 int max_ancestors = [&]()->int {
559 switch(effrule) {
560 case erule::pwc: return 0;
561 case erule::localp: return level - 2;
562 case erule::semilocalp:
563 case erule::localpb: return level - 1;
564 default: return level;
565 };
566 }();
567
568 if (max_order > 0) max_ancestors = std::min(max_ancestors, max_order - 2);
569
570 // This lambda captures most_turns and phantom_distance by reference, uses those as internal state variables, and on each
571 // call it returns the next normalized ancestor node.
572 int most_turns = 1;
573 double phantom_distance = 1.0;
574 auto update_and_get_next_node = [&]() {
575 most_turns *= 2;
576 phantom_distance = 2.0 * phantom_distance + 1.0;
577 int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns);
578 return (turns < most_turns / 2) ?
579 (phantom_distance - 2.0 * ((double) turns)) :
580 (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns)));
581 };
582
583 // This lambda is the inverse transform of update_and_get_next_node() above, and returns the previous descendant node.
584 auto rollback_and_get_prev_node = [&]() {
585 most_turns /= 2;
586 phantom_distance = 0.5 * (phantom_distance - 1.0);
587 int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns);
588 return (turns < most_turns / 2) ?
589 (phantom_distance - 2.0 * ((double) turns)) :
590 (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns)));
591 };
592
593 // Does not include the additional factor (1-x) * (1+x) and the Lagrange coefficient.
594 std::vector<double> left_prods(max_ancestors);
595 left_prods[0] = 1.0;
596 double node = update_and_get_next_node();
597 double coeff = 1.0 / (-node);
598 for(int j=1; j<max_ancestors; j++) {
599 left_prods[j] = left_prods[j-1] * (x - node);
600 node = update_and_get_next_node();
601 coeff *= 1.0 / (-node);
602 }
603 double right_prod = 1.0;
604 double derivative = left_prods[max_ancestors-1];
605 for (int j=max_ancestors-2; j>=0; j--) {
606 right_prod *= x - node;
607 derivative += right_prod * left_prods[j];
608 node = rollback_and_get_prev_node();
609 }
610
611 // Adjust for the additional factor (1-x) * (1+x) and the Lagrange coefficient.
612 derivative = derivative * (1.0 - x) * (1.0 + x) + right_prod * (x - node) * (-2.0) * x;
613 derivative *= coeff;
614
615 return derivative;
616 }
617
618 template<erule effrule>
619 double diffRaw(int max_order, int point, double x) {
620 switch(effrule) {
621 case erule::pwc: return 0.0;
622 case erule::localp: {
623 if (point == 0) return 0.0;
624 double xn = scaleX<effrule>(point, x);
625 double an = scaleDiffX<effrule>(point);
626 switch(max_order) {
627 case 1: return ((xn >= 0 ? -1.0 : 1.0) * an);
628 case 2: return an * diffPWQuadratic<effrule>(point, xn);
629 case 3: return an * diffPWCubic<effrule>(point, xn);
630 default:
631 return an * diffPWPower<effrule>(point, xn);
632 };
633 }
634 case erule::semilocalp:
635 switch(point) {
636 case 0: return 0.0;
637 case 1: return x - 0.5;
638 case 2: return x + 0.5;
639 default: {
640 double xn = scaleX<effrule>(point, x);
641 double an = scaleDiffX<effrule>(point);
642 switch(max_order) {
643 case 2: return an * diffPWQuadratic<effrule>(point, xn);
644 case 3: return an * diffPWCubic<effrule>(point, xn);
645 default:
646 return an * diffPWPower<effrule>(point, xn);
647 };
648 }
649 };
650 case erule::localp0: {
651 double xn = scaleX<effrule>(point, x);
652 double an = scaleDiffX<effrule>(point);
653 switch(max_order) {
654 case 1: return (x == 1.0 and point == 0) ? -1.0 : ((xn >= 0 ? -1.0 : 1.0) * an);
655 case 2: return an * diffPWQuadratic<effrule>(point, xn);
656 case 3: return an * diffPWCubic<effrule>(point, xn);
657 default:
658 return an * diffPWPower<effrule>(point, xn);
659 };
660 }
661 default: { // case erule::localpb:
662 double xn = scaleX<effrule>(point, x);
663 double an = scaleDiffX<effrule>(point);
664 switch(max_order) {
665 case 1: return ((xn >= 0 ? -1.0 : 1.0) * an);
666 case 2: return an * diffPWQuadratic<effrule>(point, xn);
667 case 3: return an * diffPWCubic<effrule>(point, xn);
668 default:
669 return an * diffPWPower<effrule>(point, xn);
670 };
671 }
672 };
673 }
674
675 template<erule effrule>
676 double diffSupport(int max_order, int point, double x, bool &isSupported) {
677 switch(effrule) {
678 case erule::pwc:
679 isSupported = false;
680 return 0.0;
681 case erule::localp: {
682 if (point == 0) { isSupported = true; return 0.0; }
683 double xn = scaleX<effrule>(point, x);
684 isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0);
685 if (isSupported) {
686 double an = scaleDiffX<effrule>(point);
687 switch(max_order) {
688 case 1: return (x == 1.0 and point == 2) ? an : (((xn >= 0 ? -1.0 : 1.0) * an));
689 case 2: return an * diffPWQuadratic<effrule>(point, xn);
690 case 3: return an * diffPWCubic<effrule>(point, xn);
691 default:
692 return an * diffPWPower<effrule>(max_order, point, xn);
693 };
694 } else {
695 return 0.0;
696 }
697 }
698 case erule::semilocalp:
699 switch(point) {
700 case 0: { isSupported = true; return 0.0; }
701 case 1: { isSupported = true; return x - 0.5; }
702 case 2: { isSupported = true; return x + 0.5; }
703 default: {
704 double xn = scaleX<effrule>(point, x);
705 isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0);
706 if (isSupported) {
707 double an = scaleDiffX<effrule>(point);
708 switch(max_order) {
709 case 2: return an * diffPWQuadratic<effrule>(point, xn);
710 case 3: return an * diffPWCubic<effrule>(point, xn);
711 default:
712 return an * diffPWPower<effrule>(max_order, point, xn);
713 };
714 } else {
715 return 0.0;
716 }
717 }
718 };
719 case erule::localp0: {
720 double xn = scaleX<effrule>(point, x);
721 isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0);
722 if (isSupported) {
723 double an = scaleDiffX<effrule>(point);
724 switch(max_order) {
725 case 1: return (x == 1.0 and point == 0) ? -1.0 : ((xn >= 0 ? -1.0 : 1.0) * an);
726 case 2: return an * diffPWQuadratic<effrule>(point, xn);
727 case 3: return an * diffPWCubic<effrule>(point, xn);
728 default:
729 return an * diffPWPower<effrule>(max_order, point, xn);
730 };
731 } else {
732 return 0.0;
733 }
734 }
735 default: { // case erule::localpb:
736 double xn = scaleX<effrule>(point, x);
737 isSupported = (-1.0 <= xn and xn < 1.0) or (x == 1.0 and xn == 1.0);
738 if (isSupported) {
739 double an = scaleDiffX<effrule>(point);
740 switch(max_order) {
741 case 1: return ((xn >= 0 ? -1.0 : 1.0) * an);
742 case 2: return an * diffPWQuadratic<effrule>(point, xn);
743 case 3: return an * diffPWCubic<effrule>(point, xn);
744 default:
745 return an * diffPWPower<effrule>(max_order, point, xn);
746 };
747 } else {
748 return 0.0;
749 }
750 }
751 };
752 }
753
754 template<erule effrule>
755 double getArea(int max_order, int point, std::vector<double> const &w, std::vector<double> const &x) {
756 switch(effrule) {
757 case erule::pwc:
758 return 2.0 * getSupport<effrule>(point);
759 case erule::localp:
760 switch(point) {
761 case 0: return 2.0;
762 case 1:
763 case 2: return 0.5;
764 default:
765 switch(max_order) {
766 case 1: return getSupport<effrule>(point);
767 case 2:
768 case 3: return (4.0/3.0) * getSupport<effrule>(point);
769 default:
770 if (point <= 8) return (4.0/3.0) * getSupport<effrule>(point);
771 break;
772 };
773 break;
774 };
775 break;
776 case erule::semilocalp:
777 switch(point) {
778 case 0: return 2.0;
779 case 1:
780 case 2: return 1.0/3.0;
781 default:
782 switch(max_order) {
783 case 2:
784 case 3: return (4.0/3.0) * getSupport<effrule>(point);
785 default:
786 if (point <= 4) return (4.0/3.0) * getSupport<effrule>(point);
787 break;
788 };
789 break;
790 };
791 break;
792 case erule::localp0:
793 switch(max_order) {
794 case 1: return getSupport<effrule>(point);
795 case 2:
796 case 3: return (4.0/3.0) * getSupport<effrule>(point);
797 default:
798 if (point <= 2) return (4.0/3.0) * getSupport<effrule>(point);
799 break;
800 };
801 break;
802 default: // case erule::localpb:
803 if (point <= 1) {
804 return 1.0;
805 } else {
806 switch(max_order) {
807 case 1: return getSupport<effrule>(point);
808 case 2:
809 case 3: return (4.0/3.0) * getSupport<effrule>(point);
810 default:
811 if (point <= 4) return (4.0/3.0) * getSupport<effrule>(point);
812 break;
813 };
814 }
815 break;
816 };
817 double sum = 0.0;
818 for(size_t i=0; i<w.size(); i++) sum += w[i] * evalPWPower<effrule>(max_order, point, x[i]);
819 return sum * getSupport<effrule>(point);
820 }
821
822 template<erule effrule>
823 void van_matrix(int max_order, int num_rows, std::vector<int> &pntr, std::vector<int> &indx, std::vector<double> &vals) {
824 int max_level = getLevel<effrule>(num_rows);
825
826 if (effrule == erule::pwc) {
827 switch(num_rows) {
828 case 0:
829 indx = {0,};
830 vals = {1.0,};
831 pntr = {0, 1};
832 break;
833 case 1:
834 indx = {0, 0, 1};
835 vals = {1.0, 1.0, 1.0};
836 pntr = {0, 1, 3};
837 break;
838 case 2:
839 indx = {0, 0, 1, 0, 2};
840 vals = {1.0, 1.0, 1.0, 1.0, 1.0};
841 pntr = {0, 1, 3, 5};
842 break;
843 default:
844 indx.clear();
845 indx.reserve(num_rows * max_level);
846 vals.clear();
847 vals.reserve(num_rows * max_level);
848 pntr = std::vector<int>(num_rows + 1, 0);
849 for(auto i : std::array<int, 5>{0, 0, 1, 0, 2})
850 indx.push_back(i);
851 for(size_t i=0; i<5; i++)
852 vals.push_back(1.0);
853 pntr[1] = 1;
854 pntr[2] = 3;
855 {
856 std::vector<int> ancestors;
857 ancestors.reserve(max_level);
858 for(int r=3; r<num_rows; r++) {
859 pntr[r] = static_cast<int>(indx.size());
860 ancestors.clear();
861 int kid = r;
862 int dad = kid / 3;
863 while(dad != 0) {
864 if (dad % 2 == 0) {
865 if (kid != 3 * dad)
866 ancestors.push_back(dad);
867 } else {
868 if (kid != 3 * dad + 2)
869 ancestors.push_back(dad);
870 }
871 kid = dad;
872 dad = kid / 3;
873 }
874 indx.push_back(0);
875 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
876 indx.push_back(r);
877 for(size_t i=0; i<ancestors.size() + 2; i++)
878 vals.push_back(1.0);
879 }
880 pntr.back() = static_cast<int>(indx.size());
881 }
882 break;
883 };
884 } else if (effrule == erule::localp) {
885 switch(num_rows) {
886 case 0:
887 indx = {0,};
888 vals = {1.0,};
889 pntr = {0, 1};
890 break;
891 default:
892 indx.clear();
893 indx.reserve(num_rows * max_level);
894 vals.clear();
895 vals.reserve(num_rows * max_level);
896 pntr = std::vector<int>(num_rows + 1, 0);
897 indx.push_back(0);
898 vals.push_back(1.0);
899 {
900 std::vector<int> ancestors;
901 std::vector<double> ancestors_vals;
902 ancestors.reserve(max_level);
903 ancestors_vals.reserve(max_level);
904 for(int r=1; r<num_rows; r++) {
905 double x = getNode<effrule>(r);
906 pntr[r] = static_cast<int>(indx.size());
907 ancestors.clear();
908 ancestors_vals.clear();
909 int kid = r;
910 int dad = (kid + 1) / 2;
911 if (kid < 4) dad--;
912 while(dad != 0) {
913 ancestors.push_back(dad);
914 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
915 kid = dad;
916 dad = (kid + 1) / 2;
917 if (kid < 4) dad--;
918 }
919 indx.push_back(0);
920 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
921 indx.push_back(r);
922 vals.push_back(1.0);
923 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
924 vals.push_back(1.0);
925 }
926 pntr.back() = static_cast<int>(indx.size());
927 }
928 break;
929 };
930 } else if (effrule == erule::semilocalp) {
931 switch(num_rows) {
932 case 0:
933 indx = {0,};
934 vals = {1.0,};
935 pntr = {0, 1};
936 break;
937 case 1:
938 indx = {0, 0, 1};
939 vals = {1.0, 1.0, 1.0};
940 pntr = {0, 1, 3};
941 break;
942 case 2:
943 indx = {0, 0, 1, 0, 2};
944 vals = {1.0, 1.0, 1.0, 1.0, 1.0};
945 pntr = {0, 1, 3, 5};
946 break;
947 default:
948 indx.clear();
949 indx.reserve(num_rows * max_level);
950 vals.clear();
951 vals.reserve(num_rows * max_level);
952 pntr = std::vector<int>(num_rows + 1, 0);
953 for(auto i : std::array<int, 5>{0, 0, 1, 0, 2})
954 indx.push_back(i);
955 for(int i=0; i<5; i++)
956 vals.push_back(1.0);
957 pntr[1] = 1;
958 pntr[2] = 3;
959 {
960 std::vector<int> ancestors;
961 std::vector<double> ancestors_vals;
962 ancestors.reserve(max_level);
963 ancestors_vals.reserve(max_level);
964 for(int r=3; r<num_rows; r++) {
965 pntr[r] = static_cast<int>(indx.size());
966 double x = getNode<effrule>(r);
967 ancestors.clear();
968 ancestors_vals.clear();
969 int kid = r;
970 int dad = (kid + 1) / 2;
971 while(dad > 2) {
972 ancestors.push_back(dad);
973 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
974 kid = dad;
975 dad = (kid + 1) / 2;
976 }
977 indx.push_back(0);
978 indx.push_back(1);
979 indx.push_back(2);
980 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
981 indx.push_back(r);
982 vals.push_back(1.0);
983 vals.push_back( evalRaw<effrule>(max_order, 1, x) );
984 vals.push_back( evalRaw<effrule>(max_order, 2, x) );
985 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
986 vals.push_back(1.0);
987 }
988 pntr.back() = static_cast<int>(indx.size());
989 }
990 break;
991 };
992 } else if (effrule == erule::localp0) {
993 switch(num_rows) {
994 case 0:
995 indx = {0,};
996 vals = {1.0,};
997 pntr = {0, 1};
998 break;
999 default:
1000 indx.clear();
1001 indx.reserve(num_rows * max_level);
1002 vals.clear();
1003 vals.reserve(num_rows * max_level);
1004 pntr = std::vector<int>(num_rows + 1, 0);
1005 indx.push_back(0);
1006 vals.push_back(1.0);
1007 {
1008 std::vector<int> ancestors;
1009 std::vector<double> ancestors_vals;
1010 ancestors.reserve(max_level);
1011 ancestors_vals.reserve(max_level);
1012 for(int r=1; r<num_rows; r++) {
1013 double x = getNode<effrule>(r);
1014 pntr[r] = static_cast<int>(indx.size());
1015 ancestors.clear();
1016 ancestors_vals.clear();
1017 int kid = r;
1018 int dad = (kid - 1) / 2;
1019 while(dad != 0) {
1020 ancestors.push_back(dad);
1021 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
1022 kid = dad;
1023 dad = (kid - 1) / 2;
1024 }
1025 indx.push_back(0);
1026 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
1027 indx.push_back(r);
1028 vals.push_back( evalRaw<effrule>(max_order, 0, x) );
1029 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
1030 vals.push_back(1.0);
1031 }
1032 pntr.back() = static_cast<int>(indx.size());
1033 }
1034 break;
1035 };
1036 } else { // if (effrule == erule::localpb) {
1037 switch(num_rows) {
1038 case 0:
1039 indx = {0,};
1040 vals = {1.0,};
1041 pntr = {0, 1};
1042 break;
1043 case 1:
1044 indx = {0, 1};
1045 vals = {1.0, 1.0};
1046 pntr = {0, 1, 2};
1047 break;
1048 default:
1049 indx.clear();
1050 indx.reserve(num_rows * max_level);
1051 vals.clear();
1052 vals.reserve(num_rows * max_level);
1053 pntr = std::vector<int>(num_rows + 1, 0);
1054 indx.push_back(0);
1055 indx.push_back(1);
1056 vals.push_back(1.0);
1057 vals.push_back(1.0);
1058 pntr[1] = 1;
1059 {
1060 std::vector<int> ancestors;
1061 std::vector<double> ancestors_vals;
1062 ancestors.reserve(max_level);
1063 ancestors_vals.reserve(max_level);
1064 for(int r=2; r<num_rows; r++) {
1065 pntr[r] = static_cast<int>(indx.size());
1066 double x = getNode<effrule>(r);
1067 ancestors.clear();
1068 ancestors_vals.clear();
1069 int kid = r;
1070 int dad = (kid + 1) / 2;
1071 while(dad > 1) {
1072 ancestors.push_back(dad);
1073 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
1074 kid = dad;
1075 dad = (kid + 1) / 2;
1076 }
1077 indx.push_back(0);
1078 indx.push_back(1);
1079 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
1080 indx.push_back(r);
1081 vals.push_back( evalRaw<effrule>(max_order, 0, x) );
1082 vals.push_back( evalRaw<effrule>(max_order, 1, x) );
1083 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
1084 vals.push_back(1.0);
1085 }
1086 pntr.back() = static_cast<int>(indx.size());
1087 }
1088 break;
1089 };
1090 }
1091 }
1092
1093} // namespace RuleLocal
1094
1095#endif // __TASMANIAN_DOXYGEN_SKIP
1096
1097}
1098
1099#endif
TasGrid::TypeOneDRule
TypeOneDRule
Used to specify the one dimensional family of rules that induces the sparse grid.
Definition tsgEnumerates.hpp:285
TasGrid::rule_localp
@ rule_localp
Nested rule with a hierarchy of uniformly distributed nodes and functions with compact support.
Definition tsgEnumerates.hpp:362
TasGrid::rule_localpb
@ rule_localpb
Variation of rule_localp focusing nodes on the boundary instead of the interior.
Definition tsgEnumerates.hpp:368
TasGrid::rule_localp0
@ rule_localp0
Variation of rule_localp assuming the model is zero at the domain boundary.
Definition tsgEnumerates.hpp:364
TasGrid::rule_semilocalp
@ rule_semilocalp
Variation of rule_localp using increased support in exchange for higher order basis (better for smoot...
Definition tsgEnumerates.hpp:366
TasGrid
Encapsulates the Tasmanian Sparse Grid module.
Definition TasmanianSparseGrid.hpp:68