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/ExponentialMatrixCovarianceFunction.h>
26
27 namespace QUESO {
28
29 // Default constructor -----------------------------
30 template<class P_V, class P_M, class Q_V, class Q_M>
ExponentialMatrixCovarianceFunction(const char * prefix,const VectorSet<P_V,P_M> & basicDomainSet,const VectorSet<Q_V,Q_M> & imageSet,const Q_M & sigmas,const Q_M & as)31 ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::ExponentialMatrixCovarianceFunction(
32 const char* prefix,
33 const VectorSet<P_V,P_M>& basicDomainSet,
34 const VectorSet<Q_V,Q_M>& imageSet,
35 const Q_M& sigmas,
36 const Q_M& as)
37 :
38 BaseMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>(prefix,basicDomainSet,imageSet),
39 m_sigmas(NULL),
40 m_as (NULL)
41 {
42 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
43 *m_env.subDisplayFile() << "Entering ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::constructor()"
44 << ": prefix = " << m_prefix
45 << std::endl;
46 }
47
48 m_sigmas = new Q_M(sigmas);
49 m_as = new Q_M(as);
50
51 unsigned int matrixOrder = m_imageSet.vectorSpace().dimLocal();
52
53 queso_require_equal_to_msg(m_sigmas->numRowsLocal(), matrixOrder, "m_sigmas has invalid number of rows");
54
55 queso_require_equal_to_msg(m_sigmas->numCols(), matrixOrder, "m_sigmas has invalid number of columns");
56
57 queso_require_equal_to_msg(m_as->numRowsLocal(), matrixOrder, "m_as has invalid number of rows");
58
59 queso_require_equal_to_msg(m_as->numCols(), matrixOrder, "m_as has invalid number of columns");
60
61 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
62 *m_env.subDisplayFile() << "Leaving ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::constructor()"
63 << ": prefix = " << m_prefix
64 << std::endl;
65 }
66 }
67 // Destructor ---------------------------------------
68 template<class P_V, class P_M, class Q_V, class Q_M>
~ExponentialMatrixCovarianceFunction()69 ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::~ExponentialMatrixCovarianceFunction()
70 {
71 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
72 *m_env.subDisplayFile() << "Entering ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::destructor()"
73 << ": prefix = " << m_prefix
74 << std::endl;
75 }
76
77 delete m_as;
78 delete m_sigmas;
79
80 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
81 *m_env.subDisplayFile() << "Leaving ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::destructor()"
82 << ": prefix = " << m_prefix
83 << std::endl;
84 }
85 }
86 // Math methods -------------------------------------
87 template<class P_V, class P_M, class Q_V, class Q_M>
88 void
covMatrix(const P_V & domainVector1,const P_V & domainVector2,Q_M & imageMatrix)89 ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::covMatrix(const P_V& domainVector1, const P_V& domainVector2, Q_M& imageMatrix) const
90 {
91 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
92 *m_env.subDisplayFile() << "Entering ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::covMatrix()"
93 << std::endl;
94 }
95
96 unsigned int matrixOrder = m_imageSet.vectorSpace().dimLocal();
97
98 queso_require_equal_to_msg(imageMatrix.numRowsLocal(), matrixOrder, "imageMatrix has invalid number of rows");
99
100 queso_require_equal_to_msg(imageMatrix.numCols(), matrixOrder, "imageMatrix has invalid number of columns");
101
102 double tmpSq = -(domainVector1 - domainVector2).norm2Sq();
103
104 for (unsigned int i = 0; i < matrixOrder; ++i) {
105 for (unsigned int j = 0; j < matrixOrder; ++j) {
106 double tmp = tmpSq/( (*m_sigmas)(i,j) * (*m_sigmas)(i,j) );
107 imageMatrix(i,j) = (*m_as)(i,j) * std::exp(tmp);
108 }
109 }
110
111 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
112 *m_env.subDisplayFile() << "Leaving ExponentialMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>::covMatrix()"
113 << std::endl;
114 }
115
116 return;
117 }
118
119 } // End namespace QUESO
120