C++ UQ_FATAL_TEST_MACRO函数代码示例

GslVector::GslVector(const BaseEnvironment& env, const Map& map, double value)
  :
  Vector(env,map),
  m_vec        (gsl_vector_calloc(map.NumGlobalElements()))
{
  //std::cout << "Entering GslVector::constructor(2)" << std::endl;

  UQ_FATAL_TEST_MACRO((m_vec == NULL),
                      m_env.worldRank(),
                      "GslVector::constructor(2)",
                      "null vector generated");

  UQ_FATAL_TEST_MACRO(m_vec->size != (unsigned int) map.NumMyElements(),
                      m_env.worldRank(),
                      "GslVector::constructor(2)",
                      "incompatible local vec size");

  UQ_FATAL_TEST_MACRO(m_vec->size != (unsigned int) map.NumGlobalElements(),
                      m_env.worldRank(),
                      "GslVector::constructor(2)",
                      "incompatible global vec size");

  this->cwSet(value);

  UQ_FATAL_TEST_MACRO(m_vec->size != (unsigned int) m_map.NumMyElements(),
                      m_env.worldRank(),
                      "GslVector::constructor(2)",
                      "incompatible own vec size");

  //std::cout << "Leaving GslVector::constructor(2)" << std::endl;
}

// QUESO MpiComm MPI Constructor ------------------
MpiComm::MpiComm(const BaseEnvironment& env, RawType_MPI_Comm inputRawComm)
  :
  m_env          (env),
#ifdef QUESO_HAS_TRILINOS
  m_epetraMpiComm( new Epetra_MpiComm(inputRawComm) ),
#endif
  m_rawComm      (inputRawComm),
  m_worldRank    (-1),
  m_myPid        (-1),
  m_numProc      (-1)
{
  int mpiRC = MPI_Comm_rank(inputRawComm,&m_worldRank);
  UQ_FATAL_TEST_MACRO(mpiRC != MPI_SUCCESS,
                      UQ_UNAVAILABLE_RANK,
                      "MpiComm::constructor()",
                      "failed MPI_Comm_rank() on full rank");

  mpiRC = MPI_Comm_rank(inputRawComm,&m_myPid);
  UQ_FATAL_TEST_MACRO(mpiRC != MPI_SUCCESS,
                      m_worldRank,
                      "MpiComm::constructor()",
                      "failed MPI_Comm_rank() on inputRawComm");

  mpiRC = MPI_Comm_size(inputRawComm,&m_numProc);
  UQ_FATAL_TEST_MACRO(mpiRC != MPI_SUCCESS,
                      m_worldRank,
                      "MpiComm::constructor()",
                      "failed MPI_Comm_size() on inputRawComm");
}

double                     
uqLagrangeBasis1D1DFunctionClass::deriv(double domainValue) const
{
  double value = 0.;

  if ((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)) {
    std::cerr << "In uqLagrangeBasis1D1DFunctionClass::deriv()"
              << ": requested x ("            << domainValue
              << ") is out of the interval (" << m_minDomainValue
              << ", "                         << m_maxDomainValue
              << ")"
              << std::endl;
  }

  UQ_FATAL_TEST_MACRO(((domainValue < m_minDomainValue) || (domainValue > m_maxDomainValue)),
                      UQ_UNAVAILABLE_RANK,
                      "uqLagrangeBasis1D1DFunctionClass::deriv()",
                      "x out of range");

  UQ_FATAL_TEST_MACRO(true,
                      UQ_UNAVAILABLE_RANK,
                      "uqLagrangeBasis1D1DFunctionClass::deriv()",
                      "not implemented yet");

  return value;
}

