1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24
25 #include <queso/1D1DFunction.h>
26 #include <queso/1DQuadrature.h>
27
28 namespace QUESO {
29
30 //*****************************************************
31 // Base 1D->1D class
32 //*****************************************************
Base1D1DFunction(double minDomainValue,double maxDomainValue)33 Base1D1DFunction::Base1D1DFunction(
34 double minDomainValue,
35 double maxDomainValue)
36 :
37 m_minDomainValue(minDomainValue),
38 m_maxDomainValue(maxDomainValue)
39 {
40 queso_require_less_msg(m_minDomainValue, m_maxDomainValue, "min >= max");
41 }
42
~Base1D1DFunction()43 Base1D1DFunction::~Base1D1DFunction()
44 {
45 }
46
47 double
minDomainValue()48 Base1D1DFunction::minDomainValue() const
49 {
50 return m_minDomainValue;
51 }
52
53 double
maxDomainValue()54 Base1D1DFunction::maxDomainValue() const
55 {
56 return m_maxDomainValue;
57 }
58
59 double
multiplyAndIntegrate(const Base1D1DFunction & func,unsigned int quadratureOrder,double * resultWithMultiplicationByTAsWell)60 Base1D1DFunction::multiplyAndIntegrate(const Base1D1DFunction& func, unsigned int quadratureOrder, double* resultWithMultiplicationByTAsWell) const
61 {
62 double value = 0.;
63
64 queso_not_implemented();
65
66 if (resultWithMultiplicationByTAsWell) { // Just to eliminate INTEL compiler warnings
67 func.value((double) quadratureOrder);
68 }
69
70 return value;
71 }
72
73 //*****************************************************
74 // Generic 1D->1D class
75 //*****************************************************
Generic1D1DFunction(double minDomainValue,double maxDomainValue,double (* valueRoutinePtr)(double domainValue,const void * routinesDataPtr),double (* derivRoutinePtr)(double domainValue,const void * routinesDataPtr),const void * routinesDataPtr)76 Generic1D1DFunction::Generic1D1DFunction(
77 double minDomainValue,
78 double maxDomainValue,
79 double (*valueRoutinePtr)(double domainValue, const void* routinesDataPtr),
80 double (*derivRoutinePtr)(double domainValue, const void* routinesDataPtr),
81 const void* routinesDataPtr)
82 :
83 Base1D1DFunction(minDomainValue,maxDomainValue),
84 m_valueRoutinePtr (valueRoutinePtr),
85 m_derivRoutinePtr (derivRoutinePtr),
86 m_routinesDataPtr (routinesDataPtr)
87 {
88 }
89
~Generic1D1DFunction()90 Generic1D1DFunction::~Generic1D1DFunction()
91 {
92 }
93
94 double
value(double domainValue)95 Generic1D1DFunction::value(double domainValue) const
96 {
97 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
98 std::cerr << "In Generic1D1DFunction::value()"
99 << ": requested x (" << domainValue
100 << ") is out of the interval (" << m_minDomainValue
101 << ", " << m_maxDomainValue
102 << ")"
103 << std::endl;
104 }
105
106 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
107
108 return (*m_valueRoutinePtr)(domainValue,m_routinesDataPtr);
109 }
110
111 double
deriv(double domainValue)112 Generic1D1DFunction::deriv(double domainValue) const
113 {
114 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
115 std::cerr << "In Generic1D1DFunction::deriv()"
116 << ": requested x (" << domainValue
117 << ") is out of the interval (" << m_minDomainValue
118 << ", " << m_maxDomainValue
119 << ")"
120 << std::endl;
121 }
122
123 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
124
125 return (*m_derivRoutinePtr)(domainValue,m_routinesDataPtr);
126 }
127
128 //*****************************************************
129 // Constant 1D->1D class
130 //*****************************************************
Constant1D1DFunction(double minDomainValue,double maxDomainValue,double constantValue)131 Constant1D1DFunction::Constant1D1DFunction(
132 double minDomainValue,
133 double maxDomainValue,
134 double constantValue)
135 :
136 Base1D1DFunction(minDomainValue,maxDomainValue),
137 m_constantValue (constantValue)
138 {
139 }
140
~Constant1D1DFunction()141 Constant1D1DFunction::~Constant1D1DFunction()
142 {
143 }
144
145 double
value(double domainValue)146 Constant1D1DFunction::value(double domainValue) const
147 {
148 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
149 std::cerr << "In Constant1D1DFunction::value()"
150 << ": requested x (" << domainValue
151 << ") is out of the interval (" << m_minDomainValue
152 << ", " << m_maxDomainValue
153 << ")"
154 << std::endl;
155 }
156
157 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
158
159 return m_constantValue;
160 }
161
162 double
deriv(double domainValue)163 Constant1D1DFunction::deriv(double domainValue) const
164 {
165 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
166 std::cerr << "In Constant1D1DFunction::deriv()"
167 << ": requested x (" << domainValue
168 << ") is out of the interval (" << m_minDomainValue
169 << ", " << m_maxDomainValue
170 << ")"
171 << std::endl;
172 }
173
174 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
175
176 return 0.;
177 }
178
179 //*****************************************************
180 // Linear 1D->1D class
181 //*****************************************************
Linear1D1DFunction(double minDomainValue,double maxDomainValue,double referenceDomainValue,double referenceImageValue,double rateValue)182 Linear1D1DFunction::Linear1D1DFunction(
183 double minDomainValue,
184 double maxDomainValue,
185 double referenceDomainValue,
186 double referenceImageValue,
187 double rateValue)
188 :
189 Base1D1DFunction(minDomainValue,maxDomainValue),
190 m_referenceDomainValue (referenceDomainValue),
191 m_referenceImageValue (referenceImageValue),
192 m_rateValue (rateValue)
193 {
194 }
195
~Linear1D1DFunction()196 Linear1D1DFunction::~Linear1D1DFunction()
197 {
198 }
199
200 double
value(double domainValue)201 Linear1D1DFunction::value(double domainValue) const
202 {
203 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
204 std::cerr << "In Linear1D1DFunction::value()"
205 << ": requested x (" << domainValue
206 << ") is out of the interval (" << m_minDomainValue
207 << ", " << m_maxDomainValue
208 << ")"
209 << std::endl;
210 }
211
212 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
213
214 double imageValue = m_referenceImageValue + m_rateValue*(domainValue - m_referenceDomainValue);
215
216 return imageValue;
217 }
218
219 double
deriv(double domainValue)220 Linear1D1DFunction::deriv(double domainValue) const
221 {
222 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
223 std::cerr << "In Linear1D1DFunction::deriv()"
224 << ": requested x (" << domainValue
225 << ") is out of the interval (" << m_minDomainValue
226 << ", " << m_maxDomainValue
227 << ")"
228 << std::endl;
229 }
230
231 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
232
233 return m_rateValue;
234 }
235
236 //*****************************************************
237 // PiecewiseLinear 1D->1D class
238 //*****************************************************
PiecewiseLinear1D1DFunction(double minDomainValue,double maxDomainValue,const std::vector<double> & referenceDomainValues,double referenceImageValue0,const std::vector<double> & rateValues)239 PiecewiseLinear1D1DFunction::PiecewiseLinear1D1DFunction(
240 double minDomainValue,
241 double maxDomainValue,
242 const std::vector<double>& referenceDomainValues,
243 double referenceImageValue0,
244 const std::vector<double>& rateValues)
245 :
246 Base1D1DFunction(minDomainValue,maxDomainValue),
247 m_numRefValues (referenceDomainValues.size()),
248 m_referenceDomainValues(referenceDomainValues),
249 m_rateValues (rateValues)
250 {
251 queso_require_not_equal_to_msg(m_numRefValues, 0, "num ref values = 0");
252
253 queso_require_equal_to_msg(m_numRefValues, rateValues.size(), "num rate values is inconsistent");
254
255 for (unsigned int i = 1; i < m_numRefValues; ++i) { // Yes, from '1'
256 queso_require_greater_msg(m_referenceDomainValues[i], m_referenceDomainValues[i-1], "reference domain values are inconsistent");
257 }
258
259 m_referenceImageValues.clear();
260 m_referenceImageValues.resize(m_numRefValues,0.);
261 m_referenceImageValues[0] = referenceImageValue0;
262 for (unsigned int i = 1; i < m_numRefValues; ++i) { // Yes, from '1'
263 m_referenceImageValues[i] = m_referenceImageValues[i-1] + m_rateValues[i-1]*(m_referenceDomainValues[i] - m_referenceDomainValues[i-1]);
264 }
265
266 if (false) { // For debug only
267 std::cout << "In PiecewiseLinear1D1DFunction::constructor():"
268 << std::endl;
269 for (unsigned int i = 0; i < m_numRefValues; ++i) {
270 std::cout << "i = " << i
271 << ", m_referenceDomainValues[i] = " << m_referenceDomainValues[i]
272 << ", m_referenceImageValues[i] = " << m_referenceImageValues[i]
273 << ", m_rateValues[i] = " << m_rateValues[i]
274 << std::endl;
275 }
276 }
277 }
278
~PiecewiseLinear1D1DFunction()279 PiecewiseLinear1D1DFunction::~PiecewiseLinear1D1DFunction()
280 {
281 m_rateValues.clear();
282 m_referenceImageValues.clear();
283 m_referenceDomainValues.clear();
284 }
285
286 double
value(double domainValue)287 PiecewiseLinear1D1DFunction::value(double domainValue) const
288 {
289 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
290 std::cerr << "In PiecewiseLinear1D1DFunction::value()"
291 << ": requested x (" << domainValue
292 << ") is out of the interval (" << m_minDomainValue
293 << ", " << m_maxDomainValue
294 << ")"
295 << std::endl;
296 }
297
298 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
299
300 unsigned int i = 0;
301 if (m_numRefValues == 1) {
302 // Nothing else to do
303 }
304 else {
305 bool referenceDomainValueFound = false;
306 while (referenceDomainValueFound == false) {
307 if (domainValue < m_referenceDomainValues[i+1]) {
308 referenceDomainValueFound = true;
309 }
310 else {
311 ++i;
312 if (i == (m_numRefValues-1)) {
313 referenceDomainValueFound = true;
314 }
315 }
316 }
317 }
318 double imageValue = m_referenceImageValues[i] + m_rateValues[i]*(domainValue - m_referenceDomainValues[i]);
319 if (false) { // For debug only
320 std::cout << "In PiecewiseLinear1D1DFunction::value()"
321 << ": domainValue = " << domainValue
322 << ", i = " << i
323 << ", m_referenceDomainValues[i] = " << m_referenceDomainValues[i]
324 << ", m_referenceImageValues[i] = " << m_referenceImageValues[i]
325 << ", m_rateValues[i] = " << m_rateValues[i]
326 << ", imageValue = " << imageValue
327 << std::endl;
328 }
329
330 return imageValue;
331 }
332
333 double
deriv(double domainValue)334 PiecewiseLinear1D1DFunction::deriv(double domainValue) const
335 {
336 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
337 std::cerr << "In PiecewiseLinear1D1DFunction::deriv()"
338 << ": requested x (" << domainValue
339 << ") is out of the interval (" << m_minDomainValue
340 << ", " << m_maxDomainValue
341 << ")"
342 << std::endl;
343 }
344
345 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
346
347 unsigned int i = 0;
348 if (m_numRefValues == 1) {
349 // Nothing else to do
350 }
351 else {
352 bool referenceDomainValueFound = false;
353 while (referenceDomainValueFound == false) {
354 if (domainValue < m_referenceDomainValues[i+1]) {
355 referenceDomainValueFound = true;
356 }
357 else {
358 ++i;
359 queso_require_less_equal_msg(i, m_numRefValues, "too big 'i'");
360 }
361 }
362 }
363
364 return m_rateValues[i];
365 }
366
367 //*****************************************************
368 // Quadratic 1D->1D class
369 //*****************************************************
Quadratic1D1DFunction(double minDomainValue,double maxDomainValue,double a,double b,double c)370 Quadratic1D1DFunction::Quadratic1D1DFunction(
371 double minDomainValue,
372 double maxDomainValue,
373 double a,
374 double b,
375 double c)
376 :
377 Base1D1DFunction(minDomainValue,maxDomainValue),
378 m_a (a),
379 m_b (b),
380 m_c (c)
381 {
382 }
383
~Quadratic1D1DFunction()384 Quadratic1D1DFunction::~Quadratic1D1DFunction()
385 {
386 }
387
388 double
value(double domainValue)389 Quadratic1D1DFunction::value(double domainValue) const
390 {
391 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
392 std::cerr << "In Quadratic1D1DFunction::value()"
393 << ": requested x (" << domainValue
394 << ") is out of the interval (" << m_minDomainValue
395 << ", " << m_maxDomainValue
396 << ")"
397 << std::endl;
398 }
399
400 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
401
402 double imageValue = m_a*domainValue*domainValue + m_b*domainValue + m_c;
403
404 return imageValue;
405 }
406
407 double
deriv(double domainValue)408 Quadratic1D1DFunction::deriv(double domainValue) const
409 {
410 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
411 std::cerr << "In Quadratic1D1DFunction::deriv()"
412 << ": requested x (" << domainValue
413 << ") is out of the interval (" << m_minDomainValue
414 << ", " << m_maxDomainValue
415 << ")"
416 << std::endl;
417 }
418
419 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
420
421 return 2.*m_a*domainValue + m_b;
422 }
423
424 //*****************************************************
425 // Sampled 1D->1D class
426 //*****************************************************
Sampled1D1DFunction()427 Sampled1D1DFunction::Sampled1D1DFunction()
428 :
429 Base1D1DFunction(-INFINITY,INFINITY)
430 {
431
432 // UQ_UNAVAILABLE_RANK,
433 // "SampledD1DFunction::deriv()",
434 // "invalid constructor");
435 }
436
Sampled1D1DFunction(const std::vector<double> & domainValues,const std::vector<double> & imageValues)437 Sampled1D1DFunction::Sampled1D1DFunction(
438 const std::vector<double>& domainValues,
439 const std::vector<double>& imageValues)
440 :
441 Base1D1DFunction(domainValues[0],domainValues[domainValues.size()-1]),
442 m_domainValues (domainValues.size(),0.),
443 m_imageValues (imageValues.size(), 0.)
444 {
445 unsigned int tmpSize = m_domainValues.size();
446 for (unsigned int i = 0; i < tmpSize; ++i) {
447 m_domainValues[i] = domainValues[i];
448 m_imageValues [i] = imageValues [i];
449 }
450 }
451
~Sampled1D1DFunction()452 Sampled1D1DFunction::~Sampled1D1DFunction()
453 {
454 }
455
456 double
value(double domainValue)457 Sampled1D1DFunction::value(double domainValue) const
458 {
459 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
460 std::cerr << "In Sampled1D1DFunction::value()"
461 << ": requested x (" << domainValue
462 << ") is out of the interval (" << m_minDomainValue
463 << ", " << m_maxDomainValue
464 << ")"
465 << std::endl;
466 }
467
468 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
469
470 double returnValue = 0.;
471
472 unsigned int tmpSize = m_domainValues.size();
473 //std::cout << "In Sampled1D1DFunction::value()"
474 // << ": domainValue = " << domainValue
475 // << ", tmpSize = " << tmpSize
476 // << ", m_domainValues[0] = " << m_domainValues[0]
477 // << ", m_domainValues[max] = " << m_domainValues[tmpSize-1]
478 // << std::endl;
479
480 queso_require_not_equal_to_msg(tmpSize, 0, "m_domainValues.size() = 0");
481
482 queso_require_greater_equal_msg(domainValue, m_domainValues[0], "domainValue < m_domainValues[0]");
483
484 queso_require_greater_equal_msg(m_domainValues[tmpSize-1], domainValue, "m_domainValues[max] < domainValue");
485
486 unsigned int i = 0;
487 for (i = 0; i < tmpSize; ++i) {
488 if (domainValue <= m_domainValues[i]) break;
489 }
490
491 if (domainValue == m_domainValues[i]) {
492 //if (domainValueWasMatchedExactly) *domainValueWasMatchedExactly = true;
493 returnValue = m_imageValues[i];
494 }
495 else {
496 //if (domainValueWasMatchedExactly) *domainValueWasMatchedExactly = false;
497 double ratio = (domainValue - m_domainValues[i-1])/(m_domainValues[i]-m_domainValues[i-1]);
498 returnValue = m_imageValues[i-1] + ratio * (m_imageValues[i]-m_imageValues[i-1]);
499 }
500
501 return returnValue;
502 }
503
504 double
deriv(double domainValue)505 Sampled1D1DFunction::deriv(double domainValue) const
506 {
507 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
508 std::cerr << "In Sampled1D1DFunction::deriv()"
509 << ": requested x (" << domainValue
510 << ") is out of the interval (" << m_minDomainValue
511 << ", " << m_maxDomainValue
512 << ")"
513 << std::endl;
514 }
515
516 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
517
518 queso_error_msg("this function makes no sense for this class");
519 return 0.;
520 }
521
522 const std::vector<double>&
domainValues()523 Sampled1D1DFunction::domainValues() const
524 {
525 return m_domainValues;
526 }
527
528 const std::vector<double>&
imageValues()529 Sampled1D1DFunction::imageValues() const
530 {
531 return m_imageValues;
532 }
533
534 bool
domainValueMatchesExactly(double domainValue)535 Sampled1D1DFunction::domainValueMatchesExactly(double domainValue) const
536 {
537 bool result = false;
538
539 unsigned int tmpSize = m_domainValues.size();
540 for (unsigned int i = 0; i < tmpSize; ++i) {
541 if (domainValue <= m_domainValues[i]) {
542 result = (domainValue == m_domainValues[i]);
543 break;
544 }
545 }
546
547 return result;
548 }
549
550 void
set(const std::vector<double> & domainValues,const std::vector<double> & imageValues)551 Sampled1D1DFunction::set(
552 const std::vector<double>& domainValues,
553 const std::vector<double>& imageValues)
554 {
555 m_domainValues.clear();
556 m_imageValues.clear();
557
558 unsigned int tmpSize = domainValues.size();
559 m_minDomainValue = domainValues[0];
560 m_maxDomainValue = domainValues[tmpSize-1];
561
562 m_domainValues.resize(tmpSize,0.);
563 m_imageValues.resize(tmpSize,0.);
564 for (unsigned int i = 0; i < tmpSize; ++i) {
565 m_domainValues[i] = domainValues[i];
566 m_imageValues [i] = imageValues [i];
567 }
568
569 return;
570 }
571
572 void
printForMatlab(const BaseEnvironment & env,std::ofstream & ofsvar,const std::string & prefixName)573 Sampled1D1DFunction::printForMatlab(
574 const BaseEnvironment& env,
575 std::ofstream& ofsvar,
576 const std::string& prefixName) const
577 {
578 unsigned int tmpSize = m_domainValues.size();
579 if (tmpSize == 0) {
580 tmpSize = 1;
581 ofsvar << "\n" << prefixName << "Time_sub" << env.subIdString() << " = zeros(" << tmpSize << ",1);"
582 << "\n" << prefixName << "Value_sub" << env.subIdString() << " = zeros(" << tmpSize << ",1);";
583 }
584 else {
585 ofsvar << "\n" << prefixName << "Time_sub" << env.subIdString() << " = zeros(" << tmpSize << ",1);"
586 << "\n" << prefixName << "Value_sub" << env.subIdString() << " = zeros(" << tmpSize << ",1);";
587 for (unsigned int i = 0; i < tmpSize; ++i) {
588 ofsvar << "\n" << prefixName << "Time_sub" << env.subIdString() << "(" << i+1 << ",1) = " << m_domainValues[i] << ";"
589 << "\n" << prefixName << "Value_sub" << env.subIdString() << "(" << i+1 << ",1) = " << m_imageValues[i] << ";";
590 }
591 }
592
593 return;
594 }
595
596 //*****************************************************
597 // 'ScalarTimesFunc' 1D->1D class
598 //*****************************************************
ScalarTimesFunc1D1DFunction(double scalar,const Base1D1DFunction & func)599 ScalarTimesFunc1D1DFunction::ScalarTimesFunc1D1DFunction(
600 double scalar,
601 const Base1D1DFunction& func)
602 :
603 Base1D1DFunction(func.minDomainValue(),func.maxDomainValue()),
604 m_scalar(scalar),
605 m_func (func)
606 {
607 }
608
~ScalarTimesFunc1D1DFunction()609 ScalarTimesFunc1D1DFunction::~ScalarTimesFunc1D1DFunction()
610 {
611 }
612
613 double
value(double domainValue)614 ScalarTimesFunc1D1DFunction::value(double domainValue) const
615 {
616 double value = 0.;
617
618 value += m_scalar*m_func.value(domainValue);
619
620 return value;
621 }
622
623 double
deriv(double domainValue)624 ScalarTimesFunc1D1DFunction::deriv(double domainValue) const
625 {
626 double value = 0.;
627
628 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
629 std::cerr << "In ScalarTimes1D1DFunction::deriv()"
630 << ": requested x (" << domainValue
631 << ") is out of the interval (" << m_minDomainValue
632 << ", " << m_maxDomainValue
633 << ")"
634 << std::endl;
635 }
636
637 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
638
639 queso_not_implemented();
640
641 return value;
642 }
643
644 //*****************************************************
645 // 'FuncTimesFunc' 1D->1D class
646 //*****************************************************
FuncTimesFunc1D1DFunction(const Base1D1DFunction & func1,const Base1D1DFunction & func2)647 FuncTimesFunc1D1DFunction::FuncTimesFunc1D1DFunction(
648 const Base1D1DFunction& func1,
649 const Base1D1DFunction& func2)
650 :
651 Base1D1DFunction(std::max(func1.minDomainValue(),func2.minDomainValue()),std::min(func1.maxDomainValue(),func2.maxDomainValue())),
652 m_func1(func1),
653 m_func2(func2)
654 {
655 }
656
~FuncTimesFunc1D1DFunction()657 FuncTimesFunc1D1DFunction::~FuncTimesFunc1D1DFunction()
658 {
659 }
660
661 double
value(double domainValue)662 FuncTimesFunc1D1DFunction::value(double domainValue) const
663 {
664 double value = 0.;
665
666 value += m_func1.value(domainValue)*m_func2.value(domainValue);
667
668 return value;
669 }
670
671 double
deriv(double domainValue)672 FuncTimesFunc1D1DFunction::deriv(double domainValue) const
673 {
674 double value = 0.;
675
676 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
677 std::cerr << "In FuncTimes1D1DFunction::deriv()"
678 << ": requested x (" << domainValue
679 << ") is out of the interval (" << m_minDomainValue
680 << ", " << m_maxDomainValue
681 << ")"
682 << std::endl;
683 }
684
685 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
686
687 queso_not_implemented();
688
689 return value;
690 }
691
692 //*****************************************************
693 // 'FuncPlusFunc' 1D->1D class
694 //*****************************************************
FuncPlusFunc1D1DFunction(const Base1D1DFunction & func1,const Base1D1DFunction & func2)695 FuncPlusFunc1D1DFunction::FuncPlusFunc1D1DFunction(
696 const Base1D1DFunction& func1,
697 const Base1D1DFunction& func2)
698 :
699 Base1D1DFunction(std::max(func1.minDomainValue(),func2.minDomainValue()),std::min(func1.maxDomainValue(),func2.maxDomainValue())),
700 m_func1(func1),
701 m_func2(func2)
702 {
703 }
704
~FuncPlusFunc1D1DFunction()705 FuncPlusFunc1D1DFunction::~FuncPlusFunc1D1DFunction()
706 {
707 }
708
709 double
value(double domainValue)710 FuncPlusFunc1D1DFunction::value(double domainValue) const
711 {
712 double value = 0.;
713
714 value += m_func1.value(domainValue) + m_func2.value(domainValue);
715
716 return value;
717 }
718
719 double
deriv(double domainValue)720 FuncPlusFunc1D1DFunction::deriv(double domainValue) const
721 {
722 double value = 0.;
723
724 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
725 std::cerr << "In FuncPlus1D1DFunction::deriv()"
726 << ": requested x (" << domainValue
727 << ") is out of the interval (" << m_minDomainValue
728 << ", " << m_maxDomainValue
729 << ")"
730 << std::endl;
731 }
732
733 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
734
735 queso_not_implemented();
736
737 return value;
738 }
739
740 //*****************************************************
741 // Lagrange Polynomial 1D->1D class
742 //*****************************************************
LagrangePolynomial1D1DFunction(const std::vector<double> & positionValues,const std::vector<double> * functionValues)743 LagrangePolynomial1D1DFunction::LagrangePolynomial1D1DFunction(
744 const std::vector<double>& positionValues,
745 const std::vector<double>* functionValues)
746 :
747 Base1D1DFunction(-INFINITY,INFINITY),
748 m_positionValues(positionValues),
749 m_functionValues(positionValues.size(),1.)
750 {
751 if (functionValues) {
752 queso_require_equal_to_msg(m_positionValues.size(), functionValues->size(), "invalid input");
753 m_functionValues = *functionValues;
754 }
755 }
756
~LagrangePolynomial1D1DFunction()757 LagrangePolynomial1D1DFunction::~LagrangePolynomial1D1DFunction()
758 {
759 }
760
761 double
value(double domainValue)762 LagrangePolynomial1D1DFunction::value(double domainValue) const
763 {
764 double value = 0.;
765
766 for (unsigned int k = 0; k < m_positionValues.size(); ++k) {
767 double scaleFactor = 1.;
768 double posK = m_positionValues[k];
769 for (unsigned int j = 0; j < m_positionValues.size(); ++j) {
770 if (j != k) {
771 double posJ = m_positionValues[j];
772 scaleFactor *= (domainValue-posJ)/(posK-posJ);
773 }
774 }
775
776 //if (m_env.subDisplayFile()) {
777 // *m_env.subDisplayFile() << "In sddTG<K_V,K_M>::lagrange()"
778 // << ": k = " << k
779 // << ", scaleFactor = " << scaleFactor
780 // << std::endl;
781 //}
782
783 value += scaleFactor * m_functionValues[k];
784
785 //if (m_env.subDisplayFile()) {
786 // *m_env.subDisplayFile() << "In sddTG<K_V,K_M>::lagrange()"
787 // << ": k = " << k
788 // << ", value = " << value
789 // << std::endl;
790 //}
791 }
792
793 return value;
794 }
795
796 double
deriv(double domainValue)797 LagrangePolynomial1D1DFunction::deriv(double domainValue) const
798 {
799 double value = 0.;
800
801 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
802 std::cerr << "In LagrangePolynomial1D1DFunction::deriv()"
803 << ": requested x (" << domainValue
804 << ") is out of the interval (" << m_minDomainValue
805 << ", " << m_maxDomainValue
806 << ")"
807 << std::endl;
808 }
809
810 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
811
812 queso_not_implemented();
813
814 return value;
815 }
816
817 //*****************************************************
818 // Lagrange Basis 1D->1D class
819 //*****************************************************
LagrangeBasis1D1DFunction(const std::vector<double> & positionValues,unsigned int basisIndex)820 LagrangeBasis1D1DFunction::LagrangeBasis1D1DFunction(
821 const std::vector<double>& positionValues,
822 unsigned int basisIndex)
823 :
824 Base1D1DFunction(-INFINITY,INFINITY),
825 m_positionValues(positionValues),
826 m_basisIndex (basisIndex)
827 {
828 queso_require_less_msg(m_basisIndex, m_positionValues.size(), "invalid input");
829 }
830
~LagrangeBasis1D1DFunction()831 LagrangeBasis1D1DFunction::~LagrangeBasis1D1DFunction()
832 {
833 }
834
835 double
value(double domainValue)836 LagrangeBasis1D1DFunction::value(double domainValue) const
837 {
838 double scaleFactor = 1.;
839
840 unsigned int k = m_basisIndex;
841 double posK = m_positionValues[k];
842 for (unsigned int j = 0; j < m_positionValues.size(); ++j) {
843 if (j != k) {
844 double posJ = m_positionValues[j];
845 scaleFactor *= (domainValue-posJ)/(posK-posJ);
846 }
847 }
848
849 return scaleFactor;
850 }
851
852 double
deriv(double domainValue)853 LagrangeBasis1D1DFunction::deriv(double domainValue) const
854 {
855 double value = 0.;
856
857 if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
858 std::cerr << "In LagrangeBasis1D1DFunction::deriv()"
859 << ": requested x (" << domainValue
860 << ") is out of the interval (" << m_minDomainValue
861 << ", " << m_maxDomainValue
862 << ")"
863 << std::endl;
864 }
865
866 queso_require_msg(!((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)), "x out of range");
867
868 queso_not_implemented();
869
870 return value;
871 }
872
873 template <class T>
874 double
SubF1F2Gaussian2dKdeIntegral(const ScalarSequence<T> & scalarSeq1,const ScalarSequence<T> & scalarSeq2,unsigned int initialPos,double scaleValue1,double scaleValue2,const Base1D1DFunction & func1,const Base1D1DFunction & func2,unsigned int quadratureOrder)875 SubF1F2Gaussian2dKdeIntegral(const ScalarSequence<T>& scalarSeq1,
876 const ScalarSequence<T>& scalarSeq2,
877 unsigned int initialPos,
878 double scaleValue1,
879 double scaleValue2,
880 const Base1D1DFunction& func1,
881 const Base1D1DFunction& func2,
882 unsigned int quadratureOrder)
883 {
884 double resultValue = 0.;
885
886 queso_require_equal_to_msg(initialPos, 0, "not implemented yet for initialPos != 0");
887 queso_require_equal_to_msg(scalarSeq1.subSequenceSize(), scalarSeq2.subSequenceSize(), "different sizes");
888
889 GaussianHermite1DQuadrature quadObj(0.,1.,quadratureOrder);
890 const std::vector<double>& quadPositions = quadObj.positions();
891 const std::vector<double>& quadWeights = quadObj.weights ();
892 queso_require_equal_to_msg(quadPositions.size(), quadWeights.size(), "quadObj has invalid state");
893
894 unsigned int numQuadraturePositions = quadPositions.size();
895 unsigned int dataSize = scalarSeq1.subSequenceSize();
896 for (unsigned int k = 0; k < dataSize; ++k) {
897 double value1 = 0.;
898 double value2 = 0.;
899 double x1k = scalarSeq1[k];
900 double x2k = scalarSeq2[k];
901 for (unsigned int j = 0; j < numQuadraturePositions; ++j) {
902 value1 += func1.value(scaleValue1*quadPositions[j]+x1k)*quadWeights[j];
903 value2 += func2.value(scaleValue2*quadPositions[j]+x2k)*quadWeights[j];
904 }
905 resultValue += value1*value2;
906 }
907 resultValue *= 1./(2.*M_PI)/((double) dataSize);
908
909 return resultValue;
910 }
911
912 } // End namespace QUESO
913