stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removed obsolete Kalman filter code

time-shift
Sébastien Villemot 2010-09-24 13:01:52 +02:00
parent da1c031195
commit 5504811b1a
87 changed files with 0 additions and 11345 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
mex/sources/kalman/cc/Makefile
View File

@ -1,64 +0,0 @@
# $Id: Makefile 531 2005-11-30 13:49:48Z kamenik $
# Copyright 2005, Ondra Kamenik
#DEBUG = yes
LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
CC_FLAGS := -DMATLAB -DWINDOWS -DNO_BLAS_H -DNO_LAPACK_H \
-fpic -Wall -I../sylv/cc \
-I$(MATLAB_PATH)/extern/include #-pg
ifeq ($(DEBUG),yes)
CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
else
CC_FLAGS := $(CC_FLAGS) -O3
endif
# Added by GP
# LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
LDFLAGS := -Wl,--library-path $(LD_LIBRARY_PATH) \
-Wl,-L$(MATLAB_PATH)/extern/lib/win32/microsoft/ \
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
-lg2c -lmingw32 -lstdc++ $(LDFLAGS)
# -Wl,-L'f:/CygWin/usr/local/atlas/lib' \
# -Wl,-L'f:/CygWin/lib' \
# $(LDFLAGS)
LD_LIBS :=$(LDFLAGS)
# end add
matrix_interface := GeneralMatrix Vector SylvException
matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
matcppsource := $(patsubst %, ../sylv/cc/%.cpp, $(matrix_interface))
cwebsource := $(wildcard *.cweb)
# cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource))
cppsource := $(wildcard *.cpp)
hsource := $(wildcard *.h)
objects := $(patsubst %.cpp,%.o,$(cppsource))
hwebsource := $(wildcard *.hweb)
dummy.ch:
touch dummy.ch
# %.cpp: %.cweb dummy.ch
# ctangle -bhp $*.cweb dummy.ch $*.cpp
# %.h: %.hweb dummy.ch
# ctangle -bhp $*.hweb dummy.ch $*.h
%.o: %.cpp $(hsource) $(mathsource)
c++ $(CC_FLAGS) -c $*.cpp
all: $(objects) # $(cppsource) $(hsource)
doc: main.web $(hwebsource) $(cwebsource)
cweave -bhp main.web
pdftex main
mv main.pdf ts.pdf
clear:
rm -f *.o
rm -f *.{pdf,dvi,log,scn,idx,toc}
# rm -f *.cpp
# rm -f *.h

102
mex/sources/kalman/cc/disclyap_fast.cpp
View File

@ -1,102 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/****************************************************************
% function X=disclyap_fast(G,V,tol,ch)
%
% Solve the discrete Lyapunov Equation
% X=G*X*G'+V
% Using the Doubling Algorithm
%
% INPUT:
% G and V - square General matrices of same size
% tol - double tollerance level
% flag_ch - integer flag: if 1 check if the result is positive
% definite and generate an error message if it is not
% OUTPUT:
% on completion V - square General matrice contains solution
%
% based on work of Joe Pearlman and Alejandro Justiniano
% 3/5/2005
% C++ version 28/07/09 by Dynare team
****************************************************************/
#include "ts_exception.h"
#include "cppblas.h"
#include "GeneralMatrix.h"
//#include "Vector.h"
#include "SylvException.h"
#include "utils.h"
#include "mex.h"
void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol = 1e-16, int flag_ch=0)
{
/**
if nargin == 2 | isempty( ch ) == 1
flag_ch = 0;
else
flag_ch = 1;
end
**/
//P0=V;
GeneralMatrix P0(V);
//A0=G;
GeneralMatrix A0(G);
//n=size(G,1);
int n= A0.numCols();
const double alpha=1.0;
const double half=0.5;
const double neg_alpha=-1.0;
const double omega=0.0;
GeneralMatrix A1(n,n);
GeneralMatrix Ptmp(n,n);
GeneralMatrix P1(P0);
GeneralMatrix I(n,n);
I.unit();
bool matd=true;
while (matd ) // matrix diff > tol
{
//P1=P0+A0*P0*A0';
// first step Ptmp=P0*A0';
// DGEMM: C := alpha*op( A )*op( B ) + beta*C,
BLAS_dgemm("N", "T", &n, &n, &n, &alpha, P0.base(), &n,
A0.base(), &n, &omega, Ptmp.base(), &n);
// P1=P0+A0*Ptmp;
BLAS_dgemm("N", "N", &n, &n, &n, &alpha, A0.base(), &n,
Ptmp.base(), &n, &alpha, P1.base(), &n);
// A1=A0*A0;
// A0=A1 (=A0*A0);
// A0.multRight(A0);
BLAS_dgemm("N", "N", &n, &n, &n, &alpha, A0.base(), &n,
A0.base(), &n, &omega, A1.base(), &n);
// check if max( max( abs( P1 - P0 ) ) )>tol
matd=P0.isDiffSym(P1, tol);
P0=P1;
A0=A1;
}//end while
// X=P0=(P0+P0')/2;
BLAS_dgemm("T", "N", &n, &n, &n, &half, P1.base(), &n,
I.base(), &n, &half, P0.base(), &n);
X=P0;
// Check that X is positive definite
if (flag_ch==1)
NormCholesky chol(P0);
}

34
mex/sources/kalman/cc/disclyap_fast.h
View File

@ -1,34 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/****************************************************************
% function X=disclyap_fast(G,V,ch)
%
% Solve the discrete Lyapunov Equation
% X=G*X*G'+V
% Using the Doubling Algorithm
%
% If ch is defined then the code will check if the resulting X
% is positive definite and generate an error message if it is not
%
% based on work of Joe Pearlman and Alejandro Justiniano
% 3/5/2005
% C++ version 28/07/09 by Dynare team
****************************************************************/
#include "GeneralMatrix.h"
void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol, int ch);

2481
mex/sources/kalman/cc/kalman.cpp
View File

File diff suppressed because it is too large Load Diff

339
mex/sources/kalman/cc/kalman.h
View File

@ -1,339 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
/*************
This file defines two most important classes: |KalmanTask| and
|KalmanUniTask|. Both define a linear filtering and smoothing problem;
one for multivariate observations, the other for univariate
observations. The Kalman task consists of three things: the state
space, observed data, and initialization.
In addition to these two classes, we have also an abstraction for
filter results, for diffuse filter results and smoother results. These
are only containers of series of matrices and vectors. The filter
results are used as an input for smoothing.
The important convention is that whenever a parameter $t$ stands for a
time period (which is always), it starts from 1 (not zero as in
C/C++).
In this file, we use the same naming scheme as in Durbin \& Koopman.
*************/
#ifndef KALMAN_H
#define KALMAN_H
#include "state_init.h"
#include "ssf.h"
#include "ssf_uni.h"
#include "utils.h"
#include <vector>
/*************
This is an output of the filtering and input for the smoothing. This
class is inherited by |DiffuseFilterResults| which enriches the
|FilterResults| with additional information coming from the diffuse
periods.
The object is initialized for a given number of periods, and maintains
a number of periods for which the filter results were set. If this
number, |maxt|, is less than the number of overall periods, it means
that the filter has not finished.
**************/
class FilterResults{
protected:
std::vector<PLUFact*> Finv;
std::vector<Vector*> v;
std::vector<GeneralMatrix*> L;
std::vector<Vector*> a;
std::vector<GeneralMatrix*> P;
std::vector<double> loglik;
int maxt;
public:
FilterResults(int num)
:Finv(num,(PLUFact*)0),v(num,(Vector*)0),L(num,(GeneralMatrix*)0),
a(num,(Vector*)0),P(num,(GeneralMatrix*)0),loglik(num,0.0),maxt(0)
{}
virtual~FilterResults();
void set(int t,const PLUFact&FFinv,const Vector&vv,
const GeneralMatrix&LL,const Vector&aa,
const GeneralMatrix&P,double ll);
int getMaxT()const
{return maxt;}
bool hasFinished()const
{return maxt==(int)Finv.size();}
const PLUFact&getFInverse(int t)const;
const Vector&getV(int t)const;
const GeneralMatrix&getL(int t)const;
const Vector&getA(int t)const;
const GeneralMatrix&getP(int)const;
double getLogLikelihood()const;
double getLogLikelihood(int start)const;
double getLogLikelihood(std::vector<double> *vloglik)const;
double getLogLikelihood(int start,std::vector<double> *vloglik)const;
};
class DiffuseFilterResults:public FilterResults{
protected:
std::vector<GeneralMatrix*> L_1;
std::vector<GeneralMatrix*> Pinf;
std::vector<GeneralMatrix*> F_2;
std::vector<bool> Finf_reg;
std::vector<bool> Pinf_zero;
public:
DiffuseFilterResults(int num)
:FilterResults(num),
L_1(num,(GeneralMatrix*)0),
Pinf(num,(GeneralMatrix*)0),
F_2(num,(GeneralMatrix*)0),
Finf_reg(num,true),
Pinf_zero(num,true)
{}
virtual~DiffuseFilterResults();
void set(int t,const PLUFact&FFinfinv,const GeneralMatrix&FF_2,
const Vector&vv,const GeneralMatrix&LL_0,
const GeneralMatrix&LL_1,const Vector&aa,
const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
double ll);
void set(int t,const PLUFact&FFstarinv,const Vector&vv,
const GeneralMatrix&LL_0,const Vector&aa,
const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
double ll);
int getDiffusePeriods()const;
bool isFinfRegular(int t)const;
bool isPinfZero(int t)const;
const PLUFact&getFinfInverse(int t)const;
const PLUFact&getFstarInverse(int t)const;
const GeneralMatrix&getF2(int t)const;
const GeneralMatrix&getL0(int t)const
{return getL(t);}
const GeneralMatrix&getL1(int t)const;
const GeneralMatrix&getPstar(int t)const
{return getP(t);}
const GeneralMatrix&getPinf(int t)const;
};
class SmootherResults{
protected:
std::vector<Vector*> alpha;
std::vector<Vector*> eta;
std::vector<GeneralMatrix*> V;
int mint;
public:
SmootherResults(int num)
:alpha(num,(Vector*)0),eta(num,(Vector*)0),
V(num,(GeneralMatrix*)0),
mint(num+1){}
virtual~SmootherResults();
void set(int t,const Vector&aalpha,const Vector&eeta,
const GeneralMatrix&VV);
void import(const SmootherResults&sres,int period);
void exportAlpha(GeneralMatrix&out)const;
void exportEta(GeneralMatrix&out)const;
void exportV(GeneralMatrix&out)const;
};
class BasicKalmanTask{
// friend class KalmanUniTask;
// SSForm ssf;
const GeneralMatrix &data;
const ConstGeneralMatrix &Zt;
const ConstGeneralMatrix &Ht;
const ConstGeneralMatrix &Tt;
const ConstGeneralMatrix &Rt;
const ConstGeneralMatrix &Qt;
const StateInit&init;
const double riccatiTol;
public:
BasicKalmanTask(const GeneralMatrix&d,const GeneralMatrix&ZZ,
const GeneralMatrix&HH,const GeneralMatrix&TT,
const GeneralMatrix&RR,const GeneralMatrix&QQ,
const StateInit&init_state, const double riccatiTol);
// BasicKalmanTask(const GeneralMatrix&d,const TMatrix&Z,
// const TMatrix&H,const TMatrix&T,
// const TMatrix&R,const TMatrix&Q,
// const StateInit&init_state);
BasicKalmanTask(const GeneralMatrix&d,const ConstGeneralMatrix&ZZ,
const ConstGeneralMatrix&HH,const ConstGeneralMatrix&TT,
const ConstGeneralMatrix&RR,const ConstGeneralMatrix&QQ,
const StateInit&init_state, const double riccatiTol);
virtual ~BasicKalmanTask();
// double filter(int&per,int&d)const;
// double filter(int&per,int&d, int start, std::vector<double>* vll)const;
double filter(int&per,int&d,int start, std::vector<double>* vll)const;
// double filter_and_smooth(SmootherResults&sres,int&per,int&d)const;
protected:
double filterNonDiffuse(const Vector&a,const GeneralMatrix&Pstar,
int start, std::vector<double>* vll) const; //int first,FilterResults&fres)const;
// void filterDiffuse(const Vector&a,const GeneralMatrix&Pstar,
// const GeneralMatrix&Pinf,int first,
// DiffuseFilterResults&fres)const;
// void smootherNonDiffuse(const FilterResults&fres,SmootherResults&sres)const;
// void smootherDiffuse(const DiffuseFilterResults&fres,SmootherResults&sres)const;
// void smootherNonDiffuseStep(int t,const FilterResults&fres,
// Vector&rt,GeneralMatrix&Nt,
// Vector&alphat,GeneralMatrix&Vt,
// Vector&etat)const;
static double calcStepLogLik(const PLUFact&Finv,const Vector&v);
};
class KalmanUniTask;
class KalmanTask{
friend class KalmanUniTask;
SSForm ssf;
ConstGeneralMatrix data;
const StateInit&init;
public:
KalmanTask(const GeneralMatrix&d,const GeneralMatrix&Z,
const GeneralMatrix&H,const GeneralMatrix&T,
const GeneralMatrix&R,const GeneralMatrix&Q,
const StateInit&init_state);
KalmanTask(const GeneralMatrix&d,const TMatrix&Z,
const TMatrix&H,const TMatrix&T,
const TMatrix&R,const TMatrix&Q,
const StateInit&init_state);
double filter(int&per,int&d)const;
double filter(int&per,int&d, int start, std::vector<double>* vll)const;
double filter_and_smooth(SmootherResults&sres,int&per,int&d)const;
protected:
void filterNonDiffuse(const Vector&a,const GeneralMatrix&Pstar,
int first,FilterResults&fres)const;
void filterDiffuse(const Vector&a,const GeneralMatrix&Pstar,
const GeneralMatrix&Pinf,int first,
DiffuseFilterResults&fres)const;
void smootherNonDiffuse(const FilterResults&fres,SmootherResults&sres)const;
void smootherDiffuse(const DiffuseFilterResults&fres,SmootherResults&sres)const;
void smootherNonDiffuseStep(int t,const FilterResults&fres,
Vector&rt,GeneralMatrix&Nt,
Vector&alphat,GeneralMatrix&Vt,
Vector&etat)const;
static double calcStepLogLik(const PLUFact&Finv,const Vector&v);
};
class FilterUniResults{
protected:
std::vector<double> F;
std::vector<double> v;
std::vector<GeneralMatrix*> L;
std::vector<Vector*> a;
std::vector<GeneralMatrix*> P;
std::vector<double> loglik;
int maxt;
public:
FilterUniResults(int num)
:F(num,0.0),v(num,0.0),L(num,(GeneralMatrix*)0),
a(num,(Vector*)0),P(num,(GeneralMatrix*)0),loglik(num,0.0),maxt(0)
{}
virtual~FilterUniResults();
void set(int t,double F,double vv,
const GeneralMatrix&LL,const Vector&aa,
const GeneralMatrix&P,double ll);
int getMaxT()const
{return maxt;}
bool hasFinished()const
{return maxt==(int)F.size();}
double getF(int t)const;
double getV(int t)const;
const GeneralMatrix&getL(int t)const;
const Vector&getA(int t)const;
const GeneralMatrix&getP(int)const;
double getLogLikelihood()const;
double getLogLikelihood(int start)const;
double getLogLikelihood(std::vector<double>* vloglik)const;
double getLogLikelihood(int start,std::vector<double>* vloglik)const;
};
class DiffuseFilterUniResults:public FilterUniResults{
protected:
std::vector<GeneralMatrix*> L_1;
std::vector<GeneralMatrix*> Pinf;
std::vector<double> F_2;
std::vector<bool> Finf_reg;
std::vector<bool> Pinf_zero;
public:
DiffuseFilterUniResults(int num)
:FilterUniResults(num),
L_1(num,(GeneralMatrix*)0),
Pinf(num,(GeneralMatrix*)0),
F_2(num,0.0),
Finf_reg(num,true),
Pinf_zero(num,true)
{}
virtual~DiffuseFilterUniResults();
void set(int t,double FFinf,double FF_2,
double vv,const GeneralMatrix&LL_0,
const GeneralMatrix&LL_1,const Vector&aa,
const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
double ll);
void set(int t,double FFstar,double vv,
const GeneralMatrix&LL_0,const Vector&aa,
const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
double ll);
int getDiffusePeriods()const;
bool isFinfRegular(int t)const;
bool isPinfZero(int t)const;
double getFinf(int t)const;
double getFstar(int t)const;
double getF2(int t)const;
const GeneralMatrix&getL0(int t)const
{return getL(t);}
const GeneralMatrix&getL1(int t)const;
const GeneralMatrix&getPstar(int t)const
{return getP(t);}
const GeneralMatrix&getPinf(int t)const;
};
class KalmanUniTask{
private:
MesEquation me;
protected:
SSFormUni ssf;
ConstGeneralMatrix data;
const StateInit&init;
public:
KalmanUniTask(const KalmanTask&kt);
double filter(int&per,int&d)const;
double filter(int&per,int&d, int start, std::vector<double>* vll)const;
double filter_and_smooth(SmootherResults&sres,int&per,int&d)const;
protected:
void filterNonDiffuse(const Vector&a,const GeneralMatrix&Pstar,
int first,FilterUniResults&fres)const;
void filterDiffuse(const Vector&a,const GeneralMatrix&Pstar,
const GeneralMatrix&Pinf,int first,
DiffuseFilterUniResults&fres)const;
void smootherNonDiffuse(const FilterUniResults&fres,SmootherResults&sres)const;
void smootherDiffuse(const DiffuseFilterUniResults&fres,SmootherResults&sres)const;
void smootherNonDiffuseStep(int t,const FilterUniResults&fres,
Vector&rt,GeneralMatrix&Nt,
Vector&alphat,GeneralMatrix&Vt,
Vector&etat)const;
static double calcStepLogLik(double F,double v);
};
#endif

323
mex/sources/kalman/cc/ssf.cpp
View File

@ -1,323 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#include "ssf.h"
#include "ts_exception.h"
#include "utils.h"
#include <vector>
TMatrixCycle::TMatrixCycle(int n,int nr,int nc)
:matrices(new GeneralMatrix*[n]),num(n),nrows(nr),ncols(nc)
{
for(int i= 0;i<num;i++)
matrices[i]= NULL;
}
;
TMatrixCycle::TMatrixCycle(const TMatrixCycle&m)
:matrices(new GeneralMatrix*[m.num]),num(m.num),
nrows(m.nrows),ncols(m.ncols)
{
for(int i= 0;i<num;i++)
if(m.matrices[i])
matrices[i]= new GeneralMatrix(*(m.matrices[i]));
else
matrices[i]= NULL;
}
;
TMatrixCycle::TMatrixCycle(const GeneralMatrix&m)
:matrices(new GeneralMatrix*[m.numRows()]),num(m.numRows()),
nrows(1),ncols(m.numCols())
{
for(int i= 0;i<num;i++)
matrices[i]= new GeneralMatrix(m,i,0,1,ncols);
}
;
TMatrixCycle::TMatrixCycle(const TMatrix&m,const char*dummy)
:matrices(new GeneralMatrix*[m.numRows()*m.period()]),
num(m.numRows()*m.period()),nrows(1),ncols(m.numCols())
{
TS_RAISE_IF(m.period()==0,
"Infinite period in TMatrixCycle constructor");
for(int i= 0;i<m.period();i++)
for(int j= 0;j<m.numRows();j++)
matrices[i*m.numRows()+j]
= new GeneralMatrix(m[i],j,0,1,ncols);
}
;
TMatrixCycle::~TMatrixCycle()
{
for(int i= 0;i<num;i++)
delete matrices[i];
delete[]matrices;
}
;
const GeneralMatrix&TMatrixCycle::operator[](int t)const
{
int i= (t-1)%num;
TS_RAISE_IF(matrices[i]==NULL,
"The matrix has not ben set in TMatrixCycle::operator[]");
return*(matrices[i]);
}
GeneralMatrix&TMatrixCycle::operator[](int t)
{
int i= (t-1)%num;
TS_RAISE_IF(matrices[i]==NULL,
"The matrix has not ben set in TMatrixCycle::operator[]");
return*(matrices[i]);
}
;
void TMatrixCycle::set(int t,const GeneralMatrix&m)
{
TS_RAISE_IF(m.numRows()!=numRows()||m.numCols()!=numCols(),
"Wrong matrix dimensions for TMatrixCycle::set");
int i= (t-1)%num;
if(matrices[i])
delete matrices[i];
matrices[i]= new GeneralMatrix(m);
}
;
TMatrixPadUnit::TMatrixPadUnit(const TMatrix&m,int s)
:tmat(m.clone()),skip(s),unit(m.numRows(),m.numRows())
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"TMatrix not square in TMatrixPadUnit constructor");
unit.zeros();
for(int i= 0;i<numRows();i++)
unit.get(i,i)= 1.0;
}
;
const GeneralMatrix&TMatrixPadUnit::operator[](int t)const
{
if(isUnit(t))
return unit;
else
return(*tmat)[t/skip];
}
GeneralMatrix&TMatrixPadUnit::operator[](int t)
{
TS_RAISE_IF(isUnit(t),
"Attempt to return non-const unit in TMatrixPadUnit::operator[]");
return(*tmat)[t/skip];
}
TMatrixPadZero::TMatrixPadZero(const TMatrix&m,int s)
:tmat(m.clone()),skip(s),zero(m.numRows(),m.numCols())
{
zero.zeros();
}
const GeneralMatrix&TMatrixPadZero::operator[](int t)const
{
if(isZero(t))
return zero;
else
return(*tmat)[t/skip];
}
GeneralMatrix&TMatrixPadZero::operator[](int t)
{
TS_RAISE_IF(isZero(t),
"Attempt to return non-const zero in TMatrixPadZero::operator[]");
return(*tmat)[t/skip];
}
SSForm::SSForm(const TMatrix&zz,const TMatrix&hh,const TMatrix&tt,
const TMatrix&rr,const TMatrix&qq)
:Z(zz.clone()),
H(hh.clone()),
T(tt.clone()),
R(rr.clone()),
Q(qq.clone()),
p(zz.numRows()),m(zz.numCols()),r(qq.numRows())
{
TS_RAISE_IF(T->numRows()!=m||T->numCols()!=m||
H->numRows()!=p||H->numCols()!=p||
R->numRows()!=m||R->numCols()!=r||
Q->numCols()!=r,
"Wrong TMatrix dimension in SSForm constructor");
}
SSForm::SSForm(const GeneralMatrix&zz,const GeneralMatrix&hh,
const GeneralMatrix&tt,const GeneralMatrix&rr,
const GeneralMatrix&qq)
:Z(new TMatrixInvariant(zz)),
H(new TMatrixInvariant(hh)),
T(new TMatrixInvariant(tt)),
R(new TMatrixInvariant(rr)),
Q(new TMatrixInvariant(qq)),
p(zz.numRows()),m(zz.numCols()),r(qq.numRows())
{
TS_RAISE_IF(T->numRows()!=m||T->numCols()!=m||
H->numRows()!=p||H->numCols()!=p||
R->numRows()!=m||R->numCols()!=r||
Q->numCols()!=r,
"Wrong TMatrix dimension in SSForm constructor");
}
SSForm::SSForm(const SSForm&f)
:Z(f.Z->clone()),
H(f.H->clone()),
T(f.T->clone()),
R(f.R->clone()),
Q(f.Q->clone()),
p(f.p),m(f.m),r(f.r)
{}
SSForm::~SSForm()
{
delete Z;
delete H;
delete T;
delete R;
delete Q;
}
MesEquation::MesEquation(const GeneralMatrix&data,const TMatrix&zz,
const TMatrix&hh)
:y(data),
Z((zz.period()*hh.period()==1)?(TMatrix*)new TMatrixInvariant(zz[1]):
(zz.period()*hh.period()==0)?(TMatrix*)new TMatrixCycle(y.numCols(),
zz.numRows(),zz.numCols())
:(TMatrix*)new TMatrixCycle(zz.period()*hh.period(),zz.numRows(),zz.numCols())),
H((zz.period()*hh.period()==1)?(TMatrix*)new TMatrixInvariant(hh[1]):
(zz.period()*hh.period()==0)?(TMatrix*)new TMatrixCycle(y.numCols(),
hh.numRows(),hh.numCols())
:(TMatrix*)new TMatrixCycle(zz.period()*hh.period(),hh.numRows(),hh.numCols()))
{
TS_RAISE_IF(y.numRows()!=Z->numRows()||y.numRows()!=H->numRows()||
y.numRows()!=H->numCols(),
"Incompatible dimension in MesEquation constructor");
int mper= zz.period()*hh.period();
if(mper==1)
{
construct_invariant();
}
else
{
std::vector<NormCholesky*> chols;
int per= (mper==0)?y.numCols():mper;
for(int t= 1;t<=per;t++)
{
GeneralMatrix ycol(y,0,t-1,y.numRows(),1);
int hi= t;
if(hh.period()> 0)
hi= (t-1)%hh.period()+1;
NormCholesky*ch;
if(hh.period()==0)
{
ch= new NormCholesky(hh[t]);
}
else if(hi-1>=(int)chols.size())
{
ch= new NormCholesky(hh[t]);
chols.push_back(ch);
}
else
{
ch= chols[hi-1];
}
ch->getL().multInvLeftUnit(ycol);
if(t-1<mper)
{
GeneralMatrix Zt(zz[t]);
ch->getL().multInvLeftUnit(Zt);
((TMatrixCycle*)Z)->set(t,Zt);
GeneralMatrix Ht(hh.numRows(),hh.numRows());
Ht.zeros();
for(int i= 0;i<Ht.numRows();i++)
Ht.get(i,i)= ch->getD()[i];
((TMatrixCycle*)H)->set(t,Ht);
}
if(hh.period()==0)
delete ch;
}
for(unsigned int i= 0;i<chols.size();i++)
delete chols[i];
}
}
MesEquation::MesEquation(const GeneralMatrix&data,const GeneralMatrix&zz,
const GeneralMatrix&hh)
:y(data),
Z(new TMatrixInvariant(zz)),
H(new TMatrixInvariant(hh))
{
TS_RAISE_IF(y.numRows()!=Z->numRows()||y.numRows()!=H->numRows()||
y.numRows()!=H->numCols(),
"Incompatible dimension in MesEquation constructor");
construct_invariant();
}
MesEquation::~MesEquation()
{
delete Z;
delete H;
}
void MesEquation::construct_invariant()
{
if(!TSUtils::isDiagonal((*H)[1]))
{
NormCholesky chol((*H)[1]);
chol.getL().multInvLeftUnit(y);
chol.getL().multInvLeftUnit((*Z)[1]);
(*H)[1].zeros();
for(int i= 0;i<H->numRows();i++)
(*H)[1].get(i,i)= chol.getD()[i];
}
}
;

193
mex/sources/kalman/cc/ssf.h
View File