void
GslVector::matlabDiff(
  unsigned int      firstPositionToStoreDiff,
  double            valueForRemainderPosition,
  GslVector& outputVec) const
{
  unsigned int size = this->sizeLocal();

  UQ_FATAL_TEST_MACRO(firstPositionToStoreDiff > 1,
                      m_env.worldRank(),
                      "GslVector::matlabDiff()",
                      "invalid firstPositionToStoreDiff");

  UQ_FATAL_TEST_MACRO(size != outputVec.sizeLocal(),
                      m_env.worldRank(),
                      "GslVector::matlabDiff()",
                      "invalid size of outputVecs");

  for (unsigned int i = 0; i < (size-1); ++i) {
    outputVec[firstPositionToStoreDiff+i] = (*this)[i+1]-(*this)[i];
  }
  if (firstPositionToStoreDiff == 0) {
    outputVec[size-1] = valueForRemainderPosition;
  }
  else {
    outputVec[0] = valueForRemainderPosition;
  }

  return;
}

void
qoiRoutine(
  const QUESO::GslVector&                    paramValues,
  const QUESO::GslVector*                    paramDirection,
  const void*                                functionDataPtr,
        QUESO::GslVector&                    qoiValues,
        QUESO::DistArray<QUESO::GslVector*>* gradVectors,
        QUESO::DistArray<QUESO::GslMatrix*>* hessianMatrices,
        QUESO::DistArray<QUESO::GslVector*>* hessianEffects)
{
  // Logic just to avoid warnings from INTEL compiler
  const QUESO::GslVector* aux1 = paramDirection;
  if (aux1) {};
  QUESO::DistArray<QUESO::GslVector*>* aux2 = gradVectors;
  if (aux2) {};
  aux2 = hessianEffects;
  QUESO::DistArray<QUESO::GslMatrix*>* aux3 = hessianMatrices;
  if (aux3) {};

  // Just checking: the user, at the application level, expects
  // vector 'paramValues' to have size 2 and
  // vector 'qoiValues' to have size 1.
  UQ_FATAL_TEST_MACRO(paramValues.sizeGlobal() != 2,
                      QUESO::UQ_UNAVAILABLE_RANK,
                      "qoiRoutine()",
                      "paramValues vector does not have size 2");

  UQ_FATAL_TEST_MACRO(qoiValues.sizeGlobal() != 1,
                      QUESO::UQ_UNAVAILABLE_RANK,
                      "qoiRoutine()",
                      "qoiValues vector does not have size 1");

  // Actual code
  //
  // This code exemplifies multiple Monte Carlo solvers, each calling this qoi routine. 
  // In this simple example, only node 0 in each sub-environment does the job even though 
  // there might be more than one node per sub-environment.
  // In a more realistic situation, if the user is asking for multiple nodes per sub-
  // environment, then the model code in the qoi routine might really demand more than one
  // node. Here we use 'env.subRank()' only. A realistic application might want to use 
  // either 'env.subComm()' or 'env.subComm().Comm()'.
  
  const QUESO::BaseEnvironment& env = paramValues.env();
  if (env.subRank() == 0) {
    double coef1 = ((qoiRoutine_DataType *) functionDataPtr)->coef1;
    double coef2 = ((qoiRoutine_DataType *) functionDataPtr)->coef2;
    qoiValues[0] = (coef1*paramValues[0] + coef2*paramValues[1]);
  }
  else {
    qoiValues[0] = 0.;
  }
  
  return;
}

