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/GslVector.h>
26 #include <queso/GslMatrix.h>
27 #include <queso/VectorSet.h>
28 #include <queso/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h>
29
30 namespace QUESO {
31
32 template<class V, class M>
GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(const char * prefix,const VectorSet<V,M> & domainSet,const V & observations,const GslBlockMatrix & covariance)33 GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients<V, M>::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(
34 const char * prefix, const VectorSet<V, M> & domainSet,
35 const V & observations, const GslBlockMatrix & covariance)
36 : LikelihoodBase<V, M>(prefix, domainSet, observations),
37 m_covariance(covariance)
38 {
39 unsigned int totalDim = 0;
40
41 for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) {
42 totalDim += this->m_covariance.getBlock(i).numRowsLocal();
43 }
44
45 if (totalDim != observations.sizeLocal()) {
46 queso_error_msg("Covariance matrix not same dimension as observation vector");
47 }
48 }
49
50 template<class V, class M>
~GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients()51 GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients<V, M>::~GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients()
52 {
53 }
54
55 template<class V, class M>
56 double
lnValue(const V & domainVector)57 GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients<V, M>::lnValue(const V & domainVector) const
58 {
59 V modelOutput(this->m_observations, 0, 0); // At least it's not a copy
60 V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy
61
62 this->evaluateModel(domainVector, modelOutput);
63
64 // Compute misfit G(x) - y
65 modelOutput -= this->m_observations;
66
67 // Solve \Sigma u = G(x) - y for u
68 this->m_covariance.invertMultiply(modelOutput, weightedMisfit);
69
70 // Deal with the multiplicative coefficients for each of the blocks
71 unsigned int numBlocks = this->m_covariance.numBlocks();
72 unsigned int offset = 0;
73
74 // For each block...
75 double cov_norm_factor = 0.0;
76 for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) {
77
78 // ...find the right hyperparameter
79 unsigned int index = domainVector.sizeLocal() + (i - numBlocks);
80 double coefficient = domainVector[index];
81
82 // ...divide the appropriate parts of the solution by the coefficient
83 unsigned int blockDim = this->m_covariance.getBlock(i).numRowsLocal();
84 for (unsigned int j = 0; j < blockDim; j++) {
85 // 'coefficient' is a variance, so we divide by it
86 modelOutput[offset+j] /= coefficient;
87 }
88
89 // Keep track of the part of the covariance matrix that appears in the
90 // normalising constant because of the hyperparameter
91 double cov_determinant = this->m_covariance.getBlock(i).determinant();
92 cov_determinant = std::sqrt(cov_determinant);
93
94 coefficient = std::sqrt(coefficient);
95 cov_norm_factor += std::log(std::pow(coefficient, blockDim) * cov_determinant);
96
97 offset += blockDim;
98 }
99
100 // Compute (G(x) - y)^T \Sigma^{-1} (G(x) - y)
101 modelOutput *= weightedMisfit;
102
103 double norm2_squared = modelOutput.sumOfComponents(); // This is square of 2-norm
104
105 return -0.5 * norm2_squared - cov_norm_factor;
106 }
107
108 } // End namespace QUESO
109
110 template class QUESO::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients<QUESO::GslVector, QUESO::GslMatrix>;
111