@ -1,193 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef SSF_H
#define SSF_H
#include "GeneralMatrix.h"
class TMatrix{
public:
virtual const GeneralMatrix&operator[](int t)const= 0;
virtual GeneralMatrix&operator[](int t)= 0;
virtual int numRows()const= 0;
virtual int numCols()const= 0;
virtual int period()const= 0;
virtual bool isZero(int t)const= 0;
virtual bool isUnit(int t)const= 0;
virtual~TMatrix(){}
virtual TMatrix*clone()const= 0;
};
class TMatrixInvariant:public TMatrix,public GeneralMatrix{
public:
TMatrixInvariant(const GeneralMatrix&m)
:GeneralMatrix(m){}
TMatrixInvariant(const TMatrixInvariant&m)
:GeneralMatrix(m){}
const GeneralMatrix&operator[](int t)const
{return*this;}
GeneralMatrix&operator[](int t)
{return*this;}
int numRows()const
{return GeneralMatrix::numRows();}
int numCols()const
{return GeneralMatrix::numCols();}
int period()const
{return 1;}
bool isZero(int t)const
{return false;}
bool isUnit(int t)const
{return false;}
TMatrix*clone()const
{return new TMatrixInvariant(*this);}
};
class TMatrixCycle:public TMatrix{
protected:
GeneralMatrix**const matrices;
int num;
int nrows;
int ncols;
public:
TMatrixCycle(int n,int nr,int nc);
TMatrixCycle(const TMatrixCycle&m);
TMatrixCycle(const GeneralMatrix&m);
TMatrixCycle(const TMatrix&m,const char*dummy);
~TMatrixCycle();
const GeneralMatrix&operator[](int t)const;
GeneralMatrix&operator[](int t);
int numRows()const
{return nrows;}
int numCols()const
{return ncols;}
int period()const
{return num;}
bool isZero(int t)const
{return false;}
bool isUnit(int t)const
{return false;}
TMatrix*clone()const
{return new TMatrixCycle(*this);}
void set(int t,const GeneralMatrix&m);
};
class TMatrixPadUnit:public TMatrix{
TMatrix*const tmat;
int skip;
GeneralMatrix unit;
public:
TMatrixPadUnit(const TMatrix&m,int s);
TMatrixPadUnit(const TMatrixPadUnit&m)
:tmat(m.tmat->clone()),skip(m.skip),unit(m.unit){}
~TMatrixPadUnit()
{delete tmat;}
const GeneralMatrix&operator[](int t)const;
GeneralMatrix&operator[](int t);
int numRows()const
{return tmat->numRows();}
int numCols()const
{return tmat->numCols();}
int period()const
{return skip*tmat->period();}
bool isZero(int t)const
{return false;}
bool isUnit(int t)const
{return(t/skip)*skip!=t;}
TMatrix*clone()const
{return new TMatrixPadUnit(*this);}
};
class TMatrixPadZero:public TMatrix{
TMatrix*const tmat;
int skip;
GeneralMatrix zero;
public:
TMatrixPadZero(const TMatrix&m,int s);
TMatrixPadZero(const TMatrixPadZero&m)
:tmat(m.tmat->clone()),skip(m.skip),zero(m.zero){}
~TMatrixPadZero()
{delete tmat;}
const GeneralMatrix&operator[](int t)const;
GeneralMatrix&operator[](int t);
int numRows()const
{return tmat->numRows();}
int numCols()const
{return tmat->numCols();}
int period()const
{return skip*tmat->period();}
bool isUnit(int t)const
{return false;}
bool isZero(int t)const
{return(t/skip)*skip!=t;}
TMatrix*clone()const
{return new TMatrixPadZero(*this);}
};
struct SSForm{
TMatrix*const Z;
TMatrix*const H;
TMatrix*const T;
TMatrix*const R;
TMatrix*const Q;
const int p;
const int m;
const int r;
SSForm(const TMatrix&zz,const TMatrix&hh,const TMatrix&tt,
const TMatrix&rr,const TMatrix&qq);
SSForm(const GeneralMatrix&zz,const GeneralMatrix&hh,
const GeneralMatrix&tt,const GeneralMatrix&rr,
const GeneralMatrix&qq);
SSForm(const SSForm&f);
~SSForm();
};
struct MesEquation{
GeneralMatrix y;
TMatrix*const Z;
TMatrix*const H;
MesEquation(const GeneralMatrix&data,const GeneralMatrix&zz,
const GeneralMatrix&hh);
MesEquation(const GeneralMatrix&data,const TMatrix&zz,
const TMatrix&hh);
~MesEquation();
protected:
void construct_invariant();
};
;
#endif

166
mex/sources/kalman/cc/ssf_uni.cpp
View File

@ -1,166 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#include "ssf_uni.h"
#include "ts_exception.h"
TScalarCycle::TScalarCycle(int n)
:ss(new double[n]),flags(new bool[n]),num(n)
{
for(int i= 0;i<num;i++)
flags[i]= false;
}
;
TScalarCycle::TScalarCycle(const TScalarCycle&c)
:ss(new double[c.num]),flags(new bool[c.num]),num(c.num)
{
for(int i= 0;i<num;i++){
flags[i]= c.flags[i];
ss[i]= c.ss[i];
}
}
;
TScalarCycle::TScalarCycle(const GeneralMatrix&m)
:ss(new double[m.numRows()]),flags(new bool[m.numRows()]),
num(m.numRows())
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"Matrix is not diagonal in TScalarCycle diagonal constructor");
for(int i= 0;i<m.numRows();i++){
ss[i]= m.get(i,i);
flags[i]= true;
}
}
;
TScalarCycle::TScalarCycle(const TMatrix&m)
:ss(new double[m.numRows()*m.period()]),
flags(new bool[m.numRows()*m.period()]),
num(m.numRows()*m.period())
{
TS_RAISE_IF(m.period()==0,
"Infinite period in TScalarCycle diagonal constructor");
TS_RAISE_IF(m.numRows()!=m.numCols(),
"TMatrix is not diagonal in TScalarCycle diagonal constructor");
for(int i= 0;i<m.period();i++)
for(int j= 0;j<m.numRows();j++){
ss[i*m.numRows()+j]= m[i].get(j,j);
flags[i*m.numRows()+j]= true;
}
}
;
TScalarCycle::~TScalarCycle()
{
delete[]flags;
delete[]ss;
}
;
const double&TScalarCycle::operator[](int t)const
{
int i= (t-1)%num;
TS_RAISE_IF(!flags[i],
"The scalar has not been set in TScalarCycle::operator[]");
return ss[i];
}
;
void TScalarCycle::set(int t,double s)
{
int i= (t-1)%num;
flags[i]= true;
ss[i]= s;
}
;
SSFormUni::SSFormUni(const TMatrix&zz,const TScalar&hh,const TMatrix&tt,
const TMatrix&rr,const TMatrix&qq)
:Z(zz.clone()),
H(hh.clone()),
T(tt.clone()),
R(rr.clone()),
Q(qq.clone()),
m(zz.numCols()),r(qq.numRows())
{
TS_RAISE_IF(T->numRows()!=m||T->numCols()!=m||
R->numRows()!=m||R->numCols()!=r||
Q->numCols()!=r,
"Wrong TMatrix dimension in SSFormUni constructor");
TS_RAISE_IF(Z->numRows()!=1,
"Z is not univariate in SSFormUni constructor");
}
;
SSFormUni::SSFormUni(const GeneralMatrix&zz,double hh,
const GeneralMatrix&tt,const GeneralMatrix&rr,
const GeneralMatrix&qq)
:Z(new TMatrixInvariant(zz)),
H(new TScalarInvariant(hh)),
T(new TMatrixInvariant(tt)),
R(new TMatrixInvariant(rr)),
Q(new TMatrixInvariant(qq)),
m(zz.numCols()),r(qq.numRows())
{
TS_RAISE_IF(T->numRows()!=m||T->numCols()!=m||
R->numRows()!=m||R->numCols()!=r||
Q->numCols()!=r,
"Wrong TMatrix dimension in SSFormUni constructor");
TS_RAISE_IF(Z->numRows()!=1,
"Z is not univariate in SSFormUni constructor");
}
SSFormUni::SSFormUni(const SSFormUni&ssfu)
:Z(ssfu.Z->clone()),
H(ssfu.H->clone()),
T(ssfu.T->clone()),
R(ssfu.R->clone()),
Q(ssfu.Q->clone()),
m(ssfu.m),r(ssfu.r)
{}
SSFormUni::~SSFormUni()
{
delete Z;
delete H;
delete T;
delete R;
delete Q;
}
;

101
mex/sources/kalman/cc/ssf_uni.h
View File

@ -1,101 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef SSF_UNI_H
#define SSF_UNI_H
#include "ssf.h"
class TScalar{
public:
virtual const double&operator[](int t)const= 0;
virtual~TScalar(){}
virtual int period()const= 0;
virtual TScalar*clone()const= 0;
};
;
class TScalarInvariant:public TScalar{
protected:
double s;
public:
TScalarInvariant(double ss)
:s(ss){}
TScalarInvariant(const TScalarInvariant&c)
:s(c.s){}
const double&operator[](int t)const
{return s;}
int period()const
{return 1;}
TScalar*clone()const
{return new TScalarInvariant(*this);}
};
;
class TScalarCycle:public TScalar{
protected:
double*const ss;
bool*const flags;
int num;
public:
TScalarCycle(int n);
TScalarCycle(const TScalarCycle&c);
TScalarCycle(const GeneralMatrix&m);
TScalarCycle(const TMatrix&m);
~TScalarCycle();
const double&operator[](int t)const;
int period()const
{return num;}
TScalar*clone()const
{return new TScalarCycle(*this);}
void set(int t,double s);
};
;
struct SSFormUni{
TMatrix*const Z;
TScalar*const H;
TMatrix*const T;
TMatrix*const R;
TMatrix*const Q;
const int m;
const int r;
SSFormUni(const TMatrix&zz,const TScalar&hh,const TMatrix&tt,
const TMatrix&rr,const TMatrix&qq);
SSFormUni(const GeneralMatrix&zz,double hh,
const GeneralMatrix&tt,const GeneralMatrix&rr,
const GeneralMatrix&qq);
SSFormUni(const SSFormUni&ssfu);
~SSFormUni();
};
;
#endif

59
mex/sources/kalman/cc/state_init.cpp
View File

@ -1,59 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#include "state_init.h"
#include "ts_exception.h"
#include "utils.h"
StateInit::StateInit(const GeneralMatrix&PPstar,const Vector&aa)
:m(PPstar.numRows()),ndiffuse(0),Pstar(PPstar),
Pinf(m,m),a(aa)
{
TS_RAISE_IF(Pstar.numRows()!=Pstar.numCols(),
"Pstar not square in StateInit non-diffuse constructor");
TS_RAISE_IF(m!=a.length(),
"Bad length of initial state vector in StateInit non-diffuse constructor");
Pinf.zeros();
}
;
StateInit::StateInit(const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
const Vector&aa)
:m(PPstar.numRows()),ndiffuse(0),Pstar(PPstar),
Pinf(PPinf),a(aa)
{
TS_RAISE_IF(m!=Pstar.numCols()||m!=Pinf.numRows()||
m!=Pinf.numCols()||m!=a.length(),
"Wrong dimensions for StateInit diffuse constructor");
TS_RAISE_IF(!TSUtils::isDiagonal(Pinf),
"Pinf is not diagonal in StateInit diffuse constructor");
for(int i= 0;i<m;i++)
if(Pinf.get(i,i)!=0.0)
ndiffuse++;
}
;

60
mex/sources/kalman/cc/state_init.h
View File

@ -1,60 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef STATE_INIT_H
#define STATE_INIT_H
#include "GeneralMatrix.h"
class StateInit{
const int m;
int ndiffuse;
GeneralMatrix Pstar;
GeneralMatrix Pinf;
Vector a;
public:
StateInit(const GeneralMatrix&PPstar,const Vector&aa);
StateInit(const GeneralMatrix&PPstar,const GeneralMatrix&PPinf,
const Vector&aa);
StateInit(const StateInit&init)
:m(init.m),ndiffuse(init.ndiffuse),Pstar(init.Pstar),
Pinf(init.Pinf),a(init.a){}
virtual~StateInit(){}
int getM()const
{return m;}
bool isDiffuse()const
{return ndiffuse> 0;}
const Vector&getA()const
{return a;}
const GeneralMatrix&getPstar()const
{return Pstar;}
const GeneralMatrix&getPinf()const
{return Pinf;}
int getNDiff()const
{return ndiffuse;}
};
;
#endif

64
mex/sources/kalman/cc/ts_exception.h
View File

@ -1,64 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef TS_EXCEPTION_H
#define TS_EXCEPTION_H
#include <stdio.h>
#include <string.h>
#ifdef MATLAB
#include "mex.h"
#endif
#define TS_RAISE(mes) \
throw TSException(__FILE__, __LINE__, mes);
#define TS_RAISE_IF(expr, mes) \
if (expr) throw TSException(__FILE__, __LINE__, mes);
class TSException{
char fname[50];
int lnum;
char message[500];
public:
TSException(const char*f,int l,const char*mes)
{
strncpy(fname,f,50);fname[49]= '\0';
strncpy(message,mes,500);message[499]= '\0';
lnum= l;
}
virtual void print()const
{
printf("At %s:%d:%s\n",fname,lnum,message);
#ifdef MATLAB
mexPrintf("At %s:%d:%s\n",fname,lnum,message);
#endif
}
virtual const char*getMessage()const
{return message;}
};
;
#endif

411
mex/sources/kalman/cc/utils.cpp
View File

@ -1,411 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#include "utils.h"
#include "ts_exception.h"
#include "cppblas.h"
#include "cpplapack.h"
#include <math.h>
#include <cmath>
#include <float.h>
LowerTriangle::LowerTriangle(const GeneralMatrix&m)
:GeneralMatrix(m)
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"The matrix is not square in LowerTriangle constructor");
}
void LowerTriangle::multInvLeft(GeneralMatrix&m)const
{
TS_RAISE_IF(numCols()!=m.numRows(),
"Wrong dimensions of the matrix for LowerTriangle::multInvLeft");
int mrows= m.numRows();
int mcols= m.numCols();
double alpha= 1.0;
int ld= getLD();
int ldm= m.getLD();
BLAS_dtrsm("L","L","N","N",&mrows,&mcols,&alpha,getData().base(),
&ld,m.getData().base(),&ldm);
}
;
void LowerTriangle::multInvLeftUnit(GeneralMatrix&m)const
{
TS_RAISE_IF(numCols()!=m.numRows(),
"Wrong dimensions of the matrix for LowerTriangle::multInvLeftUnit");
int mrows= m.numRows();
int mcols= m.numCols();
double alpha= 1.0;
int ld= getLD();
int ldm= m.getLD();
BLAS_dtrsm("L","L","N","U",&mrows,&mcols,&alpha,getData().base(),
&ld,m.getData().base(),&ldm);
}
;
NormCholesky::NormCholesky(const GeneralMatrix&a)
:L(a),D(a.numRows())
{
TS_RAISE_IF(a.numRows()!=a.numCols(),
"The matrix is not square in NormCholesky constructor");
int lrows= L.numRows();
int ldl= L.getLD();
int info;
LAPACK_dpotrf("L",&lrows,L.getData().base(),&ldl,&info);
TS_RAISE_IF(info> 0,
"The matrix is not positive definite in NormCholesky constructor");
TS_RAISE_IF(info<0,
"Internal error in NormCholesky constructor");
for(int i= 0;i<L.numRows();i++)
for(int j= i+1;j<L.numCols();j++)
L.get(i,j)= 0.0;
for(int j= 0;j<L.numCols();j++){
double d= L.get(j,j);
Vector Lj(L,j);
Lj.mult(1.0/d);
D[j]= d*d;
}
}
;
PLUFact::PLUFact(const GeneralMatrix&m)
:inv(m.numRows()*m.numCols()),ipiv(new int[m.numRows()]),
rows(m.numRows())
{
TS_RAISE_IF(!m.isFinite(),
"Matrix is not finite in PLUFact constructor");
TS_RAISE_IF(m.numRows()!=m.numCols(),
"Matrix not square in PLUFact constructor");
inv= m.getData();
LAPACK_dgetrf(&rows,&rows,inv.base(),&rows,ipiv,&info);
TS_RAISE_IF(info<0,
"Internal error in PLUFact constructor");
double mnorm= m.getNormInf();
double*work= new double[4*rows];
int*iwork= new int[rows];
int infotmp;
LAPACK_dgecon("I",&rows,inv.base(),&rows,&mnorm,&rcond,work,
iwork,&infotmp);
delete[]iwork;
delete[]work;
TS_RAISE_IF(infotmp<0,
"Internal error in PLUFact constructor");
;
calcDetSign();
}
PLUFact::PLUFact(const PLUFact&fact)
:inv(fact.inv),ipiv(new int[fact.rows]),
rows(fact.rows),rcond(fact.rcond),detsign(fact.detsign),info(fact.info)
{
memcpy(ipiv,fact.ipiv,rows*sizeof(int));
}
PLUFact::PLUFact(const int nc,const int nr )
:inv(nr*nc),ipiv(new int[nr]),rows(nr)
{
TS_RAISE_IF(nr!=nc,
"Matrix not square in PLUFact constructor");
}
const PLUFact&
PLUFact::operator = (const GeneralMatrix&m)
{
TS_RAISE_IF(!m.isFinite(),
"Matrix is not finite in PLUFact assignement");
TS_RAISE_IF(m.numRows()!=m.numCols(),
"Matrix not square in PLUFact assignement");
TS_RAISE_IF(m.numRows()!=rows,
"Matrix not matching PLUFact size for assignement");
inv= m.getData();
LAPACK_dgetrf(&rows,&rows,inv.base(),&rows,ipiv,&info);
TS_RAISE_IF(info<0,
"Internal error in PLUFact assignement");
double mnorm= m.getNormInf();
double*work= new double[4*rows];
int*iwork= new int[rows];
int infotmp;
LAPACK_dgecon("I",&rows,inv.base(),&rows,&mnorm,&rcond,work,
iwork,&infotmp);
delete[]iwork;
delete[]work;
TS_RAISE_IF(infotmp<0,
"Internal error in PLUFact assignement");
calcDetSign();
return *this;
}
;
void PLUFact::PL_dgetrs(const char*trans,double*b,int ldb,int bcols)const
{
if(rows> 0){
int info;
LAPACK_dgetrs(trans,&rows,&bcols,inv.base(),&rows,ipiv,b,&ldb,&info);
TS_RAISE_IF(info<0,
"Internal error in PLUFact::dgetrs");
}
}
;
void PLUFact::multInvLeft(GeneralMatrix&a)const
{
TS_RAISE_IF(rows!=a.numRows(),
"Wrong dimension of the matrix in PLUFact::multInvLeft");
PL_dgetrs("N",a.getData().base(),a.getLD(),a.numCols());
}
;
void PLUFact::multInvRight(GeneralMatrix&a)const
{
GeneralMatrix atrans(a,"trans");
TS_RAISE_IF(rows!=atrans.numRows(),
"Wrong dimension of the matrix in PLUFact::multInvRight");
PL_dgetrs("T",atrans.getData().base(),atrans.getLD(),atrans.numCols());
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
a.get(i,j)= atrans.get(j,i);
}
// pass also a temporary GM space for atrans to avoid matrix construction:
void PLUFact::multInvRight(GeneralMatrix&a, GeneralMatrix&atrans)const
{
TS_RAISE_IF(rows!=atrans.numRows(),
"Wrong dimension of the matrix in PLUFact::multInvRight");
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
atrans.get(j,i)= a.get(i,j);
PL_dgetrs("T",atrans.getData().base(),atrans.getLD(),atrans.numCols());
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
a.get(i,j)= atrans.get(j,i);
}
;
void PLUFact::multInvLeft(Vector&a)const
{
TS_RAISE_IF(rows!=a.length(),
"Wrong dimension of the vector in PLUFact::multInvLeft");
TS_RAISE_IF(a.skip()!=1,
"Not implemented error in PLUFact::multInvLeft");
PL_dgetrs("N",a.base(),a.length(),1);
}
;
void PLUFact::multInvRight(Vector&a)const
{
TS_RAISE_IF(rows!=a.length(),
"Wrong dimension of the vector in PLUFact::multInvLeft");
TS_RAISE_IF(a.skip()!=1,
"Not implemented error in PLUFact::multInvLeft");
PL_dgetrs("T",a.base(),a.length(),1);
}
;
double PLUFact::getDeterminant()const
{
double res= 1;
for(int i= 0;i<rows;i++)
res*= std::abs(inv[(rows+1)*i]);
return detsign*res;
}
;
double PLUFact::getLogDeterminant()const
{
double res= 0;
for(int i= 0;i<rows;i++)
res+= log(std::abs(inv[(rows+1)*i]));
TS_RAISE_IF(detsign==-1,
"Negative determinant in PLUFact::getLogDeterminant");
return res;
}
;
void PLUFact::calcDetSign()
{
detsign= 1;
for(int i= 0;i<rows;i++)
if(ipiv[i]!=i+1)
detsign*= -1;
for(int i= 0;i<rows;i++)
if(inv[i*(rows+1)]<0)
detsign*= -1;
}
;
void PLUFact::print()const
{
for(int i= 0;i<rows;i++)
printf(" %d",ipiv[i]);
printf("\n");
for(int i= 0;i<rows;i++){
for(int j= 0;j<rows;j++)
printf(" %15.12g",inv[j*rows+i]);
printf("\n");
}
}
;
VDVFact::VDVFact(const GeneralMatrix&m)
:V(m),D(m.numRows())
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"Matrix is not square in VDVFact constructor");
int n= m.numRows();
int lda= V.getLD();
double tmpwork;
int lwork= -1;
int info;
LAPACK_dsyev("V","U",&n,V.base(),&lda,D.base(),&tmpwork,&lwork,&info);
lwork= (int)tmpwork;
double*work= new double[lwork];
LAPACK_dsyev("V","U",&n,V.base(),&lda,D.base(),work,&lwork,&info);
delete[]work;
TS_RAISE_IF(info<0,
"Internal error in VDVFact constructor");
converged= true;
if(info)
converged= false;
}
;
bool TSUtils::isDiagonal(const ConstGeneralMatrix&m)
{
bool res= (m.numCols()==m.numRows());
for(int i= 0;i<m.numRows()&&res;i++)
for(int j= i+1;j<m.numCols()&&res;j++)
if(m.get(i,j)!=0.0||m.get(j,i)!=0.0)
res= false;
return res;
}
;
bool TSUtils::isZero(const ConstGeneralMatrix&m)
{
bool res= true;
for(int i= 0;i<m.numRows()&&res;i++)
for(int j= 0;j<m.numCols()&&res;j++)
if(m.get(i,j)!=0.0)
res= false;
return res;
}
;
bool TSUtils::hasNegativeDiagonal(const ConstGeneralMatrix&m)
{
int r= m.numRows()<m.numCols()?m.numRows():m.numCols();
bool res= false;
for(int i= 0;i<r&&!res;i++)
res= m.get(i,i)<0.0;
return res;
}
;
bool TSUtils::isSymDiagDominant(const ConstGeneralMatrix&m)
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"The matrix is not square in TSUtils::isSymDiagDominant");
bool res= true;
for(int i= 0;i<m.numRows()&&res;i++)
for(int j= i+1;j<m.numCols()&&res;j++)
res= 2*std::abs(m.get(i,j))<=
std::abs(m.get(i,i))+std::abs(m.get(j,j));
return res;
}
;
double TSUtils::correctDefinitness(GeneralMatrix&m)
{
VDVFact f(m);
if(!f.hasConverged())
return-1;
Vector d(f.getD());
double correct= 0;
int i= 0;
while(i<d.length()&&d[i]<2*DBL_EPSILON){
correct+= d[i]*d[i];
d[i]= 0.0;
i++;
}
m= f.getV();
for(int i= 0;i<d.length();i++){
Vector mi(m,i);
mi.mult(d[i]);
}
m.multRightTrans(f.getV());
return sqrt(correct);
}
;
void TSUtils::correctSymmetricity(GeneralMatrix&m)
{
TS_RAISE_IF(m.numRows()!=m.numCols(),
"Matrix is not square in TSUtils::correctSymmetricity");
GeneralMatrix tmp((const GeneralMatrix&)m,"trans");
m.add(1.0,tmp);
m.mult(0.5);
}
;

121
mex/sources/kalman/cc/utils.h
View File

@ -1,121 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef UTILS_H
#define UTILS_H
#include "GeneralMatrix.h"
class LowerTriangle:public GeneralMatrix{
public:
LowerTriangle(const GeneralMatrix&m);
LowerTriangle(const LowerTriangle&t)
:GeneralMatrix(t){}
void multInvLeft(GeneralMatrix&m)const;
void multInvLeftUnit(GeneralMatrix&m)const;
};
;
class NormCholesky{
LowerTriangle L;
Vector D;
public:
NormCholesky(const GeneralMatrix&m);
NormCholesky(const NormCholesky&chol)
:L(chol.L),D(chol.D){}
const LowerTriangle&getL()const
{return L;}
const Vector&getD()const
{return D;}
};
;
class PLUFact{
Vector inv;
int*ipiv;
int rows;
double rcond;
int detsign;
int info;
public:
PLUFact(const GeneralMatrix&m);
PLUFact(const PLUFact&plu);
PLUFact(const int nc,const int nr );
virtual~PLUFact()
{delete[]ipiv;}
const PLUFact& operator = (const GeneralMatrix&m);
void multInvLeft(GeneralMatrix&a)const;
void multInvRight(GeneralMatrix&a)const;
// pass temporary GM space for atrans to avoid matrix construction:
void multInvRight(GeneralMatrix&a, GeneralMatrix&atrans)const;
void multInvLeft(Vector&a)const;
void multInvRight(Vector&a)const;
bool isRegular()const
{return info==0;}
double getDeterminant()const;
double getLogDeterminant()const;
int getDetSign()const
{return detsign;}
int numRows()const
{return rows;}
double getRcond()const
{return rcond;}
void print()const;
private:
void PL_dgetrs(const char*trans,double*b,int ldb,int bcols)const;
void calcDetSign();
};
;
class VDVFact{
GeneralMatrix V;
Vector D;
bool converged;
public:
VDVFact(const GeneralMatrix&m);
const GeneralMatrix&getV()const
{return V;}
const Vector&getD()const
{return D;}
bool hasConverged()const
{return converged;}
};
;
struct TSUtils{
static bool isDiagonal(const ConstGeneralMatrix&m);
static bool isZero(const ConstGeneralMatrix&m);
static bool hasNegativeDiagonal(const ConstGeneralMatrix&m);
static bool isSymDiagDominant(const ConstGeneralMatrix&m);
static double correctDefinitness(GeneralMatrix&m);
static void correctSymmetricity(GeneralMatrix&m);
};
;
#endif

134
mex/sources/kalman/matlab/Makefile
View File

