queso-0.53.0
GcmSimulationTildeInfo.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-2015 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/GcmSimulationTildeInfo.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 template <class S_V,class S_M,class P_V,class P_M,class Q_V,class Q_M>
32 GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::GcmSimulationTildeInfo(
33  const GpmsaComputerModelOptions& gcmOptionsObj,
34  const GcmSimulationInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>& s)
35  :
36  m_env (s.m_env),
37  m_Kmat_tilde (m_env,s.m_eta_space.map(),s.m_Kmat_rank),
38  m_w_tilde_space (m_env, "w_tilde_", s.m_Kmat_rank, NULL), // rr0 check
39  m_Lkmat (m_env,m_w_tilde_space.map(),s.m_Kmat.numCols()), // rr0 check
40  m_Ktildet_Ktilde (m_w_tilde_space.zeroVector()),
41  m_Ktildet_Ktilde_inv (m_w_tilde_space.zeroVector()),
42  m_Zvec_tilde_hat_w (m_w_tilde_space.zeroVector()),
43  m_a_eta_modifier_tilde(0),
44  m_b_eta_modifier_tilde(0)
45 {
46  std::set<unsigned int> tmpSet;
47  tmpSet.insert(m_env.subId());
48 
49  //******************************************************************************
50  // Tilde situation: form 'm_Kmat_tilde'
51  // Tilde situation: form 'm_w_tilde_space'
52  // Tilde situation: form 'm_Lkmat'
53  //******************************************************************************
54  if (s.m_Kmat.numRowsGlobal() >= s.m_Kmat.numCols()) {
55  Q_M matkU(s.m_Kmat.svdMatU());
56  unsigned int kuMatRank = matkU.rank(0.,1.e-8 ); // todo: should be an option
57  unsigned int kuMatRank14 = matkU.rank(0.,1.e-14);
58  if (m_env.subDisplayFile()) {
59  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
60  << ": matkU.numRowsLocal() = " << matkU.numRowsLocal()
61  << ", matkU.numCols() = " << matkU.numCols()
62  << ", matkU.rank(0.,1.e-8) = " << kuMatRank
63  << ", matkU.rank(0.,1.e-14) = " << kuMatRank14
64  << std::endl;
65  }
66 
67  if (m_env.checkingLevel() >= 1) {
68  Q_M matkUcheck(s.m_w_space.zeroVector());
69  Q_V vecI(s.m_eta_space.zeroVector());
70  Q_V vecJ(s.m_eta_space.zeroVector());
71  for (unsigned int i = 0; i < matkU.numCols(); ++i) {
72  matkU.getColumn(i,vecI);
73  for (unsigned int j = i; j < matkU.numCols(); ++j) {
74  matkU.getColumn(j,vecJ);
75  matkUcheck(i,j) = scalarProduct(vecI,vecJ);
76  }
77  }
78  matkUcheck.setPrintHorizontally(false);
79  if (m_env.subDisplayFile()) {
80  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
81  << ": matkUcheck.numRowsLocal() = " << matkUcheck.numRowsLocal()
82  << ", matkUcheck.numCols() = " << matkUcheck.numCols()
83  << ", m_Kmat_rank = " << s.m_Kmat_rank
84  << ", matkUcheck =\n" << matkUcheck
85  << std::endl;
86  }
87  }
88 
89  Q_V veckJ(s.m_eta_space.zeroVector());
90  for (unsigned int j = 0; j < s.m_Kmat_rank; ++j) {
91  matkU.getColumn(j,veckJ);
92  m_Kmat_tilde.setColumn(j,veckJ);
93  }
94  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
95  m_Kmat_tilde.subWriteContents("Ktilde",
96  "Ktilde1",
97  "m",
98  tmpSet);
99  }
100  if (m_env.subDisplayFile()) {
101  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
102  << ": m_Kmat_tilde computed (1)"
103  << std::endl;
104  }
105  }
106  else {
107  queso_error_msg("code for horizontal 'm_Kmat' is not ready yet");
108  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
109  m_Kmat_tilde.subWriteContents("Ktilde",
110  "Ktilde2",
111  "m",
112  tmpSet);
113  }
114  if (m_env.subDisplayFile()) {
115  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
116  << ": m_Kmat_tilde computed (2)"
117  << std::endl;
118  }
119  }
120 
121  if (m_env.subDisplayFile()) {
122  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
123  << ": finished forming 'm_Kmat_tilde'"
124  << ", m_Kmat_tilde.numRowsLocal() = " << m_Kmat_tilde.numRowsLocal()
125  << ", m_Kmat_tilde.numCols() = " << m_Kmat_tilde.numCols()
126  << std::endl;
127  }
128 
129  m_Kmat_tilde.svdSolve(s.m_Kmat,m_Lkmat);
130  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
131  m_Lkmat.subWriteContents("Lkmat",
132  "Lkmat",
133  "m",
134  tmpSet);
135  }
136  if (m_env.subDisplayFile()) {
137  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
138  << ": m_Lkmat_tilde computed"
139  << std::endl;
140  }
141 
142  //********************************************************************************
143  // Form 'Ktilde^T' matrix
144  //********************************************************************************
145  Q_M Ktildet(m_env,m_w_tilde_space.map(),s.m_eta_size);
146  Ktildet.fillWithTranspose(0,0,s.m_Kmat,true,true);
147 
148  //********************************************************************************
149  // Compute 'Ktilde' rank
150  //********************************************************************************
151  unsigned int kTildeRank = m_Kmat_tilde.rank(0.,1.e-8 ); // todo: should be an option
152  unsigned int kRank14 = m_Kmat_tilde.rank(0.,1.e-14);
153  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
154  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
155  << ": m_Kmat_tilde.numRowsLocal() = " << m_Kmat_tilde.numRowsLocal()
156  << ", m_Kmat_tilde.numCols() = " << m_Kmat_tilde.numCols()
157  << ", m_Kmat_tilde.rank(0.,1.e-8) = " << kTildeRank
158  << ", m_Kmat_tilde.rank(0.,1.e-14) = " << kRank14
159  << std::endl;
160  }
161 
162  queso_require_equal_to_msg(kTildeRank, std::min(m_Kmat_tilde.numRowsGlobal(),m_Kmat_tilde.numCols()), "'m_Kmat_tilde' does not have a proper rank");
163 
164  //******************************************************************************
165  // Tilde situation: compute 'Ktilde^T Ktilde' matrix, and its inverse
166  //******************************************************************************
167  m_Ktildet_Ktilde = Ktildet * m_Kmat_tilde;
168  if (m_env.subDisplayFile()) {
169  m_Ktildet_Ktilde.setPrintHorizontally(false);
170  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
171  << ": finished computing 'm_Ktildet_Ktilde'"
172  //<< "\n m_Ktildet_Ktilde =\n" << m_Ktildet_Ktilde
173  << std::endl;
174  }
175  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
176  m_Ktildet_Ktilde.subWriteContents("Ktildet_Ktilde",
177  "Ktildet_Ktilde",
178  "m",
179  tmpSet);
180  }
181  if (m_env.subDisplayFile()) {
182  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
183  << ": m_Ktildet_Ktilde computed"
184  << std::endl;
185  }
186 
187  double ktildetKtildeLnDeterminant = m_Ktildet_Ktilde.lnDeterminant();
188  unsigned int ktildetKtildeRank = m_Ktildet_Ktilde.rank(0.,1.e-8 ); // todo: should be an option
189  unsigned int ktildetKtildeRank14 = m_Ktildet_Ktilde.rank(0.,1.e-14);
190  if (m_env.subDisplayFile()) {
191  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
192  << ": m_Ktildet_Ktilde.numRowsLocal() = " << m_Ktildet_Ktilde.numRowsLocal()
193  << ", m_Ktildet_Ktilde.numCols() = " << m_Ktildet_Ktilde.numCols()
194  << ", m_Ktildet_Ktilde.lnDeterminant() = " << ktildetKtildeLnDeterminant
195  << ", m_Ktildet_Ktilde.rank(0.,1.e-8) = " << ktildetKtildeRank
196  << ", m_Ktildet_Ktilde.rank(0.,1.e-14) = " << ktildetKtildeRank14
197  << std::endl;
198  }
199 
200  m_Ktildet_Ktilde_inv = m_Ktildet_Ktilde.inverse();
201  if (m_env.subDisplayFile()) {
202  m_Ktildet_Ktilde_inv.setPrintHorizontally(false);
203  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
204  << ": finished computing 'm_Ktildet_Ktilde_inv'"
205  //<< "\n m_Ktildet_Ktilde_inv =\n" << m_Ktildet_Ktilde_inv
206  << std::endl;
207  }
208  if (gcmOptionsObj.m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != gcmOptionsObj.m_ov.m_dataOutputAllowedSet.end()) {
209  m_Ktildet_Ktilde_inv.subWriteContents("Ktildet_Ktilde_inv",
210  "Ktildet_Ktilde_inv",
211  "m",
212  tmpSet);
213  }
214  if (m_env.subDisplayFile()) {
215  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
216  << ": m_Ktildet_Ktilde_inv computed"
217  << std::endl;
218  }
219 
220  double ktildetKtildeInvLnDeterminant = m_Ktildet_Ktilde_inv.lnDeterminant();
221  unsigned int ktildetKtildeInvRank = m_Ktildet_Ktilde_inv.rank(0.,1.e-8 ); // todo: should be an option
222  unsigned int ktildetKtildeInvRank14 = m_Ktildet_Ktilde_inv.rank(0.,1.e-14);
223  if (m_env.subDisplayFile()) {
224  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
225  << ": m_Ktildet_Ktilde_inv.numRowsLocal() = " << m_Ktildet_Ktilde_inv.numRowsLocal()
226  << ", m_Ktildet_Ktilde_inv.numCols() = " << m_Ktildet_Ktilde_inv.numCols()
227  << ", m_Ktildet_Ktilde_inv.lnDeterminant() = " << ktildetKtildeInvLnDeterminant
228  << ", m_Ktildet_Ktilde_inv.rank(0.,1.e-8) = " << ktildetKtildeInvRank
229  << ", m_Ktildet_Ktilde_inv.rank(0.,1.e-14) = " << ktildetKtildeInvRank14
230  << std::endl;
231  }
232 
233  //********************************************************************************
234  // Compute 'tilde' exponent modifiers
235  //********************************************************************************
236  m_a_eta_modifier_tilde = ((double) s.m_paper_m) * ((double) (s.m_paper_n_eta - s.m_paper_p_eta)) / 2.;
237  if (m_env.subDisplayFile()) {
238  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
239  << ": m_a_eta_modifier_tilde = " << m_a_eta_modifier_tilde
240  << std::endl;
241  }
242 
243  Q_V etaVec_transformed(s.m_simulationModel.etaVec_transformed("Gp.h.003"));
244  if (m_env.subDisplayFile()) {
245  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
246  << ": m_Zvec_tilde_hat_w.sizeLocal() = " << m_Zvec_tilde_hat_w.sizeLocal()
247  << ", m_Ktildet_Ktilde_inv.numRowsLocal() = " << m_Ktildet_Ktilde_inv.numRowsLocal()
248  << ", m_Ktildet_Ktilde_inv.numCols() = " << m_Ktildet_Ktilde_inv.numCols()
249  << ", Ktildet.numRowsLocal() = " << Ktildet.numRowsLocal()
250  << ", Ktildet.numCols() = " << Ktildet.numCols()
251  << ", etaVec_transformed.sizeLocal() = " << etaVec_transformed.sizeLocal()
252  << std::endl;
253  }
254  m_Zvec_tilde_hat_w = m_Ktildet_Ktilde_inv * (Ktildet * etaVec_transformed);
255  Q_V tmpVec1(etaVec_transformed - (m_Kmat_tilde * m_Zvec_tilde_hat_w));
256  m_b_eta_modifier_tilde = scalarProduct(etaVec_transformed,tmpVec1) / 2.;
257 
258  if (m_env.subDisplayFile()) {
259  *m_env.subDisplayFile() << "In GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::constructor()"
260  << ": m_b_eta_modifier_tilde = " << m_b_eta_modifier_tilde
261  << std::endl;
262  }
263 }
264 
265 template <class S_V,class S_M,class P_V,class P_M,class Q_V,class Q_M>
266 GcmSimulationTildeInfo<S_V,S_M,P_V,P_M,Q_V,Q_M>::~GcmSimulationTildeInfo()
267 {
268 }
269 
270 } // End namespace QUESO
271 
272 template class QUESO::GcmSimulationTildeInfo<QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix, QUESO::GslVector, QUESO::GslMatrix>;
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:396
QUESO::GcmSimulationTildeInfo::GcmSimulationTildeInfo
GcmSimulationTildeInfo(const GpmsaComputerModelOptions &gcmOptionsObj, const GcmSimulationInfo< S_V, S_M, P_V, P_M, Q_V, Q_M > &s)
Definition: GcmSimulationTildeInfo.C:32
QUESO::VectorSpace::map
const Map & map() const
Map.
Definition: VectorSpace.C:157
queso_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
QUESO::GcmOptionsValues::m_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
Definition: GpmsaComputerModelOptions.h:79
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
QUESO::scalarProduct
double scalarProduct(const GslVector &x, const GslVector &y)
Definition: GslVector.C:1150
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:274
QUESO::GcmSimulationInfo::m_w_space
VectorSpace< Q_V, Q_M > m_w_space
Definition: GcmSimulationInfo.h:111
QUESO::GcmSimulationInfo::m_eta_space
VectorSpace< Q_V, Q_M > m_eta_space
Definition: GcmSimulationInfo.h:108
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:295
QUESO::GcmSimulationInfo::m_simulationModel
const SimulationModel< S_V, S_M, P_V, P_M, Q_V, Q_M > & m_simulationModel
Definition: GcmSimulationInfo.h:52
QUESO::BaseEnvironment::checkingLevel
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
Definition: Environment.C:410
QUESO::GcmSimulationTildeInfo::m_w_tilde_space
VectorSpace< Q_V, Q_M > m_w_tilde_space
Definition: GcmSimulationTildeInfo.h:49
Generated on Thu Jun 11 2015 13:52:32 for queso-0.53.0 by  doxygen 1.8.5