double
LogNormalJointPdf<V,M>::actualValue(
  const V& domainVector,
  const V* domainDirection,
        V* gradVector,
        M* hessianMatrix,
        V* hessianEffect) const
{
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
    *m_env.subDisplayFile() << "Entering LogNormalJointPdf<V,M>::actualValue()"
                            << ", meanVector = "               << *m_lawExpVector
                            << ": domainVector = "             << domainVector
                            << ", domainVector.sizeLocal() = " << domainVector.sizeLocal()
                            << ", this->m_domainSet.vectorSpace().dimLocal() = " << this->m_domainSet.vectorSpace().dimLocal()
                            << std::endl;
  }

  UQ_FATAL_TEST_MACRO(domainVector.sizeLocal() != this->m_domainSet.vectorSpace().dimLocal(),
                      m_env.worldRank(),
                      "LogNormalJointPdf<V,M>::actualValue()",
                      "invalid input");

  UQ_FATAL_TEST_MACRO((gradVector || hessianMatrix || hessianEffect),
                      m_env.worldRank(),
                      "LogNormalJointPdf<V,M>::actualValue()",
                      "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");

  double returnValue = 0.;

  V zeroVector(domainVector);
  zeroVector.cwSet(0.);
  if (domainVector.atLeastOneComponentSmallerOrEqualThan(zeroVector)) {
    returnValue = 0.;
  }
  else if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04
    returnValue = 0.;
  }
  else {
    returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
  }
  //returnValue *= exp(m_logOfNormalizationFactor); // No need, because 'lnValue()' is called right above // [PDF-10]

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
    *m_env.subDisplayFile() << "Leaving LogNormalJointPdf<V,M>::actualValue()"
                            << ", meanVector = "   << *m_lawExpVector
                            << ": domainVector = " << domainVector
                            << ", returnValue = "  << returnValue
                            << std::endl;
  }

  return returnValue;
}

const D_V&
ExperimentStorage<S_V,S_M,D_V,D_M>::dataVec_transformed(unsigned int experimentId) const
{
  UQ_FATAL_TEST_MACRO(experimentId >= m_dataVecs_transformed.size(),
                      m_env.worldRank(),
                      "ExperimentStorage<S_V,S_M,P_V,P_M,Q_V,Q_M>::dataVec_transformed()",
                      "experimentId is too large");

  UQ_FATAL_TEST_MACRO(m_dataVecs_transformed[experimentId] == NULL,
                      m_env.worldRank(),
                      "ExperimentStorage<S_V,S_M,P_V,P_M,Q_V,Q_M>::dataVec_transformed()",
                      "vector is NULL");

  return *(m_dataVecs_transformed[experimentId]);
}

const V&
BaseTKGroup<V,M>::preComputingPosition(unsigned int stageId) const
{
  UQ_FATAL_TEST_MACRO(m_preComputingPositions.size() <= stageId,
                      m_env.worldRank(),
                      "BaseTKGroup<V,M>::preComputingPosition()",
                      "m_preComputingPositions.size() <= stageId");

  UQ_FATAL_TEST_MACRO(m_preComputingPositions[stageId] == NULL,
                      m_env.worldRank(),
                      "BaseTKGroup<V,M>::preComputingPosition()",
                      "m_preComputingPositions[stageId] == NULL");

  return *m_preComputingPositions[stageId];
}

//------------------------------------------------------
/// The actual (user-defined) qoi routine
//------------------------------------------------------
void
qoiRoutine(
  const uqGslVectorClass&                    paramValues,
  const uqGslVectorClass*                    paramDirection,
  const void*                                functionDataPtr,
        uqGslVectorClass&                    qoiValues,
        uqDistArrayClass<uqGslVectorClass*>* gradVectors,
        uqDistArrayClass<uqGslMatrixClass*>* hessianMatrices,
        uqDistArrayClass<uqGslVectorClass*>* hessianEffects)
{
  const uqBaseEnvironmentClass& env = paramValues.env();

  if (paramDirection && 
      gradVectors    &&
      hessianEffects &&
      hessianMatrices) {
    // Logic just to avoid warnings from INTEL compiler
  }

  // The user, at the application level, should have set
  // the vector 'paramValues' to have size 1 and
  // the vector 'qoiValues' to have size 1.
  UQ_FATAL_TEST_MACRO(paramValues.sizeGlobal() != 1,
                      env.fullRank(),
                      "qoiRoutine()",
                      "paramValues vector does not have size 1");

  UQ_FATAL_TEST_MACRO(qoiValues.sizeGlobal() != 1,
                      env.fullRank(),
                      "qoiRoutine()",
                      "qoiValues vector does not have size 1");
  
  // Compute qoi(s)
  double g = paramValues[0]; // Sample of the RV 'gravity acceleration'
  double distanceTraveled = 0.;
  if (env.subRank() == 0) {
    double velocity = ((qoiRoutine_DataClass *) functionDataPtr)->m_initialVelocity;
    double heights  = ((qoiRoutine_DataClass *) functionDataPtr)->m_initialHeight;
    double alpha    = ((qoiRoutine_DataClass *) functionDataPtr)->m_angle;
    
    double aux       = velocity * sin(alpha);
    distanceTraveled = (velocity * cos(alpha) / g) * ( aux + sqrt(pow(aux,2) + 2.*g*heights) );
  }

  qoiValues[0] = distanceTraveled;
    
  return;  
}