@ -1,134 +0,0 @@
# $Id: Makefile 531 2005-11-30 13:49:48Z kamenik $
# Copyright 2005, Ondra Kamenik
#DEBUG = yes
#LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
#CC_FLAGS := -DWINDOWS -mno-cygwin
CC_FLAGS := -DMATLAB -DNO_BLAS_H -DNO_LAPACK_H \
-Wall -fpic -I../qt/cc -I../sylv/cc -I../cc \
-I$(MATLAB_PATH)/extern/include #-pg
ifeq ($(DEBUG),yes)
CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
# CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g #-pg #-Wl,-pg
KALMANLIB := kalmanlib_dbg.a
else
# CC_FLAGS := $(CC_FLAGS) -O3
CC_FLAGS := -DTIMING_LOOP $(CC_FLAGS) -O3
KALMANLIB := kalmanlib.a
endif
# Added by GP
# LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
#LDFLAG := -Wl,--library-path $(LD_LIBRARY_PATH)
#-Wl,-L'C:/MinGW/lib/gcc-lib/i686-pc-mingw32/4.0.4'
#-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack -lf95
#-lg2c -lmingw32 kalmanlib.def -Wl,-lmwm_ir
#-Wl,-L'f:/CygWin/lib'
#LD_LIBS := -Wl,--library-path
# -Wl,-L'/usr/lib'
LD_LIBS := -Wl,-rpath-link,$(MATLAB_PATH)/bin/glnxa64 \
-Wl,-L$(MATLAB_PATH)/bin/glnxa64 \
-Wl,-lmex -lmx -lmwlapack -lmwblas -lmat -lm \
-Wl,-lstdc++ $(LDFLAGS)
#-Wl,-L'/usr/lib'
# -Wl,-L'f:/CygWin/usr/local/atlas/lib'
# -Wl,-L'f:/CygWin/lib'
# -Wl,-L'f:/MinGW/lib'
# $(LDFLAGS)
# LD_LIBS :=$(LDFLAGS)
# end add
matrix_interface := GeneralMatrix Vector SylvException
matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
matcppsource := $(patsubst %, ../sylv/cc/%.cpp, $(matrix_interface))
qtf90source := $(wildcard ../qt/f90/*.f90)
qtobjs := $(patsubst %.f90,%.o,$(qtf90source))
# cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource))
kalmancppsource := $(wildcard ../cc/*.cpp)
kalmanhsource := $(wildcard ../cc/*.h)
kalmanobjects := $(patsubst %.cpp,%.o,$(kalmancppsource))
cppsource := $(wildcard *.cpp)
hsource := $(wildcard *.h)
objects := $(patsubst %.cpp,%.o,$(cppsource))
hwebsource := $(wildcard *.hweb)
cwebsource := $(wildcard *.cweb)
dummy.ch:
touch dummy.ch
# %.cpp: %.cweb dummy.ch
# ctangle -bhp $*.cweb dummy.ch $*.cpp
# %.h: %.hweb dummy.ch
# ctangle -bhp $*.hweb dummy.ch $*.h
#%.o: %.cpp $(hsource) $(cppsource) #$(kalmanhsource) $(mathsource)
%.o: %.cpp $(hsource) $(cppsource)
c++ $(CC_FLAGS) -c $*.cpp
$(KALMANLIB): $(objects) # $(matobjs) $(qtobjs) #$(kalmanobjects)
ar cr $(KALMANLIB) $(kalmanobjects) $(matobjs) $(qtobjs)
ranlib $(KALMANLIB)
kalman_smoother_dll.dll: kalman_smoother.o $(KALMANLIB) #$(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o kalman_smoother_dll.dll kalman_smoother.o \
kalmanlib.a $(LD_LIBS)
minv.dll: minv.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o minv.dll minv.o \
kalmanlib.a $(LD_LIBS)
gmvm.dll: gmvm.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o gmvm.dll gmvm.o \
kalmanlib.a $(LD_LIBS)
qtamvm.dll: qtamvm.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o qtamvm.dll qtamvm.o \
kalmanlib.a $(LD_LIBS)
qtmvm.dll: qtmvm.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o qtmvm.dll qtmvm.o \
kalmanlib.a $(LD_LIBS)
disclyap_fast_dll.dll: disclyap_fast_dll.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o disclyap_fast_dll.dll disclyap_fast_dll.o \
$(KALMANLIB) $(LD_LIBS) kalmanlib.def
kalman_filter_dll.dll: kalman_filters.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o kalman_filter_dll.dll kalman_filters.o \
$(KALMANLIB) $(LD_LIBS)
kalman_filters_testx.exe: kalman_filters_testx.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc $(CC_FLAGS) -pg -o kalman_filters_testx.exe kalman_filters_testx.o \
$(KALMANLIB) $(LD_LIBS)
all: $(objects) $(KALMANLIB) kalman_smoother_dll.dll kalman_filter_dll.dll # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
#kalman_filter_loop.o: kalman_filters.cpp
# c++ -DTIMING_LOOP $(CC_FLAGS) -o kalman_filter_loop.o kalman_filters.cpp
#kalman_filter_loop.dll: kalman_filter_loop.o kalmanlib.a # $(hsource) $(cppsource)
# gcc -shared -DTIMING_LOOP $(CC_FLAGS) -o kalman_filter_loop.dll kalman_filter_loop.o \
# kalmanlib.a $(LD_LIBS)
doc: main.web $(hwebsource) $(cwebsource)
cweave -bhp main.web
pdftex main
mv main.pdf ts.pdf
clear:
rm -f *.o
rm -f *.{pdf,dvi,log,scn,idx,toc}
# rm -f *.cpp
# rm -f *.h

80
mex/sources/kalman/matlab/Makefile_mex
View File

@ -1,80 +0,0 @@
# $Id: Makefile 532 2005-11-30 13:51:33Z kamenik $
# Copyright 2005, Ondra Kamenik
CC = gcc
CC_FLAGS = -Wall -I../../sylv/cc/
LDFLAGS = -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++
ifeq ($(DEBUG),yes)
CC_FLAGS := $(CC_FLAGS) -g
else
CC_FLAGS := $(CC_FLAGS) -O3
endif
ifeq ($(OS),Windows_NT)
CC_FLAGS := -mno-cygwin $(CC_FLAGS)
LDFLAGS := -mno-cygwin $(LDFLAGS)
ARCH := w32
MEX_SUFFIX = dll
else
LDFLAGS := $(LDFLAGS)
ARCH := linux
MEX_SUFFIX = mexglx
endif
matrix_interface := GeneralMatrix Vector SylvException
matobjs := $(patsubst %, ../../sylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../../sylv/cc/%.h, $(matrix_interface))
matcppsource := $(patsubst %, ../../sylv/cc/%.cpp, $(matrix_interface))
tscwebsource := $(wildcard ../cc/*.cweb)
tscppsource := $(patsubst %.cweb,%.cpp,$(tscwebsource))
tshwebsource := $(wildcard ../cc/*.hweb)
tshsource := $(patsubst %.hweb,%.h,$(tshwebsource))
tsobjects := $(patsubst %.cweb,%.o,$(tscwebsource))
cppsource := $(wildcard *.cpp)
mexobjects := $(patsubst %.cpp,%_.$(MEX_SUFFIX),$(cppsource))
all: $(mexobjects)
../cc/dummy.ch:
make -C ../cc dummy.ch
../cc/%.cpp: ../cc/%.cweb ../cc/dummy.ch
make -C ../cc $*.cpp
../cc/%.h: ../cc/%.hweb ../cc/dummy.ch
make -C ../cc $*.h
../cc/%.o: ../cc/%.cpp $(tshsource)
make -C ../cc $*.o
tslib.a: $(tshwebsource) $(tscwebsoure) $(tshsource) $(tscppsource) \
$(mathsource) $(matcppsource) \
$(tsobjects) $(matobjs)
ar cr tslib.a $(tsobjects) $(matobjs)
ranlib tslib.a
# to compile mex objects for Windows do:
# 1. install gnumex
# 2. create mexopts.bat via gnumex in this directory, specify MinGW compilation, and dll output
# 3. change -Ic:/MATLAB7/extern/include according your Matlab setup
# 4. pray it works
%_.$(MEX_SUFFIX): %.cpp $(tshwebsource) $(tscwebsoure) $(tshsource) $(tscppsource) \
$(mathsource) $(matcppsource) \
tslib.a
ifeq ($(OS),Windows_NT)
mex.bat -I../../sylv/cc/ -I../cc COMPFLAGS\#"-c -DMATLAB_MEX_FILE" OPTIMFLAGS\#"-O3 -fexceptions -Ic:/MATLAB7/extern/include" GM_ADD_LIBS\#"$(LDFLAGS)" $*.cpp tslib.a
else
mex -cxx -I../../sylv/cc/ -I../cc $*.cpp CFLAGS="$(CC_FLAGS) -fexceptions" tslib.a $(LDFLAGS)
endif
mv $*.$(MEX_SUFFIX) $*_.$(MEX_SUFFIX)
clear:
rm -f tslib.a
rm -f *.mexglx
rm -f *.dll
make -C ../testing clear
make -C ../cc clear

153
mex/sources/kalman/matlab/diffuse_kalman_filter.cpp
View File

@ -1,153 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [LIK, lik] = kalman_filter_dll(T,R,Q,H,P,Y,start,Z/mf[, kalman_tol,riccati_tol])
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix with stationary variables
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z or mf: [double] Z: pp*mm matrix mapping state to pp observations
% Alternative parameters
% mf [integer] pp*1 vector of indices - alternative to Z matrix.
% Additional optional parameters
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 8 || nrhs > 10)
mexErrMsgTxt("Must have 8, 9, or 10 input parameters.\n");
if (nlhs < 1 || nlhs > 3)
mexErrMsgTxt("Must have 1, 2, or 3 output parameters.\n");
int start = 1; // default start of likelihood calculation
try
{
// make input matrices
GeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix R(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix Q(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix H(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix Pinf(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
GeneralMatrix P(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]));
GeneralMatrix Y(mxGetPr(prhs[6]), mxGetM(prhs[6]), mxGetN(prhs[6]));
if (nrhs>6) start = (int)mxGetScalar(prhs[7]);
GeneralMatrix Z(mxGetPr(prhs[8]), mxGetM(prhs[8]), mxGetN(prhs[8]));
int nper = mxGetN(prhs[5]); // no of periods
GeneralMatrix a( mxGetN(prhs[0]), 1);// initiate inital state to 0s
a.zeros();
#ifdef DEBUG
mexPrintf("kalman_filter: periods = %d ", nper);
#endif
// make storage for output
double loglik;
int per;
int d;
// output for full log-likelihood array
std::vector<double>* vll=new std::vector<double> (nper);
// create state init
StateInit* init = NULL;
init = new StateInit(P, Pinf, a.getData());
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
// developement of the output.
#ifdef DEBUG
mexPrintf("kalman_filter: running and filling outputs.\n");
#endif
loglik = kt.filter(per, d, (start-1), vll);
// destroy init
delete init;
// create output and upload output data
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
plhs[1] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[1]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
if (nlhs >= 3)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = per;
}
if (nlhs == 4)
{
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[3]))) = d;
}
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
}
};

89
mex/sources/kalman/matlab/disclyap_fast_dll.cpp
View File

@ -1,89 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/****************************************************************
% entry Matlab DLL for function X=disclyap_fast(G,V,ch)
%
% which solve the discrete Lyapunov Equation
% X=G*X*G'+V
% Using the Doubling Algorithm
%
% If ch is defined then the code will check if the resulting X
% is positive definite and generate an error message if it is not
%
****************************************************************/
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "SylvException.h"
#include "mex.h"
#include "disclyap_fast.h"
//void disclyap_fast(GeneralMatrix &G,GeneralMatrix & V, double tol= 1e-16, int ch=0);
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 2 || nrhs > 4)
mexErrMsgTxt("Must have 2, 3 or 4 input parameters.\n");
if (nlhs != 1 )
mexErrMsgTxt("Must have 1 output parameters.\n");
int cholCheck = 0;
double LyapTol=1e-06;
try
{
// make input matrices
int s = mxGetM(prhs[0]);
GeneralMatrix G(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix V(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
// create output
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]),mxGetN(prhs[0]), mxREAL);
GeneralMatrix X(mxGetPr(plhs[0]),mxGetM(plhs[0]),mxGetN(plhs[0]));
if (nrhs > 2)
LyapTol = (double)mxGetScalar(prhs[2]);
if (nrhs > 3)
cholCheck = (int)mxGetScalar(prhs[3]);
#ifdef TIMING_LOOP
for (int tt=0;tt<1000;++tt)
{
#endif
disclyap_fast(G, V, X, LyapTol, cholCheck);
#ifdef TIMING_LOOP
}
mexPrintf("disclyap_fast: finished 1000 loops");
#endif
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

30
mex/sources/kalman/matlab/gendata.m
View File

@ -1,30 +0,0 @@
% [y,epsilon,alpha,eta] = gendata(T, ssf, a0)
%
% generates random data of the length T for the given state space form
% and initial state
% $Id: gendata.m 532 2005-11-30 13:51:33Z kamenik $
% Copyright 2005, Ondra Kamenik
function [y,epsilon,alpha,eta] = gendata(T, ssf, a0)
m = size(ssf.T, 1);
p = size(ssf.Z, 1);
r = size(ssf.R, 2);
cholH = chol(ssf.H);
cholQ = chol(ssf.Q);
epsilon = cholH*randn(p,T);
eta = cholQ*randn(r, T);
y = zeros(p, T);
alpha = zeros(m,T);
alpha(:,1) = a0;
for t = 1:T
y(:,t) = ssf.Z*alpha(:,t) + epsilon(:,t);
if t ~= T
alpha(:,t+1) = ssf.T*alpha(:,t) + ssf.R*eta(:,t);
end
end

122
mex/sources/kalman/matlab/gmvm.cpp
View File

@ -1,122 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtmvm(QT,a)
%
% 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
% 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
%
% OUTPUTS
% Tinverse [double] mm*mm transition matrix of the state equation.
**********************************************************/
#include "qt.h"
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 2 )
mexErrMsgTxt("Must have min 2 input parameters.\n");
if (nlhs < 1 )
mexErrMsgTxt("Must have min 1 output parameters.\n")
;
//int start = 1; // default start of likelihood calculation
// test for univariate case
try
{
// make input matrices
int n=mxGetM(prhs[0]);
ConstGeneralMatrix QT(mxGetPr(prhs[0]), n, mxGetN(prhs[0]));
// ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
ConstVector a (mxGetPr(prhs[1]), n);
// create output and upload output data
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
// double * mxinv= mxGetPr(plhs[0]);
// GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
Vector Ta (mxGetPr(plhs[0]), n);
// Tinv.unit();
// Ta.zeros();
// make storage for output
#ifdef TIMING_LOOP
int loops=1;//000;
if (nrhs >2 )
loops = (int)mxGetScalar(prhs[2]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
QT.multVec(0.0, Ta, 1.0, a);
#ifdef DEBUG
Ta.print();
#endif
#ifdef TIMING_LOOP
}
mexPrintf("gmvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
/* if (nlhs >= 1)
{
plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
double * mxinv= mxGetPr(plhs[0]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxinv[j]=(*vll)[j];
}
*/
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

37
mex/sources/kalman/matlab/kalman_filter.m
View File

@ -1,37 +0,0 @@
%
% SYNOPSIS
%
% [loglik,per,d] = kalman_filter(Z,H,T,R,Q,Y,a,P)
% [loglik,per,d] = kalman_filter(Z,H,T,R,Q,Y,a,P,flag)
% [loglik,per,d] = kalman_filter(Z,H,T,R,Q,Y,a,Pstar,Pinf)
% [loglik,per,d] = kalman_filter(Z,H,T,R,Q,Y,a,Pstar,Pinf,flag)
%
% SEMANTICS
%
% The first two commands run a Kalman filter for non-diffuse
% initial conditions, the other two for diffuse initial conditions.
%
% Input:
% Z,H,T,R,Q gives a state space form
% Y observed data (columns correspond to periods)
% a mean of initial state
% P covariance of initial non-diffuse state
% Pstar finite part of covariance of initial diffuse state
% Pinf infinite part of covariance of initial diffuse state
% flag string starting with 'u', or 'U' runs a univariate
% form of the filter; if omitted, a multivariate version
% is run by default
%
% Output:
% loglik data log likelihood
% per number of succesfully filtered periods; if no error
% then per equals to the number of columns of Y
% d number of initial periods for which the state is
% still diffuse (d is always 0 for non-diffuse case)
%
% Copyright 2005, Ondra Kamenik
%
function [loglik, per, d] = kalman_filter(varargin)
[loglik, per, d] = kalman_filter_(varargin{:});

167
mex/sources/kalman/matlab/kalman_filter_Z.cpp
View File

@ -1,167 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [LIK, lik] = kalman_filter_dll(T,R,Q,H,P,Y,start,Z/mf[, kalman_tol,riccati_tol])
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix with stationary variables
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z or mf: [double] Z: pp*mm matrix mapping state to pp observations
% Alternative parameters
% mf [integer] pp*1 vector of indices - alternative to Z matrix.
% Additional optional parameters
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 8 || nrhs > 10)
mexErrMsgTxt("Must have 8, 9, or 10 input parameters.\n");
if (nlhs < 1 || nlhs > 3)
mexErrMsgTxt("Must have 1, 2, or 3 output parameters.\n");
// test for univariate case
bool uni = false;
const mxArray* const last = prhs[nrhs-1];
if (mxIsChar(last)
&& ((*mxGetChars(last)) == 'u' || (*mxGetChars(last)) == 'U'))
uni = true;
// test for diffuse case
bool diffuse = false;
if ((mxIsChar(last) && nrhs == 10) ||
(!mxIsChar(last) && nrhs == 9))
diffuse = true;
int start = 1; // default start of likelihood calculation
try
{
// make input matrices
GeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix R(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix Q(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix H(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix P(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
GeneralMatrix Y(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]));
if (nrhs>6) start = (int)mxGetScalar(prhs[6]);
GeneralMatrix Z(mxGetPr(prhs[7]), mxGetM(prhs[7]), mxGetN(prhs[7]));
int nper = mxGetN(prhs[5]); // no of periods
GeneralMatrix a( mxGetN(prhs[0]), 1);// initiate inital state to 0s
a.zeros();
#ifdef DEBUG
mexPrintf("kalman_filter: periods = %d ", nper);
#endif
// make storage for output
double loglik;
int per;
int d;
// output for full log-likelihood array
std::vector<double>* vll=new std::vector<double> (nper);
// create state init
StateInit* init = NULL;
init = new StateInit(P, a.getData());
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
// developement of the output.
#ifdef DEBUG
mexPrintf("kalman_filter: running and filling outputs.\n");
#endif
loglik = kt.filter(per, d, (start-1), vll);
// destroy init
delete init;
// create output and upload output data
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
plhs[1] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[1]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
if (nlhs >= 3)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = per;
}
if (nlhs == 4)
{
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[3]))) = d;
}
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
}
};

209
mex/sources/kalman/matlab/kalman_filters.cpp
View File

