Commit cc2a9d90 authored by Sébastien Villemot's avatar Sébastien Villemot

Global reindentation of MEX source files

parent f4557cb1
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
......@@ -47,8 +47,8 @@ class Interpreter : public SparseMatrix
private:
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
int EQN_lag1, EQN_lag2, EQN_lag3;
mxArray* GlobalTemporaryTerms;
protected:
mxArray *GlobalTemporaryTerms;
protected:
double pow1(double a, double b);
double divide(double a, double b);
double log1(double a);
......@@ -57,10 +57,10 @@ private:
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0, int block = -1);
int simulate_a_block(const int size, const int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
void print_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
const int Block_List_Max_Lead, const int u_count_int, int block);
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
const int Block_List_Max_Lead, const int u_count_int, int block);
vector<Block_contain_type> Block_Contain;
code_liste_type code_liste;
it_code_type it_code;
......@@ -81,14 +81,38 @@ public:
double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, mxArray* GlobalTemporaryTerms_arg);
bool global_temporary_terms_arg, bool print_arg, mxArray *GlobalTemporaryTerms_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, int block, int &nb_blocks);
inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
inline vector<double> get_residual() {return residual;};
inline mxArray* get_Temporary_Terms() {return GlobalTemporaryTerms;};
inline mxArray *
get_jacob(int block_num)
{
return jacobian_block[block_num];
};
inline mxArray *
get_jacob_exo(int block_num)
{
return jacobian_exo_block[block_num];
};
inline mxArray *
get_jacob_exo_det(int block_num)
{
return jacobian_det_exo_block[block_num];
};
inline mxArray *
get_jacob_other_endo(int block_num)
{
return jacobian_other_endo_block[block_num];
};
inline vector<double>
get_residual()
{
return residual;
};
inline mxArray *
get_Temporary_Terms()
{
return GlobalTemporaryTerms;
};
};
#endif
......@@ -85,9 +85,9 @@ Mem_Mngr::mxMalloc_NZE()
mexEvalString("drawnow;");
}
if (NZE_Mem_add)
NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/
NZE_Mem_add = (NonZeroElem **) mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to redefine the size of pointer on the memory*/
else
NZE_Mem_add = (NonZeroElem **) mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to define the size of pointer on the memory*/
NZE_Mem_add = (NonZeroElem **) mxMalloc(CHUNK_SIZE*sizeof(NonZeroElem *)); /*We have to define the size of pointer on the memory*/
if (!NZE_Mem_add)
{
mexPrintf("Not enough memory available\n");
......@@ -95,7 +95,7 @@ Mem_Mngr::mxMalloc_NZE()
}
for (i = CHUNK_heap_pos; i < CHUNK_SIZE; i++)
{
NZE_Mem_add[i] = (NonZeroElem *)(NZE_Mem+(i-CHUNK_heap_pos));
NZE_Mem_add[i] = (NonZeroElem *) (NZE_Mem+(i-CHUNK_heap_pos));
}
i = CHUNK_heap_pos++;
return (NZE_Mem_add[i]);
......@@ -109,7 +109,7 @@ Mem_Mngr::mxFree_NZE(void *pos)
size_t gap;
for (i = 0; i < Nb_CHUNK; i++)
{
gap = ((size_t)(pos)-(size_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
gap = ((size_t) (pos)-(size_t) (NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
if ((gap < CHUNK_BLCK_SIZE) && (gap >= 0))
break;
}
......
This diff is collapsed.
......@@ -27,12 +27,11 @@
#include <ctime>
#ifdef OCTAVE_MEX_FILE
#define CHAR_LENGTH 1
# define CHAR_LENGTH 1
#else
#define CHAR_LENGTH 2
# define CHAR_LENGTH 2
#endif
#include "Mem_Mngr.hh"
#include "ErrorHandling.hh"
#define NEW_ALLOC
......@@ -73,16 +72,16 @@ public:
private:
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray* x0_m);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray* x0_m);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, mxArray *x0_m);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m, bool &zero_solution, mxArray *x0_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM, bool &zero_solution);
void End_GE(int Size);
void Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bool symbolic, int Block_number);
void Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, bool steady_state, int it_);
void Solve_Matlab_Relaxation(mxArray* A_m, mxArray* b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state, mxArray* x0_m);
void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray* x0_m, bool steady_state);
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, bool steady_state, mxArray *x0_m);
void Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double slowc, int block, bool is_two_boundaries, int it_, mxArray *x0_m, bool steady_state);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
, long int *ndiv, long int *nsub
......@@ -111,14 +110,13 @@ private:
#endif
);
double simple_bksub(int it_, int Size, double slowc_l);
mxArray* Sparse_transpose(mxArray* A_m);
mxArray* Sparse_mult_SAT_SB(mxArray* A_m, mxArray* B_m);
mxArray* Sparse_mult_SAT_B(mxArray* A_m, mxArray* B_m);
mxArray* mult_SAT_B(mxArray* A_m, mxArray* B_m);
mxArray* Sparse_substract_SA_SB(mxArray* A_m, mxArray* B_m);
mxArray* Sparse_substract_A_SB(mxArray* A_m, mxArray* B_m);
mxArray* substract_A_B(mxArray* A_m, mxArray* B_m);
mxArray *Sparse_transpose(mxArray *A_m);
mxArray *Sparse_mult_SAT_SB(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_mult_SAT_B(mxArray *A_m, mxArray *B_m);
mxArray *mult_SAT_B(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_substract_SA_SB(mxArray *A_m, mxArray *B_m);
mxArray *Sparse_substract_A_SB(mxArray *A_m, mxArray *B_m);
mxArray *substract_A_B(mxArray *A_m, mxArray *B_m);
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
......@@ -153,7 +151,7 @@ private:
protected:
vector<double> residual;
int u_count_alloc, u_count_alloc_save;
vector<double*> jac;
vector<double *> jac;
double *jcb;
double res1, res2, max_res;
int max_res_idx;
......
This diff is collapsed.
......@@ -30,10 +30,10 @@
using namespace std;
#if !defined(DYN_MEX_FUNC_ERR_MSG_TXT)
#define DYN_MEX_FUNC_ERR_MSG_TXT(str) \
# define DYN_MEX_FUNC_ERR_MSG_TXT(str) \
do { \
mexPrintf("%s\n", str); \
} while(0)
} while (0)
#endif
typedef unsigned int mwIndex;
......@@ -64,66 +64,111 @@ enum mxArray_type
mxUINT64_CLASS = 14
};
class mxArray_tag {
public:
class mxArray_tag
{
public:
unsigned int size_1, size_2;
mwIndex *Ir, *Jc;
int Nzmax;
double *data;
mxArray_type type;
vector<string> field_name;
vector<vector<mxArray_tag*> > field_array;
vector<vector<mxArray_tag *> > field_array;
mxArray_tag();
};
typedef mxArray_tag mxArray;
void load_global(char* file_name);
void load_global(char *file_name);
/*Matlab clone function*/
int mexPrintf(const char *str, ...);
void mexErrMsgTxt(const string str);
void mexWarnMsgTxt(const string str);
void* mxMalloc(unsigned int amount);
void* mxRealloc(void *to_extend, unsigned int amount);
void* mxCalloc(unsigned int nb_elements, unsigned int amount_per_element);
void *mxMalloc(unsigned int amount);
void *mxRealloc(void *to_extend, unsigned int amount);
void *mxCalloc(unsigned int nb_elements, unsigned int amount_per_element);
void mxFree(void *to_release);
void mexEvalString(const string str);
double* mxGetPr(const mxArray *b_m);
inline mwIndex* mxGetIr(const mxArray *A_m) {return (mwIndex*)A_m->Ir;};
inline mwIndex* mxGetJc(const mxArray *A_m) {return (mwIndex*)A_m->Jc;};
mxArray* mxSetNzmax(mxArray* A_m, mwSize nz_max);
inline int mxGetNzmax(const mxArray *A_m) {return A_m->Nzmax;};
inline int mxGetM(const mxArray *A_m) {return A_m->size_1;};
inline int mxGetN(const mxArray *A_m) {return A_m->size_2;};
inline bool mxIsChar(const mxArray *A) {return A->type == mxCHAR_CLASS;};
inline bool mxIsFloat(const mxArray *A) {return A->type == mxDOUBLE_CLASS || A->type == mxSINGLE_CLASS;};
inline bool mxIsStruct(const mxArray *A) {return A->type == mxSTRUCT_CLASS;};
mxArray* mxCreateDoubleMatrix(unsigned int rows,unsigned int cols, mxData_type mx_type);
mxArray* mxCreateSparse(unsigned int rows, unsigned int cols, unsigned int nz_max, mxData_type mx_type);
mxArray* mxCreateCharArray(unsigned int rows, unsigned int cols, mxData_type mx_type);
mxArray* mxCreateDoubleScalar(double value);
mxArray* mxCreateStructMatrix(unsigned int rows, unsigned int cols, unsigned int nfields, const string &fieldnames);
inline mxArray* mxCreateStructArray(unsigned int rows, mwSize* cols, int nfields, const string &fieldnames) {return mxCreateStructMatrix(rows, *cols, nfields, fieldnames);};
mxArray* mxCreatNULLMatrix();
void mexCallMATLAB(unsigned int n_lhs, mxArray* lhs[], unsigned int n_rhs, mxArray* rhs[], const char* function);
void mxDestroyArray(mxArray* A_m);
mxArray* read_struct(FILE *fid);
mxArray* read_Array(FILE *fid);
mxArray* read_double_array(FILE *fid);
mxArray* mexGetVariable(const char* space_name, const char* matrix_name);
int mxGetFieldNumber(const mxArray* Struct, const char* field_name);
mxArray* mxGetFieldByNumber(mxArray* Struct, unsigned int pos, unsigned int field_number);
double *mxGetPr(const mxArray *b_m);
inline mwIndex *
mxGetIr(const mxArray *A_m)
{
return (mwIndex *) A_m->Ir;
};
inline mwIndex *
mxGetJc(const mxArray *A_m)
{
return (mwIndex *) A_m->Jc;
};
mxArray *mxSetNzmax(mxArray *A_m, mwSize nz_max);
inline int
mxGetNzmax(const mxArray *A_m)
{
return A_m->Nzmax;
};
inline int
mxGetM(const mxArray *A_m)
{
return A_m->size_1;
};
inline int
mxGetN(const mxArray *A_m)
{
return A_m->size_2;
};
inline bool
mxIsChar(const mxArray *A)
{
return A->type == mxCHAR_CLASS;
};
inline bool
mxIsFloat(const mxArray *A)
{
return A->type == mxDOUBLE_CLASS || A->type == mxSINGLE_CLASS;
};
inline bool
mxIsStruct(const mxArray *A)
{
return A->type == mxSTRUCT_CLASS;
};
mxArray *mxCreateDoubleMatrix(unsigned int rows, unsigned int cols, mxData_type mx_type);
mxArray *mxCreateSparse(unsigned int rows, unsigned int cols, unsigned int nz_max, mxData_type mx_type);
mxArray *mxCreateCharArray(unsigned int rows, unsigned int cols, mxData_type mx_type);
mxArray *mxCreateDoubleScalar(double value);
mxArray *mxCreateStructMatrix(unsigned int rows, unsigned int cols, unsigned int nfields, const string &fieldnames);
inline mxArray *
mxCreateStructArray(unsigned int rows, mwSize *cols, int nfields, const string &fieldnames)
{
return mxCreateStructMatrix(rows, *cols, nfields, fieldnames);
};
mxArray *mxCreatNULLMatrix();
void mexCallMATLAB(unsigned int n_lhs, mxArray* lhs[], unsigned int n_rhs, mxArray* rhs[], const char *function);
void mxDestroyArray(mxArray *A_m);
mxArray *read_struct(FILE *fid);
mxArray *read_Array(FILE *fid);
mxArray *read_double_array(FILE *fid);
mxArray *mexGetVariable(const char *space_name, const char *matrix_name);
int mxGetFieldNumber(const mxArray *Struct, const char *field_name);
mxArray *mxGetFieldByNumber(mxArray *Struct, unsigned int pos, unsigned int field_number);
void mxSetFieldByNumber(mxArray *Struct, mwIndex index, unsigned int field_number, mxArray *pvalue);
int mxAddField(mxArray* Struct, const char* field_name);
int mxAddField(mxArray *Struct, const char *field_name);
void mxSetFieldByNumber(mxArray *Struct, mwIndex index, unsigned int fieldnumber, mxArray *pvalue);
inline unsigned int mxGetNumberOfElements(const mxArray* A) {return A->size_1 * A->size_2;};
inline unsigned int
mxGetNumberOfElements(const mxArray *A)
{
return A->size_1 * A->size_2;
};
int mxGetString(const mxArray *array, char *buf, unsigned int buflen);
inline double mxGetScalar(const mxArray *array) {if (!mxIsFloat(array)) mexErrMsgTxt("not a float array\n");return array->data[0];};
mxArray* mxDuplicateArray(const mxArray *array);
void Free_simple_array(mxArray* array);
void Free_struct_array(mxArray* array);
void Free_Array(mxArray* array);
inline double
mxGetScalar(const mxArray *array)
{
if (!mxIsFloat(array))
mexErrMsgTxt("not a float array\n");
return array->data[0];
};
mxArray *mxDuplicateArray(const mxArray *array);
void Free_simple_array(mxArray *array);
void Free_struct_array(mxArray *array);
void Free_Array(mxArray *array);
void Free_global();
#endif
......@@ -82,7 +82,7 @@ extern "C" {
CONST_BLDOU beta, BLDOU y, CONST_BLINT incy);
#define dsymv FORTRAN_WRAPPER(dsymv)
void dsymv(BLCHAR uplo, CONST_BLINT m, CONST_BLDOU alpha, CONST_BLDOU a,
void dsymv(BLCHAR uplo, CONST_BLINT m, CONST_BLDOU alpha, CONST_BLDOU a,
CONST_BLINT lda, CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
BLDOU c, CONST_BLINT ldc);
......
......@@ -106,7 +106,7 @@ extern "C" {
#define dppsv FORTRAN_WRAPPER(dppsv)
void dppsv(LACHAR uplo, CONST_LAINT n, CONST_LAINT m, LADOU a, LADOU b, CONST_LAINT ldb,
LAINT info);
LAINT info);
#define dpotri FORTRAN_WRAPPER(dpotri)
void dpotri(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
......
......@@ -36,14 +36,14 @@ typedef int mwSize;
* Fix for trac ticket Ticket #137
*/
#if !defined(DYN_MEX_FUNC_ERR_MSG_TXT)
#define DYN_MEX_FUNC_ERR_MSG_TXT(str) \
# define DYN_MEX_FUNC_ERR_MSG_TXT(str) \
do { \
mexPrintf("%s\n", str); \
int i; \
for (i = 0; i < nlhs; i++) \
plhs[i] = mxCreateDoubleScalar(1); \
return; \
} while(0)
} while (0)
#endif
#endif
......@@ -70,7 +70,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
back_inserter(zeta_dynamic));
// Compute beta_back and pi_back
for(size_t i = 0; i < n_back_mixed; i++)
for (size_t i = 0; i < n_back_mixed; i++)
if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_back_mixed[i])
== zeta_mixed.end())
pi_back.push_back(i);
......@@ -78,7 +78,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
beta_back.push_back(i);
// Compute beta_fwrd and pi_fwrd
for(size_t i = 0; i < n_fwrd_mixed; i++)
for (size_t i = 0; i < n_fwrd_mixed; i++)
if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_fwrd_mixed[i])
== zeta_mixed.end())
pi_fwrd.push_back(i);
......@@ -95,7 +95,7 @@ DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (
assert(g_u.getRows() == n && g_u.getCols() == p);
// Construct S, perform QR decomposition and get A = Q*jacobian
for(size_t i = 0; i < n_static; i++)
for (size_t i = 0; i < n_static; i++)
mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, i);
A = MatrixConstView(jacobian, 0, 0, n, n_back_mixed + n + n_fwrd_mixed);
......
......@@ -27,7 +27,7 @@
EstimatedParameter::EstimatedParameter(const EstimatedParameter::pType type_arg,
size_t ID1_arg, size_t ID2_arg, const std::vector<size_t> &subSampleIDs_arg,
double lower_bound_arg, double upper_bound_arg, Prior* prior_arg) :
double lower_bound_arg, double upper_bound_arg, Prior *prior_arg) :
ptype(type_arg), ID1(ID1_arg), ID2(ID2_arg),
lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), prior(prior_arg),
subSampleIDs(subSampleIDs_arg)
......
......@@ -43,7 +43,7 @@ public:
EstimatedParameter(const EstimatedParameter::pType type,
size_t ID1, size_t ID2, const std::vector<size_t> &subSampleIDs,
double lower_bound, double upper_bound, Prior* prior
double lower_bound, double upper_bound, Prior *prior
);
virtual ~EstimatedParameter();
......@@ -52,7 +52,7 @@ public:
size_t ID2;
double lower_bound;
double upper_bound;
Prior* prior;
Prior *prior;
std::vector<size_t> subSampleIDs;
};
......
......@@ -55,7 +55,8 @@
*
* Time indices follow C convention: first period has index 0.
*/
class EstimationSubsample {
class EstimationSubsample
{
public:
EstimationSubsample(size_t startPeriod, size_t endPeriod);
virtual ~EstimationSubsample();
......
......@@ -44,7 +44,7 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile
g_x(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()),
g_u(n_endo_arg, n_exo_arg),
Rt(n_exo_arg, zeta_varobs_back_mixed.size()),
RQ(zeta_varobs_back_mixed.size(), n_exo_arg)
RQ(zeta_varobs_back_mixed.size(), n_exo_arg)
{
std::vector<size_t> zeta_back_mixed;
set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
......@@ -53,12 +53,12 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &dynamicDllFile
for (size_t i = 0; i < zeta_back_mixed.size(); i++)
pi_bm_vbm.push_back(find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(),
zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin());
}
// initialise parameter dependent KF matrices only but not Ps
void
InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
const Matrix &Q, Matrix &RQRt, Matrix &T,
const Matrix &Q, Matrix &RQRt, Matrix &T,
double &penalty, const MatrixConstView &dataView,
MatrixView &detrendedDataView, int &info)
{
......@@ -80,7 +80,6 @@ InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepPa
setPstar(Pstar, Pinf, T, RQRt, info);
}
void
InitializeKalmanFilter::setT(Matrix &T, int &info)
{
......@@ -104,27 +103,27 @@ InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, c
{
try
{
// disclyap_fast(T, RQR, Pstar, lyapunov_tol, 0 or 1 to check chol)
discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0);
{
// disclyap_fast(T, RQR, Pstar, lyapunov_tol, 0 or 1 to check chol)
discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0);
Pinf.setAll(0.0);
}
catch(const DiscLyapFast::DLPException &e)
{
if (e.info > 0) // The matrix is not positive definite in NormCholesky calculator
{
printf(e.message.c_str());
info = -1; //likelihood = penalty;
return;
}
else if (e.info < 0)
{
printf("Caugth unhandled TS exception with Pstar matrix: ");
printf(e.message.c_str());
info = -1; //likelihood = penalty;
throw;
}
}
Pinf.setAll(0.0);
}
catch (const DiscLyapFast::DLPException &e)
{
if (e.info > 0) // The matrix is not positive definite in NormCholesky calculator
{
printf(e.message.c_str());
info = -1; //likelihood = penalty;
return;
}
else if (e.info < 0)
{
printf("Caugth unhandled TS exception with Pstar matrix: ");
printf(e.message.c_str());
info = -1; //likelihood = penalty;
throw;
}
}
}
......@@ -43,7 +43,7 @@ public:
/*!
\param[in] zeta_varobs_back_mixed_arg The union of indices of observed, backward and mixed variables
*/
InitializeKalmanFilter(const std::string& dynamicDllFile, size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg,
InitializeKalmanFilter(const std::string &dynamicDllFile, size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
const std::vector<size_t> &zeta_varobs_back_mixed_arg,
double qz_criterium_arg, double lyapunov_tol_arg, int &info);
......
......@@ -37,8 +37,8 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz
double qz_criterium_arg, const std::vector<size_t> &varobs_arg,
double riccati_tol_arg, double lyapunov_tol_arg, int &info) :
zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)),
Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), Zt(Z.getCols(),Z.getRows()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo),
Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()),
Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), Zt(Z.getCols(), Z.getRows()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo),
Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()),
RQRt(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Ptmp(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), F(varobs_arg.size(), varobs_arg.size()),
Finv(varobs_arg.size(), varobs_arg.size()), K(zeta_varobs_back_mixed.size(), varobs_arg.size()), KFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()),
oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), a_init(zeta_varobs_back_mixed.size()),
......@@ -54,7 +54,7 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz
size_t j = find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(),
varobs_arg[i]) - zeta_varobs_back_mixed.begin();
Z(i, j) = 1.0;
Zt(j, i) = 1.0;
Zt(j, i) = 1.0;
}
}
......@@ -79,27 +79,27 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
VectorView &vll, MatrixView &detrendedDataView,
size_t start, size_t period, double &penalty, int &info)
{
double lik=INFINITY;
double lik = INFINITY;
try
{
if(period==0) // initialise all KF matrices
if (period == 0) // initialise all KF matrices
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
penalty, dataView, detrendedDataView, info);
else // initialise parameter dependent KF matrices only but not Ps
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
penalty, dataView, detrendedDataView, info);
penalty, dataView, detrendedDataView, info);
else // initialise parameter dependent KF matrices only but not Ps
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
penalty, dataView, detrendedDataView, info);
lik= filter(detrendedDataView, H, vll, start, info);
lik = filter(detrendedDataView, H, vll, start, info);
}
catch (const DecisionRules::BlanchardKahnException &bke)
{
info =22;
info = 22;
return penalty;
}
if (info != 0)
return penalty;
else
else
return lik;
};
......@@ -110,7 +110,7 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
double
KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info)
{
double loglik=0.0, ll, logFdet, Fdet, dvtFinvVt;
double loglik = 0.0, ll, logFdet, Fdet, dvtFinvVt;
size_t p = Finv.getRows();
bool nonstationary = true;
a_init.setAll(0.0);
......@@ -130,25 +130,25 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
// Finv=inv(F)
mat::set_identity(Finv);
// Pack F upper trinagle as vector
for (size_t i=1;i<=p;++i)
for(size_t j=i;j<=p;++j)
FUTP(i + (j-1)*j/2 -1) = F(i-1,j-1);
for (size_t i = 1; i <= p; ++i)
for (size_t j = i; j <= p; ++j)
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);
info=lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
if (info<0)
info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
if (info < 0)
{
std::cout << "Info:" << info << std::endl;
std::cout << "t:" << t << std::endl;
return 0;
}
if (info>0)
if (info > 0)
{
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
mat::transpose(Ptmp, Pstar);
mat::add(Pstar,Ptmp);
for (size_t i=0;i<Pstar.getCols();++i)
for(size_t j=0;j<Pstar.getCols();++j)
Pstar(i,j)*=0.5;
mat::add(Pstar, Ptmp);
for (size_t i = 0; i < Pstar.getCols(); ++i)
for (size_t j = 0; j < Pstar.getCols(); ++j)
Pstar(i, j) *= 0.5;
// K=PZ'
//blas::gemm("N", "T", 1.0, Pstar, Z, 0.0, K);
......@@ -161,12 +161,12 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
// Finv=inv(F)
mat::set_identity(Finv);
// Pack F upper trinagle as vector
for (size_t i=1;i<=p;++i)
for(size_t j=i;j<=p;++j)
FUTP(i + (j-1)*j/2 -1) = F(i-1,j-1);
for (size_t i = 1; i <= p; ++i)
for (size_t j = i; j <= p; ++j)
FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1);