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 <example_likelihood.h>
26 #include <cmath>
27
28 static unsigned int likelihoodCounter = 0;
29
likelihoodRoutine(const QUESO::GslVector & paramValues,const QUESO::GslVector * paramDirection,const void * functionDataPtr,QUESO::GslVector * gradVector,QUESO::GslMatrix * hessianMatrix,QUESO::GslVector * hessianEffect)30 double likelihoodRoutine(
31 const QUESO::GslVector& paramValues,
32 const QUESO::GslVector* paramDirection,
33 const void* functionDataPtr,
34 QUESO::GslVector* gradVector,
35 QUESO::GslMatrix* hessianMatrix,
36 QUESO::GslVector* hessianEffect)
37 {
38 //const uqGslVector& meanVector =
39 // *((likelihoodRoutine_DataType *) functionDataPtr)->meanVector;
40 //const uqGslMatrix& covMatrix =
41 // *((likelihoodRoutine_DataType *) functionDataPtr)->covMatrix;
42
43 //uqGslVector diffVec(paramValues - meanVector);
44
45 //return scalarProduct(diffVec, covMatrix.invertMultiply(diffVec));
46
47 likelihoodCounter++;
48
49 if (paramDirection ||
50 functionDataPtr ||
51 gradVector ||
52 hessianMatrix ||
53 hessianEffect) {}; // just to remove compiler warning
54
55 double x = paramValues[0];
56
57 double mean1 = 10.;
58 double sigma1 = 1.;
59 double y1 = (x-mean1)*(x-mean1)/(2.*sigma1*sigma1);
60 //double z1 = (1./sigma1/sqrt(2*M_PI))*exp(-y1);
61 double ln_z1 = log(1./sigma1/M_PI/sqrt(2*M_PI)) - y1;
62
63 double mean2 = 100.;
64 double sigma2 = 5.;
65 double y2 = (x-mean2)*(x-mean2)/(2.*sigma2*sigma2);
66 //double z2 = (1./sigma2/sqrt(2*M_PI))*exp(-y2);
67 double ln_z2 = log(1./sigma2/M_PI/sqrt(2*M_PI)) - y2;
68
69
70 //double resultValue = log((z1+2.*z2)/3.);
71 double ln_zMin = std::min(log(1./3.)+ln_z1, log(2./3.)+ln_z2);
72 double ln_zMax = std::max(log(1./3.)+ln_z1, log(2./3.)+ln_z2);
73 double resultValue = log(1. + exp(ln_zMin - ln_zMax)) + ln_zMax;
74
75 if (resultValue == -INFINITY) {
76 #if 1
77 std::cerr << "WARNING In likelihoodRoutine"
78 << ", likelihoodCounter =" << likelihoodCounter
79 << ", fullRank " << paramValues.env().fullRank()
80 << ", subEnvironment " << paramValues.env().subId()
81 << ", subRank " << paramValues.env().subRank()
82 << ", inter0Rank " << paramValues.env().inter0Rank()
83 << ": x = " << x
84 << ", ln_z1 = " << ln_z1
85 << ", ln_z2 = " << ln_z2
86 << ", resultValue = " << resultValue
87 << std::endl;
88 #endif
89 resultValue = -520.;
90 }
91
92 double returnValue = resultValue;
93
94 if (paramValues.env().exceptionalCircumstance()) {
95 if ((paramValues.env().subDisplayFile() ) &&
96 (paramValues.env().displayVerbosity() > 0)) { // detailed output debug
97 *paramValues.env().subDisplayFile() << "Leaving likelihood function"
98 << ": paramValues = " << paramValues
99 << ", returnValue = " << returnValue
100 << std::endl;
101 }
102 }
103
104 return returnValue;
105 }
106