25 #include <queso/VectorGaussianRandomField.h>
30 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
38 m_env (indexSet.env()),
39 m_prefix ((std::string)(prefix)+
"grf_"),
40 m_indexSet (indexSet),
41 m_imageSetPerIndex (imageSetPerIndex),
42 m_meanFunction (meanFunction),
43 m_covarianceFunction (covarianceFunction),
44 m_savedRvImageSpace (NULL),
45 m_savedRvLawExpVector(NULL),
46 m_savedRvLawCovMatrix(NULL),
52 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
57 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
64 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
68 return m_meanFunction;
71 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
75 return m_covarianceFunction;
78 template <
class P_V,
class P_M,
class Q_V,
class Q_M>
82 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
83 *m_env.subDisplayFile() <<
"Entering VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
89 "VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()",
90 "input data is not multiple of each other");
92 unsigned int numberOfImageValuesPerIndex = sampleValues.sizeLocal()/fieldPositions.size();
94 UQ_FATAL_TEST_MACRO(numberOfImageValuesPerIndex != m_imageSetPerIndex.vectorSpace().dimLocal(),
96 "VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()",
97 "invalid input data dimension");
99 if ((m_savedPositions.size() == 0 ) &&
100 (m_savedRvImageSpace == NULL) &&
101 (m_savedRvLawExpVector == NULL) &&
102 (m_savedRvLawCovMatrix == NULL) &&
103 (m_savedRv == NULL)) {
106 else if ((m_savedPositions.size() != 0 ) &&
107 (m_savedRvImageSpace != NULL) &&
108 (m_savedRvLawExpVector != NULL) &&
109 (m_savedRvLawCovMatrix != NULL) &&
110 (m_savedRv != NULL)) {
116 "VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()",
117 "invalid combination of pointer values");
120 unsigned int numberOfPositions = fieldPositions.size();
121 bool instantiate =
true;
122 if (m_savedPositions.size() == numberOfPositions) {
123 bool allPositionsAreEqual =
true;
124 for (
unsigned int i = 0; i < numberOfPositions; ++i) {
127 "VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()",
128 "m_savedPositions[i] should not be NULL");
129 if ((m_savedPositions[i]->sizeLocal() == fieldPositions[i]->sizeLocal()) &&
130 (*(m_savedPositions[i]) == *(fieldPositions[i]) )) {
134 allPositionsAreEqual =
false;
138 instantiate = !allPositionsAreEqual;
141 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
142 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
143 <<
": numberOfPositions = " << numberOfPositions
144 <<
", instantiate = " << instantiate
150 delete m_savedRvLawCovMatrix;
151 delete m_savedRvLawExpVector;
152 delete m_savedRvImageSpace;
153 for (
unsigned int i = 0; i < m_savedPositions.size(); ++i) {
154 delete m_savedPositions[i];
156 m_savedPositions.clear();
159 m_savedPositions.resize(numberOfPositions,NULL);
160 for (
unsigned int i = 0; i < m_savedPositions.size(); ++i) {
161 m_savedPositions[i] =
new P_V(*(fieldPositions[i]));
165 m_savedRvImageSpace =
new VectorSpace<Q_V,Q_M>(m_env,
"grf_", numberOfPositions*numberOfImageValuesPerIndex, NULL);
168 Q_V tmpVec(m_imageSetPerIndex.vectorSpace().zeroVector());
169 m_savedRvLawExpVector =
new Q_V(m_savedRvImageSpace->zeroVector());
170 for (
unsigned int i = 0; i < numberOfPositions; ++i) {
171 m_meanFunction.compute(*(fieldPositions[i]),NULL,tmpVec,NULL,NULL,NULL);
172 for (
unsigned int j = 0; j < numberOfImageValuesPerIndex; ++j) {
173 (*m_savedRvLawExpVector)[i*numberOfImageValuesPerIndex + j] = tmpVec[j];
178 Q_M tmpMat(m_imageSetPerIndex.vectorSpace().zeroVector());
179 m_savedRvLawCovMatrix =
new Q_M(m_savedRvImageSpace->zeroVector());
180 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
181 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
182 <<
": m_savedRvLawCovMatrix order = " << m_savedRvLawCovMatrix->numCols()
183 <<
", numberOfPositions = " << numberOfPositions
184 <<
", tmpMat order = " << tmpMat.numCols()
185 <<
", numberOfImageValuesPerIndex = " << numberOfImageValuesPerIndex
188 for (
unsigned int i = 0; i < numberOfPositions; ++i) {
189 for (
unsigned int j = 0; j < numberOfPositions; ++j) {
190 m_covarianceFunction.covMatrix(*(fieldPositions[i]),*(fieldPositions[j]),tmpMat);
193 if (testMat.chol() != 0) {
194 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
197 <<
", *(fieldPositions[i]) = " << *(fieldPositions[i])
198 <<
", *(fieldPositions[j]) = " << *(fieldPositions[j])
199 <<
", tmpMat = " << tmpMat
200 <<
", testMat = " << testMat
201 <<
", tmpMat is not positive definite"
205 "VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()",
206 "tmpMat is not positive definite");
209 for (
unsigned int k1 = 0; k1 < numberOfImageValuesPerIndex; ++k1) {
210 for (
unsigned int k2 = 0; k2 < numberOfImageValuesPerIndex; ++k2) {
211 unsigned int tmpI = i*numberOfImageValuesPerIndex + k1;
212 unsigned int tmpJ = j*numberOfImageValuesPerIndex + k2;
213 (*m_savedRvLawCovMatrix)(tmpI,tmpJ) = tmpMat(k1,k2);
214 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
215 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
220 <<
", tmpI = " << tmpI
221 <<
", tmpJ = " << tmpJ
222 <<
", *(fieldPositions[i]) = " << *(fieldPositions[i])
223 <<
", *(fieldPositions[j]) = " << *(fieldPositions[j])
224 <<
", (*m_savedRvLawCovMatrix)(tmpI,tmpJ) = " << (*m_savedRvLawCovMatrix)(tmpI,tmpJ)
234 *m_savedRvImageSpace,
235 *m_savedRvLawExpVector,
236 *m_savedRvLawCovMatrix);
238 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
239 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
240 <<
": just instantiated Gaussian RV"
241 <<
"\n *m_savedRvLawExpVector = " << *m_savedRvLawExpVector
242 <<
"\n *m_savedRvLawCovMatrix = " << *m_savedRvLawCovMatrix
244 for (
unsigned int i = 0; i < numberOfPositions; ++i) {
245 *m_env.subDisplayFile() <<
" *(m_savedPositions[" << i
246 <<
"]) = " << *(m_savedPositions[i])
253 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
254 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
255 <<
": about to realize sample values"
258 m_savedRv->realizer().realization(sampleValues);
259 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
260 *m_env.subDisplayFile() <<
"In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
261 <<
": just realized sample values"
265 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
266 *m_env.subDisplayFile() <<
"Leaving VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
VectorGaussianRandomField(const char *prefix, const VectorSet< P_V, P_M > &indexSet, const VectorSet< Q_V, Q_M > &imageSetPerIndex, const BaseVectorFunction< P_V, P_M, Q_V, Q_M > &meanFunction, const BaseMatrixCovarianceFunction< P_V, P_M, Q_V, Q_M > &covarianceFunction)
Constructor.
const BaseMatrixCovarianceFunction< P_V, P_M, Q_V, Q_M > & covarianceFunction() const
Covariance function; access to protected attribute m_covarianceFunction.
std::vector< P_V * > m_savedPositions
Saved positions.
const VectorSet< P_V, P_M > & indexSet() const
Index set; access to protected attribute m_indexSet.
A templated (base) class to accommodate covariance matrix of (random) vector functions.
~VectorGaussianRandomField()
Destructor.
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
const BaseVectorFunction< P_V, P_M, Q_V, Q_M > & meanFunction() const
Mean function; access to protected attribute m_meanFunction.
void sampleFunction(const std::vector< P_V * > &fieldPositions, Q_V &sampleValues)
Function that samples from a Gaussian PDF.
A templated (base) class for handling vector functions.