@ -1,209 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [loglik per d vll] = kalman_filter_dll(T,R,Q,H,Y,start,a, Z, P. [Pinf | u/U flag]
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z [double] pp*mm matrix mapping state to pp observations
% a [vector] mm vector of initial state, usually of 0s
% P [double] mm*mm variance-covariance matrix with stationary variables
% Pinf [optional] [double] mm*mm variance-covariance matrix with stationary variables
% u/U [optional] [char] u/U univariate
% OUTPUTS
% loglik [double] scalar, total likelihood
% per [int] number of succesfully filtered periods; if no error
% [int] then per equals to the number of columns of Y
% d number of initial periods for which the state is
% still diffuse (d is always 0 for non-diffuse case)
% vll [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "vll" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 9 || nrhs > 11)
mexErrMsgTxt("Must have 9, 10 or 11 input parameters.\n");
if (nlhs < 1 || nlhs > 4)
mexErrMsgTxt("Must have 1, 2, 3 or 4 output parameters.\n");
//int start = 1; // default start of likelihood calculation
// test for univariate case
bool uni = false;
double riccatiTol=0.000001;
const mxArray* const last = prhs[nrhs-1];
if (mxIsChar(last)
&& ((*mxGetChars(last)) == 'u' || (*mxGetChars(last)) == 'U'))
uni = true;
// test for diffuse case
bool diffuse = false;
if ((mxIsChar(last) && nrhs == 11) ||
(!mxIsChar(last) && nrhs == 10))
diffuse = true;
try {
// make input matrices
GeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix R(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix Q(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix H(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix Y(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
int start = (int)mxGetScalar(prhs[5]);
GeneralMatrix Z(mxGetPr(prhs[6]), mxGetM(prhs[6]), mxGetN(prhs[6]));
GeneralMatrix a(mxGetPr(prhs[7]), mxGetM(prhs[7]), mxGetN(prhs[7]));
GeneralMatrix P(mxGetPr(prhs[8]), mxGetM(prhs[8]), mxGetN(prhs[8]));
int nper=Y.numCols();
#ifdef DEBUG
mexPrintf("kalman_filter: periods=%d start=%d, a.length=%d, uni=%d diffuse=%d\n", nper, start,a.numRows(), uni, diffuse);
#endif
// make storage for output
double loglik;
int per;
int d;
// create state init
StateInit* init = NULL;
std::vector<double>* vll=new std::vector<double> (nper);
if (diffuse||uni)
{
if (diffuse)
{
GeneralMatrix Pinf(mxGetPr(prhs[9]), mxGetM(prhs[9]), mxGetN(prhs[9]));
init = new StateInit(P, Pinf, a.getData());
}
else
{
init = new StateInit(P, a.getData());
}
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
if (uni)
{
KalmanUniTask kut(kt);
loglik = kut.filter(per, d, (start-1), vll);
per = per / Y.numRows();
d = d / Y.numRows();
}
else
{
#ifdef TIMING_LOOP
for (int tt=0;tt<1000;++tt)
{
#endif
loglik = kt.filter(per, d, (start-1), vll);
#ifdef TIMING_LOOP
}
mexPrintf("kalman_filter: finished 1000 loops");
#endif
}
}
else // basic Kalman
{
init = new StateInit(P, a.getData());
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init, riccatiTol);
#ifdef TIMING_LOOP
for (int tt=0;tt<1000;++tt)
{
#endif
loglik = bkt.filter( per, d, (start-1), vll);
#ifdef DEBUG
mexPrintf("Basickalman_filter: loglik=%f \n", loglik);
#endif
#ifdef TIMING_LOOP
}
mexPrintf("Basickalman_filter: finished 1000 loops");
#endif
}
// destroy init
delete init;
// create output and upload output data
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2) {
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[1]))) = per;
}
if (nlhs >= 3) {
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = d;
}
if (nlhs >= 4)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
plhs[3] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[3]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
delete vll;
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

61
mex/sources/kalman/matlab/kalman_filters.m
View File

@ -1,61 +0,0 @@
% function [LIK per d lik] = kalman_filters(T,R,Q,H,Y,start,Z,a,P,[Pinf |u/U flag])
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z [double] pp*mm matrix mapping state to pp observations
% a [vector] mm vector of mean initial state, usually of 0s
% P [double] mm*mm variance-covariance matrix with stationary variables
% Pinf [optional] [double] mm*mm variance-covariance matrix with stationary variables
% [u/U]flag [optional] [char] u/U univariate falg
%
% SYNOPSIS
%
% [LIK,per,d,lik] = kalman_filters(T,R,Q,H,Y,start,a, Z,P)
% [LIK,per,d,lik] = kalman_filters(T,R,Q,H,Y,start,a, Z,P,flag)
% [LIK,per,d,lik] = kalman_filters(T,R,Q,H,Y,start,a, Z,Pstar,Pinf)
% [LIK,per,d,lik] = kalman_filters(T,R,Q,H,Y,start,a, Z, Pstar, Pinf, flag)
%
% SEMANTICS
%
% The first two commands run a Kalman filter for non-diffuse initial conditions,
% univariate or multivariate, the other two for diffuse initial conditions.
%
%
% Output:
% LIK data log likelihood
% per number of succesfully filtered periods; if no error
% then per equals to the number of columns of Y
% d number of initial periods for which the state is
% still diffuse (d is always 0 for non-diffuse case)
%
% Copyright 2005, Ondra Kamenik
%
%function [LIK per d lik] = kalman_filters(varargin)
function [LIK per d lik] = kalman_filters(T,R,Q,H,Y,start,Z,a,P,varargin)
if isempty(H)
H=zeros(size(Y,1), size(Y,1))
elseif H==0
H=zeros(size(Y,1), size(Y,1))
end
if isempty(a)
a=zeros(size(T,1),1)
elseif a==0
a=zeros(size(T,1),1)
end
if size(Z,1)== 1 && size(Z,2)==size(Y,1) && size(Y,1)> 1
ZM=zeros(size(Y,1), size(T,2))
for i = 1:size(Y,1)
ZM(i,Z(i))=1
end
else
ZM=Z;
end
% run basic multivariate kalman filter
[LIK per d lik] = kalman_filter_dll(T,R,Q,H,Y,start,ZM,a, P);

324
mex/sources/kalman/matlab/kalman_filters_testx.cpp
View File

@ -1,324 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [loglik per d vll] = kalman_filter_dll(T,R,Q,H,Y,start,a, Z, P. [Pinf | u/U flag]
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z [double] pp*mm matrix mapping state to pp observations
% a [vector] mm vector of initial state, usually of 0s
% P [double] mm*mm variance-covariance matrix with stationary variables
% Pinf [optional] [double] mm*mm variance-covariance matrix with stationary variables
% u/U [optional] [char] u/U univariate
% OUTPUTS
% loglik [double] scalar, total likelihood
% per [int] number of succesfully filtered periods; if no error
% [int] then per equals to the number of columns of Y
% d number of initial periods for which the state is
% still diffuse (d is always 0 for non-diffuse case)
% vll [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "vll" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include <iostream>
using namespace std;
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
/*************************************
* This main() is for testing kalman DLL entry point by linking to
* the kalman library statically and passing its hard-coded data:
* parameters, covar,
***************************************/
int main(int argc, char* argv[])
//main (int nrhs, mxArray* plhs[],
{
int nrhs=9;
int nlhs=4;
if (nrhs < 9 || nrhs > 11)
mexErrMsgTxt("Must have 9, 10 or 11 input parameters.\n");
if (nlhs < 1 || nlhs > 4)
mexErrMsgTxt("Must have 1, 2, 3 or 4 output parameters.\n");
//int start = 1; // default start of likelihood calculation
// test for univariate case
bool uni = false;
// const mxArray* const last = prhs[nrhs-1];
// if (mxIsChar(last)
// && ((*mxGetChars(last)) == 'u' || (*mxGetChars(last)) == 'U'))
// uni = true;
// test for diffuse case
bool diffuse = false;
// if ((mxIsChar(last) && nrhs == 11) ||
// (!mxIsChar(last) && nrhs == 10))
// diffuse = true;
double Tmat[]={// need to pass transposed matrices!!??
0, 0, 0, 0, 0, 0, 0, 0,
-0.0013, 0.5000, 0.0000, -0.0000, 0.0188, -0.0013, 0.1182, -0.0017,
0.2158, 0.0000, 0.9502, -0.0000, 0.0127, 0.2158, 0.0438, -0.0088,
0.0273, -0.0000, -0.0000, 0.8522, -0.1260, -0.8249, -0.4720, 0.0356,
-0.0716, -0.0000, 0.0000, 0.0000, 0.5491, -0.0716, -0.9573, -0.0935,
-0.0000, -0.0000, 0.0000, -0.0000, -0.0000, -0.0000, 0.0000, -0.0000,
0, 0, 0, 0, 0, 0, 0, 0,
0.6464, 0.0000, -0.0000, -0.0000, 0.0573, 0.6464, 0.2126, 0.8441
// 0 ,-0.001294119891461, 0.21578807493606 ,0.027263201686985, -0.071633450625617, -0, 0, 0.646379181371765,
// 0, 0.5, 0, -0, -0, -0, 0, 0,
// 0, 0, 0.9502, -0, 0, 0, 0, -0,
// 0, -0, -0, 0.8522, 0, -0, 0, -0,
// 0, 0.018758765757513, 0.012692095232426, -0.126035674083997, 0.549074256326045, -0, 0, 0.05730910985981,
// 0, -0.001294119891461, 0.21578807493606, -0.824936798313015, -0.071633450625617, -0, 0, 0.646379181371766,
// 0, 0.118192240459753, 0.04380066554165, -0.471963836695487, -0.957255289691476, 0, 0, 0.212592467520726,
// 0, -0.00168993250228, -0.008835241183444, 0.035601779209991, -0.093542875943306, -0, 0, 0.844077271823789
};
double Rmat[]={// need to pass transposed matrices!!??
0.2271, 0, 1.0000, 0, 0.0134, 0.2271, 0.0461, -0.0093,
0.0320, 0, 0, 1.0000, -0.1479, -0.9680, -0.5538, 0.0418,
-0.0026, 1.0000, 0, 0, 0.0375, -0.0026, 0.2364, -0.0034,
-0.0895, 0, 0, 0, 0.6863, -0.0895, -1.1966, -0.1169
// 0.2271, 0.0320, -0.0026, -0.0895,
// 0, 0, 1.0000, 0,
// 1.0000, 0, 0, 0,
// 0, 1.0000, 0, 0,
// 0.0134, -0.1479, 0.0375, 0.6863,
// 0.2271, -0.9680, -0.0026, -0.0895,
// 0.0461, -0.5538, 0.2364, -1.1966,
// -0.0093, 0.0418, -0.0034, -0.1169
};
double Qmat[]={
0.0931, 0, 0, 0,
0, 0.1849, 0, 0,
0, 0, 0.0931, 0,
0, 0, 0, 0.0100
};
double Zmat[]={ // need to pass transposed matrices!!??
0, 0, 1, 0,
0, 0, 0, 0,
0, 0, 0, 0,
0, 0, 0, 0,
0, 0, 0 , 1,
0, 0, 0, 0,
1, 0, 0, 0,
0, 1, 0, 0
// 0, 0, 0, 0, 0, 0, 1, 0,
// 0, 0, 0, 0, 0, 0, 0, 1,
// 1, 0, 0, 0, 0, 0, 0, 0,
// 0, 0, 0, 0, 1, 0, 0, 0
};
double Ymat[]={
-0.4073, 0.2674, 0.2896, 0.0669, 0.1166, -0.1699, -0.2518, -0.0562, -0.3269,-0.0703,-0.1046, -0.4888 ,-0.3524, -0.2485 ,-0.587, -0.4546, -0.397, -0.2353, -0.0352 -0.2171, -0.3754, -0.4322, -0.4572, -0.4903, -0.4518, -0.6435, -0.6304 ,-0.4148, -0.2892, -0.4318, -0.601, -0.4148, -0.4315, -0.3531, -0.8053, -0.468, -0.4263,
3.1739, 3.738 , 3.8285, 3.3342, 3.7447, 3.783, 3.1039, 2.8413, 3.0338, 0.3669, 0.0847 ,0.0104, 0.2115, -0.6649, -0.9625, -0.733, -0.8664, -1.4441, -1.0179, -1.2729 ,-1.9539, -1.4427, -2.0371, -1.9764, -2.5654, -2.857, -2.5842, -3.0427, -2.8312, -2.332, -2.2768, -2.1816, -2.1043, -1.8969, -2.2388, -2.1679, -2.1172,
3.2174, 3.1903, 3.3396, 3.1358, 2.8625, 3.3546, 2.4609, 1.9534, 0.9962, -0.7904,-1.1672, -1.2586, -1.3593, -1.3443 ,-0.9413, -0.6023, -0.4516, -0.5129, -0.8741, -1.0784, -1.4091, -1.3627, -1.5731, -1.6037 -1.8814, -2.1482 ,-1.3597, -1.1855, -1.1122, -0.8424, -0.9747, -1.1385, -1.4548, -1.4284, -1.4633, -1.0621, -0.7871,
0.8635, 0.9058, 0.7656, 0.7936, 0.8631, 0.9074, 0.9547, 1.2045, 1.085, 0.9178, 0.5242, 0.3178 ,0.1472, 0.0227, -0.0799, -0.0611, -0.014, 0.1132, 0.1774, 0.0782, 0.0436, -0.1596, -0.2691, -0.2895, -0.3791, -0.402, -0.4166 ,-0.4037, -0.3636, -0.4075, -0.4311, -0.447, -0.5111, -0.6274, -0.7261, -0.6974, -0.5012
};
try {
// make input matrices
GeneralMatrix T(Tmat, 8, 8);
GeneralMatrix R(Rmat, 8, 4);
GeneralMatrix Q(Qmat, 4, 4);
GeneralMatrix H(4, 4);
H.zeros();
/*********use simlated data for time being *********/
GeneralMatrix Y( 4, 109);
Y.zeros();
for (int i=0;i<4;++i)
{
for (int j=0;j<109;++j)
{
Y.get(i,j)= ((double) ( rand() % 10 -5.0))/2.0;
#ifdef DEBUG
mexPrintf("Y [%d %d] =%f, \n", i, j,Y.get(i,j));
#endif
}
}
/***********
GeneralMatrix Y(Ymat, 4, 109);
for (int i=0;i<4;++i)
{
for (int j=0;j<109;++j)
{
#ifdef DEBUG
mexPrintf("Y [%d %d] =%f, \n", i, j,Y.get(i,j));
#endif
}
}
***********/
double riccatiTol=0.000000000001;
int start = 1;
GeneralMatrix Z(Zmat, 4, 8);
GeneralMatrix a(8, 1);
a.zeros();
GeneralMatrix P( 8, 8);
P.zeros();
for (int i=0;i<8;++i)
P.get(i,i)=10.0;
int nper=Y.numCols();
#ifdef DEBUG
mexPrintf("kalman_filter: periods=%d start=%d, a.length=%d, uni=%d diffuse=%d\n", nper, start,a.numRows(), uni, diffuse);
#endif
// make storage for output
double loglik=-1.1111;
int per;
int d;
// create state init
StateInit* init = NULL;
std::vector<double>* vll=new std::vector<double> (nper);
bool basicKF=true;//false;
if (diffuse||uni||basicKF==false)
{
if (diffuse)
{
GeneralMatrix Pinf(P.numRows(),P.numCols());
Pinf.zeros();
init = new StateInit(P, Pinf, a.getData());
}
else
{
init = new StateInit(P, a.getData());
}
// fork, create objects and do filtering
#ifdef TIMING_LOOP
for (int tt=0;tt<10000;++tt)
{
#endif
KalmanTask kt(Y, Z, H, T, R, Q, *init);
if (uni)
{
KalmanUniTask kut(kt);
loglik = kut.filter(per, d, (start-1), vll);
per = per / Y.numRows();
d = d / Y.numRows();
}
else
{
loglik = kt.filter(per, d, (start-1), vll);
}
#ifdef TIMING_LOOP
// mexPrintf("kalman_filter: finished %d loops", tt);
}
mexPrintf("kalman_filter: finished 10,000 loops");
#endif
}
else // basic Kalman
{
init = new StateInit(P, a.getData());
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init, riccatiTol);
#ifdef TIMING_LOOP
for (int tt=0;tt<10000;++tt)
{
#endif
loglik = bkt.filter( per, d, (start-1), vll);
#ifdef DEBUG
// mexPrintf("Basickalman_filter: loglik=%f \n", loglik);
// cout << "loglik " << loglik << "\n";
#endif
#ifdef TIMING_LOOP
}
mexPrintf("Basickalman_filter: finished 10,000 loops");
#endif
}
// destroy init
delete init;
mexPrintf("logLik = %f \n", loglik);
delete vll;
// create output and upload output data
/************
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2) {
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[1]))) = per;
}
if (nlhs >= 3) {
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = d;
}
if (nlhs >= 4)
{
// output full log-likelihood array
// Set the output pointer to the array of log likelihood.
plhs[3] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[3]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
******************************/
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
// } // mexFunction
}; // main extern 'C'

141
mex/sources/kalman/matlab/kalman_smoother.cpp
View File

@ -1,141 +0,0 @@
/* $Id: kalman_smoother.cpp 532 2005-11-30 13:51:33Z kamenik $
*
*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nhls, mxArray* plhs[],
int nhrs, const mxArray* prhs[])
{
if (nhrs < 8 || nhrs > 10)
mexErrMsgTxt("Must have 8, 9, or 10 input parameters.\n");
if (nhls < 1 || nhls > 6)
mexErrMsgTxt("Must have 1, 2,.. or 6 output parameters.\n");
// test for univariate case
bool uni = false;
const mxArray* const last = prhs[nhrs-1];
if (mxIsChar(last)
&& ((*mxGetChars(last)) == 'u' || (*mxGetChars(last)) == 'U'))
uni = true;
// test for diffuse case
bool diffuse = false;
if ((mxIsChar(last) && nhrs == 10) ||
(!mxIsChar(last) && nhrs == 9))
diffuse = true;
try {
// make input matrices
GeneralMatrix Z(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix H(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix T(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix R(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix Q(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
GeneralMatrix Y(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]));
GeneralMatrix a(mxGetPr(prhs[6]), mxGetM(prhs[6]), mxGetN(prhs[6]));
GeneralMatrix P(mxGetPr(prhs[7]), mxGetM(prhs[7]), mxGetN(prhs[7]));
// make storage for output
double loglik;
int per;
int d;
SmootherResults sres(Y.numCols());
// create state init
StateInit* init = NULL;
if (diffuse) {
GeneralMatrix Pinf(mxGetPr(prhs[8]), mxGetM(prhs[8]), mxGetN(prhs[8]));
init = new StateInit(P, Pinf, a.getData());
} else {
init = new StateInit(P, a.getData());
}
// fork, create objects and do filtering and smoothing
KalmanTask kt(Y, Z, H, T, R, Q, *init);
if (uni) {
KalmanUniTask kut(kt);
SmootherResults sres_uni(Y.numRows()*Y.numCols());
loglik = kut.filter_and_smooth(sres_uni, per, d);
per = per / Y.numRows();
d = d / Y.numRows();
sres.import(sres_uni, Y.numRows());
} else {
loglik = kt.filter_and_smooth(sres, per, d);
// loglik = kt.filter(per, d);
}
// destroy init
delete init;
// create output and upload output data
if (nhls >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nhls >= 2) {
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[1]))) = per;
}
if (nhls >= 3) {
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = d;
}
if (nhls >= 4) {
plhs[3] = mxCreateNumericMatrix(T.numRows(), Y.numCols(), mxDOUBLE_CLASS, mxREAL);
if (per == Y.numCols()) {
GeneralMatrix tmp(mxGetPr(plhs[3]), T.numRows(), Y.numCols());
sres.exportAlpha(tmp);
}
}
if (nhls >= 5) {
plhs[4] = mxCreateNumericMatrix(R.numCols(), Y.numCols(), mxDOUBLE_CLASS, mxREAL);
if (per == Y.numCols()) {
GeneralMatrix tmp(mxGetPr(plhs[4]), R.numCols(), Y.numCols());
sres.exportEta(tmp);
}
}
if (nhls >= 6) {
int dims[3]; dims[0] = T.numRows();
dims[1] = T.numRows(); dims[2] = Y.numCols();
plhs[5] = mxCreateNumericMatrix(3, dims[0], mxDOUBLE_CLASS, mxREAL);
if (per == Y.numCols()) {
GeneralMatrix tmp(mxGetPr(plhs[5]), T.numRows(),
T.numRows()*Y.numCols());
sres.exportV(tmp);
}
}
} catch (const TSException& e) {
mexErrMsgTxt(e.getMessage());
} catch (SylvException& e) {
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
}
};

41
mex/sources/kalman/matlab/kalman_smoother.m
View File

@ -1,41 +0,0 @@
%
% SYNOPSIS
%
% [loglik,per,d,alpha,eta,V] = kalman_smoother(Z,H,T,R,Q,Y,a,P)
% [loglik,per,d,alpha,eta,V] = kalman_smoother(Z,H,T,R,Q,Y,a,P,flag)
% [loglik,per,d,alpha,eta,V] = kalman_smoother(Z,H,T,R,Q,Y,a,Pstar,Pinf)
% [loglik,per,d,alpha,eta,V] = kalman_smoother(Z,H,T,R,Q,Y,a,Pstar,Pinf,flag)
%
% SEMANTICS
%
% The first two commands run a Kalman filter and smoother for non-diffuse initial
% conditions, the other two for diffuse initial conditions.
%
% Input:
% Z,H,T,R,Q gives a state space form
% Y observed data (columns correspond to periods)
% a mean of initial state
% P covariance of initial non-diffuse state
% Pstar finite part of covariance of initial diffuse state
% Pinf infinite part of covariance of initial diffuse state
% flag string starting with 'u', or 'U' runs a univariate
% form of the filter; if omitted, a multivariate version
% is run by default
%
% Output:
% loglik data log likelihood
% per number of succesfully filtered periods; if no error
% then per equals to the number of columns of Y
% d number of initial periods for which the state is
% still diffuse (d is always 0 for non-diffuse case)
% alpha matrix of smoothed states; columns are periods
% eta matrix of smoothed shocks; columns are periods
% V 3D array of smoothed state variances; V(:,:,t) is
% smoothed state variance-covariance at time t
%
% Copyright 2005, Ondra Kamenik
%
function [loglik,per,d,alpha,eta,V] = kalman_smoother(varargin)
[loglik,per,d,alpha,eta,V] = kalman_smoother_(varargin{:});

23
mex/sources/kalman/matlab/kalmandll_test.m
View File

@ -1,23 +0,0 @@
function [LIKDLL loglik]=kalmandll_test(T,mf,R,Q,H,Pstar,Pinf,data,start)
if isempty(H)
H=zeros(size(data,1), size(data,1))
elseif H==0
H=zeros(size(data,1), size(data,1))
end
Z=zeros(size(data,1), size(T,2))
for i = 1:size(data,1)
Z(i,mf(i))=1
end
%LIKDLL= kalman_filter_dll8(T,Z,R,Q,H,Pstar,data,start)
%[loglik,per,d] = kalman_filter_dll(Z,H,T,R,Q,data,a,Pstar)
[LIK2 lik2] = kalman_filter(T,R,Q,H,Pstar,data,start,mf,options_.kalman_tol,options_.riccati_tol)
if isempty(Pinf)
Pinf=zeros(size(T));
elseif Pinf==0
Pinf=zeros(size(T));
end
% test DiffuseLikelihoodH1
%[loglikd1,per,d] = kalman_filter_dll(Z,H,T,R,Q,data,a,Pstar,Pinf)
[LIKdlikd ] = diffuse_kalman_filter(T,R,Q,H,Pinf,Pstar,data,start,Z,options_.kalman_tol,options_.riccati_tol)
%loglikd2 = dynare_filter(Z,H,T,R,Q,data,Pstar,Pinf)

2
mex/sources/kalman/matlab/kalmanlib.def
View File

@ -1,2 +0,0 @@
EXPORTS
mexFunction

58
mex/sources/kalman/matlab/mexopts.bat
View File

@ -1,58 +0,0 @@
@echo off
rem C:\ondra\tmp\dynare++\extern\matlab\mexopts.bat
rem Generated by gnumex.m script in c:\fs\gnumex-1.11
rem gnumex version: 1.11
rem Compile and link options used for building MEX etc files with
rem the Mingw/Cygwin tools. Options here are:
rem Mingw linking
rem Mex (*.dll) creation
rem Safe linking to temporary libraries
rem Language: C / C++
rem Compiling for pentium and above
rem Matlab version 7
rem
set MATLAB=C:\MATLAB7
set GM_PERLPATH=C:\MATLAB7\sys\perl\win32\bin\perl.exe
set GM_UTIL_PATH=c:\fs\gnumex-1.11
set PATH=c:\fs\mingw\bin;%PATH%
rem
rem Added libraries for linking
set GM_ADD_LIBS=
rem
rem Type of file to compile (mex or engine)
set GM_MEXTYPE=mex
rem
rem Language for compilation
set GM_MEXLANG=c
rem
rem def files to be converted to libs
set GM_DEFS2LINK=libmx.def;libmex.def;libmat.def;
rem
rem dlltool command line
set GM_DLLTOOL=dlltool
rem
rem compiler options; add compiler flags to compflags as desired
set NAME_OBJECT=-o
set COMPILER=gcc
set COMPFLAGS=-c -DMATLAB_MEX_FILE
set OPTIMFLAGS=-O3 -malign-double -fno-exceptions -mcpu=pentium
set DEBUGFLAGS=-g
set CPPCOMPFLAGS=%COMPFLAGS% -x c++
set CPPOPTIMFLAGS=%OPTIMFLAGS%
set CPPDEBUGFLAGS=%DEBUGFLAGS%
rem
rem NB Library creation commands occur in linker scripts
rem
rem Linker parameters
set LINKER=%GM_PERLPATH% %GM_UTIL_PATH%\linkmex.pl
set LINKFLAGS=
set CPPLINKFLAGS=GM_ISCPP
set LINKOPTIMFLAGS=-s
set LINKDEBUGFLAGS=-g -Wl,--image-base,0x28000000\n
set LINK_FILE=
set LINK_LIB=
set NAME_OUTPUT=-o %OUTDIR%%MEX_NAME%.dll
rem
rem Resource compiler parameters
set RC_COMPILER=%GM_PERLPATH% %GM_UTIL_PATH%\rccompile.pl -o %OUTDIR%mexversion.res
set RC_LINKER=

108
mex/sources/kalman/matlab/minv.cpp
View File

@ -1,108 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides a matrix inversion
/******************************************************
% function [Tinv] = minv(T)
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% OUTPUTS
% Tinverse [double] mm*mm transition matrix of the state equation.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 1 )
mexErrMsgTxt("Must have min 1 input parameters.\n");
if (nlhs < 1 )
mexErrMsgTxt("Must have min 1 output parameters.\n")
;
//int start = 1; // default start of likelihood calculation
// test for univariate case
try
{
// make input matrices
ConstGeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
// create output and upload output data
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxREAL);
// double * mxinv= mxGetPr(plhs[0]);
GeneralMatrix Tinv(mxGetPr(plhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
// Tinv.unit();
// Tinv.zeros();
// make storage for output
#ifdef TIMING_LOOP
int loops=1000;
if (nrhs >1 )
loops = (int)mxGetScalar(prhs[1]);
for (int tt=0;tt<loops;++tt)
{
#endif
Tinv.unit();
T.multInvLeft(Tinv);
//Tinv.print();
#ifdef TIMING_LOOP
}
mexPrintf("minv: finished: %d loops\n",loops);
#endif
// create output and upload output data
/* if (nlhs >= 1)
{
plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
double * mxinv= mxGetPr(plhs[0]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxinv[j]=(*vll)[j];
}
*/
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

148
mex/sources/kalman/matlab/qtamvm.cpp
View File

@ -1,148 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtmvm(QT,a)
%
% 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
% 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% Tinverse [double] mm*mm transition matrix of the state equation.
**********************************************************/
#include "qt.h"
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 2 )
mexErrMsgTxt("Must have min 2 input parameters.\n");
if (nlhs < 1 )
mexErrMsgTxt("Must have min 1 output parameters.\n")
;
//int start = 1; // default start of likelihood calculation
// test for univariate case
try
{
// make input matrices
int n=mxGetM(prhs[0]);
double *T1, *Ld, *TV ;
double * QT=mxGetPr(prhs[0]);
// ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
double *a =(mxGetPr(prhs[1]));
// create output and upload output data
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
// double * mxinv= mxGetPr(plhs[0]);
// GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
double * Ta =mxGetPr(plhs[0]);
// Tinv.unit();
// Ta.zeros();
T1=(double *)mxCalloc(n*n, sizeof(double));
Ld=(double *)mxCalloc(n*n,sizeof(double));
TV=(double *)mxCalloc(n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (nrhs >2 )
loops = (int)mxGetScalar(prhs[2]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
// 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
// double *T1, *Ld, *dTa;//, dT1=-7.77;
// mexPrintf("start dT1 = %f\n", dT1);
// dT1 = qt2t_(QT.base() ,&n) ;
// mexPrintf("end dT1 = %f\n", dT1);
//T1=&dT1;
qt2t_(T1, QT ,&n) ;
// T1=*T1p;
// GeneralMatrix T1gm(T1,n,n);
// Ld = qt2ld_(QT.base(),&n);
qt2ld_(Ld , QT,&n);
// 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
// dTa = ldv_(Ld,a.base() ,&n);
//Vector Ta( n);
ldv_(Ta, Ld,a ,&n);
// Ta2 = tv_(T1,a.base(),&n);
tv_(TV, T1 ,a,&n);
// Ta.add(1.0,ConstVector(Ta2.base(), n));
// Ta.add(1.0,TV);
for (int j=0; j<n;++j)
Ta[j]+=TV[j];
#ifdef TIMING_LOOP
}
mexPrintf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
/* if (nlhs >= 1)
{
plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
double * mxinv= mxGetPr(plhs[0]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxinv[j]=(*vll)[j];
}
*/
mxFree(T1);
mxFree(Ld);
mxFree(TV);
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

155
mex/sources/kalman/matlab/qtmvm.cpp
View File

@ -1,155 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtmvm(QT,a)
%
% 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
% 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% Tinverse [double] mm*mm transition matrix of the state equation.
**********************************************************/
#include "qt.h"
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 2 )
mexErrMsgTxt("Must have min 2 input parameters.\n");
if (nlhs < 1 )
mexErrMsgTxt("Must have min 1 output parameters.\n")
;
//int start = 1; // default start of likelihood calculation
// test for univariate case
try
{
// make input matrices
int n=mxGetM(prhs[0]);
ConstGeneralMatrix QT(mxGetPr(prhs[0]), n, mxGetN(prhs[0]));
// ConstGeneralMatrix a (mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
Vector a (mxGetPr(prhs[1]), n);
// create output and upload output data
plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[1]),1, mxREAL);// mxGetM(prhs[1]), mxREAL);
// double * mxinv= mxGetPr(plhs[0]);
// GeneralMatrix Ta(mxGetPr(plhs[0]), mxGetM(prhs[1]), mxGetN(prhs[1]));
Vector Ta (mxGetPr(plhs[0]), n);
// Tinv.unit();
// Ta.zeros();
GeneralMatrix T1gm(n,n);
GeneralMatrix Ld(n,n);
Vector TV( n);
// make storage for output
#ifdef TIMING_LOOP
int loops=1;//000;
if (nrhs >2 )
loops = (int)mxGetScalar(prhs[2]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
// 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
// double *T1, *Ld, *dTa;//, dT1=-7.77;
// mexPrintf("start dT1 = %f\n", dT1);
// dT1 = qt2t_(QT.base() ,&n) ;
// mexPrintf("end dT1 = %f\n", dT1);
//T1=&dT1;
qt2t_(T1gm.base(), QT.base() ,&n) ;
// T1=*T1p;
// GeneralMatrix T1gm(T1,n,n);
#ifdef DEBUG
T1gm.print();
#endif
// Ld = qt2ld_(QT.base(),&n);
qt2ld_(Ld.base() , QT.base(),&n);
#ifdef DEBUG
Ld.print();
#endif
// 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
// dTa = ldv_(Ld,a.base() ,&n);
//Vector Ta( n);
ldv_(Ta.base(), Ld.base(),a.base() ,&n);
// Ta= (ConstVector(dTa, n));
// Vector Ta(dTa, n);
#ifdef DEBUG
Ta.print();
#endif
// Ta2 = tv_(T1,a.base(),&n);
tv_(TV.base(), T1gm.base() ,a.base(),&n);
// Ta.add(1.0,ConstVector(Ta2.base(), n));
Ta.add(1.0,TV);
#ifdef DEBUG
Ta.print();
#endif
#ifdef TIMING_LOOP
}
mexPrintf("QTmvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
/* if (nlhs >= 1)
{
plhs[0] = mxCreateNumericMatrix(mxGetM(prhs[0]), mxGetM(prhs[0]), mxINT32_CLASS, mxREAL);
double * mxinv= mxGetPr(plhs[0]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxinv[j]=(*vll)[j];
}
*/
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
} // mexFunction
}; // extern 'C'

158
mex/sources/kalman/matlab/univariate_diffuse_kalman_filter.cpp
View File

@ -1,158 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [LIK, lik] = kalman_filter_dll(T,R,Q,H,P,Y,start,Z/mf[, kalman_tol,riccati_tol])
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix with stationary variables
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z or mf: [double] Z: pp*mm matrix mapping state to pp observations
% Alternative parameters
% mf [integer] pp*1 vector of indices - alternative to Z matrix.
% Additional optional parameters
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 8 || nrhs > 10)
mexErrMsgTxt("Must have 8, 9, or 10 input parameters.\n");
if (nlhs < 1 || nlhs > 3)
mexErrMsgTxt("Must have 1, 2, or 3 output parameters.\n");
int start = 1; // default start of likelihood calculation
try
{
// make input matrices
GeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix R(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix Q(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix H(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix Pinf(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
GeneralMatrix P(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]));
GeneralMatrix Y(mxGetPr(prhs[6]), mxGetM(prhs[6]), mxGetN(prhs[6]));
if (nrhs>6) start = (int)mxGetScalar(prhs[7]);
GeneralMatrix Z(mxGetPr(prhs[8]), mxGetM(prhs[8]), mxGetN(prhs[8]));
int nper = mxGetN(prhs[5]); // no of periods
GeneralMatrix a( mxGetN(prhs[0]), 1);// initiate inital state to 0s
a.zeros();
#ifdef DEBUG
mexPrintf("kalman_filter: periods = %d ", nper);
#endif
// make storage for output
double loglik;
int per;
int d;
// output for full log-likelihood array
std::vector<double>* vll=new std::vector<double> (nper);
// create state init
StateInit* init = NULL;
init = new StateInit(P, Pinf, a.getData());
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
// developement of the output.
#ifdef DEBUG
mexPrintf("kalman_filter: running and filling outputs.\n");
#endif
KalmanUniTask kut(kt);
loglik = kut.filter(per, d, (start-1), vll);
per = per / Y.numRows();
d = d / Y.numRows();
// destroy init
delete init;
// create output and upload output data
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
plhs[1] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[1]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
if (nlhs >= 3)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = per;
}
if (nlhs == 4)
{
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[3]))) = d;
}
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
}
};

156
mex/sources/kalman/matlab/univariate_kalman_filter.cpp
View File

@ -1,156 +0,0 @@
// $Id: kalman_filter.cpp 532 2005-11-30 13:51:33Z kamenik $
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
// This provides an interface to KalmanTask::filter.
/******************************************************
% kalman_filter.cpp : Defines the entry point for
% Computing the likelihood of a stationnary state space model.
% It is called from Dynare DsgeLikelihood.m,
%
% function [LIK, lik] = kalman_filter_dll(T,R,Q,H,P,Y,start,Z/mf[, kalman_tol,riccati_tol])
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% R [double] mm*rr matrix, mapping structural innovations to state variables.
% Q [double] rr*rr covariance matrix of the structural innovations.
% H [double] pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.
% P [double] mm*mm variance-covariance matrix with stationary variables
% Y [double] pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
% start [integer] scalar, likelihood evaluation starts at 'start'.
% Z or mf: [double] Z: pp*mm matrix mapping state to pp observations
% Alternative parameters
% mf [integer] pp*1 vector of indices - alternative to Z matrix.
% Additional optional parameters
% kalman_tol [double] scalar, tolerance parameter (rcond).
% riccati_tol [double] scalar, tolerance parameter (riccati iteration).
%
% OUTPUTS
% LIK [double] scalar, likelihood
% lik [double] vector, density of observations in each period.
%
% REFERENCES
% See "Filtering and Smoothing of State Vector for Diffuse State Space
% Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
% Analysis, vol. 24(1), pp. 85-98).
%
% NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood.
**********************************************************/
#include "kalman.h"
#include "ts_exception.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include "mex.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nrhs, const mxArray* prhs[])
{
if (nrhs < 8 || nrhs > 10)
mexErrMsgTxt("Must have 8, 9, or 10 input parameters.\n");
if (nlhs < 1 || nlhs > 3)
mexErrMsgTxt("Must have 1, 2, or 3 output parameters.\n");
int start = 1; // default start of likelihood calculation
try
{
// make input matrices
GeneralMatrix T(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
GeneralMatrix R(mxGetPr(prhs[1]), mxGetM(prhs[1]), mxGetN(prhs[1]));
GeneralMatrix Q(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
GeneralMatrix H(mxGetPr(prhs[3]), mxGetM(prhs[3]), mxGetN(prhs[3]));
GeneralMatrix P(mxGetPr(prhs[4]), mxGetM(prhs[4]), mxGetN(prhs[4]));
GeneralMatrix Y(mxGetPr(prhs[5]), mxGetM(prhs[5]), mxGetN(prhs[5]));
if (nrhs>6) start = (int)mxGetScalar(prhs[6]);
GeneralMatrix Z(mxGetPr(prhs[7]), mxGetM(prhs[7]), mxGetN(prhs[7]));
int nper = mxGetN(prhs[5]); // no of periods
GeneralMatrix a( mxGetN(prhs[0]), 1);// initiate inital state to 0s
a.zeros();
#ifdef DEBUG
mexPrintf("kalman_filter: periods = %d ", nper);
#endif
// make storage for output
double loglik;
int per;
int d;
// output for full log-likelihood array
std::vector<double>* vll=new std::vector<double> (nper);
// create state init
StateInit* init = NULL;
init = new StateInit(P, a.getData());
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
// developement of the output.
#ifdef DEBUG
mexPrintf("kalman_filter: running and filling outputs.\n");
#endif
KalmanUniTask kut(kt);
loglik = kut.filter(per, d, (start-1), vll);
per = per / Y.numRows();
d = d / Y.numRows();
// destroy init
delete init;
// create output and upload output data
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglik);
if (nlhs >= 2)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
plhs[1] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[1]);
// allocate likelihood array
for (int j=0;j<nper;++j)
mxll[j]=(*vll)[j];
}
if (nlhs >= 3)
{
plhs[1] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[2]))) = per;
}
if (nlhs == 4)
{
plhs[2] = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
(*((int*)mxGetData(plhs[3]))) = d;
}
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
}
};

