queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
VectorGaussianRandomField.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/VectorGaussianRandomField.h>
26 
27 namespace QUESO {
28 
29 // Default constructor -----------------------------
30 template <class P_V, class P_M, class Q_V, class Q_M>
32  const char* prefix,
33  const VectorSet<P_V,P_M>& indexSet,
34  const VectorSet<Q_V,Q_M>& imageSetPerIndex,
35  const BaseVectorFunction<P_V,P_M,Q_V,Q_M>& meanFunction,
36  const BaseMatrixCovarianceFunction<P_V,P_M,Q_V,Q_M>& covarianceFunction)
37  :
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),
47  m_savedRv (NULL)
48 {
49  m_savedPositions.clear();
50 }
51 // Destructor ---------------------------------------
52 template <class P_V, class P_M, class Q_V, class Q_M>
54 {
55 }
56 // Math methods -------------------------------------
57 template <class P_V, class P_M, class Q_V, class Q_M>
58 const VectorSet<P_V,P_M>&
60 {
61  return m_indexSet;
62 }
63 // --------------------------------------------------
64 template <class P_V, class P_M, class Q_V, class Q_M>
67 {
68  return m_meanFunction;
69 }
70 // --------------------------------------------------
71 template <class P_V, class P_M, class Q_V, class Q_M>
74 {
75  return m_covarianceFunction;
76 }
77 // --------------------------------------------------
78 template <class P_V, class P_M, class Q_V, class Q_M>
79 void
80 VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction(const std::vector<P_V*>& fieldPositions, Q_V& sampleValues)
81 {
82  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
83  *m_env.subDisplayFile() << "Entering VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
84  << std::endl;
85  }
86 
87  queso_require_equal_to_msg(( sampleValues.sizeLocal() % fieldPositions.size() ), 0, "input data is not multiple of each other");
88 
89  unsigned int numberOfImageValuesPerIndex = sampleValues.sizeLocal()/fieldPositions.size();
90 
91  queso_require_equal_to_msg(numberOfImageValuesPerIndex, m_imageSetPerIndex.vectorSpace().dimLocal(), "invalid input data dimension");
92 
93  if ((m_savedPositions.size() == 0 ) &&
94  (m_savedRvImageSpace == NULL) &&
95  (m_savedRvLawExpVector == NULL) &&
96  (m_savedRvLawCovMatrix == NULL) &&
97  (m_savedRv == NULL)) {
98  // Ok
99  }
100  else if ((m_savedPositions.size() != 0 ) &&
101  (m_savedRvImageSpace != NULL) &&
102  (m_savedRvLawExpVector != NULL) &&
103  (m_savedRvLawCovMatrix != NULL) &&
104  (m_savedRv != NULL)) {
105  // Ok
106  }
107  else {
108  queso_error_msg("invalid combination of pointer values");
109  }
110 
111  unsigned int numberOfPositions = fieldPositions.size();
112  bool instantiate = true;
113  if (m_savedPositions.size() == numberOfPositions) {
114  bool allPositionsAreEqual = true;
115  for (unsigned int i = 0; i < numberOfPositions; ++i) {
116  queso_require_msg(m_savedPositions[i], "m_savedPositions[i] should not be NULL");
117  if ((m_savedPositions[i]->sizeLocal() == fieldPositions[i]->sizeLocal()) &&
118  (*(m_savedPositions[i]) == *(fieldPositions[i]) )) {
119  // Ok
120  }
121  else {
122  allPositionsAreEqual = false;
123  break;
124  }
125  } // for i
126  instantiate = !allPositionsAreEqual;
127  }
128 
129  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
130  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
131  << ": numberOfPositions = " << numberOfPositions
132  << ", instantiate = " << instantiate
133  << std::endl;
134  }
135 
136  if (instantiate) {
137  delete m_savedRv;
138  delete m_savedRvLawCovMatrix;
139  delete m_savedRvLawExpVector;
140  delete m_savedRvImageSpace;
141  for (unsigned int i = 0; i < m_savedPositions.size(); ++i) {
142  delete m_savedPositions[i];
143  }
144  m_savedPositions.clear();
145 
146  // Set m_savedPositions
147  m_savedPositions.resize(numberOfPositions,NULL);
148  for (unsigned int i = 0; i < m_savedPositions.size(); ++i) {
149  m_savedPositions[i] = new P_V(*(fieldPositions[i]));
150  }
151 
152  // Set m_savedRvImageSpace
153  m_savedRvImageSpace = new VectorSpace<Q_V,Q_M>(m_env, "grf_", numberOfPositions*numberOfImageValuesPerIndex, NULL);
154 
155  // Set m_savedRvLawExpVector
156  Q_V tmpVec(m_imageSetPerIndex.vectorSpace().zeroVector());
157  m_savedRvLawExpVector = new Q_V(m_savedRvImageSpace->zeroVector());
158  for (unsigned int i = 0; i < numberOfPositions; ++i) {
159  m_meanFunction.compute(*(fieldPositions[i]),NULL,tmpVec,NULL,NULL,NULL);
160  for (unsigned int j = 0; j < numberOfImageValuesPerIndex; ++j) {
161  (*m_savedRvLawExpVector)[i*numberOfImageValuesPerIndex + j] = tmpVec[j];
162  }
163  }
164 
165  // Set m_savedRvLawCovMatrix
166  Q_M tmpMat(m_imageSetPerIndex.vectorSpace().zeroVector());
167  m_savedRvLawCovMatrix = new Q_M(m_savedRvImageSpace->zeroVector());
168  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
169  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
170  << ": m_savedRvLawCovMatrix order = " << m_savedRvLawCovMatrix->numCols()
171  << ", numberOfPositions = " << numberOfPositions
172  << ", tmpMat order = " << tmpMat.numCols()
173  << ", numberOfImageValuesPerIndex = " << numberOfImageValuesPerIndex
174  << std::endl;
175  }
176  for (unsigned int i = 0; i < numberOfPositions; ++i) {
177  for (unsigned int j = 0; j < numberOfPositions; ++j) {
178  m_covarianceFunction.covMatrix(*(fieldPositions[i]),*(fieldPositions[j]),tmpMat);
179 #if 1
180  Q_M testMat(tmpMat);
181  if (testMat.chol() != 0) {
182  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
183  << ": i = " << i
184  << ", j = " << j
185  << ", *(fieldPositions[i]) = " << *(fieldPositions[i])
186  << ", *(fieldPositions[j]) = " << *(fieldPositions[j])
187  << ", tmpMat = " << tmpMat
188  << ", testMat = " << testMat
189  << ", tmpMat is not positive definite"
190  << std::endl;
191  queso_error_msg("tmpMat is not positive definite");
192  }
193 #endif
194  for (unsigned int k1 = 0; k1 < numberOfImageValuesPerIndex; ++k1) {
195  for (unsigned int k2 = 0; k2 < numberOfImageValuesPerIndex; ++k2) {
196  unsigned int tmpI = i*numberOfImageValuesPerIndex + k1;
197  unsigned int tmpJ = j*numberOfImageValuesPerIndex + k2;
198  (*m_savedRvLawCovMatrix)(tmpI,tmpJ) = tmpMat(k1,k2);
199  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
200  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
201  << ": i = " << i
202  << ", j = " << j
203  << ", k1 = " << k1
204  << ", k2 = " << k2
205  << ", tmpI = " << tmpI
206  << ", tmpJ = " << tmpJ
207  << ", *(fieldPositions[i]) = " << *(fieldPositions[i])
208  << ", *(fieldPositions[j]) = " << *(fieldPositions[j])
209  << ", (*m_savedRvLawCovMatrix)(tmpI,tmpJ) = " << (*m_savedRvLawCovMatrix)(tmpI,tmpJ)
210  << std::endl;
211  }
212  }
213  }
214  }
215  }
216 
217  // Set m_savedRv
218  m_savedRv = new GaussianVectorRV<Q_V,Q_M>("grf_",
219  *m_savedRvImageSpace,
220  *m_savedRvLawExpVector,
221  *m_savedRvLawCovMatrix);
222 
223  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
224  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
225  << ": just instantiated Gaussian RV"
226  << "\n *m_savedRvLawExpVector = " << *m_savedRvLawExpVector
227  << "\n *m_savedRvLawCovMatrix = " << *m_savedRvLawCovMatrix
228  << std::endl;
229  for (unsigned int i = 0; i < numberOfPositions; ++i) {
230  *m_env.subDisplayFile() << " *(m_savedPositions[" << i
231  << "]) = " << *(m_savedPositions[i])
232  << std::endl;
233  }
234  }
235  } // if (instantiate)
236 
237  // Generate sample function
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  << ": about to realize sample values"
241  << std::endl;
242  }
243  m_savedRv->realizer().realization(sampleValues);
244  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
245  *m_env.subDisplayFile() << "In VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
246  << ": just realized sample values"
247  << std::endl;
248  }
249 
250  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
251  *m_env.subDisplayFile() << "Leaving VectorGaussianRandomField<P_V,P_M,Q_V,Q_M>::sampleFunction()"
252  << std::endl;
253  }
254 
255  return;
256 }
257 
258 } // End namespace QUESO
A templated (base) class for handling vector functions.
void sampleFunction(const std::vector< P_V * > &fieldPositions, Q_V &sampleValues)
Function that samples from a Gaussian PDF.
const VectorSet< P_V, P_M > & indexSet() const
Index set; access to protected attribute m_indexSet.
const BaseVectorFunction< P_V, P_M, Q_V, Q_M > & meanFunction() const
Mean function; access to protected attribute m_meanFunction.
const BaseMatrixCovarianceFunction< P_V, P_M, Q_V, Q_M > & covarianceFunction() const
Covariance function; access to protected attribute m_covarianceFunction.
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.
std::vector< P_V * > m_savedPositions
Saved positions.
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
A templated (base) class to accommodate covariance matrix of (random) vector functions.

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