//*****************************************************
// PiecewiseLinear 1D->1D class
//*****************************************************
uqPiecewiseLinear1D1DFunctionClass::uqPiecewiseLinear1D1DFunctionClass(
  double                     minDomainValue,
  double                     maxDomainValue,
  const std::vector<double>& referenceDomainValues,
  double                     referenceImageValue0,
  const std::vector<double>& rateValues)
  :
  uqBase1D1DFunctionClass(minDomainValue,maxDomainValue),
  m_numRefValues         (referenceDomainValues.size()),
  m_referenceDomainValues(referenceDomainValues),
  m_rateValues           (rateValues)
{
  UQ_FATAL_TEST_MACRO(m_numRefValues == 0,
                      UQ_UNAVAILABLE_RANK,
                      "uqPiecewiseLinear1D1DFunctionClass::constructor()",
                      "num ref values = 0");

  UQ_FATAL_TEST_MACRO(m_numRefValues != rateValues.size(),
                      UQ_UNAVAILABLE_RANK,
                      "uqPiecewiseLinear1D1DFunctionClass::constructor()",
                      "num rate values is inconsistent");

  for (unsigned int i = 1; i < m_numRefValues; ++i) { // Yes, from '1'
    UQ_FATAL_TEST_MACRO(m_referenceDomainValues[i] <= m_referenceDomainValues[i-1],
                        UQ_UNAVAILABLE_RANK,
                        "uqPiecewiseLinear1D1DFunctionClass::constructor()",
                        "reference domain values are inconsistent");
  }

  m_referenceImageValues.clear();
  m_referenceImageValues.resize(m_numRefValues,0.);
  m_referenceImageValues[0] = referenceImageValue0;
  for (unsigned int i = 1; i < m_numRefValues; ++i) { // Yes, from '1'
    m_referenceImageValues[i] = m_referenceImageValues[i-1] + m_rateValues[i-1]*(m_referenceDomainValues[i] - m_referenceDomainValues[i-1]);
  }

  if (false) { // For debug only
    std::cout << "In uqPiecewiseLinear1D1DFunctionClass::constructor():"
              << std::endl;
    for (unsigned int i = 0; i < m_numRefValues; ++i) {
      std::cout << "i = " << i
                << ", m_referenceDomainValues[i] = " << m_referenceDomainValues[i]
                << ", m_referenceImageValues[i] = "  << m_referenceImageValues[i]
                << ", m_rateValues[i] = "            << m_rateValues[i]
                << std::endl;
    }
  }
}

int main(int argc, char* argv[])
{
  //***********************************************************************
  // Initialize MPI
  //***********************************************************************
  MPI_Init(&argc,&argv);

  //***********************************************************************
  // Initialize QUESO environment
  //***********************************************************************
  UQ_FATAL_TEST_MACRO((argc < 2),
                      UQ_UNAVAILABLE_RANK,
                      "main()",
                      "run as <executable> 'inputFileName'");
  uqFullEnvironmentClass* env = new uqFullEnvironmentClass(MPI_COMM_WORLD,argv[1],"",NULL);

  //***********************************************************************
  // Run program
  //***********************************************************************
  solveSip(*env);

  //***********************************************************************
  // Finalize QUESO environment
  //***********************************************************************
  delete env;

  //***********************************************************************
  // Finalize MPI
  //***********************************************************************
  MPI_Finalize();

  return 0;
}