BIN
mex/sources/kalman/qt/FortranSubroutines.zip
View File

Binary file not shown.

37
mex/sources/kalman/qt/cc/qt.h
View File

@ -1,37 +0,0 @@
#ifndef QT_H
#define QT_H
#define C_CHAR const char*
#define BLINT int*
#define C_BLINT const int*
#define C_BLDOU const double*
#define BLDOU double*
extern "C"{
void ldm_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT);
void ldsld_(BLDOU, C_BLDOU, C_BLDOU, C_BLINT);
void ldv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void mtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void qt2ld_(BLDOU,C_BLDOU, C_BLINT);
void qt2t_(BLDOU, C_BLDOU, C_BLINT);
void s2d_(BLDOU,C_BLDOU, C_BLINT);
void s2u_(BLDOU,C_BLDOU, C_BLINT);
void td_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tm_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tstt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tu_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tut_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void tv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
void qtv_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
#ifdef WINDOWS
void qtv_1__(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
#else
void qtv_1_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
#endif
void qtsqtt_(BLDOU,C_BLDOU, C_BLDOU, C_BLINT);
};
#endif

24
mex/sources/kalman/qt/f90/LdM.f90
View File

@ -1,24 +0,0 @@
subroutine LdM(X,Ld,M,n)
! COMPUTATIONAL SUBROUTINE: LdM=Ld*M; Ld lower diagonal matrix; M arbitrary
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: Ld(n,n)
real(8), intent(in) :: M(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*M
do i=2,n
if ((Ld(i,i-1)/=0.0)) then
do j=1,n
X(i,j)=Ld(i,i-1)*M(i-1,j)
end do
end if
end do
return
end subroutine LdM

34
mex/sources/kalman/qt/f90/LdSLd.f90
View File

@ -1,34 +0,0 @@
subroutine LdSLd(X,Ld,S,n)
! COMPUTATIONAL SUBROUTINE: LdM=Ld*S*Ld; Ld lower diagonal matrix; S symmetric
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: Ld(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4), dimension(1,n) :: vk
integer(4) :: i,j, jj, h
h=0.0
do i=2,n
if ((Ld(i,i-1)/= 0.0)) then
h=h+1.0
vk(1,h)=i
end if
end do
if (h==0.0) then
X=0*S
return
end if
do j=1,h
do jj=1,h
X(vk(1,j),vk(1,jj))=Ld(vk(1,j),vk(1,j)-1)*S(vk(1,j)-1,vk(1,jj)-1)*Ld(vk(1,jj),vk(1,jj)-1)
end do
end do
return
end subroutine LdSLd

17
mex/sources/kalman/qt/f90/LdV.f90
View File

@ -1,17 +0,0 @@
subroutine LdV(X,Ld,v,n)
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: ld(n,n)
real(8), intent(in) :: v(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i
X=0*v
do i=2,n
X(i,1)=Ld(i,i-1)*v(i-1,1)
end do
return
end subroutine LdV

31
mex/sources/kalman/qt/f90/MTt.f90
View File

@ -1,31 +0,0 @@
subroutine MTt(X,M,Tt,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: M(n,n)
real(8), intent(in) :: Tt(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
X=0.0*M
do i=1,n
do j=1,n
stemp = 0.0
do k = j,n
stemp = stemp + M(i,k+i) * Tt(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine MTt

37
mex/sources/kalman/qt/f90/Makefile
View File

@ -1,37 +0,0 @@
# $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Makefile,v 1.4 2005/01/18 21:28:26 kamenik Exp $
# Tag $Name: $
FF=gfortran
CC_FLAGS := -Wall -Winline -fpic
#CC_FLAGS := -Wall -Winline -I../testing -I../cc -DMATLAB \
# $(CC_INCLUDE_PATH) -I$(MATLAB_PATH)/extern/include #-pg
LDFLAGS = -Wl,-L$(MATLAB_PATH)/extern/lib/win32/microsoft/ \
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
-lg2c -lmingw32 -lstdc++
LD_LIBS=$(LDFLAGS)
ifeq ($(DEBUG),yes)
# CC_FLAGS := $(CC_FLAGS) -g -DTL_DEBUG=2
CC_FLAGS := $(CC_FLAGS) -g -DPOSIX_THREADS
else
CC_FLAGS := $(CC_FLAGS) -O3
endif
ifeq ($(OS),Windows_NT)
CC_FLAGS := -mno-cygwin -mthreads $(CC_FLAGS)
endif
objects := $(patsubst %.f90,%.o,$(wildcard *.f90))
all: $(objects)
clear:
rm -f *.o
%.o : %.f90
$(FF) $(CC_FLAGS) $(EXTERN_DEFS) -c $*.f90

26
mex/sources/kalman/qt/f90/QT2Ld.f90
View File

@ -1,26 +0,0 @@
subroutine QT2Ld(X,QT,n)
! COMPUTATIONAL SUBROUTINE: extracts lower diagonal from Quasi triangular matrix
! COMPUTATIONAL SUBROUTINE: extracts lower diagonal from Quasi triangular matrix
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i
X=0.0*QT
do i=2,n
if (QT(i,i-1)/=0.0) then
X(i,i-1)=QT(i,i-1)
end if
end do
return
end subroutine QT2Ld

20
mex/sources/kalman/qt/f90/QT2T.f90
View File

@ -1,20 +0,0 @@
subroutine QT2T(X,QT,n)
! COMPUTATIONAL SUBROUTINE: extracts upper triangular from Quasi triangular matrix
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*QT
do i=1,n
do j=i,n
X(i,j)=QT(i,j)
end do
end do
return
end subroutine QT2T

50
mex/sources/kalman/qt/f90/QTSQTt.f90
View File

@ -1,50 +0,0 @@
subroutine QTSQTt(X,QT,S,n)
! COMPUTATIONAL SUBROUTINE X=QT*SQT' QT upper quasi-triangular; S symmetric
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h,k0
real(8) :: stemp
X=0*S
do i=1,n
do j=i, n
if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
h=i-1
else
h=i
end if
if (j > 1 .AND. (QT(j,j-1)/= 0.0)) then
k=j-1
k0=k
else
k=j
k0=k
end if
stemp=0
do while (h <= n)
do while (k <= n)
stemp=stemp+QT(i,h)*S(h,k)*QT(j,k)
k=k+1
end do
k=k0
h=h+1
end do
X(i,j)=stemp
if (i /= j) then
X(j,i)=stemp
end if
end do
end do
return
end subroutine QTSQTt

33
mex/sources/kalman/qt/f90/QTV.f90
View File

@ -1,33 +0,0 @@
subroutine QTV(X,QT,V,n)
! COMPUTATIONAL SUBROUTINE: X=QT*V QT upper quasi-triangular; V vector
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(in) :: V(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i,k
real(8) :: stemp
do i=1,n
stemp = 0.0
if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
stemp = QT(i,i-1)*v(i-1,1)
end if
do k = 0,n-i
stemp = stemp + QT(i,k+i) * V(k+i,1)
end do
X(i,1)=stemp
end do
return
end subroutine QTV

33
mex/sources/kalman/qt/f90/QTV_1.f90
View File

@ -1,33 +0,0 @@
subroutine QTV_1(X,QT,V,n)
! COMPUTATIONAL SUBROUTINE: X=QT*V QT upper quasi-triangular; V vector
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: QT(n,n)
real(8), intent(in) :: V(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i,k
real(8) :: stemp
X=0*V
do i=1,n
stemp = 0.0
if (i > 1 .AND. (QT(i,i-1)/= 0.0)) then
stemp = QT(i,i-1)*v(i-1,1)
end if
do k = 0,n-i
stemp = stemp + QT(i,k+i) * V(k+i,1)
end do
X(i,1)=stemp
end do
return
end subroutine QTV_1

19
mex/sources/kalman/qt/f90/S2D.f90
View File

@ -1,19 +0,0 @@
subroutine S2D(X,S,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty upper triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i
X=0.0*S
do i=1,n
X(i,i)=sqrt(S(i,i))
end do
return
end subroutine S2D

23
mex/sources/kalman/qt/f90/S2U.f90
View File

@ -1,23 +0,0 @@
subroutine S2U(X,S,n)
! COMPUTATIONAL SUBROUTINE: S2U extracts striclty upper triangular from symmetric
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*S
do i=1,n
do j=i+1,n
X(i,j)=S(i,j)
end do
end do
return
end subroutine S2U

22
mex/sources/kalman/qt/f90/TD.f90
View File

@ -1,22 +0,0 @@
subroutine TD(X,T,D,n)
! COMPUTATIONAL SUBROUTINE TUt=T*D T upper triangular; D diagonal
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: D(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j
X=0.0*T
do i=1,n
do j=i,n
X(i,j)=T(i,j)*D(j,j)
end do
end do
return
end subroutine TD

29
mex/sources/kalman/qt/f90/TM.f90
View File

@ -1,29 +0,0 @@
subroutine TM(X,T,M,n)
! COMPUTATIONAL SUBROUTINE: TM=T*M T upper triangular; M arbitrary
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: M(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
do i=1,n
do j=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T(i,k+i) * M(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TM

34
mex/sources/kalman/qt/f90/TSTt.f90
View File

@ -1,34 +0,0 @@
subroutine TSTt(X,T,S,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: S(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*S
do i=1,n
do j=1,n
stemp = 0.0
do h = i,n
do k = j,n
stemp = stemp + T(i,h) * S(h,k) * T(j,k)
end do
end do
X(i,j)=stemp
X(j,i)=stemp
end do
end do
return
end subroutine TSTt

29
mex/sources/kalman/qt/f90/TT.f90
View File

@ -1,29 +0,0 @@
subroutine TT(X, T1,T2,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T1(n,n)
real(8), intent(in) :: T2(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*T1
do i=1,n
do j=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T1(i,k+i) * T2(k+i,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TT

31
mex/sources/kalman/qt/f90/TU.f90
View File

@ -1,31 +0,0 @@
subroutine TU(X,T,U,n)
! COMPUTATIONAL SUBROUTINE: TU=T*U; T upper triangular matrix; U strictly upper triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: U(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k
real(8) :: stemp
X=0.0*T
do i=1,n-1
do j=i+1,n
stemp = 0.0
do k = i,j-1
stemp = stemp + T(i,k) * U(k,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TU

29
mex/sources/kalman/qt/f90/TUt.f90
View File

@ -1,29 +0,0 @@
subroutine TUt(X,T,Ut,n)
! COMPUTATIONAL SUBROUTINE TUt=T*Ut T upper triangular; U striclty lower triangular
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: Ut(n,n)
real(8), intent(out) :: X(n,n)
integer(4) :: i,j,k,h
real(8) :: stemp
X=0.0*T
do i=1,n
do j=1,n-1
h=max(i,j)
stemp = 0.0
do k = 0,n-h
stemp = stemp + T(i,k+h) * Ut(k+h,j)
end do
X(i,j)=stemp
end do
end do
return
end subroutine TUt

27
mex/sources/kalman/qt/f90/TV.f90
View File

@ -1,27 +0,0 @@
subroutine TV(X,T,V,n)
! COMPUTATIONAL SUBROUTINE: TV=T*V T upper triangular; V vector
implicit none
integer(4), intent(in) :: n
real(8), intent(in) :: T(n,n)
real(8), intent(in) :: V(n,1)
real(8), intent(out) :: X(n,1)
integer(4) :: i,k
real(8) :: stemp
do i=1,n
stemp = 0.0
do k = 0,n-i
stemp = stemp + T(i,k+i) * V(k+i,1)
end do
X(i,1)=stemp
end do
return
end subroutine TV

85
mex/sources/kalman/qt/test/Makefile
View File

@ -1,85 +0,0 @@
# $Id: Makefile 531 2005-11-30 13:49:48Z kamenik $
# Copyright 2005, Ondra Kamenik
#DEBUG = yes
#LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
#-mno-cygwin -DWINDOWS
CC_FLAGS := -DMATLAB -DNO_BLAS_H -DNO_LAPACK_H \
-Wall -I../cc -I../../sylv/cc -I../cc \
-I$(MATLAB_PATH)/extern/include #-pg
ifeq ($(DEBUG),yes)
CC_FLAGS := -DDEBUG $(CC_FLAGS) -g -pg
# CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g
KALMANLIB := kalmanlib_dbg.a
else
# CC_FLAGS := $(CC_FLAGS) -O3
CC_FLAGS := -DTIMING_LOOP $(CC_FLAGS) -O3
KALMANLIB := kalmanlib.a
endif
# Added by GP
# LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
#LDFLAGS := -Wl,--library-path $(LD_LIBRARY_PATH)
# -Wl,-L'f:/MinGW/lib'
#-Wl,-L'C:/MinGW/lib/gcc-lib/i686-pc-mingw32/4.0.4' -Wl,-L'C:/MinGW/lib'
LD_LIBS := -Wl,-rpath-link,$(MATLAB_PATH)/bin/glnxa64 \
-Wl,-L$(MATLAB_PATH)/bin/glnxa64 \
-Wl,-lmex -lmx -lmwlapack -lmwblas -lmat -lm \
-Wl,-lstdc++ -lgfortran $(LDFLAGS)
# -Wl,-L'f:/CygWin/usr/local/atlas/lib'
# -Wl,-L'f:/CygWin/lib'
# $(LDFLAGS)
# LD_LIBS :=$(LDFLAGS)
# end add
#matrix_interface := GeneralMatrix Vector SylvException
#matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
#mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
#matcppsource := $(patsubst %, ../sylv/cc/%.cpp, $(matrix_interface))
qtf90source := $(wildcard ../f90/*.f90)
qtobjs := $(patsubst %.f90,%.o,$(qtf90source))
cppsource := $(wildcard *.cpp)
hsource := $(wildcard *.h)
objects := $(patsubst %.cpp,%.o,$(cppsource))
dummy.ch:
touch dummy.ch
%.o: %.cpp $(hsource) $(cppsource)
c++ $(CC_FLAGS) -c $*.cpp
dgemmmtm_exe.exe: dgemmmtm_exe.o $(objects) $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o dgemmmtm_exe.exe dgemmmtm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
qtmmmtm_exe.exe: qtmmmtm_exe.o $(objects) $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o qtmmmtm_exe.exe qtmmmtm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
dgemvm_exe.exe: dgemvm_exe.o $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o dgemvm_exe.exe dgemvm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
qtamvm_exe.exe: qtamvm_exe.o $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o qtamvm_exe.exe qtamvm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
qtvmvm_exe.exe: qtvmvm_exe.o $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o qtvmvm_exe.exe qtvmvm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
qtv1mvm_exe.exe: qtv1mvm_exe.o $(qtobjs) $(hsource) $(cppsource)
gcc $(CC_FLAGS) -o qtv1mvm_exe.exe qtv1mvm_exe.o ascii_array.o \
$(qtobjs) $(LD_LIBS)
all: $(objects) dgemmmtm_exe.exe qtmmmtm_exe.exe qtv1mvm_exe.exe dgemvm_exe.exe qtvmvm_exe.exe qtamvm_exe.exe # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
clear:
rm -f *.o
rm -f *.a
rm -f *.{exe,dll}

11
mex/sources/kalman/qt/test/QT_tab
View File

@ -1,11 +0,0 @@
0.153060262389052 0.001380442914015 0.001695719778633 -0.096495788866461 -0.006326378534057 0.000968937494692 0.00940481674898 -0.003274924112445 0.001481132666149 0.004112932319104 -0.012701307804641
0 0.073952133671219 -0.001504207761519 0.057433281320695 0.009662045343528 0.006617685912409 0.013462091408161 0.007891283166287 0.00060603779672 0.003369481465008 0.018748481144843
0 0 0.016722949277241 -0.043158612323585 -0.004579527748968 0.014353708194285 -0.002590188137874 0.015728887447029 -0.000497609190313 -0.000142760660804 -0.019083813953199
0 0 0 0.265453840405543 0.00050417813553 0.002313090008356 0.027009998927894 -0.00302237841155 -0.006605759131604 -0.009638158849163 0.01228621141328
0 0 0 0 0.040159616212026 -0.035581727103737 -0.016954216473738 -0.007949691886576 0.003628861926114 -0.001291923621455 0.039842615500298
0 0 0 0 0 0.10084154296869 -0.001042311199417 -0.005382521361517 0.001696402254541 -0.000540820717391 0.023206920949216
0 0 0 0 0 0 0.054131830294643 -0.020987700402038 0.000402689298563 -0.000018817482395 0.001917827090534
0 0 0 0 0 0 0 0.14427448853052 -0.000904268210001 -0.004155316544672 0.018404810011954
0 0 0 0 0 0 0 0 0.009309814517978 0.000447330038827 0.00529992488884
0 0 0 0 0 0 0 0 0 0.005213746475416 0.012395564932675
0 0 0 0 0 0 0 0 0 0 0.026158448255444

45
mex/sources/kalman/qt/test/QTmatlabTest.m
View File

@ -1,45 +0,0 @@
size=100;
AA=rand(size);
BB=rand(size);
[QT1 QT2 U1 U2]=qz(AA,BB);
a=rand(size,1);
Loops=10000
% Calling QT without use of Sylv Vector and General Matrix
t = clock;
[AAAA]=qtamvm(QT1,a,Loops);
QTcpp_noSylv_TaInnerLoop_time=etime(clock, t)
% Calling QT using of Sylv Vector and General Matrix
t = clock;
[AAAA]=qtmvm(QT1,a,Loops);
QTcppTaInnerLoop_time=etime(clock, t)
t = clock;
[AAAA]=qtmvm_sub(QT1,a,Loops);
QTcppTaInnerLoop_time=etime(clock, t)
t = clock;
for tt=1:Loops%0
[AAAA]=qtmvm_sub(QT1,a,1);
end
QTcppTaOuterLoop_time=etime(clock, t)
t = clock;
[AAAA]=gmvm(QT1,a,Loops);
GMcppTaInnrLoop_time=etime(clock, t)
t = clock;
for tt=1:Loops%0
[AAAA]=gmvm(QT1,a,1);
end
GMcppTaOuterLoop_time=etime(clock, t)
t = clock;
for tt=1:Loops%0
MTA=QT1*a;
end
matlabTa_time=etime(clock, t)

11
mex/sources/kalman/qt/test/aa_dat
View File

@ -1,11 +0,0 @@
0.130845306698951
0.784832609671869
0.213653301284573
0.81536432324564
0.021265975953734
0.226338606586172
0.447782600002914
0.43948408700417
0.687965758422965
0.34628969510696
0.996770191062048

84
mex/sources/kalman/qt/test/ascii_array.cpp
View File

@ -1,84 +0,0 @@
// ascii_matrix.cpp
// Based on work of 2005, Ondra Kamenik
#include "ascii_array.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <fstream>
#include <string>
// if the file doesn't exist, the number array is empty
void
AsciiNumberArray::GetMX(const char* fname, int INrows, int INcols)
{
rows = 0;
cols = INcols;
std::ifstream file(fname);
std::string line;
data=(double*)calloc(INcols*INrows,sizeof(double));
while (getline(file, line))
{
rows++;
int icols = 0;
const char delims[] = " \r\n\t";
char* lineptr = strdup(line.c_str());
char* tok = strtok(lineptr, delims);
while (tok)
{
icols++;
double item;
if (1 != sscanf(tok, "%lf", &item))
{
fprintf(stderr, "Couldn't parse a token %s as double.\n", tok);
exit(1);
}
data[(rows-1)*INcols+icols-1]=item;
tok = strtok(NULL, delims);
}
free(lineptr);
if (cols)
{
if (cols != icols)
{
fprintf(stderr, "Asserted a different number of columns.\n");
exit(1);
}
}
else
{
cols = icols;
}
}
}
void
AsciiNumberArray::WriteMX()
{
std::ofstream outFile(strcat(fname, "_out"));
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
outFile << data[i*cols+j] << "\t";
}
outFile << std::endl;
}
outFile.close();
}
void WriteMX(char* fname, double* data, int rows, int cols)
{
AsciiNumberArray OutArray;
OutArray.fname=fname;
OutArray.rows=rows;
OutArray.cols=cols;
OutArray.data=data;
OutArray.WriteMX();
}

15
mex/sources/kalman/qt/test/ascii_array.h
View File

@ -1,15 +0,0 @@
// ascii_matrix.h
// Based on work of Ondra Kamenik
struct AsciiNumberArray {
int rows;
int cols;
double * data;
char* fname;
void GetMX(const char* fname, int rows, int cols);
void WriteMX();
};
void WriteMX(char* fname, double* data, int rows, int cols);

117
mex/sources/kalman/qt/test/dgemmmtm_exe.cpp
View File

@ -1,117 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to BLAS dgemv f90 library function
% to multiply Quasi trinagular matrix (T) with a vector a
%
%
% use:
% dgemmmtm_exe QTt_file SS_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
%
% INPUTS
% QT [double] mm*mm transition matrix of the state equation.
% SS [double] mm*mm state cov matrix.
%
% OUTPUTS
% TSTt update [double] mm*mm state cov matrix updated.
% as file: a_file_out
**********************************************************/
#include "cppblas.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *TSTt, *TS ;
AsciiNumberArray QT, SS;
QT.GetMX(argv[1],n,n);
SS.GetMX(argv[2],n,n);
const double alpha=1.0;
const double beta=0.0;
// create output and upload output data
TS=(double *)calloc(n*n, sizeof(double));
TSTt=(double *)calloc(n*n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
/* DGEMM performs one of the matrix-matrix operations
*
* C := alpha*op( A )*op( B ) + beta*C,
*
* where op( X ) is one of
*
* op( X ) = X or op( X ) = X',
*
* void BLAS_dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n,
* CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
* CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
* BLDOU c, CONST_BLINT ldc);
*
* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
*/
BLAS_dgemm("N", "N", &n, &n, &n, &alpha, QT.data, &n, SS.data, &n, &beta, TS, &n);
BLAS_dgemm("N", "T", &n, &n, &n, &alpha, TS, &n, QT.data, &n, &beta, TSTt, &n);
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], TSTt,n,n);
free(TSTt);
free(TS);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

104
mex/sources/kalman/qt/test/dgemvm_exe.cpp
View File

@ -1,104 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to BLAS dgemv f90 library function
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = dgevm(QT,a)
%
% use:
% dgemvm_exe QTt_file a_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
% 1. Ta = dgemv(Ld;a;n)+TV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% a update [double] mm vector of the state equation.
% as file: a_file_out
**********************************************************/
#include "cppblas.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *Ta ;
AsciiNumberArray QT, a;
QT.GetMX(argv[1],n,n);
a.GetMX(argv[2],n,1);
const double alpha=1.0;
const double beta=0.0;
int inc =1;
// create output and upload output data
Ta=(double *)calloc(n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
// void BLAS_dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
// CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
// CONST_BLDOU beta, BLDOU y, CONST_BLINT incy);
BLAS_dgemv("N", &n, &n, &alpha, QT.data, &n, a.data, &inc, &beta, Ta, &inc);
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], Ta,n,1);
free(Ta);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

115
mex/sources/kalman/qt/test/qtamvm_exe.cpp
View File

@ -1,115 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtamvm(QT,a)
%
% use:
% qtamvm_exe QTt_file a_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
% 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
% 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% a update [double] mm vector of the state equation.
% as file: a_file_out
**********************************************************/
#include "qt.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *T1, *Ld, *TV, * Ta ;
AsciiNumberArray QT, a;
QT.GetMX(argv[1],n,n);
a.GetMX(argv[2],n,1);
T1=(double *)calloc(n*n, sizeof(double));
Ld=(double *)calloc(n*n,sizeof(double));
TV=(double *)calloc(n, sizeof(double));
// create output and upload output data
Ta=(double *)calloc(n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
// 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
qt2t_(T1, QT.data ,&n) ;
qt2ld_(Ld , QT.data,&n);
// 2. Ta = LdV(Ld;a;n)+TV(T1;a;n).
ldv_(Ta, Ld,a.data ,&n);
tv_(TV, T1 ,a.data,&n);
// Ta.add(1.0,TV);
for (int j=0; j<n;++j)
Ta[j]+=TV[j];
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], Ta,n,1);
free(T1);
free(Ld);
free(TV);
free(Ta);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

100
mex/sources/kalman/qt/test/qtmmmtm_exe.cpp
View File

@ -1,100 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to BLAS dgemm library function
% to multiply Quasi trinagular matrix (T) with another matrix
%
% use:
% dgemmmtm_exe QTt_file SS_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
%
% INPUTS
% QT [double] mm*mm transition matrix of the state equation.
% SS [double] mm*mm state cov matrix.
%
% OUTPUTS
% TSTt update [double] mm*mm state cov matrix updated.
% as file: a_file_out
**********************************************************/
#include "qt.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *TSTt ;
AsciiNumberArray QT, SS;
QT.GetMX(argv[1],n,n);
SS.GetMX(argv[2],n,n);
// create output and upload output data
TSTt=(double *)calloc(n*n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
/* qtsqtt_ performs one of the matrix-matrix operations
*
* C := QT*S*QT'
*/
qtsqtt_(TSTt, QT.data, SS.data, &n);
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], TSTt,n,n);
free(TSTt);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

105
mex/sources/kalman/qt/test/qtv1mvm_exe.cpp
View File

@ -1,105 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtvmvm(QT,a)
%
% use:
% qtvmvm_exe QTt_file a_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
% 2. Ta = QTV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% a update [double] mm vector of the state equation.
% as file: a_file_out
**********************************************************/
#include "qt.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *Ta ;
AsciiNumberArray QT, a;
QT.GetMX(argv[1],n,n);
a.GetMX(argv[2],n,1);
// create output and upload output data
Ta=(double *)calloc(n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
// 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
#ifdef WINDOWS
qtv_1__(Ta, QT.data, a.data, &n) ;
#else
qtv_1_(Ta, QT.data, a.data, &n) ;
#endif
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], Ta,n,1);
free(Ta);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

101
mex/sources/kalman/qt/test/qtvmvm_exe.cpp
View File

@ -1,101 +0,0 @@
/*
* Copyright (C) 2008-2009 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/******************************************************
%
% This provides an interface to QT f90 library by Andrea Pagano
% to multiply Quasi trinagular matrix (T) with a vector a
%
% function [a] = qtvmvm(QT,a)
%
% use:
% qtvmvm_exe QTt_file a_file size [loops - if enabled]
%
% NOTE: due to fortran matrix orientation, input matrices need to be passed
% as transposed so QTt instead QT
%
% 2. Ta = QTV(T1;a;n).
%
% INPUTS
% T [double] mm*mm transition matrix of the state equation.
% a [double] mm state vector.
%
% OUTPUTS
% a update [double] mm vector of the state equation.
% as file: a_file_out
**********************************************************/
#include "qt.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ascii_array.h"
#include <iostream>
#include <stdexcept>
#include <malloc.h>
int main(int argc, char* argv[])
{
if (argc < 3 )
{
printf("Must have min 2 input parameters.\n");
exit(1);
}
try
{
// make input matrices
int n=atoi(argv[3]);
double *Ta ;
AsciiNumberArray QT, a;
QT.GetMX(argv[1],n,n);
a.GetMX(argv[2],n,1);
// create output and upload output data
Ta=(double *)calloc(n, sizeof(double));
#ifdef TIMING_LOOP
int loops=1;//000;
if (argc >3 )
loops = atoi(argv[4]);
for (int tt=0;tt<loops;++tt)
{
#endif
#ifdef DEBUG
// QT.print();
#endif
// 1. T1 = QT2T(QT;n) and Ld = QT2Ld(QT;n);
qtv_(Ta, QT.data, a.data, &n) ;
#ifdef TIMING_LOOP
}
printf("QT array mvm: finished: %d loops\n",loops);
#endif
// create output and upload output data
WriteMX(argv[2], Ta,n,1);
free(Ta);
}
catch (std::exception e)
{
std::cout <<"Error" << std::endl;
}
}; //main

825
mex/sources/kalman/sylv/cc/GeneralMatrix.cpp
View File

@ -1,825 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.cpp,v 1.4 2004/11/24 20:41:59 kamenik Exp $ */
/* Tag $Name: $ */
#include "SylvException.h"
#include "GeneralMatrix.h"
#include "cppblas.h"
#include "cpplapack.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <cmath>
#include <limits>
//vector<int>nullVec(0);
int GeneralMatrix::md_length = 32;
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m)
: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
{
copy(m);
}
GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m)
: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
{
copy(m);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, const char* dummy)
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
{
for (int i = 0; i < m.rows; i++)
for (int j = 0; j < m.cols; j++)
get(j,i) = m.get(i,j);
}
GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m, const char* dummy)
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
{
for (int i = 0; i < m.rows; i++)
for (int j = 0; j < m.cols; j++)
get(j,i) = m.get(i,j);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(nrows*ncols), rows(nrows), cols(ncols), ld(nrows)
{
copy(m, i, j);
}
GeneralMatrix::GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.base()+m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld)
{}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b)
: data(a.rows*b.cols), rows(a.rows), cols(b.cols), ld(a.rows)
{
gemm("N", a, "N", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum)
: data(a.rows*b.rows), rows(a.rows), cols(b.rows), ld(a.rows)
{
gemm("N", a, "T", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b)
: data(a.cols*b.cols), rows(a.cols), cols(b.cols), ld(a.cols)
{
gemm("T", a, "N", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2)
: data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols)
{
gemm("T", a, "T", b, 1.0, 0.0);
}
/* generate new matrix b as subset or whole of matrix a but reordered by
vrows and vcols as Matlab b=a(vrows,vcols) where vectors vrows and vcols start from 1.
It ignores non-positive elements passing zero length vector is equivalent
to Matlab operator ":" = all elements of that dimension in its order */
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const vector<int>&vrows, const vector<int>&vcols)
{
int nrows=0, ncols=0;
if (vrows.size()==0 && vcols.size()==0)
{
*this=a;
return;
}
else
{
if (vrows.size()==0)
nrows=a.numRows();
else
{
for (int i=0;i<vrows.size();++i)
{
if (vrows[i]>0)
nrows++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in construction by vector.");
}
}
if (nrows>a.numRows())
throw SYLV_MES_EXCEPTION("Wrong dimensions for construction by vector.");
if (vcols.size()==0)
ncols=a.numCols();
else
{
for (int i=0;i<vcols.size();++i)
{
if (vcols[i]>0)
ncols++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in construction by vector.");
}
}
if (ncols>a.numCols())
throw SYLV_MES_EXCEPTION("Wrong dimensions for construction by vector.");
data= *(new Vector(nrows*ncols));
rows=nrows;
cols=ncols;
if(nrows*ncols==0) return;
for (int i=0;i<nrows;++i)
{
for (int j=0;j<nrows;++j)
if (vrows.size()!=0 && vcols.size()!=0)
{
if (vrows[i]>0 && vcols[j] >0)
get(i,j)=a.get(vrows[i]-1, vcols[j]-1);
}
else if (vrows.size()!=0 && vcols.size()==0)
{
if (vrows[i]>0 )
get(i,j)=a.get(vrows[i]-1, j);
}
else if (vrows.size()==0 && vcols.size()!=0)
{
if (vcols[j] >0)
get(i,j)=a.get(i, vcols[j]-1);
}
}
}
}
GeneralMatrix::~GeneralMatrix()
{
}
/* Matlab element product: this = this .*m */
void
GeneralMatrix::multElements(const GeneralMatrix& m)
{
if(cols!=m.numCols() || rows!=m.numRows())
throw SYLV_MES_EXCEPTION("multiply Element porduct: matrices must be same dimension.");
for (int i=0;i<cols;++i)
for (int j=0;j<rows;++j)
get(i,j)*=m.get(i,j);
};
/* emulates Matlab repmat: new matrix = multv*multh*this */
GeneralMatrix&
GeneralMatrix::repmat(int multv, int multh)
{
GeneralMatrix* repMat=(new GeneralMatrix ( multv*rows, multh*cols));
for (int i=0;i<multv;++i)
for (int j=0;j<multh;++j)
(*repMat).place(*this, multv*i, multh*j);
return *repMat;
};
void GeneralMatrix::place(const ConstGeneralMatrix& m, int i, int j)
{
if (i + m.numRows() > numRows() ||
j + m.numCols() > numCols())
throw SYLV_MES_EXCEPTION("Bad submatrix placement, matrix dimensions exceeded.");
GeneralMatrix tmpsub(*this, i, j, m.numRows(), m.numCols());
tmpsub.copy(m);
}
/* this = a*b */
void GeneralMatrix::mult(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b)
{
gemm("N", a, "N", b, 1.0, 0.0);
}
/* this = this + scalar*a*b */
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
double mult)
{
gemm("N", a, "N", b, mult, 1.0);
}
/* this = this + scalar*a*b' */
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
const char* dum, double mult)
{
gemm("N", a, "T", b, mult, 1.0);
}
/* this = this + scalar*a'*b */
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum,
const ConstGeneralMatrix& b, double mult)
{
gemm("T", a, "N", b, mult, 1.0);
}
/* this = this + scalar*a'*b' */
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
const ConstGeneralMatrix& b, const char* dum2, double mult)
{
gemm("T", a, "T", b, mult, 1.0);
}
void GeneralMatrix::addOuter(const ConstVector& a, double mult)
{
if (numRows() != numCols())
throw SYLV_MES_EXCEPTION("Matrix is not square in GeneralMatrix::addOuter.");
if (numRows() != a.length())
throw SYLV_MES_EXCEPTION("Wrong length of a vector in GeneralMatrix::addOuter.");
// since BLAS dsyr (symmetric rank 1 update) assumes symmetricity, we do this manually
for (int i = 0; i < numRows(); i++)
for (int j = i; j < numRows(); j++) {
double s = mult*a[i]*a[j];
get(i,j) = get(i,j) + s;
if (i != j)
get(j,i) = get(j,i) + s;
}
}
void GeneralMatrix::multRight(const ConstGeneralMatrix& m)
{
gemm_partial_right("N", m, 1.0, 0.0);
}
void GeneralMatrix::multLeft(const ConstGeneralMatrix& m)
{
gemm_partial_left("N", m, 1.0, 0.0);
}
void GeneralMatrix::multRightTrans(const ConstGeneralMatrix& m)
{
gemm_partial_right("T", m, 1.0, 0.0);
}
void GeneralMatrix::multLeftTrans(const ConstGeneralMatrix& m)
{
gemm_partial_left("T", m, 1.0, 0.0);
}
/* this = this * A^(-1) */
void
GeneralMatrix::multInvRight( GeneralMatrix&A)
{
// check or allocate tmp space for Transpose *this
/**
if (tmpGMp)
{
if (tmpGMp->numCols()!=rows || tmpGMp->numRows()!=cols)
delete (tmpGMp);
}
if (!tmpGMp)
********/
tmpGMp= new GeneralMatrix(cols,rows); // allocate space only once if and when needed!
// tmpGMp=(*this)' i.e. Transpose (*this)
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
tmpGMp->get(j,i) = get(i,j);
// check A and this suiability
const int mcols=A.numCols();
if (A.numRows() != mcols) {
throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
}
if (cols != A.numRows()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
}
if (rows > 0)
{
/* out =tmpGMp out = inv(A')*(*this)' = inv(A')*(*tmpGMp) */
int* ipiv = new int[cols];
int info;
const int mld=A.getLD();
LAPACK_dgetrf(&cols, &cols, A.base(), &cols, ipiv, &info);
LAPACK_dgetrs("T", &cols, &mcols, A.base(), &cols, ipiv, tmpGMp->base(),
&mld, &info);
delete [] ipiv;
// *this= Transpose(tmpGMp out)
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = tmpGMp->get(j,i);
}
delete tmpGMp;
}
// here we must be careful for ld
void GeneralMatrix::zeros()
{
if (ld == rows)
data.zeros();
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = 0;
}
}
void GeneralMatrix::unit()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
if (i == j)
get(i,j) = 1.0;
else
get(i,j) = 0.0;
}
void GeneralMatrix::nans()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = std::numeric_limits<double>::quiet_NaN();
}
void GeneralMatrix::infs()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = std::numeric_limits<double>::infinity();
}
// here we must be careful for ld
void GeneralMatrix::mult(double a)
{
if (ld == rows)
data.mult(a);
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) *= a;
}
}
// here we must be careful for ld
void GeneralMatrix::add(double a, const ConstGeneralMatrix& m)
{
if (m.numRows() != rows || m.numCols() != cols)
throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::add.");
if (ld == rows && m.ld == m.rows)
data.add(a, m.data);
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) += a*m.get(i,j);
}
}
void GeneralMatrix::add(double a, const ConstGeneralMatrix& m, const char* dum)
{
if (m.numRows() != cols || m.numCols() != rows)
throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::add.");
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) += a*m.get(j,i);
}
bool GeneralMatrix::isDiff(const GeneralMatrix& m, const double tol=0.0)const
{
if (m.numRows() != rows || m.numCols() != cols)
throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::isDiff.");
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
if (fabs(get(i,j) - m.get(i,j))>tol)
return true;
return false;
}
bool GeneralMatrix::isDiffSym(const GeneralMatrix& m, const double tol=0.0)const
{
if (m.numRows() != rows || m.numCols() != cols || m.numRows() != cols || m.numCols() != rows)
throw SYLV_MES_EXCEPTION("Matrix has different size or not square in GeneralMatrix::isDiffSym.");
for (int i = 0; i < cols; i++)
for (int j = 0; i+j < cols ; j++) // traverse the upper triangle only
if (fabs(get(j,j+i) - m.get(j,j+i))>tol) // along diagonals where higher changes occur
return true;
return false;
}
/* x = scalar(a)*x + scalar(b)*this*d */
void GeneralMatrix::multVec(double a, Vector& x, double b, const ConstVector& d) const
{
if (x.length() != rows || cols != d.length()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
int mm = rows;
int nn = cols;
double alpha = b;
int lda = ld;
int incx = d.skip();
double beta = a;
int incy = x.skip();
BLAS_dgemv("N", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
void GeneralMatrix::copy(const ConstGeneralMatrix& m, int ioff, int joff)
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = m.get(i+ioff,j+joff);
}
void
GeneralMatrix::copy(const ConstGeneralMatrix& m)
{
memcpy(data.base() ,m.getData().base() ,m.numCols()*m.numRows()*sizeof(double));
};
void
GeneralMatrix::copy(const GeneralMatrix& m)
{
memcpy(data.base(),m.getData().base(),m.numCols()*m.numRows()*sizeof(double));
};
void
GeneralMatrix::copyColumns(const GeneralMatrix& m, int istart, int iend, int ito)
{
if ((rows!=m.numRows())|| istart<iend|| istart> m.numCols()-1 || iend> m.numCols()-1
|| ito> cols-1 )
throw SYLV_MES_EXCEPTION("Wrong dimensions for copying matrix columns.");
memcpy(data.base()+ito*rows*sizeof(double)
,m.getData().base()+istart*rows*sizeof(double)
,(iend-istart+1)*rows*sizeof(double));
};
/* emulates Matlab command A(a,b)=B(c,d) where a,b,c,d are vectors or ":")*/
void
GeneralMatrix::AssignByVectors(GeneralMatrix& a, const vector<int>& vToRows, const vector<int>& vToCols
, const GeneralMatrix& b, const vector<int>& vrows, const vector<int>& vcols)
{
int nrows=0, ncols=0, tonrows=0, toncols=0;
const vector<int> *vpToCols=0, *vpToRows=0, *vpRows=0, *vpCols=0;
vector<int> *tmpvpToCols=0, *tmpvpToRows=0, *tmpvpRows=0, *tmpvpCols=0;
if (vToRows.size()==0 && vToCols.size()==0 &&vrows.size()==0 && vcols.size()==0)
a=b;
else
{
if (vToRows.size()==0)
{
tonrows=a.numRows();
// vpToRows=new vector<int>(tonrows);
tmpvpToRows=new vector<int>(tonrows);
for (int i=0;i<tonrows;++i)
(*tmpvpToRows)[i]=i+1;
vpToRows=(const vector<int>*)tmpvpToRows;
}
else
{
for (int i=0;i<vToRows.size();++i)
{
if (vToRows[i]>0)
tonrows++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in assignment by vector.");
}
vpToRows=&vToRows;
}
if (tonrows>a.numRows())
throw SYLV_MES_EXCEPTION("Wrong dimensions for assignment by vector.");
if (vToCols.size()==0)
{
toncols=a.numCols();
tmpvpToCols=new vector<int>(toncols);
for (int i=0;i<toncols;++i)
(*tmpvpToCols)[i]=i+1;
vpToCols=(const vector<int>*)tmpvpToCols;
}
else
{
for (int i=0;i<vToCols.size();++i)
{
if (vToCols[i]>0)
toncols++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in assignment by vector.");
}
vpToCols=&vToCols;
}
if (toncols>a.numCols())
throw SYLV_MES_EXCEPTION("Wrong dimensions for assignment by vector.");
if (vrows.size()==0)
{
nrows=b.numRows();
tmpvpRows=new vector<int>(nrows);
for (int i=0;i<nrows;++i)
(*tmpvpToRows)[i]=i+1;
vpRows=(const vector<int>*)tmpvpRows;
}
else
{
for (int i=0;i<vrows.size();++i)
{
if (vrows[i]>0)
nrows++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in assignment by vector.");
}
vpRows=&vrows;
}
if (nrows>b.numRows())
throw SYLV_MES_EXCEPTION("Wrong dimensions for assignment by vector.");
if (vcols.size()==0)
{
ncols=b.numCols();
tmpvpCols=new vector<int>(ncols);
for (int i=0;i<ncols;++i)
(*tmpvpCols)[i]=i+1;
vpCols=(const vector<int>*)tmpvpCols;
}
else
{
for (int i=0;i<vcols.size();++i)
{
if (vcols[i]>0)
ncols++;
else
throw SYLV_MES_EXCEPTION("Non-positive indices in assignment by vector.");
}
vpCols=&vcols;
}
if (ncols>b.numCols())
throw SYLV_MES_EXCEPTION("Wrong dimensions for assignment by vector.");
if (tonrows!=nrows || toncols!=ncols)
throw SYLV_MES_EXCEPTION("Wrong indices dimensions for assignment by vector.");
if(!(nrows*ncols==0 || tonrows*toncols==0))
{
for (int i=0;i<nrows;++i)
{
for (int j=0;j<nrows;++j)
a.get((*vpToRows)[i]-1,(*vpToCols)[j]-1)=b.get((*vpRows)[i]-1, (*vpCols)[j]-1);
}
}
if (tmpvpToCols) delete(tmpvpToCols);
if (tmpvpToRows) delete(tmpvpToRows);
if (tmpvpRows) delete(tmpvpRows);
if (tmpvpCols) delete(tmpvpCols);
}
}
void GeneralMatrix::gemm(const char* transa, const ConstGeneralMatrix& a,
const char* transb, const ConstGeneralMatrix& b,
double alpha, double beta)
{
int opa_rows = a.numRows();
int opa_cols = a.numCols();
if (!strcmp(transa, "T")) {
opa_rows = a.numCols();
opa_cols = a.numRows();
}
int opb_rows = b.numRows();
int opb_cols = b.numCols();
if (!strcmp(transb, "T")) {
opb_rows = b.numCols();
opb_cols = b.numRows();
}
if (opa_rows != numRows() ||
opb_cols != numCols() ||
opa_cols != opb_rows) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
}
int m = opa_rows;
int n = opb_cols;
int k = opa_cols;
int lda = a.ld;
int ldb = b.ld;
int ldc = ld;
if (lda > 0 && ldb > 0 && ldc > 0) {
BLAS_dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda,
b.data.base(), &ldb, &beta, data.base(), &ldc);
} else if (numRows()*numCols() > 0) {
if (beta == 0.0)
zeros();
else
mult(beta);
}
}
void GeneralMatrix::gemm_partial_left(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta)
{
int icol;
for (icol = 0; icol + md_length < cols; icol += md_length) {
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, md_length);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, md_length);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
if (cols > icol) {
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, cols - icol);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, cols - icol);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
}
void GeneralMatrix::gemm_partial_right(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta)
{
int irow;
for (irow = 0; irow + md_length < rows; irow += md_length) {
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, md_length, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, md_length, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
if (rows > irow) {
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, rows - irow, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, rows - irow, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
}
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix& m)
: data(m.data), rows(m.rows), cols(m.cols), ld(m.ld) {}
double ConstGeneralMatrix::getNormInf() const
{
double norm = 0.0;
for (int i = 0; i < numRows(); i++) {
ConstVector rowi(data.base()+i, ld, cols);
double normi = rowi.getNorm1();
if (norm < normi)
norm = normi;
}
return norm;
}
double ConstGeneralMatrix::getNorm1() const
{
double norm = 0.0;
for (int j = 0; j < numCols(); j++) {
ConstVector colj(data.base()+ld*j, 1, rows);
double normj = colj.getNorm1();
if (norm < normj)
norm = normj;
}
return norm;
}
/* x = scalar(a)*x + scalar(b)*this*d */
void ConstGeneralMatrix::multVec(double a, Vector& x, double b, const ConstVector& d) const
{
if (x.length() != rows || cols != d.length()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
int mm = rows;
int nn = cols;
double alpha = b;
int lda = ld;
int incx = d.skip();
double beta = a;
int incy = x.skip();
BLAS_dgemv("N", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
void ConstGeneralMatrix::multVecTrans(double a, Vector& x, double b,
const ConstVector& d) const
{
if (x.length() != cols || rows != d.length()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
int mm = rows;
int nn = cols;
double alpha = b;
int lda = rows;
int incx = d.skip();
double beta = a;
int incy = x.skip();
BLAS_dgemv("T", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
/* m = inv(this)*m */
void ConstGeneralMatrix::multInvLeft(const char* trans, int mrows, int mcols, int mld, double* d) const
{
if (rows != cols) {
throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
}
if (cols != mrows) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
}
if (rows > 0) {
GeneralMatrix inv(*this);
int* ipiv = new int[rows];
int info;
LAPACK_dgetrf(&rows, &rows, inv.getData().base(), &rows, ipiv, &info);
LAPACK_dgetrs(trans, &rows, &mcols, inv.base(), &rows, ipiv, d,
&mld, &info);
delete [] ipiv;
}
}
/* m = inv(this)*m */
void ConstGeneralMatrix::multInvLeft(GeneralMatrix& m) const
{
multInvLeft("N", m.numRows(), m.numCols(), m.getLD(), m.getData().base());
}
/* m = inv(this')*m */
void ConstGeneralMatrix::multInvLeftTrans(GeneralMatrix& m) const
{
multInvLeft("T", m.numRows(), m.numCols(), m.getLD(), m.getData().base());
}
/* d = inv(this)*d */
void ConstGeneralMatrix::multInvLeft(Vector& d) const
{
if (d.skip() != 1) {
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
}
multInvLeft("N", d.length(), 1, d.length(), d.base());
}
/* d = inv(this')*d */
void ConstGeneralMatrix::multInvLeftTrans(Vector& d) const
{
if (d.skip() != 1) {
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
}
multInvLeft("T", d.length(), 1, d.length(), d.base());
}
bool ConstGeneralMatrix::isFinite() const
{
for (int i = 0; i < numRows(); i++)
for (int j = 0; j < numCols(); j++)
if (! std::isfinite(get(i,j)))
return false;
return true;
}
bool ConstGeneralMatrix::isZero() const
{
for (int i = 0; i < numRows(); i++)
for (int j = 0; j < numCols(); j++)
if (get(i,j) != 0.0)
return false;
return true;
}
void ConstGeneralMatrix::print() const
{
printf("rows=%d, cols=%d\n",rows, cols);
for (int i = 0; i < rows; i++) {
printf("row %d:\n",i);
for (int j = 0; j < cols; j++) {
printf("%6.3g ",get(i,j));
}
printf("\n");
}
}

355
mex/sources/kalman/sylv/cc/GeneralMatrix.h
View File

@ -1,355 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.h,v 1.3 2004/11/24 20:41:59 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef GENERAL_MATRIX_H
#define GENERAL_MATRIX_H
#include "Vector.h"
#include <vector>
using namespace std;
//#define nullVec (const vector<int>(0))
const vector<int>nullVec(0);
class GeneralMatrix;
class ConstGeneralMatrix {
friend class GeneralMatrix;
protected:
ConstVector data;
int rows;
int cols;
int ld;
public:
ConstGeneralMatrix(const double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
ConstGeneralMatrix(const GeneralMatrix& m);
ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols);
virtual ~ConstGeneralMatrix() {}
const double& get(int i, int j) const
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
const double* base() const {return data.base();}
const ConstVector& getData() const {return data;}
double getNormInf() const;
double getNorm1() const;
/* x = scalar(a)*x + scalar(b)*this*d */
void multVec(double a, Vector& x, double b, const ConstVector& d) const;
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector& x, double b, const ConstVector& d) const;
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, 1.0, d);}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, 1.0, d);}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, -1.0, d);}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, -1.0, d);}
/* m = inv(this)*m */
void multInvLeft(GeneralMatrix& m) const;
/* m = inv(this')*m */
void multInvLeftTrans(GeneralMatrix& m) const;
/* d = inv(this)*d */
void multInvLeft(Vector& d) const;
/* d = inv(this')*d */
void multInvLeftTrans(Vector& d) const;
bool isFinite() const;
/** Returns true of the matrix is exactly zero. */
bool isZero() const;
virtual void print() const;
protected:
void multInvLeft(const char* trans, int mrows, int mcols, int mld, double* d) const;
};
class GeneralMatrix {
friend class ConstGeneralMatrix;
protected:
Vector data;
int rows;
int cols;
int ld;
GeneralMatrix * tmpGMp;
public:
GeneralMatrix(int m, int n)
: data(m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(const double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(const GeneralMatrix& m);
GeneralMatrix(const ConstGeneralMatrix& m);
GeneralMatrix(const GeneralMatrix&m, const char* dummy); // transpose
GeneralMatrix(const ConstGeneralMatrix&m, const char* dummy); // transpose
GeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols);
/* this = a*b */
GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b);
/* this = a*b' */
GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum);
/* this = a'*b */
GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b);
/* this = a'*b */
GeneralMatrix(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2);
/* generate new matrix b as subset or whole of matrix a but
reordered by vrows and vcols as Matlab b=a(vrows,vcols)
ignores non-positive elements and passing zero length vector is equivalent to
Matlab operator ":" = all elements of that dimension in its order */
GeneralMatrix(const GeneralMatrix& a, const vector<int>& vrows, const vector<int>& vcols);
GeneralMatrix(const ConstGeneralMatrix& a, const vector<int> &vrows, const vector<int> &vcols)
{GeneralMatrix(GeneralMatrix(a), vrows, vcols);};
virtual ~GeneralMatrix();
const GeneralMatrix& operator=(const GeneralMatrix& m)
{
if (this!=&m)
{
if (rows==m.rows && cols==m.cols && ld==m.ld && data.base()!=m.data.base() )
copy(m);
else
{
data=m.data; rows=m.rows; cols=m.cols; ld=m.ld;
}
}
return *this;
}
/* emulate Matlab repmat: new matrix = multv*multh*this */
GeneralMatrix& repmat(int multv, int multh);
const double& get(int i, int j) const
{return data[j*ld+i];}
double& get(int i, int j)
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
double* base() {return data.base();}
const double* base() const {return data.base();}
Vector& getData() {return data;}
const Vector& getData() const {return data;}
double getNormInf() const
{return ConstGeneralMatrix(*this).getNormInf();}
double getNorm1() const
{return ConstGeneralMatrix(*this).getNorm1();}
/* place matrix m to the position (i,j) */
void place(const ConstGeneralMatrix& m, int i, int j);
void place(const GeneralMatrix& m, int i, int j)
{place(ConstGeneralMatrix(m), i, j);}
/* this = a*b */
void mult(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b);
void mult(const GeneralMatrix& a, const GeneralMatrix& b)
{mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));}
/* this = this + scalar*a*b */
void multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);}
/* this = this + scalar*a*b' */
void multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
const char* dum, double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
const char* dum, double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);}
/* this = this + scalar*a'*b */
void multAndAdd(const ConstGeneralMatrix& a, const char* dum, const ConstGeneralMatrix& b,
double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b,
double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);}
/* this = this + scalar*a'*b' */
void multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
const ConstGeneralMatrix& b, const char* dum2, double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2, double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b),dum2, mult);}
/* this = this + scalar*a*a' */
void addOuter(const ConstVector& a, double mult=1.0);
void addOuter(const Vector& a, double mult=1.0)
{addOuter(ConstVector(a), mult);}
/* this = this * m */
void multRight(const ConstGeneralMatrix& m);
void multRight(const GeneralMatrix& m)
{multRight(ConstGeneralMatrix(m));}
/* this = m * this */
void multLeft(const ConstGeneralMatrix& m);
void multLeft(const GeneralMatrix& m)
{multLeft(ConstGeneralMatrix(m));}
/* this = this * m' */
void multRightTrans(const ConstGeneralMatrix& m);
void multRightTrans(const GeneralMatrix& m)
{multRightTrans(ConstGeneralMatrix(m));}
/* this = m' * this */
void multLeftTrans(const ConstGeneralMatrix& m);
void multLeftTrans(const GeneralMatrix& m)
{multLeftTrans(ConstGeneralMatrix(m));}
/* Matlab element product: this = this .*m */
void multElements(const GeneralMatrix& m);
void multElements(const ConstGeneralMatrix& m)
{multElements(GeneralMatrix(m));}
/* this = this * m^(-1) */
void multInvRight(GeneralMatrix&m);
/* x = scalar(a)*x + scalar(b)*this*d */
void multVec(double a, Vector& x, double b, const ConstVector& d) const;
// {ConstGeneralMatrix(*this).multVec(a, x, b, d);}
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector& x, double b, const ConstVector& d) const
{ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);}
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVec(x, d);}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVecTrans(x, d);}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVec(x, d);}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVecTrans(x, d);}
/* this = zero */
void zeros();
/** this = unit (on main diagonal) */
void unit();
/* this = NaN */
void nans();
/* this = Inf */
void infs();
/* this = scalar*this */
void mult(double a);
/* this = this + scalar*m */
void add(double a, const ConstGeneralMatrix& m);
void add(double a, const GeneralMatrix& m)
{add(a, ConstGeneralMatrix(m));}
/* this = this + scalar*m' */
void add(double a, const ConstGeneralMatrix& m, const char* dum);
void add(double a, const GeneralMatrix& m, const char* dum)
{add(a, ConstGeneralMatrix(m), dum);}
/* Returns true if this and m matrices are different for more than tolerance tol */
bool isDiff(const GeneralMatrix& m, const double tol)const;
bool isDiffSym(const GeneralMatrix& m, const double tol)const;
bool isDiffUpprTriang(const GeneralMatrix& m, const double tol=0.0)const
{return isDiffSym(m, tol);}
bool isFinite() const
{return (ConstGeneralMatrix(*this)).isFinite();}
bool isZero() const
{return (ConstGeneralMatrix(*this)).isZero();}
virtual void print() const
{ConstGeneralMatrix(*this).print();}
void copyColumns(const GeneralMatrix& m, int istart, int iend, int ito);
void copyColumns(const ConstGeneralMatrix& m, int istart, int iend, int ito)
{copyColumns(GeneralMatrix( m), istart, iend, ito);};
/* emulates Matlab command A(a,b)=B(c,d) where a,b,c,d are Matlab index vectors starting from 1 or ":") */
static void AssignByVectors(GeneralMatrix& a, const vector<int>& vToRows, const vector<int>& vToCols
, const GeneralMatrix& b, const vector<int>& vrows, const vector<int>& vcols);
static void AssignByVectors(GeneralMatrix& a, const vector<int>& vToRows, const vector<int>& vToCols
, const ConstGeneralMatrix& b, const vector<int> &vrows, const vector<int> &vcols)
{AssignByVectors(a, vToRows, vToCols, GeneralMatrix(b), vrows, vcols);};
void AssignByVectors(const vector<int>& vToRows, const vector<int>& vToCols
, const GeneralMatrix& b, const vector<int>& vrows, const vector<int>& vcols)
{AssignByVectors( *this, vToRows, vToCols, b, vrows, vcols);};
void AssignByVectors( const vector<int>& vToRows, const vector<int>& vToCols
, const ConstGeneralMatrix& b, const vector<int> &vrows, const vector<int> &vcols)
{AssignByVectors(*this, vToRows, vToCols, GeneralMatrix(b), vrows, vcols);};
void AssignByVectors( const GeneralMatrix& b, const vector<int>& vrows, const vector<int>& vcols)
{AssignByVectors( *this, nullVec, nullVec, b, vrows, vcols);};
void AssignByVectors( const ConstGeneralMatrix& b, const vector<int> &vrows, const vector<int> &vcols)
{AssignByVectors(*this, nullVec, nullVec, GeneralMatrix(b), vrows, vcols);};
void AssignByVectors(const vector<int>& vToRows, const vector<int>& vToCols, const GeneralMatrix& b)
{AssignByVectors( *this, vToRows, vToCols, b, nullVec, nullVec);};
void AssignByVectors( const vector<int>& vToRows, const vector<int>& vToCols, const ConstGeneralMatrix& b)
{AssignByVectors(*this, vToRows, vToCols, GeneralMatrix(b), nullVec, nullVec);};
private:
void copy(const ConstGeneralMatrix& m, int ioff , int joff );
void copy(const GeneralMatrix& m, int ioff , int joff )
{copy(ConstGeneralMatrix(m), ioff, joff);}
void copy(const ConstGeneralMatrix& m);
void copy(const GeneralMatrix& m);
void gemm(const char* transa, const ConstGeneralMatrix& a,
const char* transb, const ConstGeneralMatrix& b,
double alpha, double beta);
void gemm(const char* transa, const GeneralMatrix& a,
const char* transb, const GeneralMatrix& b,
double alpha, double beta)
{gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b),
alpha, beta);}
/* this = this * op(m) (without whole copy of this) */
void gemm_partial_right(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta);
void gemm_partial_right(const char* trans, const GeneralMatrix& m,
double alpha, double beta)
{gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);}
/* this = op(m) *this (without whole copy of this) */
void gemm_partial_left(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta);
void gemm_partial_left(const char* trans, const GeneralMatrix& m,
double alpha, double beta)
{gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta);}
/* number of rows/columns for copy used in gemm_partial_* */
static int md_length;
};
#endif /* GENERAL_MATRIX_H */
// Local Variables:
// mode:C++
// End:

