Case-Study II: d-dimensional Gaussian distribution

To test the performance of DREAM in the presence of high-dimensionality, the second case study involves a d = 100-dimensional multivariate normal distribution

where the vector stores the mean of the normal distribution, and the matrix signifies the variance-covariance matrix. The probability density function of is given by

We assume a mean of zero, , and define the covariance matrix, , in such a way that the variance of the th variable is equal to th variable is equal to and all pairwise correlations are equivalent to 0.5 The main diagonal thus has values of . The correlation, , between two entries and of can be calculated from their respective covariance value and individual variances according to

Entry of matrix can thus be calculated using

where the variances of and are found on the main diagonal of , and thus equivalent to entries and , respectively.

The initial population is drawn from the multivariate uniform distribution, , using Latin hypercube sampling. We assume a lower bound of each dimension equal to and upper bound of . The initial value of each parameter {1,…,d} is thus drawn from this range. We use N = 100 chains with DREAM and apply default values of the algorithmic variables, with the exception that we use thinning to reduce memory storage. Indeed, only every 10th sample is stored in each Markov chain.

Implementation of plugin functions

The complete source code can be found in DREAM SDK - Examples\D1\Drm_D1_Example02\Plugin\Src_Cpp

///////////////////////////////////////////////////////////////////////////////////////////////////

// EvaluatorModel.cpp

///////////////////////////////////////////////////////////////////////////////////////////////////

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample02::GetInputParams

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParams* pIDreamInputParams // Pointer to DREAM input data to be defined

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 The purpose of this function is to define default DREAM parameters for this particular project.

 This function is called by the main program to initialize data of each newly created simulation

 case. This data can be later modified in the GUI (except parameters locked by flag

 DreamSolver::disable_edit) or can be overwritten while reading case data from a file.

 In the AppTest console application, there is no GUI and this functions actually defines all

 input parameters that are subsequently used by DREAM Solver.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Always check pIDreamInputParams pointer

 if(pIDreamInputParams == nullptr)

 {

   return E_INVALIDARG;

 }

 // Set default input parameters

 pIDreamInputParams->SetProblemDimension(100, DreamSolver::enable_edit);

 pIDreamInputParams->SetNumberOfMarkovChains(50, DreamSolver::enable_edit);

 pIDreamInputParams->SetNumberOfGenerations(10000, DreamSolver::enable_edit);

 pIDreamInputParams->SetThinnigSampleToStore(10, DreamSolver::enable_edit);

 pIDreamInputParams->SetLikelihoodChoice(eLikelihood::LIK_102, DreamSolver::disable_edit);

 pIDreamInputParams->SetInitialDistribution(eInitDistrib::eInitLatin, DreamSolver::enable_edit);

 // Example can run in parallel mode (recommended mode)

 pIDreamInputParams->SetParallelMode(eParallelModeYes, DreamSolver::enable_edit);

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample02::InitData

//-------------------------------------------------------------------------------------------------

(

 IDreamPluginInfo* pIPluginInfo            // Plugin information

, IDreamInputParViewer* pIDreamInputParams  // Current DREAM input parameters

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 This function is called to prepare:

 A/ DREAM input data before accessing them via the IDreamDataSource interface and

 B/ all other local data required by evaluator for the calculation.  

 It can be also used to read needed data from a file in directory pIPluginInfo->m_strModelDataDir

 (directory ...Drm_ExampleXX\Model\Data).

         

 Remark 1: This function prepares DREAM input data in the evaluator according to current

 DREAM parameters defined by pIPluginInfo interface. Note that some parameters can be

 different from default parameters defined in function CEvaluatorExample::GetInputParams().

 Remark 2: This function should be called before using functions

 CEvaluatorExample::GetMinMaxData, ::GetNormalData, ::GetPriorData, ::GetPriorDataCustom,

 ::GetMeasurementData and ::GetBayesData. However, function CEvaluatorExample::GetInputParams()

 can be called independently.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Check parameters

 if(pIPluginInfo == nullptr || pIDreamInputParams == nullptr)

 {

   // Invalid argument

   return E_INVALIDARG;

 }

 // Free memory from previous calculation

 FunctionsSC::DeleteMatrix(m_invC, m_invCSize);

 int d = pIDreamInputParams->GetProblemDimension();

 // Construct the d x d covariance matrix

 double** matrix1 = nullptr;

 FunctionsSC::IdentityMatrix(matrix1, d, d, 0.5);

 double** matrix2 = nullptr;

 FunctionsSC::CreateMatrix(matrix2, d, d, 0.5);

 double** A = nullptr;

 FunctionsSC::MatrixMatrixAddition(A, matrix1, d, d, matrix2, d, d);

 // Rescale to variance - covariance matrix of interest

 double** C = nullptr;

 FunctionsSC::CreateMatrix(C, d);

 for(int i = 0; i < d; i++)

 {

   for(int j = 0; j < d; j++)

   {

     C[i][j] = A[i][j] * sqrt((i + 1.0) * (j + 1.0));

   }

 }

 m_invCSize = d;

 FunctionsSC::CopyMatrix(m_invC, C, d);

 if(!FunctionsSC::Inv(m_invC, d))

 {

   return E_FAIL;

 }

 // Calculate integration constant

 double temp = pow(2 * pi<double>(), -(double)d / 2.0);

 double det = FunctionsSC::DeterminantLU(C, d);

 double temp2 = pow(det, -0.5);

 m_logF = log(temp * temp2);

 // log_F becomes - Inf for large d--> hence we set log_F to zero

 // Also need to resolve this in importance distribution as well!!

 if(d > 150)

 {

   m_logF = 0;

 }

 // Free allocated memory

 FunctionsSC::DeleteMatrix(matrix1, d);

 FunctionsSC::DeleteMatrix(matrix2, d);

 FunctionsSC::DeleteMatrix(A, d);

 FunctionsSC::DeleteMatrix(C, d);

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample02::GetMinMaxData

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParViewer* pInputPar // Interface to DREAM input parameters

, IDreamDataMinMax* pIMinMaxData  // Interface to Min/Max data

)

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Always check pIMinMaxData pointer

 if(pIMinMaxData == nullptr)

 {

   return E_INVALIDARG;

 }

 int d = pInputPar->GetProblemDimension();

 // Fill min max data

 for(int i = 0; i < d; i++)

 {

   pIMinMaxData->AddItem(-5.0, 15.0);

 }

 // Data added successfully

 return S_OK;

}

