Doxygen 1.17.0
Toolkit for Adaptive Stochastic Modeling and Non-Intrusive ApproximatioN: Tasmanian v8.3
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 // The edge test is more complicated but it guards against ambiguity when dealing with derivatives
376 // at the boundary and interior cell edges
377 // The standard interpolation is not affected since the absis smoothly goes to zero around edges,
378 // thus small perturbations do not get magnified.
379 // However, derivatives have jump discontinuities and infinitesimal change can yield O(1) perturbation.
380 template<erule effective_rule, bool use_edge_test = false>
381 double evalSupport(int max_order, int point, double x, bool &isSupported) {
382 auto test_support = [&](double xn, int p) -> bool {
383 if constexpr (use_edge_test) {
384 // the first part of the test (-1 <= xn and xn < 1) restricts the support to [left, rigth)
385 // where left/right are edges of the cells
386 // This is OK for the interior, but will fail at the right boundary, thus the second test
387 // checks if the support of the basis funciton ends at the right-most point of the domain.
388 // The test ((p & (p - 1)) == 0))) uses bit-patterns to check if p is an exact power of 2.
389 if constexpr (effective_rule == erule::localp0) {
390 int const pp = p + 2; // shift by 2 points
391 return (-1.0 <= xn and (xn < 1.0 or (pp < 3 or ((pp & (pp - 1)) == 0))));
392 } else if constexpr (effective_rule == erule::localpb) {
393 return (-1.0 <= xn and (xn < 1.0 or (p < 3 or ((p & (p - 1)) == 0))));
394 } else {
395 return (-1.0 <= xn and (xn < 1.0 or (p > 1 and ((p & (p - 1)) == 0))));
396 }
397 } else {
398 return (std::abs(xn) <= 1.0);
399 }
400 };
401
402 switch(effective_rule) {
403 case erule::pwc: {
404 double distance = std::abs(x - getNode<effective_rule>(point));
405 double support = getSupport<effective_rule>(point);
406 isSupported = (distance <= 2.0 * support);
407 return (distance <= support) ? 1.0 : 0.0;
408 }
409 case erule::localp:
410 isSupported = true;
411 if (point == 0) {
412 return 1.0;
413 } else {
414 const double xn = scaleX<effective_rule>(point, x);
415 if (test_support(xn, point)) {
416 switch(max_order) {
417 case 1: return 1.0 - std::abs(xn);
418 case 2: return evalPWQuadratic<effective_rule>(point, xn);
419 case 3: return evalPWCubic<effective_rule>(point, xn);
420 default:
421 return evalPWPower<effective_rule>(max_order, point, xn);
422 };
423 } else {
424 isSupported = false;
425 return 0.0;
426 }
427 }
428 case erule::semilocalp:
429 isSupported = true;
430 switch(point) {
431 case 0: return 1.0;
432 case 1: return 0.5 * x * (x - 1.0);
433 case 2: return 0.5 * x * (x + 1.0);
434 default: {
435 const double xn = scaleX<effective_rule>(point, x);
436 if (test_support(xn, point)) {
437 switch(max_order) {
438 case 1: return 1.0 - std::abs(xn);
439 case 2: return evalPWQuadratic<effective_rule>(point, xn);
440 case 3: return evalPWCubic<effective_rule>(point, xn);
441 default:
442 return evalPWPower<effective_rule>(max_order, point, xn);
443 };
444 } else {
445 isSupported = false;
446 return 0.0;
447 }
448 }
449 };
450 case erule::localp0:
451 default: { // case erule::localpb:
452 const double xn = scaleX<effective_rule>(point, x);
453 if (test_support(xn, point)) {
454 isSupported = true;
455 switch(max_order) {
456 case 1: return 1.0 - std::abs(xn);
457 case 2: return evalPWQuadratic<effective_rule>(point, xn);
458 case 3: return evalPWCubic<effective_rule>(point, xn);
459 default:
460 return evalPWPower<effective_rule>(max_order, point, xn);
461 };
462 } else {
463 isSupported = false;
464 return 0.0;
465 }
466 }
467 };
468 }
469
470 template<erule effective_rule>
471 double evalRaw(int max_order, int point, double x) {
472 switch(effective_rule) {
473 case erule::pwc:
474 return (std::abs(x - getNode<effective_rule>(point)) <= getSupport<effective_rule>(point)) ? 1.0 : 0.0;
475 case erule::localp:
476 if (point == 0) {
477 return 1.0;
478 } else {
479 double xn = scaleX<effective_rule>(point, x);
480 if (std::abs(xn) <= 1.0) {
481 switch(max_order) {
482 case 1: return 1.0 - std::abs(xn);
483 case 2: return evalPWQuadratic<effective_rule>(point, xn);
484 case 3: return evalPWCubic<effective_rule>(point, xn);
485 default:
486 return evalPWPower<effective_rule>(max_order, point, xn);
487 };
488 } else {
489 return 0.0;
490 }
491 }
492 case erule::semilocalp:
493 switch(point) {
494 case 0: return 1.0;
495 case 1: return 0.5 * x * (x - 1.0);
496 case 2: return 0.5 * x * (x + 1.0);
497 default: {
498 double xn = scaleX<effective_rule>(point, x);
499 if (std::abs(xn) <= 1.0) {
500 switch(max_order) {
501 case 1: return 1.0 - std::abs(xn);
502 case 2: return evalPWQuadratic<effective_rule>(point, xn);
503 case 3: return evalPWCubic<effective_rule>(point, xn);
504 default:
505 return evalPWPower<effective_rule>(max_order, point, xn);
506 };
507 } else {
508 return 0.0;
509 }
510 }
511 };
512 case erule::localp0:
513 default: { // case erule::localpb:
514 double xn = scaleX<effective_rule>(point, x);
515 if (std::abs(xn) <= 1.0) {
516 switch(max_order) {
517 case 1: return 1.0 - std::abs(xn);
518 case 2: return evalPWQuadratic<effective_rule>(point, xn);
519 case 3: return evalPWCubic<effective_rule>(point, xn);
520 default:
521 return evalPWPower<effective_rule>(max_order, point, xn);
522 };
523 } else {
524 return 0.0;
525 }
526 }
527 };
528 }
529
530 template<erule effective_rule>
531 double diffPWQuadratic(int point, double x) {
532 if (effective_rule == erule::localp) {
533 switch(point) {
534 case 1: return -1.0;
535 case 2: return 1.0;
536 default:
537 return -2.0 * x;
538 };
539 } else if (effective_rule == erule::localpb) {
540 switch(point) {
541 case 0: return -1.0;
542 case 1: return 1.0;
543 default:
544 return -2.0 * x;
545 };
546 }
547 return -2.0 * x;
548 }
549 template<erule effective_rule>
550 double diffPWCubic(int point, double x) {
551 if (effective_rule == erule::localp) {
552 switch(point) {
553 case 0: return 0.0;
554 case 1: return -1.0;
555 case 2: return 1.0;
556 case 3:
557 case 4: return -2.0 * x;
558 default:
559 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
560 };
561 } else if (effective_rule == erule::localpb) {
562 switch(point) {
563 case 0: return -1.0;
564 case 1: return 1.0;
565 case 2: return -2.0 * x;
566 default:
567 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
568 };
569 } else if (effective_rule == erule::localp0) {
570 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));
571 }
572 return (point % 2 == 0) ? 1.0 / 3.0 - x * (x + 2.0) : -1.0 / 3.0 + x * (x - 2.0);
573 }
574
575 template<erule effrule>
576 double diffPWPower(int max_order, int point, double x) {
577 // use the cubic implementation until we have enough points
578 if (effrule == erule::localp and point <= 8) return diffPWCubic<effrule>(point, x);
579 if (effrule == erule::semilocalp and point <= 4) return diffPWCubic<effrule>(point, x);
580 if (effrule == erule::localpb and point <= 4) return diffPWCubic<effrule>(point, x);
581 if (effrule == erule::localp0 and point <= 2) return diffPWCubic<effrule>(point, x);
582 int level = getLevel<effrule>(point);
583 int max_ancestors = [&]()->int {
584 switch(effrule) {
585 case erule::pwc: return 0;
586 case erule::localp: return level - 2;
587 case erule::semilocalp:
588 case erule::localpb: return level - 1;
589 default: return level;
590 };
591 }();
592
593 if (max_order > 0) max_ancestors = std::min(max_ancestors, max_order - 2);
594
595 // This lambda captures most_turns and phantom_distance by reference, uses those as internal state variables, and on each
596 // call it returns the next normalized ancestor node.
597 int most_turns = 1;
598 double phantom_distance = 1.0;
599 auto update_and_get_next_node = [&]() {
600 most_turns *= 2;
601 phantom_distance = 2.0 * phantom_distance + 1.0;
602 int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns);
603 return (turns < most_turns / 2) ?
604 (phantom_distance - 2.0 * ((double) turns)) :
605 (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns)));
606 };
607
608 // This lambda is the inverse transform of update_and_get_next_node() above, and returns the previous descendant node.
609 auto rollback_and_get_prev_node = [&]() {
610 most_turns /= 2;
611 phantom_distance = 0.5 * (phantom_distance - 1.0);
612 int turns = (effrule == erule::localp0) ? ((point+1) % most_turns) : ((point-1) % most_turns);
613 return (turns < most_turns / 2) ?
614 (phantom_distance - 2.0 * ((double) turns)) :
615 (-phantom_distance + 2.0 * ((double) (most_turns - 1 - turns)));
616 };
617
618 // Does not include the additional factor (1-x) * (1+x) and the Lagrange coefficient.
619 std::vector<double> left_prods(max_ancestors);
620 left_prods[0] = 1.0;
621 double node = update_and_get_next_node();
622 double coeff = 1.0 / (-node);
623 for(int j=1; j<max_ancestors; j++) {
624 left_prods[j] = left_prods[j-1] * (x - node);
625 node = update_and_get_next_node();
626 coeff *= 1.0 / (-node);
627 }
628 double right_prod = 1.0;
629 double derivative = left_prods[max_ancestors-1];
630 for (int j=max_ancestors-2; j>=0; j--) {
631 right_prod *= x - node;
632 derivative += right_prod * left_prods[j];
633 node = rollback_and_get_prev_node();
634 }
635
636 // Adjust for the additional factor (1-x) * (1+x) and the Lagrange coefficient.
637 derivative = derivative * (1.0 - x) * (1.0 + x) + right_prod * (x - node) * (-2.0) * x;
638 derivative *= coeff;
639
640 return derivative;
641 }
642
643 template<erule effrule>
644 double scaleXclamp(int i, double x) {
645 return std::clamp(scaleX<effrule>(i, x), -1.0, 1.0);
646 }
647
648 template<erule effrule>
649 double diffRaw(int max_order, int point, double x) {
650 switch(effrule) {
651 case erule::pwc: return 0.0;
652 case erule::localp: {
653 if (point < 3)
654 return std::array<double, 3>{0, -1, 1}[point];
655 double xn = scaleXclamp<effrule>(point, x);
656 double an = scaleDiffX<effrule>(point);
657 switch(max_order) {
658 case 1: return (xn >= 0 ? -1.0 : 1.0) * an;
659 case 2: return an * diffPWQuadratic<effrule>(point, xn);
660 case 3: return an * diffPWCubic<effrule>(point, xn);
661 default:
662 return an * diffPWPower<effrule>(point, xn);
663 };
664 }
665 case erule::semilocalp:
666 switch(point) {
667 case 0: return 0.0;
668 case 1: return x - 0.5;
669 case 2: return x + 0.5;
670 default: {
671 double xn = scaleXclamp<effrule>(point, x);
672 double an = scaleDiffX<effrule>(point);
673 switch(max_order) {
674 case 2: return an * diffPWQuadratic<effrule>(point, xn);
675 case 3: return an * diffPWCubic<effrule>(point, xn);
676 default:
677 return an * diffPWPower<effrule>(point, xn);
678 };
679 }
680 };
681 case erule::localp0: {
682 double xn = scaleXclamp<effrule>(point, x);
683 double an = scaleDiffX<effrule>(point);
684 switch(max_order) {
685 case 1: return (xn >= 0 ? -1.0 : 1.0) * an;
686 case 2: return an * diffPWQuadratic<effrule>(point, xn);
687 case 3: return an * diffPWCubic<effrule>(point, xn);
688 default:
689 return an * diffPWPower<effrule>(point, xn);
690 };
691 }
692 default: { // case erule::localpb:
693 if (point < 2)
694 return std::array<double, 2>{-0.5, 0.5}[point];
695 double xn = scaleXclamp<effrule>(point, x);
696 double an = scaleDiffX<effrule>(point);
697 switch(max_order) {
698 case 1: return (xn >= 0 ? -1.0 : 1.0) * an;
699 case 2: return an * diffPWQuadratic<effrule>(point, xn);
700 case 3: return an * diffPWCubic<effrule>(point, xn);
701 default:
702 return an * diffPWPower<effrule>(point, xn);
703 };
704 }
705 };
706 }
707
708 template<erule effrule>
709 double diffSupport(int max_order, int point, double x) {
710 switch(effrule) {
711 case erule::pwc:
712 return 0.0;
713 case erule::localp: {
714 if (point < 3)
715 return std::array<double, 3>{0, -1, 1}[point];
716 double xn = scaleXclamp<effrule>(point, x);
717 double an = scaleDiffX<effrule>(point);
718 switch(max_order) {
719 case 1: return (xn >= 0 ? -1.0 : 1.0) * an;
720 case 2: return an * diffPWQuadratic<effrule>(point, xn);
721 case 3: return an * diffPWCubic<effrule>(point, xn);
722 default:
723 return an * diffPWPower<effrule>(max_order, point, xn);
724 };
725 }
726 case erule::semilocalp:
727 switch(point) {
728 case 0: return 0.0;
729 case 1: return x - 0.5;
730 case 2: return x + 0.5;
731 default: {
732 double xn = scaleXclamp<effrule>(point, x);
733 double an = scaleDiffX<effrule>(point);
734 switch(max_order) {
735 case 2: return an * diffPWQuadratic<effrule>(point, xn);
736 case 3: return an * diffPWCubic<effrule>(point, xn);
737 default:
738 return an * diffPWPower<effrule>(max_order, point, xn);
739 };
740 }
741 };
742 case erule::localp0: {
743 double xn = scaleXclamp<effrule>(point, x);
744 double an = scaleDiffX<effrule>(point);
745 switch(max_order) {
746 case 1: return ((xn >= 0 ? -1.0 : 1.0) * an);
747 case 2: return an * diffPWQuadratic<effrule>(point, xn);
748 case 3: return an * diffPWCubic<effrule>(point, xn);
749 default:
750 return an * diffPWPower<effrule>(max_order, point, xn);
751 };
752 }
753 default: { // case erule::localpb:
754 if (point < 2)
755 return std::array<double, 2>{-0.5, 0.5}[point];
756 double xn = scaleXclamp<effrule>(point, x);
757 double an = scaleDiffX<effrule>(point);
758 switch(max_order) {
759 case 1: return ((xn >= 0 ? -1.0 : 1.0) * an);
760 case 2: return an * diffPWQuadratic<effrule>(point, xn);
761 case 3: return an * diffPWCubic<effrule>(point, xn);
762 default:
763 return an * diffPWPower<effrule>(max_order, point, xn);
764 };
765 }
766 };
767 }
768
769 template<erule effrule>
770 double getArea(int max_order, int point, std::vector<double> const &w, std::vector<double> const &x) {
771 switch(effrule) {
772 case erule::pwc:
773 return 2.0 * getSupport<effrule>(point);
774 case erule::localp:
775 switch(point) {
776 case 0: return 2.0;
777 case 1:
778 case 2: return 0.5;
779 default:
780 switch(max_order) {
781 case 1: return getSupport<effrule>(point);
782 case 2:
783 case 3: return (4.0/3.0) * getSupport<effrule>(point);
784 default:
785 if (point <= 8) return (4.0/3.0) * getSupport<effrule>(point);
786 break;
787 };
788 break;
789 };
790 break;
791 case erule::semilocalp:
792 switch(point) {
793 case 0: return 2.0;
794 case 1:
795 case 2: return 1.0/3.0;
796 default:
797 switch(max_order) {
798 case 2:
799 case 3: return (4.0/3.0) * getSupport<effrule>(point);
800 default:
801 if (point <= 4) return (4.0/3.0) * getSupport<effrule>(point);
802 break;
803 };
804 break;
805 };
806 break;
807 case erule::localp0:
808 switch(max_order) {
809 case 1: return getSupport<effrule>(point);
810 case 2:
811 case 3: return (4.0/3.0) * getSupport<effrule>(point);
812 default:
813 if (point <= 2) return (4.0/3.0) * getSupport<effrule>(point);
814 break;
815 };
816 break;
817 default: // case erule::localpb:
818 if (point <= 1) {
819 return 1.0;
820 } else {
821 switch(max_order) {
822 case 1: return getSupport<effrule>(point);
823 case 2:
824 case 3: return (4.0/3.0) * getSupport<effrule>(point);
825 default:
826 if (point <= 4) return (4.0/3.0) * getSupport<effrule>(point);
827 break;
828 };
829 }
830 break;
831 };
832 double sum = 0.0;
833 for(size_t i=0; i<w.size(); i++) sum += w[i] * evalPWPower<effrule>(max_order, point, x[i]);
834 return sum * getSupport<effrule>(point);
835 }
836
837 template<erule effrule>
838 void van_matrix(int max_order, int num_rows, std::vector<int> &pntr, std::vector<int> &indx, std::vector<double> &vals) {
839 int max_level = getLevel<effrule>(num_rows);
840
841 if (effrule == erule::pwc) {
842 switch(num_rows) {
843 case 0:
844 indx = {0,};
845 vals = {1.0,};
846 pntr = {0, 1};
847 break;
848 case 1:
849 indx = {0, 0, 1};
850 vals = {1.0, 1.0, 1.0};
851 pntr = {0, 1, 3};
852 break;
853 case 2:
854 indx = {0, 0, 1, 0, 2};
855 vals = {1.0, 1.0, 1.0, 1.0, 1.0};
856 pntr = {0, 1, 3, 5};
857 break;
858 default:
859 indx.clear();
860 indx.reserve(num_rows * max_level);
861 vals.clear();
862 vals.reserve(num_rows * max_level);
863 pntr = std::vector<int>(num_rows + 1, 0);
864 for(auto i : std::array<int, 5>{0, 0, 1, 0, 2})
865 indx.push_back(i);
866 for(size_t i=0; i<5; i++)
867 vals.push_back(1.0);
868 pntr[1] = 1;
869 pntr[2] = 3;
870 {
871 std::vector<int> ancestors;
872 ancestors.reserve(max_level);
873 for(int r=3; r<num_rows; r++) {
874 pntr[r] = static_cast<int>(indx.size());
875 ancestors.clear();
876 int kid = r;
877 int dad = kid / 3;
878 while(dad != 0) {
879 if (dad % 2 == 0) {
880 if (kid != 3 * dad)
881 ancestors.push_back(dad);
882 } else {
883 if (kid != 3 * dad + 2)
884 ancestors.push_back(dad);
885 }
886 kid = dad;
887 dad = kid / 3;
888 }
889 indx.push_back(0);
890 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
891 indx.push_back(r);
892 for(size_t i=0; i<ancestors.size() + 2; i++)
893 vals.push_back(1.0);
894 }
895 pntr.back() = static_cast<int>(indx.size());
896 }
897 break;
898 };
899 } else if (effrule == erule::localp) {
900 switch(num_rows) {
901 case 0:
902 indx = {0,};
903 vals = {1.0,};
904 pntr = {0, 1};
905 break;
906 default:
907 indx.clear();
908 indx.reserve(num_rows * max_level);
909 vals.clear();
910 vals.reserve(num_rows * max_level);
911 pntr = std::vector<int>(num_rows + 1, 0);
912 indx.push_back(0);
913 vals.push_back(1.0);
914 {
915 std::vector<int> ancestors;
916 std::vector<double> ancestors_vals;
917 ancestors.reserve(max_level);
918 ancestors_vals.reserve(max_level);
919 for(int r=1; r<num_rows; r++) {
920 double x = getNode<effrule>(r);
921 pntr[r] = static_cast<int>(indx.size());
922 ancestors.clear();
923 ancestors_vals.clear();
924 int kid = r;
925 int dad = (kid + 1) / 2;
926 if (kid < 4) dad--;
927 while(dad != 0) {
928 ancestors.push_back(dad);
929 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
930 kid = dad;
931 dad = (kid + 1) / 2;
932 if (kid < 4) dad--;
933 }
934 indx.push_back(0);
935 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
936 indx.push_back(r);
937 vals.push_back(1.0);
938 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
939 vals.push_back(1.0);
940 }
941 pntr.back() = static_cast<int>(indx.size());
942 }
943 break;
944 };
945 } else if (effrule == erule::semilocalp) {
946 switch(num_rows) {
947 case 0:
948 indx = {0,};
949 vals = {1.0,};
950 pntr = {0, 1};
951 break;
952 case 1:
953 indx = {0, 0, 1};
954 vals = {1.0, 1.0, 1.0};
955 pntr = {0, 1, 3};
956 break;
957 case 2:
958 indx = {0, 0, 1, 0, 2};
959 vals = {1.0, 1.0, 1.0, 1.0, 1.0};
960 pntr = {0, 1, 3, 5};
961 break;
962 default:
963 indx.clear();
964 indx.reserve(num_rows * max_level);
965 vals.clear();
966 vals.reserve(num_rows * max_level);
967 pntr = std::vector<int>(num_rows + 1, 0);
968 for(auto i : std::array<int, 5>{0, 0, 1, 0, 2})
969 indx.push_back(i);
970 for(int i=0; i<5; i++)
971 vals.push_back(1.0);
972 pntr[1] = 1;
973 pntr[2] = 3;
974 {
975 std::vector<int> ancestors;
976 std::vector<double> ancestors_vals;
977 ancestors.reserve(max_level);
978 ancestors_vals.reserve(max_level);
979 for(int r=3; r<num_rows; r++) {
980 pntr[r] = static_cast<int>(indx.size());
981 double x = getNode<effrule>(r);
982 ancestors.clear();
983 ancestors_vals.clear();
984 int kid = r;
985 int dad = (kid + 1) / 2;
986 while(dad > 2) {
987 ancestors.push_back(dad);
988 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
989 kid = dad;
990 dad = (kid + 1) / 2;
991 }
992 indx.push_back(0);
993 indx.push_back(1);
994 indx.push_back(2);
995 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
996 indx.push_back(r);
997 vals.push_back(1.0);
998 vals.push_back( evalRaw<effrule>(max_order, 1, x) );
999 vals.push_back( evalRaw<effrule>(max_order, 2, x) );
1000 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
1001 vals.push_back(1.0);
1002 }
1003 pntr.back() = static_cast<int>(indx.size());
1004 }
1005 break;
1006 };
1007 } else if (effrule == erule::localp0) {
1008 switch(num_rows) {
1009 case 0:
1010 indx = {0,};
1011 vals = {1.0,};
1012 pntr = {0, 1};
1013 break;
1014 default:
1015 indx.clear();
1016 indx.reserve(num_rows * max_level);
1017 vals.clear();
1018 vals.reserve(num_rows * max_level);
1019 pntr = std::vector<int>(num_rows + 1, 0);
1020 indx.push_back(0);
1021 vals.push_back(1.0);
1022 {
1023 std::vector<int> ancestors;
1024 std::vector<double> ancestors_vals;
1025 ancestors.reserve(max_level);
1026 ancestors_vals.reserve(max_level);
1027 for(int r=1; r<num_rows; r++) {
1028 double x = getNode<effrule>(r);
1029 pntr[r] = static_cast<int>(indx.size());
1030 ancestors.clear();
1031 ancestors_vals.clear();
1032 int kid = r;
1033 int dad = (kid - 1) / 2;
1034 while(dad != 0) {
1035 ancestors.push_back(dad);
1036 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
1037 kid = dad;
1038 dad = (kid - 1) / 2;
1039 }
1040 indx.push_back(0);
1041 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
1042 indx.push_back(r);
1043 vals.push_back( evalRaw<effrule>(max_order, 0, x) );
1044 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
1045 vals.push_back(1.0);
1046 }
1047 pntr.back() = static_cast<int>(indx.size());
1048 }
1049 break;
1050 };
1051 } else { // if (effrule == erule::localpb) {
1052 switch(num_rows) {
1053 case 0:
1054 indx = {0,};
1055 vals = {1.0,};
1056 pntr = {0, 1};
1057 break;
1058 case 1:
1059 indx = {0, 1};
1060 vals = {1.0, 1.0};
1061 pntr = {0, 1, 2};
1062 break;
1063 default:
1064 indx.clear();
1065 indx.reserve(num_rows * max_level);
1066 vals.clear();
1067 vals.reserve(num_rows * max_level);
1068 pntr = std::vector<int>(num_rows + 1, 0);
1069 indx.push_back(0);
1070 indx.push_back(1);
1071 vals.push_back(1.0);
1072 vals.push_back(1.0);
1073 pntr[1] = 1;
1074 {
1075 std::vector<int> ancestors;
1076 std::vector<double> ancestors_vals;
1077 ancestors.reserve(max_level);
1078 ancestors_vals.reserve(max_level);
1079 for(int r=2; r<num_rows; r++) {
1080 pntr[r] = static_cast<int>(indx.size());
1081 double x = getNode<effrule>(r);
1082 ancestors.clear();
1083 ancestors_vals.clear();
1084 int kid = r;
1085 int dad = (kid + 1) / 2;
1086 while(dad > 1) {
1087 ancestors.push_back(dad);
1088 ancestors_vals.push_back( evalRaw<effrule>(max_order, dad, x) );
1089 kid = dad;
1090 dad = (kid + 1) / 2;
1091 }
1092 indx.push_back(0);
1093 indx.push_back(1);
1094 indx.insert(indx.end(), ancestors.rbegin(), ancestors.rend());
1095 indx.push_back(r);
1096 vals.push_back( evalRaw<effrule>(max_order, 0, x) );
1097 vals.push_back( evalRaw<effrule>(max_order, 1, x) );
1098 vals.insert(vals.end(), ancestors_vals.rbegin(), ancestors_vals.rend());
1099 vals.push_back(1.0);
1100 }
1101 pntr.back() = static_cast<int>(indx.size());
1102 }
1103 break;
1104 };
1105 }
1106 }
1107
1108} // namespace RuleLocal
1109
1110#endif // __TASMANIAN_DOXYGEN_SKIP
1111
1112}
1113
1114#endif
TasGrid::TypeOneDRule
TypeOneDRule
Used to specify the one dimensional family of rules that induces the sparse grid.
Definition tsgEnumerates.hpp:288
TasGrid::rule_localp
@ rule_localp
Nested rule with a hierarchy of uniformly distributed nodes and functions with compact support.
Definition tsgEnumerates.hpp:365
TasGrid::rule_localpb
@ rule_localpb
Variation of rule_localp focusing nodes on the boundary instead of the interior.
Definition tsgEnumerates.hpp:371
TasGrid::rule_localp0
@ rule_localp0
Variation of rule_localp assuming the model is zero at the domain boundary.
Definition tsgEnumerates.hpp:367
TasGrid::rule_semilocalp
@ rule_semilocalp
Variation of rule_localp using increased support in exchange for higher order basis (better for smoot...
Definition tsgEnumerates.hpp:369
TasGrid
Encapsulates the Tasmanian Sparse Grid module.
Definition TasmanianSparseGrid.hpp:68