queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C
Go to the documentation of this file.
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>
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>
52 {
53 }
54 
55 template<class V, class M>
56 double
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 
A class representing a Gaussian likelihood with block-diagonal covariance matrix. ...
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
A templated class for handling sets.
Definition: VectorSet.h:52
unsigned int numBlocks() const
Return the number of blocks in the block diagonal matrix.
Class for representing block matrices using GSL library.
GslMatrix & getBlock(unsigned int i) const
Return block i in the block diagonal matrix.
virtual double lnValue(const V &domainVector) const
Logarithm of the value of the scalar function.
Base class for canned Gaussian likelihoods.
GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(const char *prefix, const VectorSet< V, M > &domainSet, const V &observations, const GslBlockMatrix &covariance)
Default constructor.

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5