//-------------------------------------------------------------------------------------------------

 HRESULT STDMETHODCALLTYPE CEvaluatorExample02::EvaluateProposal

//-------------------------------------------------------------------------------------------------

(

 IDreamInputParViewer* pInputPar // Interface to DREAM input parameters

, int iSim              // input parameter (i-th simulation / Markov chain)

, IDreamMatrix* x       // input matrix (m_NumberOfMarkovChains x m_ProblemDimension)

, IDreamMatrix* res     // output matrix (Data_Dimension x m_NumberOfMarkovChains)

, BSTR strWorkingDir    // path to the current working directory (for current calculation thread)

, BSTR strModelDataDir  // path to the permanent directory with model data

, BSTR strModelBinDir   // path to the permanent directory with model executables

)

///////////////////////////////////////////////////////////////////////////////////////////////////

/*

 DESCRIPTION:

 The purpose of this function is to evaluate the current proposal (i.e. one single step

 of IDreamInputParams::GetNumberOfGenerations() steps) for iSim-th Markov chain. It computes

 the likelihood (or log-likelihood or the residual vector to be compared with the measurement

 - see the technical manual) for iSim-th row of input matrix "x" and returns the result

 in output matrix "res".  

 

 Parameter iSim is from interval <1, IDreamInputParams::m_NumberOfMarkovChains> except

 the case when pInputPar->GetVectorization() == eVectorization::eVectYes, when iSim = 0

 and all proposals (for all Markov chains) are evaluated in one step.

 RETURN VALUE:

 S_OK    - Operation successful

 S_FALSE - Function not implemented (program will continue)

 E_FAIL  - Operation failed (program will stop)

*/

///////////////////////////////////////////////////////////////////////////////////////////////////

{

 // Temporary vector used for calculations

 double* tmpx = nullptr;

 FunctionsSC::CopyRowFromIDreamMatrix(iSim, tmpx, x);

 double logLtmp = 0.0;

 if(x->_colCount == 1) // 1-dimensional Gaussian

 {

   // 1-dimensional twisted Gaussian probability density function

   logLtmp = m_logF - 0.5 * m_invC[0][0] * tmpx[0] * tmpx[0];

 }

 else // N-dimensional Gaussian

 {

   // N-dimensional twisted Gaussian probability density function

   // logL = m_logF - 1/2 * sum(x'.*(invC * x'));  

   // where:

   // x'  - transposition of matrix "x"

   // sum - returns the sum of the elements of input along the first array dimension

   // .*  - multiplies arrays A and B element by element

   // vector1 = invC * x'

   double* vector1 = nullptr;

   FunctionsSC::MatrixColumnVectorProduct(vector1, m_invC, m_invCSize, m_invCSize,

                                          tmpx, x->_colCount);

   // vector2 = x'.*(invC * x') = vector1 .* vector2

   double* vector2 = nullptr;

   FunctionsSC::VectorVectorElementProduct(vector2, tmpx, x->_colCount, vector1, x->_colCount);

     

   // logL = m_logF - 1/2 * sum(x'.*(invC * x')) = m_logF - 0.5 sum(vector2)

   logLtmp = m_logF - 0.5 * FunctionsSC::VectorElementSum(vector2, x->_colCount);

   // Free allocated memory

   FunctionsSC::DeleteVector(vector1);

   FunctionsSC::DeleteVector(vector2);

 }

 // Replace content of res by tmp values

 res->_data[iSim] = logLtmp;

 // Free allocated memory

 FunctionsSC::DeleteVector(tmpx);

 return S_OK;

}