41
mex/sources/kalman/sylv/cc/Makefile
View File

@ -1,41 +0,0 @@
# $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Makefile,v 1.4 2005/01/18 21:28:26 kamenik Exp $
# Tag $Name: $
CC_FLAGS := -Wall -Winline -fpic -I../testing -I../cc -DMATLAB \
$(CC_INCLUDE_PATH) -I$(MATLAB)/extern/include #-pg
LDFLAGS = -Wl,-L$(MATLAB_PATH)/extern/lib/win32/microsoft/ \
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
-lg2c -lmingw32 -lstdc++
LD_LIBS=$(LDFLAGS)
ifeq ($(DEBUG),yes)
# CC_FLAGS := $(CC_FLAGS) -g -DTL_DEBUG=2
CC_FLAGS := $(CC_FLAGS) -g -DPOSIX_THREADS
else
CC_FLAGS := $(CC_FLAGS) -O3
endif
ifeq ($(OS),Windows_NT)
CC_FLAGS := -mno-cygwin -mthreads $(CC_FLAGS)
endif
objects := $(patsubst %.cpp,%.o,$(wildcard *.cpp))
headers := $(wildcard *.h)
all: $(objects)
sylvester.a: $(objects)
ar cr sylvester.a $(objects)
ranlib sylvester.a
clear:
rm -f *.o
rm -f sylvester.a
%.o : %.cpp $(headers)
$(CC) $(CC_FLAGS) $(EXTERN_DEFS) -c $*.cpp