//*****************************************************
// WignerInverseChebyshev1st 1D quadrature class
//*****************************************************
uqWignerInverseChebyshev1st1DQuadratureClass::uqWignerInverseChebyshev1st1DQuadratureClass(
  double       minDomainValue,
  double       maxDomainValue,
  unsigned int order)
  :
  uqBase1DQuadratureClass(minDomainValue,maxDomainValue,order)
{
  m_positions.resize(m_order+1,0.); // Yes, '+1'
  m_weights.resize  (m_order+1,0.); // Yes, '+1'

  // http://en.wikipedia.org/wiki/Chebyshev-Gauss_quadrature
  switch (m_order) {
    default:
      UQ_FATAL_TEST_MACRO(true,
                          UQ_UNAVAILABLE_RANK,
                          "uqWignerInverseChebyshev1st1DQuadratureClass::constructor()",
                          "order not supported");
    break;
  }

  // Scale positions from the interval [-1, 1] to the interval [min,max]
  for (unsigned int j = 0; j < m_positions.size(); ++j) {
    m_positions[j] = .5*(m_maxDomainValue - m_minDomainValue)*m_positions[j] + .5*(m_maxDomainValue + m_minDomainValue);
    m_weights[j] *= .5*(m_maxDomainValue - m_minDomainValue);
  }
}

int main(int argc, char* argv[])
{
  // Initialize environment
  MPI_Init(&argc,&argv);

  UQ_FATAL_TEST_MACRO(argc != 3,
                      QUESO::UQ_UNAVAILABLE_RANK,
                      "main()",
                      "after executable argv[0], input file must be specified in command line as argv[1], then numModes (1 or 2) must be specified as argv[2]");

  QUESO::FullEnvironment* env =
    new QUESO::FullEnvironment(MPI_COMM_WORLD,argv[1],"",NULL);

  // Compute
  unsigned int numModes = (unsigned int) atoi(argv[2]);

  compute(*env,numModes);

  // Finalize environment
  delete env;
  MPI_Finalize();

  std::cout << std::endl << "FIM!" << std::endl << std::endl;

  return 0;
}

double
InverseGammaJointPdf<V,M>::lnValue(
  const V& domainVector,
  const V* domainDirection,
        V* gradVector,
        M* hessianMatrix,
        V* hessianEffect) const
{
  UQ_FATAL_TEST_MACRO((domainDirection || gradVector || hessianMatrix || hessianEffect),
                      m_env.worldRank(),
                      "InverseGammaJointPdf<V,M>::lnValue()",
                      "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");

  double result = 0.;
  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
    result -= (m_alpha[i]+1.)*log(domainVector[i]);
    result -= m_beta[i]/domainVector[i];
    if (m_normalizationStyle == 0) {
      // Code needs to be done yet
    }
  }
  result += m_logOfNormalizationFactor; // [PDF-07]

  return result;
}

void
SimulationModelOptions::scanOptionsValues()
{
  UQ_FATAL_TEST_MACRO(m_optionsDesc == NULL,
                      m_env.worldRank(),
                      "SimulationModelOptions::scanOptionsValues()",
                      "m_optionsDesc variable is NULL");

  defineMyOptions                (*m_optionsDesc);
  m_env.scanInputFileForMyOptions(*m_optionsDesc);
  //std::cout << "scan 000\n"
  //          << std::endl;
  getMyOptionValues              (*m_optionsDesc);
  //std::cout << "scan 001\n"
  //          << std::endl;

  if (m_env.subDisplayFile() != NULL) {
    *m_env.subDisplayFile() << "In SimulationModelOptions::scanOptionsValues()"
                            << ": after reading values of options with prefix '" << m_prefix
                            << "', state of  object is:"
                            << "\n" << *this
                            << std::endl;
  }

  return;
}