Logarithm of the value of the scalar function.
151 const unsigned int totalOutputs = totalRuns * numOutputs;
153 const unsigned int residualSize = (numOutputs == 1) ? totalOutputs :
156 double prodScenario = 1.0;
157 double prodParameter = 1.0;
158 double prodDiscrepancy = 1.0;
163 unsigned int dimSum = 3 +
164 (numOutputs > 1) * 2 +
172 const unsigned int offset1 = (numOutputs == 1) ?
176 const unsigned int offset1b = offset1 +
180 const unsigned int offset2 = (numOutputs == 1) ?
184 const MpiComm & comm = domainVector.map().Comm();
185 Map z_map(residualSize, 0, comm);
186 M covMatrix(this->
m_env, z_map, residualSize);
189 for (
unsigned int k = 0;
k < dimParameter;
k++) {
191 domainVectorParameter[
k] = domainVector[
k];
195 for (
unsigned int i = 0; i < totalRuns; i++) {
207 parameter1 = &domainVectorParameter;
216 for (
unsigned int j = 0; j < totalRuns; j++) {
221 if (j < this->m_numExperiments) {
223 parameter2 = &domainVectorParameter;
234 unsigned int emulatorCorrStrStart =
236 for (
unsigned int k = 0;
k < dimScenario;
k++) {
237 const double & emulator_corr_strength =
238 domainVector[emulatorCorrStrStart+
k];
239 prodScenario *= std::pow(emulator_corr_strength,
240 4.0 * ((*scenario1)[
k] - (*scenario2)[
k]) *
241 ((*scenario1)[k] - (*scenario2)[k]));
248 for (
unsigned int k = 0; k < dimParameter; k++) {
252 const double & emulator_corr_strength =
253 domainVector[emulatorCorrStrStart+dimScenario+
k];
254 prodParameter *= std::pow(
255 emulator_corr_strength,
256 4.0 * ((*parameter1)[k] - (*parameter2)[k]) *
257 ((*parameter1)[k] - (*parameter2)[k]));
264 for (
unsigned int basis = 0; basis !=
num_svd_terms; ++basis)
269 const double relevant_precision =
270 domainVector[dimParameter+basis+1+(numOutputs>1)];
273 const unsigned int stridei =
275 const unsigned int offseti =
277 const unsigned int stridej =
279 const unsigned int offsetj =
282 covMatrix(offseti+basis*stridei+i,
283 offsetj+basis*stridej+j) =
284 prodScenario * prodParameter / relevant_precision;
289 if (i < this->m_numExperiments && j < this->m_numExperiments) {
292 prodDiscrepancy = 1.0;
293 unsigned int discrepancyCorrStrStart = dimParameter +
298 for (
unsigned int k = 0; k < dimScenario; k++) {
299 const double & discrepancy_corr_strength =
300 domainVector[discrepancyCorrStrStart+
k];
302 std::pow(discrepancy_corr_strength, 4.0 *
303 ((*cross_scenario1)[k] - (*cross_scenario2)[k]) *
304 ((*cross_scenario1)[k] - (*cross_scenario2)[k]));
307 const double discrepancy_precision =
308 domainVector[discrepancyCorrStrStart-1];
314 const double R_v = prodDiscrepancy / discrepancy_precision;
315 for (
unsigned int disc = 0; disc != num_discrepancy_bases;
317 covMatrix(disc*m_numExperiments+i,
318 disc*m_numExperiments+j) += R_v;
325 const double experimentalError =
330 covMatrix(i,j) += experimentalError;
337 const double emulator_data_precision = domainVector[dimSum-1-(numOutputs>1)];
339 double nugget = 1.0 / emulator_data_precision;
341 for (
unsigned int disc = 0; disc != num_discrepancy_bases;
343 covMatrix(disc*m_numExperiments+i,
344 disc*m_numExperiments+i) += nugget;
351 const double lambda_y = domainVector[dimSum-1];
352 const double inv_lambda_y = 1.0/lambda_y;
355 for (
unsigned int i=0; i != BT_Wy_B_size; ++i)
356 for (
unsigned int j=0; j != BT_Wy_B_size; ++j)
359 const double emulator_precision =
360 domainVector[dimParameter+1];
361 const double inv_emulator_precision = 1.0/emulator_precision;
363 unsigned int KT_K_size =
KT_K_inv.numCols();
364 for (
unsigned int i=0; i != KT_K_size; ++i)
365 for (
unsigned int j=0; j != KT_K_size; ++j)
366 covMatrix(i+offset2,j+offset2) +=
367 KT_K_inv(i,j) * inv_emulator_precision;
373 V sol(covMatrix.invertMultiply(
residual));
376 double minus_2_log_lhd = 0.0;
378 for (
unsigned int i = 0; i < residualSize; i++) {
379 minus_2_log_lhd += sol[i] *
residual[i];
383 for (
unsigned int i = 0; i < residualSize; i++) {
385 std::cout <<
"NaN sol[" << i <<
']' << std::endl;
387 std::cout <<
"NaN residual[" << i <<
']' << std::endl;
389 std::cout <<
"Covariance Matrix:" << std::endl;
390 covMatrix.print(std::cout);
397 double cov_det = covMatrix.determinant();
401 std::cout <<
"Non-positive determinant for covMatrix = " << std::endl;
402 covMatrix.print(std::cout);
406 minus_2_log_lhd += std::log(covMatrix.determinant());
409 return -0.5 * minus_2_log_lhd;
const unsigned int m_numExperiments
#define queso_assert_greater(expr1, expr2)
#define queso_assert_greater_equal(expr1, expr2)
const std::vector< V * > & m_experimentScenarios
const M & m_experimentErrors
const VectorSpace< V, M > & m_scenarioSpace
std::vector< V > m_discrepancyBases
const std::vector< V * > & m_simulationParameters
const VectorSpace< V, M > & m_experimentOutputSpace
#define queso_assert(asserted)
unsigned int num_svd_terms
const unsigned int m_numSimulations
const BaseEnvironment & m_env
const std::vector< V * > & m_simulationScenarios
const VectorSpace< V, M > & m_parameterSpace