69
mex/sources/kalman/sylv/cc/SylvException.cpp
View File

@ -1,69 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/SylvException.cpp,v 1.2 2004/10/01 10:30:40 kamenik Exp $ */
/* Tag $Name: $ */
#include "SylvException.h"
#include <string.h>
#include <stdio.h>
SylvException::SylvException(const char* f, int l, const SylvException* s)
{
strcpy(file,f);
line = l;
source = s;
}
SylvException::~SylvException()
{
if (source != NULL) {
delete source;
}
}
void SylvException::printMessage() const
{
char mes[1500];
mes[0] = '\0';
printMessage(mes, 1499);
printf(mes);
}
int SylvException::printMessage(char* str, int maxlen) const
{
int remain = maxlen;
if (source != NULL) {
remain = source->printMessage(str, maxlen);
}
char aux[100];
sprintf(aux, "From %s:%d\n", file, line);
int newremain = remain - strlen(aux);
if (newremain < 0) {
aux[remain] = '\0';
newremain = 0;
}
strcat(str, aux);
return newremain;
}
SylvExceptionMessage::SylvExceptionMessage(const char* f, int i,
const char* mes)
: SylvException(f,i,NULL)
{
strcpy(message,mes);
}
int SylvExceptionMessage::printMessage(char* str, int maxlen) const
{
char aux[600];
sprintf(aux, "At %s:%d:%s\n", file, line, message);
int newremain = maxlen - strlen(aux);
if (newremain < 0) {
aux[maxlen] = '\0';
newremain = 0;
}
strcat(str, aux);
return newremain;
}

39
mex/sources/kalman/sylv/cc/SylvException.h
View File

@ -1,39 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/SylvException.h,v 1.1.1.1 2004/06/04 13:00:44 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef SYLV_EXCEPTION_H
#define SYLV_EXCEPTION_H
#include "SylvMemory.h"
class SylvException : public MallocAllocator {
protected:
char file[50];
int line;
const SylvException* source;
public:
SylvException(const char* f, int l, const SylvException* s);
virtual ~SylvException();
virtual int printMessage(char* str, int maxlen) const;
void printMessage() const;
};
class SylvExceptionMessage : public SylvException {
char message[500];
public:
SylvExceptionMessage(const char* f, int l, const char* mes);
virtual int printMessage(char* str, int maxlen) const;
};
// define macros:
#define SYLV_EXCEPTION(exc) (SylvException(__FILE__, __LINE__, exc))
#define SYLV_MES_EXCEPTION(mes) (SylvExceptionMessage(__FILE__, __LINE__, mes))
#endif /* SYLV_EXCEPTION_H */
// Local Variables:
// mode:C++
// End:

63
mex/sources/kalman/sylv/cc/SylvMemory.h
View File

@ -1,63 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/SylvMemory.h,v 1.1.1.1 2004/06/04 13:00:49 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef SYLV_MEMORY_H
#define SYLV_MEMORY_H
//#include "SylvParams.h"
#include <new>
class MallocAllocator {
#ifdef USE_MEMORY_POOL
public:
void* operator new(size_t size);
void* operator new[](size_t size);
void operator delete(void* p);
void operator delete[](void* p);
#endif
};
/*
#ifdef USE_MEMORY_POOL
void* operator new(size_t size);
void* operator new[](size_t size);
void operator delete(void* p);
void operator delete[](void* p);
#endif
class SylvMemoryPool {
char* base;
size_t length;
size_t allocated;
bool stack_mode;
SylvMemoryPool(const SylvMemoryPool&);
const SylvMemoryPool& operator=(const SylvMemoryPool&);
public:
SylvMemoryPool();
~SylvMemoryPool();
void init(size_t size);
void* allocate(size_t size);
void free(void* p);
void reset();
void setStackMode(bool);
};
class SylvMemoryDriver {
SylvMemoryDriver(const SylvMemoryDriver&);
const SylvMemoryDriver& operator=(const SylvMemoryDriver&);
public:
SylvMemoryDriver(int num_d, int m, int n, int order);
SylvMemoryDriver(const SylvParams& pars, int num_d, int m, int n, int order);
static void setStackMode(bool);
~SylvMemoryDriver();
protected:
void allocate(int num_d, int m, int n, int order);
};
*/
#endif /* SYLV_MEMORY_H */
// Local Variables:
// mode:C++
// End:

370
mex/sources/kalman/sylv/cc/Vector.cpp
View File

@ -1,370 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Vector.cpp,v 1.1.1.1 2004/06/04 13:01:13 kamenik Exp $ */
/* Tag $Name: $ */
#include "Vector.h"
#include "GeneralMatrix.h"
#include "SylvException.h"
#include "cppblas.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <cmath>
#include <algorithm>
#include <limits>
using namespace std;
ZeroPad zero_pad;
Vector::Vector(const ConstVector& v)
: len(v.length()), s(1), data(new double[len]), destroy(true)
{
copy(v.base(), v.skip());
}
const Vector& Vector::operator=(const Vector& v)
{
if (this == &v)
return *this;
if (v.length() != length()) {
throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
}
if (s == v.s &&
(data <= v.data && v.data < data+len*s ||
v.data <= data && data < v.data+v.len*v.s) &&
(data-v.data) % s == 0) {
printf("this destroy=%d, v destroy=%d, data-v.data=%d, len=%d\n", destroy, v.destroy, data-v.data, len);
throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
}
copy(v.base(), v.skip());
return *this;
}
const Vector& Vector::operator=(const ConstVector& v)
{
if (v.length() != length()) {
throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
}
if (v.skip() == 1 && skip() == 1 && (
(base() < v.base() + v.length() && base() >= v.base()) ||
(base() + length() < v.base() + v.length() &&
base() + length() > v.base()))) {
throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
}
copy(v.base(), v.skip());
return *this;
}
void Vector::copy(const double* d, int inc)
{
int n = length();
int incy = skip();
BLAS_dcopy(&n, d, &inc, base(), &incy);
}
Vector::Vector(Vector& v, int off, int l)
: len(l), s(v.skip()), data(v.base()+off*v.skip()), destroy(false)
{
if (off < 0 || off + length() > v.length())
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
}
Vector::Vector(const Vector& v, int off, int l)
: len(l), s(1), data(new double[len]), destroy(true)
{
if (off < 0 || off + length() > v.length())
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
copy(v.base()+off*v.skip(), v.skip());
}
Vector::Vector(GeneralMatrix& m, int col)
: len(m.numRows()), s(1), data(&(m.get(0, col))), destroy(false)
{
}
Vector::Vector(int row, GeneralMatrix& m)
: len(m.numCols()), s(m.getLD()), data(&(m.get(row, 0))), destroy(false)
{
}
bool Vector::operator==(const Vector& y) const
{
return ConstVector(*this) == y;
}
bool Vector::operator!=(const Vector& y) const
{
return ConstVector(*this) != y;
}
bool Vector::operator<(const Vector& y) const
{
return ConstVector(*this) < y;
}
bool Vector::operator<=(const Vector& y) const
{
return ConstVector(*this) <= y;
}
bool Vector::operator>(const Vector& y) const
{
return ConstVector(*this) > y;
}
bool Vector::operator>=(const Vector& y) const
{
return ConstVector(*this) >= y;
}
void Vector::zeros()
{
if (skip() == 1) {
double* p = base();
for (int i = 0; i < length()/ZeroPad::length;
i++, p += ZeroPad::length)
memcpy(p, zero_pad.getBase(), sizeof(double)*ZeroPad::length);
for ( ; p < base()+length(); p++)
*p = 0.0;
} else {
for (int i = 0; i < length(); i++)
operator[](i) = 0.0;
}
}
void Vector::nans()
{
for (int i = 0; i < length(); i++)
operator[](i) = std::numeric_limits<double>::quiet_NaN();
}
void Vector::infs()
{
for (int i = 0; i < length(); i++)
operator[](i) = std::numeric_limits<double>::infinity();
}
Vector::~Vector()
{
if (destroy) {
delete [] data;
}
}
void Vector::rotatePair(double alpha, double beta1, double beta2, int i)
{
double tmp = alpha*operator[](i) - beta1*operator[](i+1);
operator[](i+1) = alpha*operator[](i+1) - beta2*operator[](i);
operator[](i) = tmp;
}
void Vector::add(double r, const Vector& v)
{
add(r, ConstVector(v));
}
void Vector::add(double r, const ConstVector& v)
{
int n = length();
int incx = v.skip();
int incy = skip();
BLAS_daxpy(&n, &r, v.base(), &incx, base(), &incy);
}
void Vector::add(const double* z, const Vector& v)
{
add(z, ConstVector(v));
}
void Vector::add(const double* z, const ConstVector& v)
{
int n = length()/2;
int incx = v.skip();
int incy = skip();
BLAS_zaxpy(&n, z, v.base(), &incx, base(), &incy);
}
void Vector::mult(double r)
{
int n = length();
int incx = skip();
BLAS_dscal(&n, &r, base(), &incx);
}
void Vector::mult2(double alpha, double beta1, double beta2,
Vector& x1, Vector& x2,
const Vector& b1, const Vector& b2)
{
x1.zeros();
x2.zeros();
mult2a(alpha, beta1, beta2, x1, x2, b1, b2);
}
void Vector::mult2a(double alpha, double beta1, double beta2,
Vector& x1, Vector& x2,
const Vector& b1, const Vector& b2)
{
x1.add(alpha, b1);
x1.add(-beta1, b2);
x2.add(alpha, b2);
x2.add(-beta2, b1);
}
double Vector::getNorm() const
{
ConstVector v(*this);
return v.getNorm();
}
double Vector::getMax() const
{
ConstVector v(*this);
return v.getMax();
}
double Vector::getNorm1() const
{
ConstVector v(*this);
return v.getNorm1();
}
double Vector::dot(const Vector& y) const
{
return ConstVector(*this).dot(ConstVector(y));
}
bool Vector::isFinite() const
{
return (ConstVector(*this)).isFinite();
}
void Vector::print() const
{
for (int i = 0; i < length(); i++) {
printf("%d\t%8.4g\n", i, operator[](i));
#ifdef MATLAB
//#include "mex.h"
mexPrintf("%d\t%8.4g\n", i, operator[](i));
#endif
}
}
ConstVector::ConstVector(const Vector& v, int off, int l)
: BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
{
if (off < 0 || off + length() > v.length()) {
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
}
}
ConstVector::ConstVector(const ConstVector& v, int off, int l)
: BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
{
if (off < 0 || off + length() > v.length()) {
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
}
}
ConstVector::ConstVector(const double* d, int skip, int l)
: BaseConstVector(l, skip, d)
{
}
ConstVector::ConstVector(const ConstGeneralMatrix& m, int col)
: BaseConstVector(m.numRows(), 1, &(m.get(0, col)))
{
}
ConstVector::ConstVector(int row, const ConstGeneralMatrix& m)
: BaseConstVector(m.numCols(), m.getLD(), &(m.get(row, 0)))
{
}
bool ConstVector::operator==(const ConstVector& y) const
{
if (length() != y.length())
return false;
if (length() == 0)
return true;
int i = 0;
while (i < length() && operator[](i) == y[i])
i++;
return i == length();
}
bool ConstVector::operator<(const ConstVector& y) const
{
int i = std::min(length(), y.length());
int ii = 0;
while (ii < i && operator[](ii) == y[ii])
ii++;
if (ii < i)
return operator[](ii) < y[ii];
else
return length() < y.length();
}
double ConstVector::getNorm() const
{
double s = 0;
for (int i = 0; i < length(); i++) {
s+=operator[](i)*operator[](i);
}
return sqrt(s);
}
double ConstVector::getMax() const
{
double r = 0;
for (int i = 0; i < length(); i++) {
if (abs(operator[](i))>r)
r = abs(operator[](i));
}
return r;
}
double ConstVector::getNorm1() const
{
double norm = 0.0;
for (int i = 0; i < length(); i++) {
norm += abs(operator[](i));
}
return norm;
}
double ConstVector::dot(const ConstVector& y) const
{
if (length() != y.length())
throw SYLV_MES_EXCEPTION("Vector has different length in ConstVector::dot.");
int n = length();
int incx = skip();
int incy = y.skip();
return BLAS_ddot(&n, base(), &incx, y.base(), &incy);
}
bool ConstVector::isFinite() const
{
int i = 0;
while (i < length() && isfinite(operator[](i)))
i++;
return i == length();
}
void ConstVector::print() const
{
for (int i = 0; i < length(); i++) {
printf("%d\t%8.4g\n", i, operator[](i));
}
}
ZeroPad::ZeroPad()
{
for (int i = 0; i < length; i++)
pad[i] = 0.0;
}

184
mex/sources/kalman/sylv/cc/Vector.h
View File

@ -1,184 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Vector.h,v 1.1.1.1 2004/06/04 13:01:13 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef VECTOR_H
#define VECTOR_H
/* NOTE! Vector and ConstVector have not common super class in order
* to avoid running virtual method invokation mechanism. Some
* members, and methods are thus duplicated */
#ifdef MATLAB
#include "mex.h"
#define printf mexPrintf
#endif
#include <stdio.h>
class GeneralMatrix;
class ConstVector;
class Vector {
protected:
int len;
int s;
double* data;
bool destroy;
public:
Vector() : len(0), s(1), data(0), destroy(false) {}
Vector(int l) : len(l), s(1), data(new double[l]), destroy(true) {}
Vector(Vector& v) : len(v.length()), s(v.skip()), data(v.base()), destroy(false) {}
Vector(const Vector& v)
: len(v.length()), s(1), data(new double[len]), destroy(true)
{copy(v.base(), v.skip());}
Vector(const ConstVector& v);
Vector(const double* d, int l)
: len(l), s(1), data(new double[len]), destroy(true)
{copy(d, 1);}
Vector(double* d, int l)
: len(l), s(1), data(d), destroy(false) {}
Vector(double* d, int skip, int l)
: len(l), s(skip), data(d), destroy(false) {}
Vector(Vector& v, int off, int l);
Vector(const Vector& v, int off, int l);
Vector(GeneralMatrix& m, int col);
Vector(int row, GeneralMatrix& m);
const Vector& operator=(const Vector& v);
const Vector& operator=(const ConstVector& v);
double& operator[](int i)
{return data[s*i];}
const double& operator[](int i) const
{return data[s*i];}
const double* base() const
{return data;}
double* base()
{return data;}
int length() const
{return len;}
int skip() const
{return s;}
/** Exact equality. */
bool operator==(const Vector& y) const;
bool operator!=(const Vector& y) const;
/** Lexicographic ordering. */
bool operator<(const Vector& y) const;
bool operator<=(const Vector& y) const;
bool operator>(const Vector& y) const;
bool operator>=(const Vector& y) const;
virtual ~Vector();
void zeros();
void nans();
void infs();
bool toBeDestroyed() const {return destroy;}
void rotatePair(double alpha, double beta1, double beta2, int i);
void add(double r, const Vector& v);
void add(double r, const ConstVector& v);
void add(const double* z, const Vector& v);
void add(const double* z, const ConstVector& v);
void mult(double r);
double getNorm() const;
double getMax() const;
double getNorm1() const;
double dot(const Vector& y) const;
bool isFinite() const;
void print() const;
/* multiplies | alpha -beta1| |b1| |x1|
| |\otimes I .| | = | |
| -beta2 alpha| |b2| |x2|
*/
static void mult2(double alpha, double beta1, double beta2,
Vector& x1, Vector& x2,
const Vector& b1, const Vector& b2);
/* same, but adds instead of set */
static void mult2a(double alpha, double beta1, double beta2,
Vector& x1, Vector& x2,
const Vector& b1, const Vector& b2);
/* same, but subtracts instead of add */
static void mult2s(double alpha, double beta1, double beta2,
Vector& x1, Vector& x2,
const Vector& b1, const Vector& b2)
{mult2a(-alpha, -beta1, -beta2, x1, x2, b1, b2);}
private:
void copy(const double* d, int inc);
const Vector& operator=(int); // must not be used (not implemented)
const Vector& operator=(double); // must not be used (not implemented)
};
class BaseConstVector {
protected:
int len;
int s;
const double* data;
public:
BaseConstVector(int l, int si, const double* d) : len(l), s(si), data(d) {}
BaseConstVector(const BaseConstVector& v) : len(v.len), s(v.s), data(v.data) {}
const BaseConstVector& operator=(const BaseConstVector& v)
{len = v.len; s = v.s; data = v.data; return *this;}
const double& operator[](int i) const
{return data[s*i];}
const double* base() const
{return data;}
int length() const
{return len;}
int skip() const
{return s;}
};
class ConstGeneralMatrix;
class ConstVector : public BaseConstVector {
public:
ConstVector(const Vector& v) : BaseConstVector(v.length(), v.skip(), v.base()) {}
ConstVector(const ConstVector& v) : BaseConstVector(v) {}
ConstVector(const double* d, int l) : BaseConstVector(l, 1, d) {}
ConstVector(const Vector& v, int off, int l);
ConstVector(const ConstVector& v, int off, int l);
ConstVector(const double* d, int skip, int l);
ConstVector(const ConstGeneralMatrix& m, int col);
ConstVector(int row, const ConstGeneralMatrix& m);
virtual ~ConstVector() {}
/** Exact equality. */
bool operator==(const ConstVector& y) const;
bool operator!=(const ConstVector& y) const
{return ! operator==(y);}
/** Lexicographic ordering. */
bool operator<(const ConstVector& y) const;
bool operator<=(const ConstVector& y) const
{return operator<(y) || operator==(y);}
bool operator>(const ConstVector& y) const
{return ! operator<=(y);}
bool operator>=(const ConstVector& y) const
{return ! operator<(y);}
double getNorm() const;
double getMax() const;
double getNorm1() const;
double dot(const ConstVector& y) const;
bool isFinite() const;
void print() const;
};
class ZeroPad {
public:
//static const int length = 16;
enum { length = 16};
private:
double pad[16];
public:
ZeroPad();
const double* getBase() const {return pad;}
};
#endif /* VECTOR_H */
// Local Variables:
// mode:C++
// End:

68
mex/sources/kalman/sylv/cc/cppblas.h
View File

@ -1,68 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/cppblas.h,v 1.2 2004/11/24 20:42:52 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef CPPBLAS_H
#define CPPBLAS_H
#if defined(MATLAB) && !defined(__linux__)
#define BLAS_dgemm dgemm
#define BLAS_dgemv dgemv
#define BLAS_dtrsv dtrsv
#define BLAS_dtrmv dtrmv
#define BLAS_daxpy daxpy
#define BLAS_dcopy dcopy
#define BLAS_zaxpy zaxpy
#define BLAS_dscal dscal
#define BLAS_dtrsm dtrsm
#define BLAS_ddot ddot
#else /* defined(MATLAB) && !defined(__linux__) */
#define BLAS_dgemm dgemm_
#define BLAS_dgemv dgemv_
#define BLAS_dtrsv dtrsv_
#define BLAS_dtrmv dtrmv_
#define BLAS_daxpy daxpy_
#define BLAS_dcopy dcopy_
#define BLAS_zaxpy zaxpy_
#define BLAS_dscal dscal_
#define BLAS_dtrsm dtrsm_
#define BLAS_ddot ddot_
#endif /* defined(MATLAB) && !defined(__linux__) */
#define BLCHAR const char*
#define CONST_BLINT const int*
#define CONST_BLDOU const double*
#define BLDOU double*
extern "C" {
void BLAS_dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n,
CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
BLDOU c, CONST_BLINT ldc);
void BLAS_dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
CONST_BLDOU beta, BLDOU y, CONST_BLINT incy);
void BLAS_dtrsv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
void BLAS_dtrmv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
void BLAS_daxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy);
void BLAS_dcopy(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy);
void BLAS_zaxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
BLDOU y, CONST_BLINT incy);
void BLAS_dscal(CONST_BLINT n, CONST_BLDOU a, BLDOU x, CONST_BLINT incx);
void BLAS_dtrsm(BLCHAR side, BLCHAR uplo, BLCHAR transa, BLCHAR diag, CONST_BLINT m,
CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
BLDOU b, CONST_BLINT ldb);
double BLAS_ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y,
CONST_BLINT incy);
};
#endif /* CPPBLAS_H */
// Local Variables:
// mode:C++
// End:

82
mex/sources/kalman/sylv/cc/cpplapack.h
View File

@ -1,82 +0,0 @@
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/cpplapack.h,v 1.3 2004/11/24 20:43:10 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef CPPLAPACK_H
#define CPPLAPACK_H
#if defined(MATLAB) && !defined(__linux__)
#define LAPACK_dgetrs dgetrs
#define LAPACK_dgetrf dgetrf
#define LAPACK_dgees dgees
#define LAPACK_dgecon dgecon
#define LAPACK_dtrexc dtrexc
#define LAPACK_dtrsyl dtrsyl
#define LAPACK_dpotrf dpotrf
#define LAPACK_dgges dgges
#define LAPACK_dsyev dsyev
#define LAPACK_dsyevr dsyevr
#else /* MATLAB */
#define LAPACK_dgetrs dgetrs_
#define LAPACK_dgetrf dgetrf_
#define LAPACK_dgees dgees_
#define LAPACK_dgecon dgecon_
#define LAPACK_dtrexc dtrexc_
#define LAPACK_dtrsyl dtrsyl_
#define LAPACK_dpotrf dpotrf_
#define LAPACK_dgges dgges_
#define LAPACK_dsyev dsyev_
#define LAPACK_dsyevr dsyevr_
#endif /* MATLAB */
#define LACHAR const char*
#define CONST_LAINT const int*
#define LAINT int*
#define CONST_LADOU const double*
#define LADOU double*
typedef int (*DGGESCRIT)(const double*, const double*, const double*);
extern "C" {
void LAPACK_dgetrs(LACHAR trans, CONST_LAINT n, CONST_LAINT nrhs,
CONST_LADOU a, CONST_LAINT lda, CONST_LAINT ipiv,
LADOU b, CONST_LAINT ldb, LAINT info);
void LAPACK_dgetrf(CONST_LAINT m, CONST_LAINT n, LADOU a,
CONST_LAINT lda, LAINT ipiv, LAINT info);
void LAPACK_dgees(LACHAR jobvs, LACHAR sort, const void* select,
CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT sdim,
LADOU wr, LADOU wi, LADOU vs, CONST_LAINT ldvs,
LADOU work, CONST_LAINT lwork, const void* bwork, LAINT info);
void LAPACK_dgecon(LACHAR norm, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda,
CONST_LADOU anorm, LADOU rnorm, LADOU work, LAINT iwork,
LAINT info);
void LAPACK_dtrexc(LACHAR compq, CONST_LAINT n, LADOU t, CONST_LAINT ldt,
LADOU q, CONST_LAINT ldq, LAINT ifst, LAINT ilst, LADOU work,
LAINT info);
void LAPACK_dtrsyl(LACHAR trana, LACHAR tranb, CONST_LAINT isgn, CONST_LAINT m,
CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU b,
CONST_LAINT ldb, LADOU c, CONST_LAINT ldc, LADOU scale,
LAINT info);
void LAPACK_dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
LAINT info);
void LAPACK_dgges(LACHAR jobvsl, LACHAR jobvsr, LACHAR sort, DGGESCRIT delztg,
CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU b, CONST_LAINT ldb,
LAINT sdim, LADOU alphar, LADOU alphai, LADOU beta,
LADOU vsl, CONST_LAINT ldvsl, LADOU vsr, CONST_LAINT ldvsr,
LADOU work, CONST_LAINT lwork, LAINT bwork, LAINT info);
void LAPACK_dsyev(LACHAR jobz, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
LADOU w, LADOU work, CONST_LAINT lwork, LAINT info);
void LAPACK_dsyevr(LACHAR jobz, LACHAR range, LACHAR uplo, CONST_LAINT n, LADOU a,
CONST_LAINT lda, LADOU lv, LADOU vu, CONST_LAINT il, CONST_LAINT iu,
CONST_LADOU abstol, LAINT m, LADOU w, LADOU z, CONST_LAINT ldz,
LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork,
LAINT info);
};
#endif /* CPPLAPACK_H */
// Local Variables:
// mode:C++
// End:

65
mex/sources/kalman/testing/Makefile
View File

@ -1,65 +0,0 @@
# $Id: Makefile 534 2005-11-30 13:58:11Z kamenik $
# Copyright 2005, Ondra Kamenik
#LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
#CC_FLAGS := -Wall -I../cc -I../../sylv/cc
CC_FLAGS := -Wall -I../cc -I../sylv/cc
DEBUG = yes
MATLAB = 1
# Added by GP
# LDFLAGS := -llapack -lcblas -lf77blas -latlas -lg2c -lstdc++ -lmingw32
LDFLAGS := -Wl,--library-path $(LD_LIBRARY_PATH) \
-Wl,-L'f:/MinGW/lib' \
-Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ \
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
-lg2c -lmingw32 -lstdc++ $(LDFLAGS)
# -Wl,-L'f:/CygWin/usr/local/atlas/lib' \
# -Wl,-L'f:/CygWin/lib' \
# $(LDFLAGS)
LD_LIBS :=$(LDFLAGS)
# end add
ifeq ($(DEBUG),yes)
CC_FLAGS := $(CC_FLAGS) -g
else
CC_FLAGS := $(CC_FLAGS) -O2
endif
matrix_interface := GeneralMatrix Vector SylvException
matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
#cwebsource := $(wildcard ../cc/*.cweb)
cppsource := $(wildcard ../cc/*.cpp)
#cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource))
mainobjects := $(patsubst %.cpp,%.o,$(cppsource))
#objects := $(patsubst %.cpp,%.o, *.cpp)
objects := $(wildcard *.o)
#hwebsource := $(wildcard ../cc/*.hweb)
hsource := $(wildcard ../cc/*.h)
#hsource := $(patsubst %.hweb,%.h,$(hwebsource))
../cc/dummy.ch:
make -C ../cc dummy.ch
../cc/%.o: ../cc/%.cpp $(hsource) $(mathsource)
make -C ../cc $*.o
%.o: %.cpp $(hsource) ascii_matrix.h
c++ $(CC_FLAGS) -c $*.cpp
#tests: $(hmainsource) $(cppmainsource) $(hsource) $(cppsource) objects) \
tests: $(hsource) $(cppsource) $(mainobjects) \
tests.o ascii_matrix.o
c++ $(CC_FLAGS) $(mainobjects) $(matobjs) ascii_matrix.o tests.o -o tests $(LD_LIBS)
clear:
make -C ../cc clear
rm -f *.o
rm -f tests

56
mex/sources/kalman/testing/ascii_matrix.cpp
View File

@ -1,56 +0,0 @@
// $Id: ascii_matrix.cpp 534 2005-11-30 13:58:11Z kamenik $
// Copyright 2005, Ondra Kamenik
#include "ascii_matrix.h"
#include <stdio.h>
#include <string.h>
#include <fstream>
#include <string>
// if the file doesn't exist, the number array is empty
void AsciiNumberArray::parse(const char* fname)
{
rows = 0;
cols = 0;
std::ifstream file(fname);
std::string line;
while (getline(file, line)) {
rows++;
int icols = 0;
const char delims[] = " \r\n\t";
char* lineptr = strdup(line.c_str());
char* tok = strtok(lineptr, delims);
while (tok) {
icols++;
double item;
if (1 != sscanf(tok, "%lf", &item)) {
fprintf(stderr, "Couldn't parse a token %s as double.\n", tok);
exit(1);
}
data.push_back(item);
tok = strtok(NULL, delims);
}
free(lineptr);
if (cols) {
if (cols != icols) {
fprintf(stderr, "Asserted a different number of columns.\n");
exit(1);
}
} else {
cols = icols;
}
}
}
AsciiMatrix::AsciiMatrix(const AsciiNumberArray& na)
: GeneralMatrix(na.rows, na.cols)
{
for (int i = 0; i < numRows(); i++)
for (int j = 0; j < numCols(); j++)
get(i, j) = na.data[i*numCols()+j];
}

24
mex/sources/kalman/testing/ascii_matrix.h
View File

@ -1,24 +0,0 @@
// $Id: ascii_matrix.h 534 2005-11-30 13:58:11Z kamenik $
// Copyright 2005, Ondra Kamenik
#include "GeneralMatrix.h"
#include <vector>
#include <string>
struct AsciiNumberArray {
int rows;
int cols;
std::vector<double> data;
AsciiNumberArray(const char* fname)
{parse(fname);}
AsciiNumberArray(std::string fname)
{parse(fname.c_str());}
protected:
void parse(const char* fname);
};
class AsciiMatrix : public GeneralMatrix {
public:
AsciiMatrix(const AsciiNumberArray& na);
};

33
mex/sources/kalman/testing/dynare_filter.m
View File

@ -1,33 +0,0 @@
%
% loglik = dynare_filter(Z,H,T,R,Q,Y,Pstar,Pinf)
%
% This is just an interface to DiffuseLikelihoodH1 of Dynare. It
% takes state space in the form
% y_t = Z*alpha_t + epsilon_t
% alpha_{t+1} = T*alpha_t + R*eta_t
% where epsilon covariance is H, eta covariance is Q
%
% It returns log likelihood.
%
% Copyright 2005, Ondra Kamenik
% $Id: dynare_filter.m 534 2005-11-30 13:58:11Z kamenik $
function lik = dynare_filter(Z,H,T,R,Q,Y,Pstar,Pinf)
global bayestopt_ options_
pp = size(Z,1);
mm = size(T,1);
rr = size(R,2);
dT = [zeros(pp,pp) Z*T; zeros(mm,pp) T];
dR = [eye(pp) Z*R; zeros(mm,pp) R];
dQ = [zeros(pp,pp) zeros(pp,rr); zeros(rr,pp) Q];
dPinf = [zeros(pp,pp) zeros(pp,mm); zeros(mm,pp) Pinf];
dPstar = [Z*Pstar*Z' Z*Pstar; Pstar*Z' Pstar];
bayestopt_.mf = [1:pp];
options_.kalman_tol = 1e-10;
lik = DiffuseLikelihoodH1(dT,dR,dQ,H,dPinf,dPstar,Y,zeros(pp,size(Y,2)),1);

37
mex/sources/kalman/testing/dynare_smoother.m
View File

@ -1,37 +0,0 @@
%
% [alpha,epsilon,eta] = dynare_smoother(Z,H,T,R,Q,Y,Pstar,Pinf)
%
% This is just an interface to DiffuseKalmanSmootherH1 of Dynare. It
% takes state space in the form
% y_t = Z*alpha_t + epsilon_t
% alpha_{t+1} = T*alpha_t + R*eta_t
% where epsilon covariance is H, eta covariance is Q
%
% It returns smoothed alpha, epsilon and eta.
%
% Copyright 2005, Ondra Kamenik
% $Id: dynare_smoother.m 534 2005-11-30 13:58:11Z kamenik $
function [alpha,epsilon,eta] = dynare_smoother(Z,H,T,R,Q,Y,Pstar,Pinf)
global options_
pp = size(Z,1);
mm = size(T,1);
rr = size(R,2);
dT = [zeros(pp,pp) Z*T; zeros(mm,pp) T];
dR = [eye(pp) Z*R; zeros(mm,pp) R];
dQ = [zeros(pp,pp) zeros(pp,rr); zeros(rr,pp) Q];
dPinf = [Z*Pinf*Z' Z*Pinf; Pinf*Z' Pinf];
dPstar = [Z*Pstar*Z' Z*Pstar; Pstar*Z' Pstar];
mf = [1:pp];
options_.kalman_tol = 1e-10;
options_.nk = 0;
% if you want DiffuseKalmanSmootherH3, uncomment the following and set
% diffuse_d (possibly empty [])
% options_.diffuse_d = [7];
[alpha,epsilon,eta] = DiffuseKalmanSmootherH1(dT,dR,dQ,H,dPinf,dPstar,Y,zeros(pp,size(Y,2)),pp,mm+pp,size(Y,2),mf);
alpha = alpha(pp+1:end,:);
eta = eta(pp+1:end,:);

20
mex/sources/kalman/testing/kalmandll_test.m
View File

@ -1,20 +0,0 @@
function [LIKDLL loglik]=kalmandll_test(T,mf,R,Q,H,Pstar,Pinf,data,start)
if isempty(H)
H=zeros(size(data,1), size(data,1))
elseif H==0
H=zeros(size(data,1), size(data,1))
end
Z=zeros(size(data,1), size(T,2))
for i = 1:size(data,1)
Z(i,mf(i))=1
end
LIKDLL= kalman_filter_dll4(T,Z,R,Q,H,Pstar,data,start)
%Y=data;
if isempty(Pinf)
Pinf=zeros(size(T));
elseif Pinf==0
Pinf=zeros(size(T));
end
% test DiffuseLikelihoodH1
loglik = dynare_filter(Z,H,T,R,Q,data,Pstar,Pinf)

BIN
mex/sources/kalman/testing/test_data.tgz
View File

Binary file not shown.

311
mex/sources/kalman/testing/tests.asv
View File

@ -1,311 +0,0 @@
// $Id: tests.cpp 534 2005-11-30 13:58:11Z kamenik $
// Copyright 2005, Ondra Kamenik
#include "../cc/kalman.h"
#include "../cc/ts_exception.h"
#include "ascii_matrix.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include <sys/time.h>
#include <math.h>
// gettimeofday for MinGW
#ifdef __MINGW32__
#define _W32_FT_OFFSET (116444736000000000LL)
typedef struct _filetime {
unsigned long dwLowDateTime;
unsigned long dwHighDateTime;
} filetime;
extern "C" {
void __stdcall GetSystemTimeAsFileTime(filetime*);
};
typedef union {
long long ns100; // time since 1 Jan 1601 in 100ns units
filetime ft;
} w32_ftv;
void D_gettimeofday(struct timeval* p, struct timezone* tz)
{
w32_ftv _now;
GetSystemTimeAsFileTime( &(_now.ft) );
p->tv_usec=(long)((_now.ns100 / 10LL) % 1000000LL );
p->tv_sec= (long)((_now.ns100-_W32_FT_OFFSET)/10000000LL);
return;
}
#else
#define D_gettimeofday gettimeofday
#endif // gettimeofday for MinGW
struct AsciiKalmanTask {
AsciiMatrix Z;
AsciiMatrix H;
AsciiMatrix T;
AsciiMatrix R;
AsciiMatrix Q;
AsciiMatrix Pstar;
AsciiMatrix Pinf;
AsciiMatrix a;
AsciiMatrix Y;
AsciiKalmanTask(const char* prefix)
: Z(std::string(prefix) + "_Z.dat"),
H(std::string(prefix) + "_H.dat"),
T(std::string(prefix) + "_T.dat"),
R(std::string(prefix) + "_R.dat"),
Q(std::string(prefix) + "_Q.dat"),
Pstar(std::string(prefix) + "_Pstar.dat"),
Pinf(std::string(prefix) + "_Pinf.dat"),
a(std::string(prefix) + "_a.dat"),
Y(std::string(prefix) + "_Y.dat")
{}
};
// WallTimer class. Constructor saves the wall time, destructor
// cancels the current time from the saved, and prints the message
// with time information
class WallTimer {
char mes[100];
struct timeval start;
bool new_line;
public:
WallTimer(const char* m, bool nl = true)
{strcpy(mes, m);new_line = nl; D_gettimeofday(&start, NULL);}
~WallTimer()
{
struct timeval end;
D_gettimeofday(&end, NULL);
printf("%s%8.4g", mes,
end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
if (new_line)
printf("\n");
}
};
/****************************************************/
/* declaration of TestRunnable class */
/****************************************************/
class TestRunnable {
char name[100];
public:
TestRunnable(const char* n)
{strncpy(name, n, 100);}
bool test() const;
virtual bool run() const =0;
const char* getName() const
{return name;}
protected:
static bool filter_and_smoother(const char* prefix, bool diffuse_flag);
};
bool TestRunnable::test() const
{
printf("Running test <%s>\n",name);
bool passed;
{
WallTimer tim("Wall clock time ", false);
passed = run();
}
if (passed) {
printf("............................ passed\n\n");
return passed;
} else {
printf("............................ FAILED\n\n");
return passed;
}
}
/****************************************************/
/* definition of TestRunnable static methods */
/****************************************************/
bool TestRunnable::filter_and_smoother(const char* prefix, bool diffuse_flag)
{
AsciiKalmanTask akt(prefix);
StateInit* init;
if (diffuse_flag)
init = new StateInit(akt.Pstar, akt.Pinf, akt.a.getData());
else
init = new StateInit(akt.Pstar, akt.a.getData());
KalmanTask kt(akt.Y, akt.Z, akt.H, akt.T, akt.R, akt.Q, *init);
// multivariate
int per;
int d;
double ll;
GeneralMatrix alpha(akt.T.numRows(), akt.Y.numCols());
GeneralMatrix eta(akt.R.numCols(), akt.Y.numCols());
GeneralMatrix V(akt.T.numRows(), akt.T.numRows()*akt.Y.numCols());
SmootherResults sres(akt.Y.numCols());
{
WallTimer tim("\tMultivariate time ", true);
ll = kt.filter_and_smooth(sres, per, d);
printf("\t\tll=%f per=%d d=%d\n", ll, per, d);
if (per == akt.Y.numCols()) {
sres.exportAlpha(alpha);
sres.exportEta(eta);
sres.exportV(V);
} else {
printf("\t\tNot finished.\n");
}
}
// univariate
KalmanUniTask kut(kt);
int per1;
int d1;
double ll1;
GeneralMatrix alpha1(akt.T.numRows(), akt.Y.numCols());
GeneralMatrix eta1(akt.R.numCols(), akt.Y.numCols());
GeneralMatrix V1(akt.T.numRows(), akt.T.numRows()*akt.Y.numCols());
SmootherResults sres1(akt.Y.numCols()*akt.Y.numRows());
{
WallTimer tim("\tUnivariate time ", true);
int dd;
ll1 = kut.filter_and_smooth(sres1, per1, dd);
per1 /= akt.Y.numRows();
d1 = dd/akt.Y.numRows();
printf("\t\tll=%f per=%d d=%d(%d)\n", ll1, per1, d1, dd);
if (per1 == akt.Y.numCols()) {
SmootherResults sres_uni(akt.Y.numCols());
sres_uni.import(sres1, akt.Y.numRows());
sres_uni.exportAlpha(alpha1);
sres_uni.exportEta(eta1);
sres_uni.exportV(V1);
} else {
printf("\t\tNot finished.\n");
}
}
// compare
if (per == per1 && per == akt.Y.numCols()) {
WallTimer tim("\tComparison time ", true);
alpha.add(-1.0, alpha1);
eta.add(-1.0, eta1);
V.add(-1.0, V1);
int maxd = std::max(d,d1);
for (int t = 1; t <= maxd; t++) {
Vector alphat(alpha, t-1);
printf("\t\tt=%d alpha error %10.6g\n",t,alphat.getMax());
Vector etat(eta, t-1);
printf("\t\tt=%d eta error %10.6g\n",t,etat.getMax());
GeneralMatrix Vt(V, 0, (t-1)*akt.T.numRows(), akt.T.numRows(), akt.T.numRows());
printf("\t\tt=%d V error %10.6g\n",t,V.getData().getMax());
}
GeneralMatrix alpha_rest(alpha, 0, maxd, akt.T.numRows(), alpha.numCols()-maxd);
printf("\t\tt=%d.. alpha error %10.6g\n",maxd+1,alpha_rest.getData().getMax());
GeneralMatrix eta_rest(eta, 0, maxd, akt.R.numCols(), eta.numCols()-maxd);
printf("\t\tt=%d.. eta error %10.6g\n",maxd+1,eta_rest.getData().getMax());
GeneralMatrix V_rest(V, 0, maxd*akt.T.numRows(), akt.T.numRows(),
V.numCols()-maxd*akt.T.numRows());
printf("\t\tt=%d.. V error %10.6g\n",maxd+1,V_rest.getData().getMax());
}
delete init;
return true;
}
/****************************************************/
/* definition of TestRunnable subclasses */
/****************************************************/
class SmallNonDiffuse : public TestRunnable {
public:
SmallNonDiffuse()
: TestRunnable("Non-diffuse small (p=2,m=3,r=4)") {}
bool run() const
{
filter_and_smoother("small2x3x4", false);
return true;
}
};
class SmallDiffuse : public TestRunnable {
public:
SmallDiffuse()
: TestRunnable("Diffuse small (p=2,m=3,r=4)") {}
bool run() const
{
return filter_and_smoother("small2x3x4", true);
}
};
class MiddleNonDiffuse : public TestRunnable {
public:
MiddleNonDiffuse()
: TestRunnable("Non-diffuse middle (p=10,m=15,r=12)") {}
bool run() const
{
return filter_and_smoother("10x15x12", false);
}
};
class MiddleDiffuse : public TestRunnable {
public:
MiddleDiffuse()
: TestRunnable("Diffuse middle (p=10,m=15,r=12)") {}
bool run() const
{
return filter_and_smoother("10x15x12", true);
}
};
class SOEDiffuse : public TestRunnable {
public:
SOEDiffuse()
: TestRunnable("Diffuse soe (p=8,m=25,r=15)") {}
bool run() const
{
return filter_and_smoother("soe8x25x15", true);
}
};
int main()
{
TestRunnable* all_tests[50];
// fill in vector of all tests
int num_tests = 0;
all_tests[num_tests++] = new SmallNonDiffuse();
all_tests[num_tests++] = new SmallDiffuse();
all_tests[num_tests++] = new MiddleNonDiffuse();
all_tests[num_tests++] = new MiddleDiffuse();
all_tests[num_tests++] = new SOEDiffuse();
// launch the tests
int success = 0;
for (int i = 0; i < num_tests; i++) {
try {
if (all_tests[i]->test())
success++;
} catch (const TSException& e) {
printf("Caugth TS exception in <%s>:\n", all_tests[i]->getName());
e.print();
} catch (SylvException& e) {
printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
e.printMessage();
}
}
printf("There were %d tests that failed out of %d tests run.\n",
num_tests - success, num_tests);
// destroy
for (int i = 0; i < num_tests; i++) {
delete all_tests[i];
}
return 0;
}

311
mex/sources/kalman/testing/tests.cpp
View File

@ -1,311 +0,0 @@
// $Id: tests.cpp 534 2005-11-30 13:58:11Z kamenik $
// Copyright 2005, Ondra Kamenik
#include "../cc/kalman.h"
#include "../cc/ts_exception.h"
#include "ascii_matrix.h"
#include "GeneralMatrix.h"
#include "Vector.h"
#include "SylvException.h"
#include <sys/time.h>
#include <math.h>
// gettimeofday for MinGW
#ifdef __MINGW32__
#define _W32_FT_OFFSET (116444736000000000LL)
typedef struct _filetime {
unsigned long dwLowDateTime;
unsigned long dwHighDateTime;
} filetime;
extern "C" {
void __stdcall GetSystemTimeAsFileTime(filetime*);
};
typedef union {
long long ns100; // time since 1 Jan 1601 in 100ns units
filetime ft;
} w32_ftv;
void D_gettimeofday(struct timeval* p, struct timezone* tz)
{
w32_ftv _now;
GetSystemTimeAsFileTime( &(_now.ft) );
p->tv_usec=(long)((_now.ns100 / 10LL) % 1000000LL );
p->tv_sec= (long)((_now.ns100-_W32_FT_OFFSET)/10000000LL);
return;
}
#else
#define D_gettimeofday gettimeofday
#endif // gettimeofday for MinGW
struct AsciiKalmanTask {
AsciiMatrix Z;
AsciiMatrix H;
AsciiMatrix T;
AsciiMatrix R;
AsciiMatrix Q;
AsciiMatrix Pstar;
AsciiMatrix Pinf;
AsciiMatrix a;
AsciiMatrix Y;
AsciiKalmanTask(const char* prefix)
: Z(std::string(prefix) + "_Z.dat"),
H(std::string(prefix) + "_H.dat"),
T(std::string(prefix) + "_T.dat"),
R(std::string(prefix) + "_R.dat"),
Q(std::string(prefix) + "_Q.dat"),
Pstar(std::string(prefix) + "_Pstar.dat"),
Pinf(std::string(prefix) + "_Pinf.dat"),
a(std::string(prefix) + "_a.dat"),
Y(std::string(prefix) + "_Y.dat")
{}
};
// WallTimer class. Constructor saves the wall time, destructor
// cancels the current time from the saved, and prints the message
// with time information
class WallTimer {
char mes[100];
struct timeval start;
bool new_line;
public:
WallTimer(const char* m, bool nl = true)
{strcpy(mes, m);new_line = nl; D_gettimeofday(&start, NULL);}
~WallTimer()
{
struct timeval end;
D_gettimeofday(&end, NULL);
printf("%s%8.4g", mes,
end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
if (new_line)
printf("\n");
}
};
/****************************************************/
/* declaration of TestRunnable class */
/****************************************************/
class TestRunnable {
char name[100];
public:
TestRunnable(const char* n)
{strncpy(name, n, 100);}
bool test() const;
virtual bool run() const =0;
const char* getName() const
{return name;}
protected:
static bool filter_and_smoother(const char* prefix, bool diffuse_flag);
};
bool TestRunnable::test() const
{
printf("Running test <%s>\n",name);
bool passed;
{
WallTimer tim("Wall clock time ", false);
passed = run();
}
if (passed) {
printf("............................ passed\n\n");
return passed;
} else {
printf("............................ FAILED\n\n");
return passed;
}
}
/****************************************************/
/* definition of TestRunnable static methods */
/****************************************************/
bool TestRunnable::filter_and_smoother(const char* prefix, bool diffuse_flag)
{
AsciiKalmanTask akt(prefix);
StateInit* init;
if (diffuse_flag)
init = new StateInit(akt.Pstar, akt.Pinf, akt.a.getData());
else
init = new StateInit(akt.Pstar, akt.a.getData());
KalmanTask kt(akt.Y, akt.Z, akt.H, akt.T, akt.R, akt.Q, *init);
// multivariate
int per;
int d;
double ll;
GeneralMatrix alpha(akt.T.numRows(), akt.Y.numCols());
GeneralMatrix eta(akt.R.numCols(), akt.Y.numCols());
GeneralMatrix V(akt.T.numRows(), akt.T.numRows()*akt.Y.numCols());
SmootherResults sres(akt.Y.numCols());
{
WallTimer tim("\tMultivariate time ", true);
ll = kt.filter_and_smooth(sres, per, d);
printf("\t\tll=%f per=%d d=%d\n", ll, per, d);
if (per == akt.Y.numCols()) {
sres.exportAlpha(alpha);
sres.exportEta(eta);
sres.exportV(V);
} else {
printf("\t\tNot finished.\n");
}
}
// univariate
KalmanUniTask kut(kt);
int per1;
int d1;
double ll1;
GeneralMatrix alpha1(akt.T.numRows(), akt.Y.numCols());
GeneralMatrix eta1(akt.R.numCols(), akt.Y.numCols());
GeneralMatrix V1(akt.T.numRows(), akt.T.numRows()*akt.Y.numCols());
SmootherResults sres1(akt.Y.numCols()*akt.Y.numRows());
{
WallTimer tim("\tUnivariate time ", true);
int dd;
ll1 = kut.filter_and_smooth(sres1, per1, dd);
per1 /= akt.Y.numRows();
d1 = dd/akt.Y.numRows();
printf("\t\tll=%f per=%d d=%d(%d)\n", ll1, per1, d1, dd);
if (per1 == akt.Y.numCols()) {
SmootherResults sres_uni(akt.Y.numCols());
sres_uni.import(sres1, akt.Y.numRows());
sres_uni.exportAlpha(alpha1);
sres_uni.exportEta(eta1);
sres_uni.exportV(V1);
} else {
printf("\t\tNot finished.\n");
}
}
// compare
if (per == per1 && per == akt.Y.numCols()) {
WallTimer tim("\tComparison time ", true);
alpha.add(-1.0, alpha1);
eta.add(-1.0, eta1);
V.add(-1.0, V1);
int maxd = std::max(d,d1);
for (int t = 1; t <= maxd; t++) {
Vector alphat(alpha, t-1);
printf("\t\tt=%d alpha error %10.6g\n",t,alphat.getMax());
Vector etat(eta, t-1);
printf("\t\tt=%d eta error %10.6g\n",t,etat.getMax());
GeneralMatrix Vt(V, 0, (t-1)*akt.T.numRows(), akt.T.numRows(), akt.T.numRows());
printf("\t\tt=%d V error %10.6g\n",t,V.getData().getMax());
}
GeneralMatrix alpha_rest(alpha, 0, maxd, akt.T.numRows(), alpha.numCols()-maxd);
printf("\t\tt=%d.. alpha error %10.6g\n",maxd+1,alpha_rest.getData().getMax());
GeneralMatrix eta_rest(eta, 0, maxd, akt.R.numCols(), eta.numCols()-maxd);
printf("\t\tt=%d.. eta error %10.6g\n",maxd+1,eta_rest.getData().getMax());
GeneralMatrix V_rest(V, 0, maxd*akt.T.numRows(), akt.T.numRows(),
V.numCols()-maxd*akt.T.numRows());
printf("\t\tt=%d.. V error %10.6g\n",maxd+1,V_rest.getData().getMax());
}
delete init;
return true;
}
/****************************************************/
/* definition of TestRunnable subclasses */
/****************************************************/
class SmallNonDiffuse : public TestRunnable {
public:
SmallNonDiffuse()
: TestRunnable("Non-diffuse small (p=2,m=3,r=4)") {}
bool run() const
{
filter_and_smoother("small2x3x4", false);
return true;
}
};
class SmallDiffuse : public TestRunnable {
public:
SmallDiffuse()
: TestRunnable("Diffuse small (p=2,m=3,r=4)") {}
bool run() const
{
return filter_and_smoother("small2x3x4", true);
}
};
class MiddleNonDiffuse : public TestRunnable {
public:
MiddleNonDiffuse()
: TestRunnable("Non-diffuse middle (p=10,m=15,r=12)") {}
bool run() const
{
return filter_and_smoother("10x15x12", false);
}
};
class MiddleDiffuse : public TestRunnable {
public:
MiddleDiffuse()
: TestRunnable("Diffuse middle (p=10,m=15,r=12)") {}
bool run() const
{
return filter_and_smoother("10x15x12", true);
}
};
class SOEDiffuse : public TestRunnable {
public:
SOEDiffuse()
: TestRunnable("Diffuse soe (p=8,m=25,r=15)") {}
bool run() const
{
return filter_and_smoother("soe8x25x15", true);
}
};
int main()
{
TestRunnable* all_tests[50];
// fill in vector of all tests
int num_tests = 0;
all_tests[num_tests++] = new SmallNonDiffuse();
all_tests[num_tests++] = new SmallDiffuse();
all_tests[num_tests++] = new MiddleNonDiffuse();
all_tests[num_tests++] = new MiddleDiffuse();
all_tests[num_tests++] = new SOEDiffuse();
// launch the tests
int success = 0;
for (int i = 0; i < num_tests; i++) {
try {
if (all_tests[i]->test())
success++;
} catch (const TSException& e) {
printf("Caugth TS exception in <%s>:\n", all_tests[i]->getName());
e.print();
} catch (SylvException& e) {
printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
e.printMessage();
}
}
printf("There were %d tests that failed out of %d tests run.\n",
num_tests - success, num_tests);
// destroy
for (int i = 0; i < num_tests; i++) {
delete all_tests[i];
}
return 0;
}