dynare++/kord/journal.cweb

-@q Copyright (C) 2004-2011, Ondra Kamenik @>
-
-@ Start of {\tt journal.cpp} file.
-
-@c
-#include "journal.h"
-#include "kord_exception.h"
-
-#if !defined(__MINGW32__)
-# include <sys/resource.h>
-# include <sys/utsname.h>
-#endif
-#include <cstdlib>
-#include <unistd.h>
-#include <ctime>
-
-SystemResources _sysres;
-#if defined(__MINGW32__)
-@<|sysconf| Win32 implementation@>;
-#endif
-
-#if defined(__APPLE__)
-#define _SC_PHYS_PAGES 2
-#define _SC_AVPHYS_PAGES 3
-#endif
-
-@<|SystemResources| constructor code@>;
-@<|SystemResources::pageSize| code@>;
-@<|SystemResources::physicalPages| code@>;
-@<|SystemResources::onlineProcessors| code@>;
-@<|SystemResources::availableMemory| code@>;
-@<|SystemResources::getRUS| code@>;
-@<|SystemResourcesFlash| constructor code@>;
-@<|SystemResourcesFlash::diff| code@>;
-@<|JournalRecord::operator<<| symmetry code@>;
-@<|JournalRecord::writePrefix| code@>;
-@<|JournalRecord::writePrefixForEnd| code@>;
-@<|JournalRecordPair| destructor code@>;
-@<|endrec| code@>;
-@<|Journal::printHeader| code@>;
-
-@ 
-@<|SystemResources| constructor code@>=
-SystemResources::SystemResources()
-{
-	gettimeofday(&start, NULL);
-}
-
-
-@ 
-@<|SystemResources::pageSize| code@>=
-long int SystemResources::pageSize()
-{
-	return sysconf(_SC_PAGESIZE);
-}
-
-@ 
-@<|SystemResources::physicalPages| code@>=
-long int SystemResources::physicalPages()
-{
-	return sysconf(_SC_PHYS_PAGES);
-}
-
-@ 
-@<|SystemResources::onlineProcessors| code@>=
-long int SystemResources::onlineProcessors()
-{
-	return sysconf(_SC_NPROCESSORS_ONLN);
-}
-
-@ 
-@<|SystemResources::availableMemory| code@>=
-long int SystemResources::availableMemory()
-{
-	return pageSize()*sysconf(_SC_AVPHYS_PAGES);
-}
-
-@ Here we read the current values of resource usage. For MinGW, we
-implement only a number of available physical memory pages.
-
-@<|SystemResources::getRUS| code@>=
-void SystemResources::getRUS(double& load_avg, long int& pg_avail,
-							 double& utime, double& stime, double& elapsed,
-							 long int& idrss, long int& majflt)
-{
-	struct timeval now;
-	gettimeofday(&now, NULL);
-	elapsed = now.tv_sec-start.tv_sec + (now.tv_usec-start.tv_usec)*1.0e-6;
-
-#if !defined(__MINGW32__)
-	struct rusage rus;
-	getrusage(RUSAGE_SELF, &rus);
-	utime = rus.ru_utime.tv_sec+rus.ru_utime.tv_usec*1.0e-6;
-	stime = rus.ru_stime.tv_sec+rus.ru_stime.tv_usec*1.0e-6;
-	idrss = rus.ru_idrss;
-	majflt = rus.ru_majflt;
-#else
-	utime = -1.0;
-	stime = -1.0;
-	idrss = -1;
-	majflt = -1;
-#endif
-
-#if !defined(__MINGW32__) && !defined(__CYGWIN32__)
-	getloadavg(&load_avg, 1);
-#else
-	load_avg = -1.0;
-#endif
-
-	pg_avail = sysconf(_SC_AVPHYS_PAGES);
-}
-
-@ 
-@<|SystemResourcesFlash| constructor code@>=
-SystemResourcesFlash::SystemResourcesFlash()
-{
-	_sysres.getRUS(load_avg, pg_avail, utime, stime,
-				   elapsed, idrss, majflt);
-}
-
-@ 
-@<|SystemResourcesFlash::diff| code@>=
-void SystemResourcesFlash::diff(const SystemResourcesFlash& pre)
-{
-	utime -= pre.utime;
-	stime -= pre.stime;
-	elapsed -= pre.elapsed;
-	idrss -= pre.idrss;
-	majflt -= pre.majflt;
-}
-
-@ 
-@<|JournalRecord::operator<<| symmetry code@>=
-JournalRecord& JournalRecord::operator<<(const IntSequence& s)
-{
-	operator<<("[");
-	for (int i = 0; i < s.size(); i++) {
-		operator<<(s[i]);
-		if (i < s.size()-1)
-			operator<<(",");
-	}
-	operator<<("]");
-	return *this;
-}
-
-@ 
-@<|JournalRecord::writePrefix| code@>=
-void JournalRecord::writePrefix(const SystemResourcesFlash& f)
-{
-	for (int i = 0; i < MAXLEN; i++)
-		prefix[i] = ' ';
-	double mb = 1024*1024;
-	sprintf(prefix, "%07.6g", f.elapsed);
-	sprintf(prefix+7, ":%c%05d", recChar, ord);
-	sprintf(prefix+14, ":%1.1f", f.load_avg);
-	sprintf(prefix+18, ":%05.4g", f.pg_avail*_sysres.pageSize()/mb);
-	sprintf(prefix+24, "%s", ":      : ");
-	for (int i = 0; i < 2*journal.getDepth(); i++)
-		prefix[i+33]=' ';
-	prefix[2*journal.getDepth()+33]='\0';
-}
-
-@ 
-@<|JournalRecord::writePrefixForEnd| code@>=
-void JournalRecordPair::writePrefixForEnd(const SystemResourcesFlash& f)
-{
-	for (int i = 0; i < MAXLEN; i++)
-		prefix_end[i] = ' ';
-	double mb = 1024*1024;
-	SystemResourcesFlash difnow;
-	difnow.diff(f);
-	sprintf(prefix_end, "%07.6g", f.elapsed+difnow.elapsed);
-	sprintf(prefix_end+7, ":E%05d", ord);
-	sprintf(prefix_end+14, ":%1.1f", difnow.load_avg);
-	sprintf(prefix_end+18, ":%05.4g", difnow.pg_avail*_sysres.pageSize()/mb);
-	sprintf(prefix_end+24, ":%06.5g", difnow.majflt*_sysres.pageSize()/mb);
-	sprintf(prefix_end+31, "%s", ": ");
-	for (int i = 0; i < 2*journal.getDepth(); i++)
-		prefix_end[i+33]=' ';
-	prefix_end[2*journal.getDepth()+33]='\0';
-}
-
-@ 
-@<|JournalRecordPair| destructor code@>=
-JournalRecordPair::~JournalRecordPair()
-{
-	journal.decrementDepth();
-	writePrefixForEnd(flash);
-	journal << prefix_end;
-	journal << mes;
-	journal << endl;
-	journal.flush();
-}
-
-@ 
-@<|endrec| code@>=
-JournalRecord& endrec(JournalRecord& rec)
-{
-	rec.journal << rec.prefix;
-	rec.journal << rec.mes;
-	rec.journal << endl;
-	rec.journal.flush();
-	rec.journal.incrementOrd();
-	return rec;
-}
-
-@ 
-@<|Journal::printHeader| code@>=
-void Journal::printHeader()
-{
-	(*this)<< "This is Dynare++, Copyright (C) 2004-2011, Ondra Kamenik\n"
-	       << "Dynare++ comes with ABSOLUTELY NO WARRANTY and is distributed under\n"
-		     << "GPL: modules integ, tl, kord, sylv, src, extern and documentation\n"
-	       << "LGPL: modules parser, utils\n"
-	       << " for GPL  see http://www.gnu.org/licenses/gpl.html\n"
-	       << " for LGPL see http://www.gnu.org/licenses/lgpl.html\n"
-	       << "\n\n";
-
-#if !defined(__MINGW32__)
-	utsname info;
-	uname(&info);
-	(*this)<< "System info: ";
-	(*this)<< info.sysname << " " << info.release << " " << info.version << " ";
-	(*this)<< info.machine << ", processors online: " << _sysres.onlineProcessors();
-
-	(*this)<< "\n\nStart time: ";
-	char ts[100];
-	time_t curtime = time(NULL);
-	tm loctime;
-	localtime_r(&curtime, &loctime);
-	asctime_r(&loctime, ts);
-	(*this)<< ts << "\n";
-#else
-	(*this) << "System info: (not implemented for MINGW)\n";
-	(*this) << "Start time:  (not implemented for MINGW)\n\n";
-#endif
-	
-	(*this)<< "  ------ elapsed time (seconds)                     \n";
-	(*this)<< "  |       ------ record unique identifier           \n";
-	(*this)<< "  |       |     ------ load average                 \n";
-	(*this)<< "  |       |     |    ------ available memory (MB)   \n";
-	(*this)<< "  |       |     |    |     ------  major faults (MB)\n";
-	(*this)<< "  |       |     |    |     |                        \n";
-	(*this)<< "  V       V     V    V     V                        \n";
-	(*this)<< "\n";
-}
-
-
-@ Here we implement |sysconf| for MinGW. We implement only page size,
-number of physial pages, and a number of available physical pages. The
-pagesize is set to 1024 bytes, real pagesize can differ but it is not
-important. We can do this since Windows kernel32 |GlobalMemoryStatus|
-call returns number of bytes.
-
-Number of online processors is not implemented and returns -1, since
-Windows kernel32 |GetSystemInfo| call is too complicated.
-
-@<|sysconf| Win32 implementation@>=
-#ifndef NOMINMAX
-# define NOMINMAX // Do not define "min" and "max" macros
-#endif
-#include <windows.h>
-
-#define _SC_PAGESIZE 1
-#define _SC_PHYS_PAGES 2
-#define _SC_AVPHYS_PAGES 3
-#define _SC_NPROCESSORS_ONLN 4
-@#
-long sysconf(int name)
-{
-	switch (name) {
-	case _SC_PAGESIZE:@;
-		return 1024;
-	case _SC_PHYS_PAGES:@;
-		{
-			MEMORYSTATUS memstat;
-			GlobalMemoryStatus(&memstat);
-			return memstat.dwTotalPhys/1024;
-		}
-	case _SC_AVPHYS_PAGES:@;
-		{
-			MEMORYSTATUS memstat;
-			GlobalMemoryStatus(&memstat);
-			return memstat.dwAvailPhys/1024;
-		}
-	case _SC_NPROCESSORS_ONLN:@;
-		return -1;
-	default:@; 
-		KORD_RAISE("Not implemented in Win32 sysconf.");
-		return -1;
-	}
-}
-
-@ End of {\tt journal.cpp} file.
\ No newline at end of file

dynare++/kord/journal.hweb

-@q $Id: journal.hweb 417 2005-08-16 15:04:24Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@*2 Resource usage journal. Start of {\tt journal.h} file.
-
-@s timeval int
-@s rusage int
-@s SystemResources int
-@s SystemResourcesFlash int
-@s Journal int
-@s JournalRecord int
-@s JournalRecordPair int
-
-@c
-#ifndef JOURNAL_H
-#define JOURNAL_H
-
-#include "int_sequence.h"
-
-#include <sys/time.h>
-#include <cstdio>
-#include <iostream>
-#include <fstream>
-
-@<|SystemResources| class declaration@>;
-@<|SystemResourcesFlash| struct declaration@>;
-@<|Journal| class declaration@>;
-@<|JournalRecord| class declaration@>;
-@<|JournalRecordPair| class declaration@>;
-
-#endif
-
-@ 
-@<|SystemResources| class declaration@>=
-class SystemResources {
-	timeval start;
-public:@;
-	SystemResources();
-	static long int pageSize();
-	static long int physicalPages();
-	static long int onlineProcessors();
-	static long int availableMemory();
-	void getRUS(double& load_avg, long int& pg_avail, double& utime,
-				double& stime, double& elapsed, long int& idrss,
-				long int& majflt);
-};
-
-@ 
-@<|SystemResourcesFlash| struct declaration@>=
-struct SystemResourcesFlash {
-	double load_avg;
-	long int pg_avail;
-	double utime;
-	double stime;
-	double elapsed;
-	long int idrss;
-	long int majflt;
-	SystemResourcesFlash();
-	void diff(const SystemResourcesFlash& pre);
-};
-
-
-@ 
-@s stringstream int
-@d MAXLEN 1000
-
-@<|JournalRecord| class declaration@>=
-class JournalRecord;
-JournalRecord& endrec(JournalRecord&);
-
-class JournalRecord {
-protected:@;
-	char recChar;
-	int ord;
-public:@;
-	Journal& journal;
-	char prefix[MAXLEN];
-	char mes[MAXLEN];
-	SystemResourcesFlash flash;
-	typedef JournalRecord& (*_Tfunc)(JournalRecord&);
-
-	JournalRecord(Journal& jr, char rc = 'M')
-		: recChar(rc), ord(jr.getOrd()), journal(jr) 
-		{@+ prefix[0]='\0';mes[0]='\0';writePrefix(flash); @+}
-	virtual ~JournalRecord() @+{}
-	JournalRecord& operator<<(const IntSequence& s);
-	JournalRecord& operator<<(_Tfunc f)
-		{@+ (*f)(*this); return *this;@+}
-	JournalRecord& operator<<(const char* s)
-		{@+ strcat(mes, s); return *this; @+}
-	JournalRecord& operator<<(int i)
-		{@+ sprintf(mes+strlen(mes), "%d", i); return *this;@+}
-	JournalRecord& operator<<(double d)
-		{@+ sprintf(mes+strlen(mes), "%f", d); return *this;@+}
-protected:@;
-	void writePrefix(const SystemResourcesFlash& f);
-};
-
-@ 
-@<|JournalRecordPair| class declaration@>=
-class JournalRecordPair : public JournalRecord {
-	char prefix_end[MAXLEN];
-public:@;
-	JournalRecordPair(Journal& jr)
-		: JournalRecord(jr, 'S')
-		{@+ prefix_end[0] = '\0'; journal.incrementDepth(); @+}
-	~JournalRecordPair();
-private:@;
-	void writePrefixForEnd(const SystemResourcesFlash& f);
-};
-
-@ 
-@<|Journal| class declaration@>=
-class Journal : public ofstream {
-	int ord;
-	int depth;
-public:@;
-	Journal(const char* fname)
-		: ofstream(fname), ord(0), depth(0)
-		{@+ printHeader();@+}
-	~Journal()
-		{@+ flush();@+}
-	void printHeader();
-	void incrementOrd()
-		{@+ ord++; @+}
-	int getOrd() const
-		{@+ return ord; @+}
-	void incrementDepth()
-		{@+ depth++; @+}
-	void decrementDepth()
-		{@+ depth--; @+}
-	int getDepth() const
-		{return depth;}
-};
-
-
-@ End of {\tt journal.h} file.
\ No newline at end of file

dynare++/kord/kord_exception.hweb

-@q $Id: kord_exception.hweb 1452 2007-11-21 11:33:30Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Exception. Start of {\tt kord\_exception.h} file.
-
-This is a simple code defining an exception and two convenience macros.
-
-@s KordException int
-@c
-#ifndef KORD_EXCEPTION_H
-#define KORD_EXCEPTION_H
-
-#include <cstring>
-#include <cstdio>
-
-#define KORD_RAISE(mes) \
-throw KordException(__FILE__, __LINE__, mes);
-
-#define KORD_RAISE_IF(expr, mes) \
-if (expr) throw KordException(__FILE__, __LINE__, mes);
-
-#define KORD_RAISE_X(mes, c) \
-throw KordException(__FILE__, __LINE__, mes, c);
-
-#define KORD_RAISE_IF_X(expr, mes, c) \
-if (expr) throw KordException(__FILE__, __LINE__, mes, c);
-
-@<|KordException| class definition@>;
-@<|KordException| error code definitions@>;
- 
-#endif
-
-@ 
-@<|KordException| class definition@>=
-class KordException {
-protected:@;
-	char fname[50];
-	int lnum;
-	char message[500];
-	int cd;
-public:@;
-	KordException(const char* f, int l, const char* mes, int c=255)
-		{
-			strncpy(fname, f, 50);@+ fname[49] = '\0';
-			strncpy(message, mes, 500);@+ message[499] = '\0';
-			lnum = l;
-			cd = c;
-		}
-	virtual ~KordException()@+ {}
-	virtual void print() const
-		{@+ printf("At %s:%d:(%d):%s\n", fname, lnum, cd, message);@+}
-	virtual int code() const
-		{@+ return cd; @+}
-	const char* get_message() const
-		{@+ return message; @+}
-};
-
-@ 
-@<|KordException| error code definitions@>=
-#define KORD_FP_NOT_CONV 254
-#define KORD_FP_NOT_FINITE 253
-#define KORD_MD_NOT_STABLE 252
-
-@ End of {\tt kord\_exception.h} file.

dynare++/kord/korder.cweb

-@q $Id: korder.cweb 1831 2008-05-18 20:13:42Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@ Start of {\tt korder.cpp} file.
-
-@c
-
-#include "kord_exception.h"
-#include "korder.h"
-
-@<|PLUMatrix| copy constructor@>;
-@<|PLUMatrix::calcPLU| code@>;
-@<|PLUMatrix::multInv| code@>;
-@<|MatrixA| constructor code@>;
-@<|MatrixS| constructor code@>;
-@<|KOrder| member access method specializations@>;
-@<|KOrder::sylvesterSolve| unfolded specialization@>;
-@<|KOrder::sylvesterSolve| folded specialization@>;
-@<|KOrder::switchToFolded| code@>;
-@<|KOrder| constructor code@>;
-
-@ 
-@<|PLUMatrix| copy constructor@>=
-PLUMatrix::PLUMatrix(const PLUMatrix& plu)
-	: TwoDMatrix(plu), inv(plu.inv), ipiv(new lapack_int[nrows()])
-{
-	memcpy(ipiv, plu.ipiv, nrows()*sizeof(lapack_int));
-}
-
-
-@ Here we set |ipiv| and |inv| members of the |PLUMatrix| depending on
-its content. It is assumed that subclasses will call this method at
-the end of their constructors.
-
-@<|PLUMatrix::calcPLU| code@>=
-void PLUMatrix::calcPLU()
-{
-	lapack_int info;
-	lapack_int rows = nrows();
-	inv = (const Vector&)getData();
-	dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
-}
-
-@ Here we just call the LAPACK machinery to multiply by the inverse.
-
-@<|PLUMatrix::multInv| code@>=
-void PLUMatrix::multInv(TwoDMatrix& m) const
-{
-	KORD_RAISE_IF(m.nrows() != ncols(),
-				  "The matrix is not square in PLUMatrix::multInv");
-	lapack_int info;
-	lapack_int mcols = m.ncols();
-	lapack_int mrows = m.nrows();
-	double* mbase = m.getData().base();
-	dgetrs("N", &mrows, &mcols, inv.base(), &mrows, ipiv,
-				  mbase, &mrows, &info);
-	KORD_RAISE_IF(info != 0,
-				  "Info!=0 in PLUMatrix::multInv");
-}
-
-@ Here we construct the matrix $A$. Its dimension is |ny|, and it is
-$$A=\left[f_{y}\right]+
-\left[0 \left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right] 0\right]$$,
-where the first zero spans |nstat| columns, and last zero spans
-|nforw| columns.
-
-@<|MatrixA| constructor code@>=
-MatrixA::MatrixA(const FSSparseTensor& f, const IntSequence& ss,
-				 const TwoDMatrix& gy, const PartitionY& ypart)
-	: PLUMatrix(ypart.ny())
-{
-	zeros();
-
-	IntSequence c(1); c[0] = 1;
-	FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
-	add(1.0, f_y);
-
-	ConstTwoDMatrix gss_ys(ypart.nstat+ypart.npred, ypart.nyss(), gy);
-	c[0] = 0;
-	FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
-	TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
-	sub.multAndAdd(ConstTwoDMatrix(f_yss), gss_ys);
-
-	calcPLU();
-}
-
-@ Here we construct the matrix $S$. Its dimension is |ny|, and it is
-$$S=\left[f_{y}\right]+
-\left[0\quad\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right]\quad
-0\right]+ \left[0\quad 0\quad\left[f_{y^{**}_+}\right]\right]$$
-It is, in fact, the matrix $A$ plus the third summand. The first zero
-in the summand spans |nstat| columns, the second zero spans |npred|
-columns.
-
-@<|MatrixS| constructor code@>=
-MatrixS::MatrixS(const FSSparseTensor& f, const IntSequence& ss,
-				 const TwoDMatrix& gy, const PartitionY& ypart)
-	: PLUMatrix(ypart.ny())
-{
-	zeros();
-
-	IntSequence c(1); c[0] = 1;
-	FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
-	add(1.0, f_y);
-
-	ConstTwoDMatrix gss_ys(ypart.nstat+ypart.npred, ypart.nyss(), gy);
-	c[0] = 0;
-	FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
-	TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
-	sub.multAndAdd(ConstTwoDMatrix(f_yss), gss_ys);
-
-	TwoDMatrix sub2(*this, ypart.nstat+ypart.npred, ypart.nyss());
-	sub2.add(1.0, f_yss);
-
-	calcPLU();
-}
-
-
-@ Here is the constructor of the |KOrder| class. We pass what we have
-to. The partitioning of the $y$ vector, a sparse container with model
-derivatives, then the first order approximation, these are $g_y$ and
-$g_u$ matrices, and covariance matrix of exogenous shocks |v|.
-
-We build the members, it is nothing difficult. Note that we do not make
-a physical copy of sparse tensors, so during running the class, the
-outer world must not change them.
-
-In the body, we have to set |nvs| array, and initialize $g$ and $G$
-containers to comply to preconditions of |performStep|.
-
-@<|KOrder| constructor code@>=
-KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
-			   const TensorContainer<FSSparseTensor>& fcont,
-			   const TwoDMatrix& gy, const TwoDMatrix& gu, const TwoDMatrix& v,
-			   Journal& jr)
-	: ypart(num_stat, num_pred, num_both, num_forw),@/
-	  ny(ypart.ny()), nu(gu.ncols()), maxk(fcont.getMaxDim()),@/
-	  nvs(4),@/
-	  _ug(4), _fg(4), _ugs(4), _fgs(4), _ugss(4), _fgss(4), @/
-	  _uG(4), _fG(4),@/
-	  _uZstack(&_uG, ypart.nyss(), &_ug, ny, ypart.nys(), nu),@/
-	  _fZstack(&_fG, ypart.nyss(), &_fg, ny, ypart.nys(), nu),@/
-	  _uGstack(&_ugs, ypart.nys(), nu),@/
-	  _fGstack(&_fgs, ypart.nys(), nu),@/
-	  _um(maxk, v), _fm(_um), f(fcont),@/
-	  matA(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),@/
-	  matS(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),@/
-	  matB(*(f.get(Symmetry(1))), _uZstack.getStackSizes()),@/
-	  journal(jr)@/
-{
-	KORD_RAISE_IF(gy.ncols() != ypart.nys(),
-				  "Wrong number of columns in gy in KOrder constructor");
-	KORD_RAISE_IF(v.ncols() != nu,
-				  "Wrong number of columns of Vcov in KOrder constructor");
-	KORD_RAISE_IF(nu != v.nrows(),
-				  "Wrong number of rows of Vcov in KOrder constructor");
-	KORD_RAISE_IF(maxk < 2,
-				  "Order of approximation must be at least 2 in KOrder constructor");
-	KORD_RAISE_IF(gy.nrows() != ypart.ny(),
-				  "Wrong number of rows in gy in KOrder constructor");
-	KORD_RAISE_IF(gu.nrows() != ypart.ny(),
-				  "Wrong number of rows in gu in KOrder constuctor");
-	KORD_RAISE_IF(gu.ncols() != nu,
-				  "Wrong number of columns in gu in KOrder constuctor");
-
-	// set nvs:
-	nvs[0] = ypart.nys(); nvs[1] = nu; nvs[2] = nu; nvs[3] = 1;
-
-	@<put $g_y$ and $g_u$ to the container@>;
-	@<put $G_y$, $G_u$ and $G_{u'}$ to the container@>;@q'@>
-}
-
-@ Note that $g_\sigma$ is zero by the nature and we do not insert it to
-the container. We insert a new physical copies.
-
-@<put $g_y$ and $g_u$ to the container@>=
-	UGSTensor* tgy = new UGSTensor(ny, TensorDimens(Symmetry(1,0,0,0), nvs));
-	tgy->getData() = gy.getData();
-	insertDerivative<unfold>(tgy);
-	UGSTensor* tgu = new UGSTensor(ny, TensorDimens(Symmetry(0,1,0,0), nvs));
-	tgu->getData() = gu.getData();
-	insertDerivative<unfold>(tgu);
-
-@ Also note that since $g_\sigma$ is zero, so $G_\sigma$.
-@<put $G_y$, $G_u$ and $G_{u'}$ to the container@>=
-	UGSTensor* tGy = faaDiBrunoG<unfold>(Symmetry(1,0,0,0));
-	G<unfold>().insert(tGy);
-	UGSTensor* tGu = faaDiBrunoG<unfold>(Symmetry(0,1,0,0));
-	G<unfold>().insert(tGu);
-	UGSTensor* tGup = faaDiBrunoG<unfold>(Symmetry(0,0,1,0));
-	G<unfold>().insert(tGup);
-
-
-
-@ Here we have an unfolded specialization of |sylvesterSolve|. We
-simply create the sylvester object and solve it. Note that the $g^*_y$
-is not continuous in memory as assumed by the sylvester code, so we
-make a temporary copy and pass it as matrix $C$.
-
-If the $B$ matrix is empty, in other words there are now forward
-looking variables, then the system becomes $AX=D$ which is solved by
-simple |matA.multInv()|.
-
-If one wants to display the diagnostic messages from the Sylvester
-module, then after the |sylv.solve()| one needs to call
-|sylv.getParams().print("")|.
-
-
-@<|KOrder::sylvesterSolve| unfolded specialization@>=
-template<>@/
-void KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor& der) const
-{
-	JournalRecordPair pa(journal);
-	pa << "Sylvester equation for dimension = " << der.getSym()[0] << endrec;
-	if (ypart.nys() > 0 && ypart.nyss() > 0) {
-		KORD_RAISE_IF(! der.isFinite(),
-					  "RHS of Sylverster is not finite");
-		TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry(1,0,0,0))));
-		GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
-							  ypart.nstat+ypart.npred,
-							  matA.getData().base(), matB.getData().base(),
-							  gs_y.getData().base(), der.getData().base());
-		sylv.solve();
-	} else if (ypart.nys() > 0 && ypart.nyss() == 0) {
-		matA.multInv(der);
-	}
-}
-
-@ Here is the folded specialization of sylvester. We unfold the right
-hand side. Then we solve it by the unfolded version of
-|sylvesterSolve|, and fold it back and copy to output vector.
-
-@<|KOrder::sylvesterSolve| folded specialization@>=
-template<>@/
-void KOrder::sylvesterSolve<KOrder::fold>(ctraits<fold>::Ttensor& der) const
-{
-	ctraits<unfold>::Ttensor tmp(der);
-	sylvesterSolve<unfold>(tmp);
-	ctraits<fold>::Ttensor ftmp(tmp);
-	der.getData() = (const Vector&)(ftmp.getData());
-}
-
-@ 
-@<|KOrder::switchToFolded| code@>=
-void KOrder::switchToFolded()
-{
-	JournalRecordPair pa(journal);
-	pa << "Switching from unfolded to folded" << endrec;
-
-	int maxdim = g<unfold>().getMaxDim();
-	for (int dim = 1; dim <= maxdim; dim++) {
-		SymmetrySet ss(dim, 4);
-		for (symiterator si(ss); !si.isEnd(); ++si) {
-			if ((*si)[2] == 0 && g<unfold>().check(*si)) {
-				FGSTensor* ft = new FGSTensor(*(g<unfold>().get(*si)));
-				insertDerivative<fold>(ft);
-				if (dim > 1) {
-					gss<unfold>().remove(*si);
-					gs<unfold>().remove(*si);
-					g<unfold>().remove(*si);
-				}
-			}
-			if (G<unfold>().check(*si)) {
-				FGSTensor* ft = new FGSTensor(*(G<unfold>().get(*si)));
-				G<fold>().insert(ft);
-				if (dim > 1) {
-					G<fold>().remove(*si);
-				}
-			}
-		}
-	}
-}
-
-
-
-@ These are the specializations of container access methods. Nothing
-interesting here.
-
-@<|KOrder| member access method specializations@>=
-	template<> ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>()
-		{@+ return _ug;@+}
-	template<>@; const ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>()@+const@;
-		{@+ return _ug;@+}
-	template<> ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>()
-		{@+ return _fg;@+}
-	template<> const ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>()@+const@;
-		{@+ return _fg;@+}
-	template<> ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>()
-		{@+ return _ugs;@+}
-	template<> const ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>()@+const@;
-		{@+ return _ugs;@+}
-	template<> ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>()
-		{@+ return _fgs;@+}
-	template<> const ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>()@+const@;
-		{@+ return _fgs;@+}
-	template<> ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>()
-		{@+ return _ugss;@+}
-	template<> const ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>()@+const@;
-		{@+ return _ugss;@+}
-	template<> ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>()
-		{@+ return _fgss;@+}
-	template<> const ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>()@+const@;
-		{@+ return _fgss;@+}
-	template<> ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>()
-		{@+ return _uG;@+}
-	template<> const ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>()@+const@;
-		{@+ return _uG;@+}
-	template<> ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>()
-		{@+ return _fG;@+}
-	template<> const ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>()@+const@;
-		{@+ return _fG;@+}
-	template<> ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>()
-		{@+ return _uZstack;@+}
-	template<> const ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>()@+const@;
-		{@+ return _uZstack;@+}
-	template<> ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>()
-		{@+ return _fZstack;@+}
-	template<> const ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>()@+const@;
-		{@+ return _fZstack;@+}
-	template<> ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>()
-		{@+ return _uGstack;@+}
-	template<> const ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>()@+const@;
-		{@+ return _uGstack;@+}
-	template<> ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>()
-		{@+ return _fGstack;@+}
-	template<> const ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>()@+const@;
-		{@+ return _fGstack;@+}
-	template<> ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>()
-		{@+ return _um;@+}
-	template<> const ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>()@+const@;
-		{@+ return _um;@+}
-	template<> ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>()
-		{@+ return _fm;@+}
-	template<> const ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>()@+const@;
-		{@+ return _fm;@+}
-
-
-@ End of {\tt korder.cpp} file.

dynare++/kord/korder.hweb

-@q $Id: korder.hweb 2332 2009-01-14 10:26:54Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@*2 Higher order at deterministic steady. Start of {\tt korder.h} file.
-
-The main purpose of this file is to implement a perturbation method
-algorithm for an SDGE model for higher order approximations. The input
-of the algorithm are sparse tensors as derivatives of the dynamic
-system, then dimensions of vector variables, then the first order
-approximation to the decision rule and finally a covariance matrix of
-exogenous shocks. The output are higher order derivatives of decision
-rule $y_t=g(y^*_{t-1},u_t,\sigma)$. The class provides also a method
-for checking a size of residuals of the solved equations.
-
-The algorithm is implemented in |KOrder| class. The class contains
-both unfolded and folded containers to allow for switching (usually
-from unfold to fold) during the calculations. The algorithm is
-implemented in a few templated methods. To do this, we need some
-container type traits, which are in |ctraits| struct. Also, the
-|KOrder| class contains some information encapsulated in other
-classes, which are defined here. These include: |PartitionY|,
-|MatrixA|, |MatrixS| and |MatrixB|.
-
-@s KOrder int
-@s ctraits int
-@s PartitionY int
-@s MatrixA int
-@s MatrixB int
-@s MatrixS int
-@s PLUMatrix int
-@s FGSTensor int
-@s UGSTensor int
-@s FGSContainer int
-@s UGSContainer int
-@s FSSparseTensor int
-@s TensorContainer int
-@s UNormalMoments int
-@s FNormalMoments int
-@s FoldedZContainer int
-@s UnfoldedZContainer int
-@s FoldedGContainer int
-@s UnfoldedGContainer int
-@s FoldedZXContainer int
-@s UnfoldedZXContainer int
-@s FoldedGXContainer int
-@s UnfoldedGXContainer int
-@s TwoDMatrix int
-@s ConstTwoDMatrix int
-@s IntSequence int
-@s Symmetry int
-@s SymmetrySet int
-@s symiterator int
-@s TensorDimens int
-@s Vector int
-@s ConstVector int
-@s UTensorPolynomial int
-@s FTensorPolynomial int
-@s UFSTensor int
-@s FFSTensor int
-@s GeneralSylvester int
-
-@c
-#ifndef KORDER_H
-#define KORDER_H
-
-#include "int_sequence.h"
-#include "fs_tensor.h"
-#include "gs_tensor.h"
-#include "t_container.h"
-#include "stack_container.h"
-#include "normal_moments.h"
-#include "t_polynomial.h"
-#include "faa_di_bruno.h"
-#include "journal.h"
-
-#include "kord_exception.h"
-#include "GeneralSylvester.h"
-
-#include <dynlapack.h>
-
-#include <cmath>
-
-#define TYPENAME typename
-
-@<|ctraits| type traits declaration@>;
-@<|PartitionY| struct declaration@>;
-@<|PLUMatrix| class declaration@>;
-@<|MatrixA| class declaration@>;
-@<|MatrixS| class declaration@>;
-@<|MatrixB| class declaration@>;
-@<|KOrder| class declaration@>;
-
-
-#endif
-
-@ Here we use a classical IF template, and in |ctraits| we define a
-number of types. We have a type for tensor |Ttensor|, and types for
-each pair of folded/unfolded containers used as a member in |KOrder|.
-
-Note that we have enumeration |fold| and |unfold|. These must have the
-same value as the same enumeration in |KOrder|.
-
-@s IF int
-@s Then int
-@s Else int
-@s RET int
-
-@<|ctraits| type traits declaration@>=
-class FoldedZXContainer;
-class UnfoldedZXContainer;
-class FoldedGXContainer;
-class UnfoldedGXContainer;
-
-template<bool condition, class Then, class Else>
-struct IF {
-	typedef Then RET;
-};
-
-template<class Then, class Else>
-struct IF<false, Then, Else> {
-	typedef Else RET;
-};
-
-template <int type>
-class ctraits {
-public:@;
-	enum {@+ fold, unfold@+};
-	typedef TYPENAME IF<type==fold, FGSTensor, UGSTensor>::RET Ttensor;
-	typedef TYPENAME IF<type==fold, FFSTensor, UFSTensor>::RET Ttensym;
-	typedef TYPENAME IF<type==fold, FGSContainer, UGSContainer>::RET Tg;
-	typedef TYPENAME IF<type==fold, FGSContainer, UGSContainer>::RET Tgs;
-	typedef TYPENAME IF<type==fold, FGSContainer, UGSContainer>::RET Tgss;
-	typedef TYPENAME IF<type==fold, FGSContainer, UGSContainer>::RET TG;
-	typedef TYPENAME IF<type==fold, FoldedZContainer, UnfoldedZContainer>::RET TZstack;
-	typedef TYPENAME IF<type==fold, FoldedGContainer, UnfoldedGContainer>::RET TGstack;
-	typedef TYPENAME IF<type==fold, FNormalMoments, UNormalMoments>::RET Tm;
-	typedef TYPENAME IF<type==fold, FTensorPolynomial, UTensorPolynomial>::RET Tpol;
-	typedef TYPENAME IF<type==fold, FoldedZXContainer, UnfoldedZXContainer>::RET TZXstack;
-	typedef TYPENAME IF<type==fold, FoldedGXContainer, UnfoldedGXContainer>::RET TGXstack;
-};
-
-
-@ The |PartitionY| class defines the partitioning of state variables
-$y$. The vector $y$, and subvector $y^*$, and $y^{**}$ are defined.
-$$y=\left[\matrix{\hbox{static}\cr\hbox{predeter}\cr\hbox{both}\cr
-				  \hbox{forward}}\right],\quad
-  y^*=\left[\matrix{\hbox{predeter}\cr\hbox{both}}\right],\quad
-  y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right],$$
-where ``static'' means variables appearing only at time $t$,
-``predeter'' means variables appearing at time $t-1$, but not at
-$t+1$, ``both'' means variables appearing both at $t-1$ and $t+1$
-(regardless appearance at $t$), and ``forward'' means variables
-appearing at $t+1$, but not at $t-1$.
-
-The class maintains the four lengths, and returns the whole length,
-length of $y^s$, and length of $y^{**}$.
-
-@<|PartitionY| struct declaration@>=
-struct PartitionY {
-	const int nstat;
-	const int npred;
-	const int nboth;
-	const int nforw;
-	PartitionY(int num_stat, int num_pred,
-			   int num_both, int num_forw)
-		: nstat(num_stat), npred(num_pred),
-		  nboth(num_both), nforw(num_forw)
-		{}
-	int ny() const
-		{@+ return nstat+npred+nboth+nforw;@+}
-	int nys() const
-		{@+ return npred+nboth;@+}
-	int nyss() const
-		{@+ return nboth+nforw;@+}
-};
-
-
-@ This is an abstraction for a square matrix with attached PLU
-factorization. It can calculate the PLU factorization and apply the
-inverse with some given matrix.
-
-We use LAPACK $PLU$ decomposition for the inverse. We store the $L$
-and $U$ in the |inv| array and |ipiv| is the permutation $P$.
-
-@<|PLUMatrix| class declaration@>=
-class PLUMatrix : public TwoDMatrix {
-public:@;
-	PLUMatrix(int n)
-		: TwoDMatrix(n,n),
-		  inv(nrows()*ncols()),
-		  ipiv(new lapack_int[nrows()]) {}
-	PLUMatrix(const PLUMatrix& plu);
-	virtual ~PLUMatrix()
-		{delete [] ipiv;}
-	void multInv(TwoDMatrix& m) const;
-private:@;
-	Vector inv;
-	lapack_int* ipiv;
-protected:@;
-	void calcPLU();
-};
-
-@ The class |MatrixA| is used for matrix $\left[f_{y}\right]+ \left[0
-\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right] 0\right]$,
-which is central for the perturbation method step.
-
-@<|MatrixA| class declaration@>=
-class MatrixA : public PLUMatrix {
-public:@;
-	MatrixA(const FSSparseTensor& f, const IntSequence& ss,
-			const TwoDMatrix& gy, const PartitionY& ypart);
-};
-
-@ The class |MatrixS| slightly differs from |MatrixA|. It is used for
-matrix $$\left[f_{y}\right]+ \left[0
-\quad\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right]\quad
-0\right]+\left[0\quad 0\quad\left[f_{y^{**}_+}\right]\right]$$, which is
-needed when recovering $g_{\sigma^k}$.
-
-@<|MatrixS| class declaration@>=
-class MatrixS : public PLUMatrix {
-public:@;
-	MatrixS(const FSSparseTensor& f, const IntSequence& ss,
-			const TwoDMatrix& gy, const PartitionY& ypart);
-};
-
-
-@ The $B$ matrix is equal to $\left[f_{y^{**}_+}\right]$. We have just
-a constructor.
-
-@<|MatrixB| class declaration@>=
-class MatrixB : public TwoDMatrix {
-public:@;
-	MatrixB(const FSSparseTensor& f, const IntSequence& ss)
-		: TwoDMatrix(FGSTensor(f, ss, IntSequence(1,0),
-							   TensorDimens(ss, IntSequence(1,0))))
-		{}
-};
-
-@ Here we have the class for the higher order approximations. It
-contains the following data:
-
-\halign{\kern\parindent\vrule height12pt width0pt
-\vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
-variable sizes ypart& |PartitionY| struct maintaining partitions of
-$y$, see |@<|PartitionY| struct declaration@>|\cr
-tensor variable dimension |nvs|& variable sizes of all tensors in
-containers, sizes of $y^*$, $u$, $u'$@q'@> and $\sigma$\cr
-tensor containers & folded and unfolded containers for $g$, $g_{y^*}$,
-$g_{y^**}$ (the latter two collect appropriate subtensors of $g$, they
-do not allocate any new space), $G$, $G$ stack, $Z$ stack\cr
-dynamic model derivatives & just a reference to the container of
-sparse tensors of the system derivatives, lives outside the class\cr
-moments & both folded and unfolded normal moment containers, both are
-calculated at initialization\cr
-matrices & matrix $A$, matrix $S$, and matrix $B$, see |@<|MatrixA| class
- declaration@>| and |@<|MatrixB| class declaration@>|\cr
-}
-
-\kern 0.4cm
-
-The methods are the following:
-\halign{\kern\parindent\vrule height12pt width0pt
-\vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
-member access & we declare template methods for accessing containers
-depending on |fold| and |unfold| flag, we implement their
-specializations\cr
-|performStep| & this performs $k$-order step provided that $k=2$ or
- the $k-1$-th step has been run, this is the core method\cr
-|check| & this calculates residuals of all solved equations for
- $k$-order and reports their sizes, it is runnable after $k$-order
- |performStep| has been run\cr
-|insertDerivative| & inserts a $g$ derivative to the $g$ container and
- also creates subtensors and insert them to $g_{y^*}$ and $g_{y^{**}}$
- containers\cr
-|sylvesterSolve| & solve the sylvester equation (templated fold, and
- unfold)\cr
-|faaDiBrunoZ| & calculates derivatives of $F$ by Faa Di Bruno for the
-sparse container of system derivatives and $Z$ stack container\cr
-|faaDiBrunoG| & calculates derivatives of $G$ by Faa Di Bruno for the
- dense container $g^{**}$ and $G$ stack\cr
-|recover_y| & recovers $g_{y^{*i}}$\cr
-|recover_yu| & recovers $g_{y^{*i}u^j}$\cr
-|recover_ys| & recovers $g_{y^{*i}\sigma^j}$\cr
-|recover_yus| & recovers $g_{y^{*i}u^j\sigma^k}$\cr
-|recover_s| & recovers $g_{\sigma^i}$\cr
-|fillG| & calculates specified derivatives of $G$ and inserts them to
-the container\cr
-|calcE_ijk|& calculates $E_{ijk}$\cr
-|calcD_ijk|& calculates $D_{ijk}$\cr
- }
-
-\kern 0.3cm
-
-Most of the code is templated, and template types are calculated in
-|ctraits|. So all templated methods get a template argument |T|, which
-can be either |fold|, or |unfold|. To shorten a reference to a type
-calculated by |ctraits| for a particular |t|, we define the following
-macros.
-
-@s _Ttensor int
-@s _Ttensym int
-@s _Tg int
-@s _Tgs int
-@s _Tgss int
-@s _TG int
-@s _TZstack int
-@s _TGstack int
-@s _TZXstack int
-@s _TGXstack int
-@s __Tm int
-@s _Tpol int
-@d _Ttensor TYPENAME ctraits<t>::Ttensor 
-@d _Ttensym TYPENAME ctraits<t>::Ttensym 
-@d _Tg TYPENAME ctraits<t>::Tg
-@d _Tgs TYPENAME ctraits<t>::Tgs
-@d _Tgss TYPENAME ctraits<t>::Tgss
-@d _TG TYPENAME ctraits<t>::TG
-@d _TZstack TYPENAME ctraits<t>::TZstack
-@d _TGstack TYPENAME ctraits<t>::TGstack
-@d _TZXstack TYPENAME ctraits<t>::TZXstack
-@d _TGXstack TYPENAME ctraits<t>::TGXstack
-@d __Tm TYPENAME ctraits<t>::Tm
-@d _Tpol TYPENAME ctraits<t>::Tpol
-
-
-@<|KOrder| class declaration@>=
-class KOrder {
-protected:@;
-	const PartitionY ypart;
-	const int ny;
-	const int nu;
-	const int maxk;
-	IntSequence nvs;
-	@<|KOrder| container members@>;
-	const MatrixA matA;
-	const MatrixS matS;
-	const MatrixB matB;
-	@<|KOrder| member access method declarations@>;
-	Journal& journal;
-public:@;
-	KOrder(int num_stat, int num_pred, int num_both, int num_forw,
-		   const TensorContainer<FSSparseTensor>& fcont,
-		   const TwoDMatrix& gy, const TwoDMatrix& gu, const TwoDMatrix& v,
-		   Journal& jr);
-	enum {@+ fold, unfold@+ };
-	@<|KOrder::performStep| templated code@>;
-	@<|KOrder::check| templated code@>;
-	@<|KOrder::calcStochShift| templated code@>;
-	void switchToFolded();
-	const PartitionY& getPartY() const
-		{@+ return ypart;@+}
-	const FGSContainer& getFoldDers() const
-		{@+ return _fg;@+}
-	const UGSContainer& getUnfoldDers() const
-		{@+ return _ug;@+}
-	static bool is_even(int i)
-		{@+ return (i/2)*2 == i;@+}
-protected:@;
-	@<|KOrder::insertDerivative| templated code@>;
-	template<int t>
-	void sylvesterSolve(_Ttensor& der) const;
-
-	@<|KOrder::faaDiBrunoZ| templated code@>;
-	@<|KOrder::faaDiBrunoG| templated code@>;
-
-	@<|KOrder::recover_y| templated code@>;
-	@<|KOrder::recover_yu| templated code@>;
-	@<|KOrder::recover_ys| templated code@>;
-	@<|KOrder::recover_yus| templated code@>;
-	@<|KOrder::recover_s| templated code@>;
-	@<|KOrder::fillG| templated code@>;
-
-	@<|KOrder::calcD_ijk| templated code@>;
-	@<|KOrder::calcD_ik| templated code@>;
-	@<|KOrder::calcD_k| templated code@>;
-
-	@<|KOrder::calcE_ijk| templated code@>;
-	@<|KOrder::calcE_ik| templated code@>;
-	@<|KOrder::calcE_k| templated code@>;
-};
-
-
-
-@ Here we insert the result to the container. Along the insertion, we
-also create subtensors and insert as well.
-
-@<|KOrder::insertDerivative| templated code@>=
-template <int t>
-void insertDerivative(_Ttensor* der)
-{
-	g<t>().insert(der);
-	gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *der));
-	gss<t>().insert(new _Ttensor(ypart.nstat+ypart.npred,
-								 ypart.nyss(), *der));
-}
-
-
-@ Here we implement Faa Di Bruno formula
-$$\sum_{l=1}^k\left[f_{z^l}\right]_{\gamma_1\ldots\gamma_l}
-\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s(c_m)}\right]^{\gamma_m},
-$$
-where $s$ is a given outer symmetry and $k$ is the dimension of the
-symmetry.
-
-@<|KOrder::faaDiBrunoZ| templated code@>=
-template <int t>
-_Ttensor* faaDiBrunoZ(const Symmetry& sym) const
-{
-	JournalRecordPair pa(journal);
-	pa << "Faa Di Bruno Z container for " << sym << endrec;
-	_Ttensor* res = new _Ttensor(ny, TensorDimens(sym, nvs));
-	FaaDiBruno bruno(journal);
-	bruno.calculate(Zstack<t>(), f, *res);
-	return res;
-}
-
-@ The same as |@<|KOrder::faaDiBrunoZ| templated code@>|, but for
-$g^{**}$ and $G$ stack.
-
-@<|KOrder::faaDiBrunoG| templated code@>=
-template <int t>
-_Ttensor* faaDiBrunoG(const Symmetry& sym) const
-{
-	JournalRecordPair pa(journal);
-	pa << "Faa Di Bruno G container for " << sym << endrec;
-	TensorDimens tdims(sym, nvs);
-	_Ttensor* res = new _Ttensor(ypart.nyss(), tdims);
-	FaaDiBruno bruno(journal);
-	bruno.calculate(Gstack<t>(), gss<t>(), *res);
-	return res;
-}
-
-@ Here we solve $\left[F_{y^i}\right]=0$. First we calculate
-conditional $G_{y^i}$ (it misses $l=1$ and $l=i$ since $g_{y^i}$ does
-not exist yet). Then calculate conditional $F_{y^i}$ and we have the
-right hand side of equation. Since we miss two orders, we solve by
-Sylvester, and insert the solution as the derivative $g_{y^i}$. Then
-we need to update $G_{y^i}$ running |multAndAdd| for both dimensions
-$1$ and $i$.
-
-{\bf Requires:} everything at order $\leq i-1$
-
-{\bf Provides:} $g_{y^i}$, and $G_{y^i}$
-
-@<|KOrder::recover_y| templated code@>=
-template<int t>
-void recover_y(int i)
-{
-	Symmetry sym(i,0,0,0);
-	JournalRecordPair pa(journal);
-	pa << "Recovering symmetry " << sym << endrec;
-
-	_Ttensor* G_yi = faaDiBrunoG<t>(sym);
-	G<t>().insert(G_yi);
-
-	_Ttensor* g_yi = faaDiBrunoZ<t>(sym);
-	g_yi->mult(-1.0);
-
-	sylvesterSolve<t>(*g_yi);
-
-	insertDerivative<t>(g_yi);
-
-	_Ttensor* gss_y = gss<t>().get(Symmetry(1,0,0,0));
-	gs<t>().multAndAdd(*gss_y, *G_yi);
-	_Ttensor* gss_yi = gss<t>().get(sym);
-	gs<t>().multAndAdd(*gss_yi, *G_yi);
-}
-
-
-@ Here we solve $\left[F_{y^iu^j}\right]=0$ to obtain $g_{y^iu^j}$ for
-$j>0$. We calculate conditional $G_{y^iu^j}$ (this misses only $l=1$)
-and calculate conditional $F_{y^iu^j}$ and we have the right hand
-side. It is solved by multiplication of inversion of $A$. Then we insert
-the result, and update $G_{y^iu^j}$ by |multAndAdd| for $l=1$.
-
-{\bf Requires:} everything at order $\leq i+j-1$, $G_{y^{i+j}}$, and
-$g_{y^{i+j}}$.
-
-{\bf Provides:} $g_{y^iu^j}$, and $G_{y^iu^j}$
-
-@<|KOrder::recover_yu| templated code@>=
-template <int t>
-void recover_yu(int i, int j)
-{
-	Symmetry sym(i,j,0,0);
-	JournalRecordPair pa(journal);
-	pa << "Recovering symmetry " << sym << endrec;
-
-	_Ttensor* G_yiuj = faaDiBrunoG<t>(sym);
-	G<t>().insert(G_yiuj);
-
-	_Ttensor* g_yiuj = faaDiBrunoZ<t>(sym);
-	g_yiuj->mult(-1.0);
-	matA.multInv(*g_yiuj);
-	insertDerivative<t>(g_yiuj);
-
-	gs<t>().multAndAdd(*(gss<t>().get(Symmetry(1,0,0,0))), *G_yiuj);
-}
-
-@ Here we solve
-$\left[F_{y^i\sigma^j}\right]+\left[D_{ij}\right]+\left[E_{ij}\right]=0$
-to obtain $g_{y^i\sigma^j}$. We calculate conditional
-$G_{y^i\sigma^j}$ (missing dimensions $1$ and $i+j$), calculate
-conditional $F_{y^i\sigma^j}$. Before we can calculate $D_{ij}$ and
-$E_{ij}$, we have to calculate $G_{y^iu'^m\sigma^{j-m}}$ for
-$m=1,\ldots,j$. Then we add the $D_{ij}$ and $E_{ij}$ to obtain the
-right hand side. Then we solve the sylvester to obtain
-$g_{y^i\sigma^j}$. Then we update $G_{y^i\sigma^j}$ for $l=1$ and
-$l=i+j$.
-
-{\bf Requires:} everything at order $\leq i+j-1$, $g_{y^{i+j}}$,
-$G_{y^iu'^j}$ and $g_{y^iu^j}$ through $D_{ij}$,
-$g_{y^iu^m\sigma^{j-m}}$ for
-$m=1,\ldots,j-1$ through $E_{ij}$.
-
-{\bf Provides:} $g_{y^i\sigma^j}$ and $G_{y^i\sigma^j}$, and finally
-$G_{y^iu'^m\sigma^{j-m}}$ for $m=1,\ldots,j$. The latter is calculated
-by |fillG| before the actual calculation.
-
-@<|KOrder::recover_ys| templated code@>=
-template <int t>
-void recover_ys(int i, int j)
-{
-	Symmetry sym(i,0,0,j);
-	JournalRecordPair pa(journal);
-	pa << "Recovering symmetry " << sym << endrec;
-
-	fillG<t>(i, 0, j);
-
-	if (is_even(j)) {
-		_Ttensor* G_yisj = faaDiBrunoG<t>(sym);
-		G<t>().insert(G_yisj);
-
-		_Ttensor* g_yisj = faaDiBrunoZ<t>(sym);
-
-		{
-			_Ttensor* D_ij = calcD_ik<t>(i, j);
-			g_yisj->add(1.0, *D_ij);
-			delete D_ij;
-		}
-
-		if (j >= 3) {
-			_Ttensor* E_ij = calcE_ik<t>(i, j);
-			g_yisj->add(1.0, *E_ij);
-			delete E_ij;
-		}
-
-		g_yisj->mult(-1.0);
-
-		sylvesterSolve<t>(*g_yisj);
-
-		insertDerivative<t>(g_yisj);
-
-		Gstack<t>().multAndAdd(1, gss<t>(), *G_yisj);
-		Gstack<t>().multAndAdd(i+j, gss<t>(), *G_yisj);
-	}
-}
-
-@ Here we solve
-$\left[F_{y^iu^j\sigma^k}\right]+\left[D_{ijk}\right]+\left[E_{ijk}\right]=0$
-to obtain $g_{y^iu^j\sigma^k}$. First we calculate conditional
-$G_{y^iu^j\sigma^k}$ (missing only for dimension $l=1$), then we
-evaluate conditional $F_{y^iu^j\sigma^k}$. Before we can calculate
-$D_{ijk}$, and $E_{ijk}$, we need to insert
-$G_{y^iu^ju'^m\sigma^{k-m}}$ for $m=1,\ldots, k$. This is done by
-|fillG|. Then we have right hand side and we multiply by $A^{-1}$ to
-obtain $g_{y^iu^j\sigma^k}$. Finally we have to update
-$G_{y^iu^j\sigma^k}$ by |multAndAdd| for dimension $l=1$.
-
-{\bf Requires:} everything at order $\leq i+j+k$, $g_{y^{i+j}\sigma^k}$
-through $G_{y^iu^j\sigma^k}$ involved in right hand side, then
-$g_{y^iu^{j+k}}$ through $D_{ijk}$, and $g_{y^iu^{j+m}\sigma^{k-m}}$
-for $m=1,\ldots,k-1$ through $E_{ijk}$.
-
-{\bf Provides:} $g_{y^iu^j\sigma^k}$, $G_{y^iu^j\sigma^k}$, and
-$G_{y^iu^ju'^m\sigma^{k-m}}$ for $m=1,\ldots, k$
-
-@<|KOrder::recover_yus| templated code@>=
-template <int t>
-void recover_yus(int i, int j, int k)
-{
-	Symmetry sym(i,j,0,k);
-	JournalRecordPair pa(journal);
-	pa << "Recovering symmetry " << sym << endrec;
-
-	fillG<t>(i, j, k);
-
-	if (is_even(k)) {
-		_Ttensor* G_yiujsk = faaDiBrunoG<t>(sym);
-		G<t>().insert(G_yiujsk);
-
-		_Ttensor* g_yiujsk = faaDiBrunoZ<t>(sym);
-
-		{
-			_Ttensor* D_ijk = calcD_ijk<t>(i,j,k);
-			g_yiujsk->add(1.0, *D_ijk);
-			delete D_ijk;
-		}
-
-		if (k >= 3) {
-			_Ttensor* E_ijk = calcE_ijk<t>(i,j,k);
-			g_yiujsk->add(1.0, *E_ijk);
-			delete E_ijk;
-		}
-
-		g_yiujsk->mult(-1.0);
-
-		matA.multInv(*g_yiujsk);
-		insertDerivative<t>(g_yiujsk);
-
-		Gstack<t>().multAndAdd(1, gss<t>(), *G_yiujsk);
-	}
-}
-
-@ Here we solve
-$\left[F_{\sigma^i}\right]+\left[D_i\right]+\left[E_i\right]=0$ to
-recover $g_{\sigma^i}$. First we calculate conditional $G_{\sigma^i}$
-(missing dimension $l=1$ and $l=i$), then we calculate conditional
-$F_{\sigma^i}$. Before we can calculate $D_i$ and $E_i$, we have to
-obtain $G_{u'm\sigma^{i-m}}$ for $m=1,\ldots,i$. Than
-adding $D_i$ and $E_i$ we have the right hand side. We solve by
-$S^{-1}$ multiplication and update $G_{\sigma^i}$ by calling
-|multAndAdd| for dimension $l=1$.
-
-Recall that the solved equation here is:
-$$
-\left[f_y\right]\left[g_{\sigma^k}\right]+
-\left[f_{y^{**}_+}\right]\left[g^{**}_{y^*}\right]\left[g^*_{\sigma^k}\right]+
-\left[f_{y^{**}_+}\right]\left[g^{**}_{\sigma^k}\right]=\hbox{RHS}
-$$
-This is a sort of deficient sylvester equation (sylvester equation for
-dimension=0), we solve it by $S^{-1}$. See |@<|MatrixS| constructor
-code@>| to see how $S$ looks like.
-
-{\bf Requires:} everything at order $\leq i-1$, $g_{y^i}$ and
-$g_{y^{i-j}\sigma^j}$, then $g_{u^k}$ through $F_{u'^k}$, and
-$g_{y^mu^j\sigma^k}$ for $j=1,\ldots,i-1$ and $m+j+k=i$ through
-$F_{u'j\sigma^{i-j}}$.
-
-{\bf Provides:} $g_{\sigma^i}$, $G_{\sigma^i}$, and
-$G_{u'^m\sigma^{i-m}}$ for $m=1,\ldots,i$
-
-@<|KOrder::recover_s| templated code@>=
-template <int t>
-void recover_s(int i)
-{
-	Symmetry sym(0,0,0,i);
-	JournalRecordPair pa(journal);
-	pa << "Recovering symmetry " << sym << endrec;
-
-	fillG<t>(0, 0, i);
-
-	if (is_even(i)) {
-		_Ttensor* G_si = faaDiBrunoG<t>(sym);
-		G<t>().insert(G_si);
-
-		_Ttensor* g_si = faaDiBrunoZ<t>(sym);
-
-		{
-			_Ttensor* D_i = calcD_k<t>(i);
-			g_si->add(1.0, *D_i);
-			delete D_i;
-		}
-
-		if (i >= 3) {
-			_Ttensor* E_i = calcE_k<t>(i);
-			g_si->add(1.0, *E_i);
-			delete E_i;
-		}
-
-		g_si->mult(-1.0);
-
-
-		matS.multInv(*g_si);
-		insertDerivative<t>(g_si);
-
-		Gstack<t>().multAndAdd(1, gss<t>(), *G_si);
-		Gstack<t>().multAndAdd(i, gss<t>(), *G_si);
-	}
-}
-
-@ Here we calculate and insert $G_{y^iu^ju'^m\sigma^{k-m}}$ for
-$m=1,\ldots, k$. The derivatives are inserted only for $k-m$ being
-even.
-
-@<|KOrder::fillG| templated code@>=
-template<int t>
-void fillG(int i, int j, int k)
-{
-	for (int m = 1; m <= k; m++) {
-		if (is_even(k-m)) {
-			_Ttensor* G_yiujupms = faaDiBrunoG<t>(Symmetry(i,j,m,k-m));
-			G<t>().insert(G_yiujupms);
-		}
-	}
-}
-
-
-@ Here we calculate
-$$\left[D_{ijk}\right]_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j}=
-\left[F_{y^iu^ju'^k}\right]
-_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j\gamma_1\ldots\gamma_k}
-\left[\Sigma\right]^{\gamma_1\ldots\gamma_k}$$
-So it is non zero only for even $k$.
-
-@<|KOrder::calcD_ijk| templated code@>=
-template <int t>
-_Ttensor* calcD_ijk(int i, int j, int k) const
-{
-	_Ttensor* res =	new _Ttensor(ny, TensorDimens(Symmetry(i,j,0,0), nvs));
-	res->zeros();
-	if (is_even(k)) {
-		_Ttensor* tmp = faaDiBrunoZ<t>(Symmetry(i,j,k,0));
-		tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(k))));
-		delete tmp;
-	}
-	return res;
-}		   
-
-
-@ Here we calculate
-$$\left[E_{ijk}\right]_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j}=
-\sum_{m=1}^{k-1}\left(\matrix{k\cr m}\right)\left[F_{y^iu^ju'^m\sigma^{k-m}}\right]
-_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j\gamma_1\ldots\gamma_m}
-\left[\Sigma\right]^{\gamma_1\ldots\gamma_m}$$
-The sum can sum only for even $m$.
-
-@<|KOrder::calcE_ijk| templated code@>=
-template <int t>
-_Ttensor* calcE_ijk(int i, int j, int k) const
-{
-	_Ttensor* res = new _Ttensor(ny, TensorDimens(Symmetry(i,j,0,0), nvs));
-	res->zeros();
-	for (int n = 2; n <= k-1; n+=2) {
-		_Ttensor* tmp = faaDiBrunoZ<t>(Symmetry(i,j,n,k-n));
-		tmp->mult((double)(Tensor::noverk(k,n)));
-		tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(n))));
-		delete tmp;
-	}
-	return res;
-}
-
-@ 
-@<|KOrder::calcD_ik| templated code@>=
-template <int t>
-_Ttensor* calcD_ik(int i, int k) const
-{
-	return calcD_ijk<t>(i, 0, k);
-}
-
-@ 
-@<|KOrder::calcD_k| templated code@>=
-template <int t>
-_Ttensor* calcD_k(int k) const
-{
-	return calcD_ijk<t>(0, 0, k);
-}
-
-@ 
-@<|KOrder::calcE_ik| templated code@>=
-template <int t>
-_Ttensor* calcE_ik(int i, int k) const
-{
-	return calcE_ijk<t>(i, 0, k);
-}
-
-@ 
-@<|KOrder::calcE_k| templated code@>=
-template <int t>
-_Ttensor* calcE_k(int k) const
-{
-	return calcE_ijk<t>(0, 0, k);
-}
-
-@ Here is the core routine. It calls methods recovering derivatives in
-the right order. Recall, that the code, namely Faa Di Bruno's formula,
-is implemented as to be run conditionally on the current contents of
-containers. So, if some call of Faa Di Bruno evaluates derivatives,
-and some derivatives are not present in the container, then it is
-considered to be zero. So, we have to be very careful to put
-everything in the right order. The order here can be derived from
-dependencies, or it is in the paper.
-
-The method recovers all the derivatives of the given |order|.
-
-The precondition of the method is that all tensors of order |order-1|,
-which are not zero, exist (including $G$). The postcondition of of the
-method is derivatives of $g$ and $G$ of order |order| are calculated
-and stored in the containers. Responsibility of precondition lays upon
-the constructor (for |order==2|), or upon the previous call of
-|performStep|.
-
-From the code, it is clear, that all $g$ are calculated. If one goes
-through all the recovering methods, he should find out that also all
-$G$ are provided.
-
-@<|KOrder::performStep| templated code@>=
-template <int t>
-void performStep(int order)
-{
-	KORD_RAISE_IF(order-1 != g<t>().getMaxDim(),
-				  "Wrong order for KOrder::performStep");
-	JournalRecordPair pa(journal);
-	pa << "Performing step for order = " << order << endrec;
-
-	recover_y<t>(order);
-
-	for (int i = 0; i < order; i++) {
-		recover_yu<t>(i, order-i);
-	}
-
-	for (int j = 1; j < order; j++) {
-		for (int i = j-1; i >= 1; i--) {
-			recover_yus<t>(order-j,i,j-i);
-		}
-   		recover_ys<t>(order-j, j);
-	}
-
-	for (int i = order-1; i >= 1; i--) {
-		recover_yus<t>(0, i, order-i);
-	}
-	recover_s<t>(order);
-}
-
-@ Here we check for residuals of all the solved equations at the given
-order. The method returns the largest residual size. Each check simply
-evaluates the equation.
-
-@<|KOrder::check| templated code@>=
-template <int t>
-double check(int dim) const
-{
-	KORD_RAISE_IF(dim > g<t>().getMaxDim(),
-				  "Wrong dimension for KOrder::check");
-	JournalRecordPair pa(journal);
-	pa << "Checking residuals for order = " << dim << endrec;
-
-	double maxerror = 0.0;
-
-	@<check for $F_{y^iu^j}=0$@>;
-	@<check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$@>;@q'@>
-	@<check for $F_{\sigma^i}+D_i+E_i=0$@>;
-
-	return maxerror;
-}
-
-
-@ 
-@<check for $F_{y^iu^j}=0$@>=
-	for (int i = 0;	i <= dim; i++) {
-		Symmetry sym(dim-i, i, 0,  0);
-		_Ttensor* r = faaDiBrunoZ<t>(sym);
-		double err = r->getData().getMax();
-		JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-		if (err > maxerror)
-			maxerror = err;
-		delete r;
-	}
-
-@ 
-@<check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$@>=
-	SymmetrySet ss(dim, 3);
-	for (symiterator si(ss); !si.isEnd(); ++si) {
-		int i = (*si)[0];
-		int j = (*si)[1];
-		int k = (*si)[2];
-		if (i+j > 0 && k > 0) {
-			Symmetry sym(i, j, 0, k);
-			_Ttensor* r = faaDiBrunoZ<t>(sym);
-			_Ttensor* D_ijk = calcD_ijk<t>(i,j,k);
-			r->add(1.0, *D_ijk);
-			delete D_ijk;
-			_Ttensor* E_ijk = calcE_ijk<t>(i,j,k);
-			r->add(1.0, *E_ijk);
-			delete E_ijk;
-			double err = r->getData().getMax();
-			JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-			delete r;
-		}
-	}
-
-
-@ 
-@<check for $F_{\sigma^i}+D_i+E_i=0$@>=
-	_Ttensor* r = faaDiBrunoZ<t>(Symmetry(0,0,0,dim));
-	_Ttensor* D_k = calcD_k<t>(dim);
-	r->add(1.0, *D_k);
-	delete D_k;
-	_Ttensor* E_k = calcE_k<t>(dim);
-	r->add(1.0, *E_k);
-	delete E_k;
-	double err = r->getData().getMax();
-	Symmetry sym(0,0,0,dim);
-	JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
-	if (err > maxerror)
-		maxerror = err;
-	delete r;
-
-@ 
-@<|KOrder::calcStochShift| templated code@>=
-template <int t>
-Vector* calcStochShift(int order, double sigma) const
-{
-	Vector* res = new Vector(ny);
-	res->zeros();
-	int jfac = 1;
-	for (int j = 1; j <= order; j++, jfac *= j)
-		if (is_even(j)) {
-			_Ttensor* ten = calcD_k<t>(j);
-			res->add(std::pow(sigma, j)/jfac, ten->getData());
-			delete ten;
-		}
-	return res;
-}
-
-
-@ These are containers. The names are not important because they do
-not appear anywhere else since we access them by template functions.
-
-@<|KOrder| container members@>=
-	UGSContainer _ug;
-	FGSContainer _fg;
-	UGSContainer _ugs;
-	FGSContainer _fgs;
-	UGSContainer _ugss;
-	FGSContainer _fgss;
-	UGSContainer _uG;
-	FGSContainer _fG;
-	UnfoldedZContainer _uZstack;
-	FoldedZContainer _fZstack;
-	UnfoldedGContainer _uGstack;
-	FoldedGContainer _fGstack;
-	UNormalMoments _um;
-	FNormalMoments _fm;
-	const TensorContainer<FSSparseTensor>& f;
-
-@ These are the declarations of the template functions accessing the
-containers.
-
-@<|KOrder| member access method declarations@>=
-	template<int t> _Tg& g();
-	template<int t> const _Tg& g() const;
-	template<int t> _Tgs& gs();
-	template<int t> const _Tgs& gs() const;
-	template<int t> _Tgss& gss();
-	template<int t> const _Tgss& gss() const;
-	template<int t> _TG& G();
-	template<int t> const _TG& G() const;
-	template<int t> _TZstack& Zstack();
-	template<int t> const _TZstack& Zstack() const;
-	template<int t> _TGstack& Gstack();
-	template<int t> const _TGstack& Gstack() const;
-	template<int t> __Tm& m();
-	template<int t> const __Tm& m() const;
-
-
-@ End of {\tt korder.h} file.

dynare++/kord/korder_stoch.cweb

-@q $Id: korder_stoch.cweb 148 2005-04-19 15:12:26Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ Start of {\tt korder\_stoch.cpp} file.
-@c
-#include "korder_stoch.h"
-
-@<|MatrixAA| constructor code@>;
-@<|KOrderStoch| folded constructor code@>;
-@<|KOrderStoch| unfolded constructor code@>;
-@<|KOrderStoch| convenience method specializations@>;
-
-@ Same as |@<|MatrixA| constructor code@>|, but the submatrix |gss_ys| is passed directly.
-@<|MatrixAA| constructor code@>=
-MatrixAA::MatrixAA(const FSSparseTensor& f, const IntSequence& ss,
-				   const TwoDMatrix& gss_ys, const PartitionY& ypart)
-	: PLUMatrix(ypart.ny())
-{
-	zeros();
-
-	IntSequence c(1); c[0] = 1;
-	FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
-	add(1.0, f_y);
-
-	c[0] = 0;
-	FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
-	TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
-	sub.multAndAdd(f_yss, gss_ys);
-
-	calcPLU();
-}
-	
-
-@ 
-@<|KOrderStoch| folded constructor code@>=
-KOrderStoch::KOrderStoch(const PartitionY& yp, int nu,
-						 const TensorContainer<FSSparseTensor>& fcont,
-						 const FGSContainer& hh, Journal& jr)
-	: nvs(4), ypart(yp), journal(jr),@/
-	  _ug(4), _fg(4), _ugs(4), _fgs(4), _uG(4), _fG(4),@/
-	  _uh(NULL), _fh(&hh),@/
-	  _uZstack(&_uG, ypart.nyss(), &_ug, ypart.ny(), ypart.nys(), nu),@/
-	  _fZstack(&_fG, ypart.nyss(), &_fg, ypart.ny(), ypart.nys(), nu),@/
-	  _uGstack(&_ugs, ypart.nys(), nu),@/
-	  _fGstack(&_fgs, ypart.nys(), nu),@/
-	  f(fcont),@/
-	  matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1,0,0,0))),
-		   ypart)
-{
-	nvs[0] = ypart.nys();
-	nvs[1] = nu;
-	nvs[2] = nu;
-	nvs[3] = 1;
-}
-
-@ 
-@<|KOrderStoch| unfolded constructor code@>=
-KOrderStoch::KOrderStoch(const PartitionY& yp, int nu,
-						 const TensorContainer<FSSparseTensor>& fcont,
-						 const UGSContainer& hh, Journal& jr)
-	: nvs(4), ypart(yp), journal(jr),@/
-	  _ug(4), _fg(4), _ugs(4), _fgs(4), _uG(4), _fG(4),@/
-	  _uh(&hh), _fh(NULL),@/
-	  _uZstack(&_uG, ypart.nyss(), &_ug, ypart.ny(), ypart.nys(), nu),@/
-	  _fZstack(&_fG, ypart.nyss(), &_fg, ypart.ny(), ypart.nys(), nu),@/
-	  _uGstack(&_ugs, ypart.nys(), nu),@/
-	  _fGstack(&_fgs, ypart.nys(), nu),@/
-	  f(fcont),@/
-	  matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1,0,0,0))),
-		   ypart)
-{
-	nvs[0] = ypart.nys();
-	nvs[1] = nu;
-	nvs[2] = nu;
-	nvs[3] = 1;
-}
-
-
-@ 
-@<|KOrderStoch| convenience method specializations@>=
-	template<> ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>()
-		{@+ return _ug;@+}
-	template<>@; const ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>()@+const@;
-		{@+ return _ug;@+}
-	template<> ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>()
-		{@+ return _fg;@+}
-	template<> const ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>()@+const@;
-		{@+ return _fg;@+}
-	template<> ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>()
-		{@+ return _ugs;@+}
-	template<> const ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>()@+const@;
-		{@+ return _ugs;@+}
-	template<> ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>()
-		{@+ return _fgs;@+}
-	template<> const ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>()@+const@;
-		{@+ return _fgs;@+}
-	template<> const ctraits<KOrder::unfold>::Tgss& KOrderStoch::h<KOrder::unfold>()@+const@;
-		{@+ return *_uh;@+}
-	template<> const ctraits<KOrder::fold>::Tgss& KOrderStoch::h<KOrder::fold>()@+const@;
-		{@+ return *_fh;@+}
-	template<> ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>()
-		{@+ return _uG;@+}
-	template<> const ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>()@+const@;
-		{@+ return _uG;@+}
-	template<> ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>()
-		{@+ return _fG;@+}
-	template<> const ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>()@+const@;
-		{@+ return _fG;@+}
-	template<> ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>()
-		{@+ return _uZstack;@+}
-	template<> const ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>()@+const@;
-		{@+ return _uZstack;@+}
-	template<> ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>()
-		{@+ return _fZstack;@+}
-	template<> const ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>()@+const@;
-		{@+ return _fZstack;@+}
-	template<> ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>()
-		{@+ return _uGstack;@+}
-	template<> const ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>()@+const@;
-		{@+ return _uGstack;@+}
-	template<> ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>()
-		{@+ return _fGstack;@+}
-	template<> const ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>()@+const@;
-		{@+ return _fGstack;@+}
-
-
-@ End of {\tt korder\_stoch.cpp} file.

dynare++/kord/korder_stoch.hweb

-@q $Id: korder_stoch.hweb 418 2005-08-16 15:10:06Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Higher order at stochastic steady. Start of {\tt korder\_stoch.h} file.
-
-This file defines a number of classes of which |KOrderStoch| is the
-main purpose. Basically, |KOrderStoch| calculates first and higher
-order Taylor expansion of a policy rule at $\sigma>0$ with explicit
-forward $g^{**}$. More formally, we have to solve a policy rule $g$
-from
-$$E_t\left[f(g^{**}(g^*(y^*_t,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y^*,u_t)\right]$$
-As an introduction in {\tt approximation.hweb} argues, $g^{**}$ at
-tine $t+1$ must be given from outside. Let the explicit
-$E_t(g^{**}(y^*,u_{t+1},\sigma)$ be equal to $h(y^*,\sigma)$. Then we
-have to solve
-$$f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u),$$
-which is much easier than fully implicit system for $\sigma=0$.
-
-Besides the class |KOrderStoch|, we declare here also classes for the
-new containers corresponding to
-$f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. Further, we
-declare |IntegDerivs| and |StochForwardDerivs| classes which basically
-calculate $h$ as an extrapolation based on an approximation to $g$ at
-lower $\sigma$.
-
-@s IntegDerivs int
-@s StochForwardDerivs int
-@s GContainer int
-@s ZContainer int
-@s GXContainer int
-@s ZXContainer int
-@s MatrixAA int
-@s KOrderStoch int
-@s StackContainer int
-@s _Ttype int
-@s _Ctype int
-@s UnfoldedStackContainer int
-@s FoldedStackContainer int
-
-@c
-#include "korder.h"
-#include "faa_di_bruno.h"
-#include "journal.h"
-
-
-@<|IntegDerivs| class declaration@>;
-@<|StochForwardDerivs| class declaration@>;
-@<|GXContainer| class declaration@>;
-@<|ZXContainer| class declaration@>;
-@<|UnfoldedGXContainer| class declaration@>;
-@<|FoldedGXContainer| class declaration@>;
-@<|UnfoldedZXContainer| class declaration@>;
-@<|FoldedZXContainer| class declaration@>;
-@<|MatrixAA| class declaration@>;
-@<|KOrderStoch| class declaration@>;
-
-@ This class is a container, which has a specialized constructor
-integrating the policy rule at given $\sigma$.
-
-@<|IntegDerivs| class declaration@>=
-template <int t>
-class IntegDerivs : public ctraits<t>::Tgss {
-public:@;
-	@<|IntegDerivs| constructor code@>;
-};
-
-@ This constuctor integrates a rule (namely its $g^{**}$ part) with
-respect to $u=\tilde\sigma\eta$, and stores to the object the
-derivatives of this integral $h$ at $(y^*,u,\sigma)=(\tilde
-y^*,0,\tilde\sigma)$. The original container of $g^{**}$, the moments of
-the stochastic shocks |mom| and the $\tilde\sigma$ are input.
-
-The code follows the following derivation
-\def\lims{\vbox{\baselineskip=0pt\lineskip=1pt
-    \setbox0=\hbox{$\scriptstyle n+k=p$}\hbox to\wd0{\hss$\scriptstyle m=0$\hss}\box0}}
-$$
-\eqalign{h(y,\sigma)&=E_t\left[g(y,u',\sigma)\right]=\cr
-&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+j+k=d}\pmatrix{d\cr i,j,k}\left[g_{y^iu^j\sigma^k}\right]
-  (y^*-\tilde y^*)^i\sigma^j\Sigma^j(\sigma-\tilde\sigma)^k\cr
-&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m+n,k}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \hat y^{*i}\Sigma^{m+n}\pmatrix{m+n\cr m,n}{\tilde\sigma}^m\hat\sigma^{k+n}\cr
-&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m,n,k}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
-&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\sum_{\lims}\pmatrix{d\cr i,m,n,k}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
-&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\pmatrix{d\cr i,p}
-  \left[\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \Sigma^{m+n}{\tilde\sigma}^m\right]\hat y^{*i}\hat\sigma^{k+n},
-}
-$$
-where $\pmatrix{a\cr b_1,\ldots, b_n}$ is a generalized combination
-number, $p=k+n$, $\hat\sigma=\sigma-\tilde\sigma$, $\hat
-y^*=y^*-\tilde y$, and we dropped writing the multidimensional indexes
-in Einstein summation.
-
-This implies that
-$$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \Sigma^{m+n}{\tilde\sigma}^m$$
-and this is exactly what the code does.
-
-@<|IntegDerivs| constructor code@>=
-IntegDerivs(int r, const IntSequence& nvs, const _Tgss& g, const __Tm& mom,
-			double at_sigma)
-	: ctraits<t>::Tgss(4)
-{
-	int maxd = g.getMaxDim();
-	for (int d = 1; d <= maxd; d++) {
-		for (int i = 0; i <= d; i++) {
-			int p = d-i;
-			Symmetry sym(i,0,0,p);
-			_Ttensor* ten = new _Ttensor(r, TensorDimens(sym, nvs));
-			@<calculate derivative $h_{y^i\sigma^p}$@>;
-			this->insert(ten);
-		}
-	}
-}
-
-@ This code calculates
-$$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
-  \left[g_{y^iu^{m+n}\sigma^k}\right]
-  \Sigma^{m+n}{\tilde\sigma}^m$$
-and stores it in |ten|.
-
-@<calculate derivative $h_{y^i\sigma^p}$@>=
-	ten->zeros();
-	for (int n = 0; n <= p; n++) {
-		int k = p-n;
-		int povern = Tensor::noverk(p,n);
-		int mfac = 1;
-		for (int m = 0; i+m+n+k <= maxd; m++, mfac*=m) {
-			double mult = (pow(at_sigma,m)*povern)/mfac;
-			Symmetry sym_mn(i,m+n,0,k);
-			if (m+n == 0 && g.check(sym_mn))
-				ten->add(mult, *(g.get(sym_mn)));
-			if (m+n > 0 && KOrder::is_even(m+n) && g.check(sym_mn)) {
-				_Ttensor gtmp(*(g.get(sym_mn)));
-				gtmp.mult(mult);
-				gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry(m+n))));
-			}
-		}
-	}
-
-@ This class calculates an extrapolation of expectation of forward
-derivatives. It is a container, all calculations are done in a
-constructor.
-
-The class calculates derivatives of $E[g(y*,u,\sigma)]$ at $(\bar
-y^*,\bar\sigma)$. The derivatives are extrapolated based on
-derivatives at $(\tilde y^*,\tilde\sigma)$.
-
-@<|StochForwardDerivs| class declaration@>=
-template <int t>
-class StochForwardDerivs : public ctraits<t>::Tgss {
-public:@;
-	@<|StochForwardDerivs| constructor code@>;
-};
-
-@ This is the constructor which performs the integration and the
-extrapolation. Its parameters are: |g| is the container of derivatives
-at $(\tilde y,\tilde\sigma)$; |m| are the moments of stochastic
-shocks; |ydelta| is a difference of the steady states $\bar y-\tilde
-y$; |sdelta| is a difference between new sigma and old sigma
-$\bar\sigma-\tilde\sigma$, and |at_sigma| is $\tilde\sigma$. There is
-no need of inputing the $\tilde y$.
-
-We do the operation in four steps:
-\orderedlist
-\li Integrate $g^{**}$, the derivatives are at $(\tilde y,\tilde\sigma)$
-\li Form the (full symmetric) polynomial from the derivatives stacking
-$\left[\matrix{y^*\cr\sigma}\right]$
-\li Centralize this polynomial about $(\bar y,\bar\sigma)$
-\li Recover general symmetry tensors from the (full symmetric) polynomial
-\endorderedlist
-
-@<|StochForwardDerivs| constructor code@>=
-StochForwardDerivs(const PartitionY& ypart, int nu,
-				   const _Tgss& g, const __Tm& m,
-				   const Vector& ydelta, double sdelta,
-				   double at_sigma)
-	: ctraits<t>::Tgss(4)
-{
-	int maxd = g.getMaxDim();
-	int r = ypart.nyss();
-
-	@<make |g_int| be integral of $g^{**}$ at $(\tilde y,\tilde\sigma)$@>;
-	@<make |g_int_sym| be full symmetric polynomial from |g_int|@>;
-	@<make |g_int_cent| the centralized polynomial about $(\bar y,\bar\sigma)$@>;
-	@<pull out general symmetry tensors from |g_int_cent|@>;
-}
-
-@ This simply constructs |IntegDerivs| class. Note that the |nvs| of
-the tensors has zero dimensions for shocks, this is because we need to
-make easily stacks of the form $\left[\matrix{y^*\cr\sigma}\right]$ in
-the next step.
-
-@<make |g_int| be integral of $g^{**}$ at $(\tilde y,\tilde\sigma)$@>=
-	IntSequence nvs(4);
-	nvs[0] = ypart.nys(); nvs[1] = 0; nvs[2] = 0; nvs[3] = 1;
-	IntegDerivs<t> g_int(r, nvs, g, m, at_sigma);
-
-@ Here we just form a polynomial whose unique variable corresponds to
-$\left[\matrix{y^*\cr\sigma}\right]$ stack.
-
-@<make |g_int_sym| be full symmetric polynomial from |g_int|@>=
-	_Tpol g_int_sym(r, ypart.nys()+1);
-	for (int d = 1; d <= maxd; d++) {
-		_Ttensym* ten = new _Ttensym(r, ypart.nys()+1, d);
-		ten->zeros();
-		for (int i = 0; i <= d; i++) {
-			int k = d-i;
-			if (g_int.check(Symmetry(i,0,0,k)))
-				ten->addSubTensor(*(g_int.get(Symmetry(i,0,0,k))));
-		}
-		g_int_sym.insert(ten);
-	}
-
-@ Here we centralize the polynomial to $(\bar y,\bar\sigma)$ knowing
-that the polynomial was centralized about $(\tilde
-y,\tilde\sigma)$. This is done by derivating and evaluating the
-derivated polynomial at $(\bar y-\tilde
-y,\bar\sigma-\tilde\sigma)$. The stack of this vector is |delta| in
-the code.
-
-@<make |g_int_cent| the centralized polynomial about $(\bar y,\bar\sigma)$@>=
-	Vector delta(ypart.nys()+1);
-	Vector dy(delta, 0, ypart.nys());
-	ConstVector dy_in(ydelta, ypart.nstat, ypart.nys());
-	dy = dy_in;
-	delta[ypart.nys()] = sdelta;
-	_Tpol g_int_cent(r, ypart.nys()+1);
-	for (int d = 1; d <= maxd; d++) {
-		g_int_sym.derivative(d-1);
-		_Ttensym* der = g_int_sym.evalPartially(d, delta);
-		g_int_cent.insert(der);
-	}
-
-@ Here we only recover the general symmetry derivatives from the full
-symmetric polynomial. Note that the derivative get the true |nvs|.
-
-@<pull out general symmetry tensors from |g_int_cent|@>=
-	IntSequence ss(4);
-	ss[0]=ypart.nys(); ss[1]=0; ss[2]=0; ss[3]=1;
-	IntSequence pp(4);
-	pp[0]=0; pp[1]=1; pp[2]=2; pp[3]=3;
-	IntSequence true_nvs(nvs);
-	true_nvs[1]=nu; true_nvs[2]=nu;
-	for (int d = 1; d <= maxd; d++) {
-		if (g_int_cent.check(Symmetry(d))) {
-			for (int i = 0;	 i <= d; i++) {
-				Symmetry sym(i, 0, 0, d-i);
-				IntSequence coor(sym, pp);
-				_Ttensor* ten = new _Ttensor(*(g_int_cent.get(Symmetry(d))), ss, coor,
-											 TensorDimens(sym, true_nvs));
-				this->insert(ten);
-			}
-		}
-	}
-
-
-@ This container corresponds to $h(g^*(y,u,\sigma),\sigma)$. Note that
-in our application, the $\sigma$ as a second argument to $h$ will be
-its fourth variable in symmetry, so we have to do four member stack
-having the second and third stack dummy.
-
-@<|GXContainer| class declaration@>=
-template <class _Ttype>
-class GXContainer : public GContainer<_Ttype> {
-public:@;
-	typedef StackContainerInterface<_Ttype> _Stype;
-	typedef typename StackContainer<_Ttype>::_Ctype _Ctype;
-	typedef typename StackContainer<_Ttype>::itype itype;
-	GXContainer(const _Ctype* gs, int ngs, int nu)
-		: GContainer<_Ttype>(gs, ngs, nu)@+ {}
-	@<|GXContainer::getType| code@>;
-};
-
-@ This routine corresponds to this stack:
-$$\left[\matrix{g^*(y,u,\sigma)\cr dummy\cr dummy\cr\sigma}\right]$$
-
-@<|GXContainer::getType| code@>=
-itype getType(int i, const Symmetry& s) const
-{
-	if (i == 0)
-		if (s[2] > 0)
-			return _Stype::zero;
-		else
-			return _Stype::matrix;
-	if (i == 1)
-		return _Stype::zero;
-	if (i == 2)
-		return _Stype::zero;
-	if (i == 3)
-		if (s == Symmetry(0,0,0,1))
-			return _Stype::unit;
-		else
-			return _Stype::zero;
-
-	KORD_RAISE("Wrong stack index in GXContainer::getType");
-}
-
-
-@ This container corresponds to $f(H(y,u,\sigma),g(y,u,sigma),y,u)$,
-where the $H$ has the size (number of rows) as $g^{**}$. Since it is
-very simmilar to |ZContainer|, we inherit form it and override only
-|getType| method.
-
-@<|ZXContainer| class declaration@>=
-template <class _Ttype>
-class ZXContainer : public ZContainer<_Ttype> {
-public:@;
-	typedef StackContainerInterface<_Ttype> _Stype;
-	typedef typename StackContainer<_Ttype>::_Ctype _Ctype;
-	typedef typename StackContainer<_Ttype>::itype itype;
-	ZXContainer(const _Ctype* gss, int ngss, const _Ctype* g, int ng, int ny, int nu)
-		: ZContainer<_Ttype>(gss, ngss, g, ng, ny, nu) @+{}
-	@<|ZXContainer::getType| code@>;
-};
-
-@ This |getType| method corresponds to this stack:
-$$\left[\matrix{H(y,u,\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$
-
-@<|ZXContainer::getType| code@>=
-itype getType(int i, const Symmetry& s) const
-{
-	if (i == 0)
-		if (s[2] > 0)
-			return _Stype::zero;
-		else
-			return _Stype::matrix;
-	if (i == 1)
-		if (s[2] > 0)
-			return _Stype::zero;
-		else
-			return _Stype::matrix;
-	if (i == 2)
-		if (s == Symmetry(1,0,0,0))
-			return _Stype::unit;
-		else
-			return _Stype::zero;
-	if (i == 3)
-		if (s == Symmetry(0,1,0,0))
-			return _Stype::unit;
-		else
-			return _Stype::zero;
-
-	KORD_RAISE("Wrong stack index in ZXContainer::getType");
-}
-
-@ 
-@<|UnfoldedGXContainer| class declaration@>=
-class UnfoldedGXContainer : public GXContainer<UGSTensor>, public UnfoldedStackContainer {
-public:@;
-	typedef TensorContainer<UGSTensor> _Ctype;
-	UnfoldedGXContainer(const _Ctype* gs, int ngs, int nu)
-		: GXContainer<UGSTensor>(gs, ngs, nu)@+ {}
-};
-
-@ 
-@<|FoldedGXContainer| class declaration@>=
-class FoldedGXContainer : public GXContainer<FGSTensor>, public FoldedStackContainer {
-public:@;
-	typedef TensorContainer<FGSTensor> _Ctype;
-	FoldedGXContainer(const _Ctype* gs, int ngs, int nu)
-		: GXContainer<FGSTensor>(gs, ngs, nu)@+ {}
-};
-
-@ 
-@<|UnfoldedZXContainer| class declaration@>=
-class UnfoldedZXContainer : public ZXContainer<UGSTensor>, public UnfoldedStackContainer {
-public:@;
-	typedef TensorContainer<UGSTensor> _Ctype;
-	UnfoldedZXContainer(const _Ctype* gss, int ngss, const _Ctype* g, int ng, int ny, int nu)
-		: ZXContainer<UGSTensor>(gss, ngss, g, ng, ny, nu)@+ {}
-};
-
-@ 
-@<|FoldedZXContainer| class declaration@>=
-class FoldedZXContainer : public ZXContainer<FGSTensor>, public FoldedStackContainer {
-public:@;
-	typedef TensorContainer<FGSTensor> _Ctype;
-	FoldedZXContainer(const _Ctype* gss, int ngss, const _Ctype* g, int ng, int ny, int nu)
-		: ZXContainer<FGSTensor>(gss, ngss, g, ng, ny, nu)@+ {}
-};
-
-@ This matrix corresponds to
-$$\left[f_{y}\right]+ \left[0
-\left[f_{y^{**}_+}\right]\cdot\left[h^{**}_{y^*}\right] 0\right]$$
-This is very the same as |MatrixA|, the only difference that the
-|MatrixA| is constructed from whole $h_{y^*}$, not only from
-$h^{**}_{y^*}$, hence the new abstraction.
-
-@<|MatrixAA| class declaration@>=
-class MatrixAA : public PLUMatrix {
-public:@;
-	MatrixAA(const FSSparseTensor& f, const IntSequence& ss,
-			 const TwoDMatrix& gyss, const PartitionY& ypart);
-};
-
-
-@ This class calculates derivatives of $g$ given implicitly by
-$f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$, where $h(y,\sigma)$
-is given from outside.
-
-Structurally, the class is very similar to |KOrder|, but calculations
-are much easier. The two constructors construct an object from sparse
-derivatives of $f$, and derivatives of $h$. The caller must ensure
-that the both derivatives are done at the same point.
-
-The calculation for order $k$ (including $k=1$) is done by a call
-|performStep(k)|. The derivatives can be retrived by |getFoldDers()|
-or |getUnfoldDers()|.
-
-@<|KOrderStoch| class declaration@>=
-class KOrderStoch {
-protected:@;
-	IntSequence nvs;
-	PartitionY ypart;
-	Journal& journal;
-	UGSContainer _ug;
-	FGSContainer _fg;
-	UGSContainer _ugs;
-	FGSContainer _fgs;
-	UGSContainer _uG;
-	FGSContainer _fG;
-	const UGSContainer* _uh;
-	const FGSContainer* _fh;
-	UnfoldedZXContainer _uZstack;
-	FoldedZXContainer _fZstack;
-	UnfoldedGXContainer _uGstack;
-	FoldedGXContainer _fGstack;
-	const TensorContainer<FSSparseTensor>& f;
-	MatrixAA matA;
-public:@;
-	KOrderStoch(const PartitionY& ypart, int nu, const TensorContainer<FSSparseTensor>& fcont,
-				const FGSContainer& hh, Journal& jr);
-	KOrderStoch(const PartitionY& ypart, int nu, const TensorContainer<FSSparseTensor>& fcont,
-				const UGSContainer& hh, Journal& jr);
-	@<|KOrderStoch::performStep| templated code@>;
-	const FGSContainer& getFoldDers() const
-		{@+ return _fg;@+}
-	const UGSContainer& getUnfoldDers() const
-		{@+ return _ug;@+}	
-protected:@;
-	@<|KOrderStoch::faaDiBrunoZ| templated code@>;
-	@<|KOrderStoch::faaDiBrunoG| templated code@>;
-	@<|KOrderStoch| convenience access methods@>;
-};
-
-@ This calculates a derivative of $f(G(y,u,\sigma),g(y,u,\sigma),y,u)$
-of a given symmetry.
-
-@<|KOrderStoch::faaDiBrunoZ| templated code@>=
-template <int t>
-_Ttensor* faaDiBrunoZ(const Symmetry& sym) const
-{
-	JournalRecordPair pa(journal);
-	pa << "Faa Di Bruno ZX container for " << sym << endrec;
-	_Ttensor* res = new _Ttensor(ypart.ny(), TensorDimens(sym, nvs));
-	FaaDiBruno bruno(journal);
-	bruno.calculate(Zstack<t>(), f, *res);
-	return res;
-}
-
-@ This calculates a derivative of
-$G(y,u,\sigma)=h(g^*(y,u,\sigma),\sigma)$ of a given symmetry.
-
-@<|KOrderStoch::faaDiBrunoG| templated code@>=
-template <int t>
-_Ttensor* faaDiBrunoG(const Symmetry& sym) const
-{
-	JournalRecordPair pa(journal);
-	pa << "Faa Di Bruno GX container for " << sym << endrec;
-	TensorDimens tdims(sym, nvs);
-	_Ttensor* res = new _Ttensor(ypart.nyss(), tdims);
-	FaaDiBruno bruno(journal);
-	bruno.calculate(Gstack<t>(), h<t>(), *res);
-	return res;
-}
-
-@ This retrives all $g$ derivatives of a given dimension from implicit
-$f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. It suppose that all
-derivatives of smaller dimensions have been retrieved.
-
-So, we go through all symmetries $s$, calculate $G_s$ conditional on
-$g_s=0$, insert the derivative to the $G$ container, then calculate
-$F_s$ conditional on $g_s=0$. This is a righthand side. The left hand
-side is $matA\cdot g_s$. The $g_s$ is retrieved as
-$$g_s=-matA^{-1}\cdot RHS.$$ Finally we have to update $G_s$ by
-calling |Gstack<t>().multAndAdd(1, h<t>(), *G_sym)|.
-
-@<|KOrderStoch::performStep| templated code@>=
-template <int t>
-void performStep(int order)
-{
-	int maxd = g<t>().getMaxDim();
-	KORD_RAISE_IF(order-1 != maxd && (order != 1 || maxd != -1),
-				  "Wrong order for KOrderStoch::performStep");
-	SymmetrySet ss(order, 4);
-	for (symiterator si(ss); !si.isEnd(); ++si) {
-		if ((*si)[2] == 0) {
-			JournalRecordPair pa(journal);
-			pa << "Recovering symmetry " << *si << endrec;
-
-			_Ttensor* G_sym = faaDiBrunoG<t>(*si);
-			G<t>().insert(G_sym);
-
-			_Ttensor* g_sym = faaDiBrunoZ<t>(*si);
-			g_sym->mult(-1.0);
-			matA.multInv(*g_sym);
-			g<t>().insert(g_sym);
-			gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *g_sym));
-
-			Gstack<t>().multAndAdd(1, h<t>(), *G_sym);
-		}
-	}
-}
-
-@ 
-@<|KOrderStoch| convenience access methods@>=
-	template<int t> _Tg& g();
-	template<int t> const _Tg& g() const;
-	template<int t> _Tgs& gs();
-	template<int t> const _Tgs& gs() const;
-	template<int t> const _Tgss& h() const;
-	template<int t> _TG& G();
-	template<int t> const _TG& G() const;
-	template<int t> _TZXstack& Zstack();
-	template<int t> const _TZXstack& Zstack() const;
-	template<int t> _TGXstack& Gstack();
-	template<int t> const _TGXstack& Gstack() const;
-
-
-@ End of {\tt korder\_stoch.h} file.

dynare++/kord/main.web

-@q $Id: main.web 2333 2009-01-14 10:32:55Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@q cwebmac.tex defines its own \ifpdf, which is incompatible with the @>
-@q \ifpdf defined by eplain, so undefine it @>
-\let\ifpdf\relax
-\input eplain
-
-@q now define \ifpdf to be always false: PDF macros of cwebmac are buggy @>
-\newif\ifpdf
-\iffalse\fi
-
-\def\title{{\mainfont Dynare++}}
-
-
-@i ../c++lib.w
-@s const_reverse_iterator int
-@s value_type int
-
-\titletrue
-\null\vfill
-\centerline{\titlefont Dynare++ DSGE solver}
-\vskip\baselineskip
-\centerline{\vtop{\hsize=10cm\leftskip=0pt plus 1fil
-  \rightskip=0pt plus 1fil\noindent
-	solves higher order approximation to a decision rule of a Dynamic Stochastic
-	General Equilibrium model about deterministic and stochastic fix point}}
-\vfill\vfill
-Copyright \copyright\ 2004, 2005, 2006, 2007, 2008, 2009 by Ondra Kamenik
-
-@*1 Utilities.
-@i kord_exception.hweb
-
-@i journal.hweb
-@i journal.cweb
-
-@i normal_conjugate.hweb
-@i normal_conjugate.cweb
-
-@i random.hweb
-@i random.cweb
-
-@i mersenne_twister.hweb
-
-@i faa_di_bruno.hweb
-@i faa_di_bruno.cweb
-
-@*1 Retrieving derivatives.
-
-@i first_order.hweb
-@i first_order.cweb
-
-@i korder.hweb
-@i korder.cweb
-
-@i korder_stoch.hweb
-@i korder_stoch.cweb
-
-@*1 Putting all together.
-
-@i dynamic_model.hweb
-@i dynamic_model.cweb
-
-@i approximation.hweb
-@i approximation.cweb
-
-@i decision_rule.hweb
-@i decision_rule.cweb
-
-@i global_check.hweb
-@i global_check.cweb
-
-@*1 Index.
\ No newline at end of file

dynare++/kord/mersenne_twister.hweb

-@q $Id: mersenne_twister.hweb 1490 2007-12-19 14:29:46Z kamenik $ @>
-@q Copyright 2007, Ondra Kamenik @>
-
-@*2 Mersenne Twister PRNG. Start of {\tt mersenne\_twister.h} file.
-
-This file provides a class for generating random numbers with
-encapsulated state. It is based on the work of Makoto Matsumoto and
-Takuji Nishimura, implementation inspired by code of Richard Wagner
-and Geoff Kuenning.
-
-@s uint32 int
-@s MersenneTwister int
-
-@c
-#ifndef MERSENNE_TWISTER_H
-#define MERSENNE_TWISTER_H
-
-#include "random.h"
-#include <cstring>
-
-@<|MersenneTwister| class declaration@>;
-@<|MersenneTwister| inline method definitions@>;
-
-#endif
-
-@ 
-@<|MersenneTwister| class declaration@>=
-class MersenneTwister : public RandomGenerator {
-protected:@;
-	typedef unsigned int uint32; 
-	enum {STATE_SIZE = 624};
-	enum {RECUR_OFFSET = 397};
-	uint32 statevec[STATE_SIZE];
-	int stateptr;
-public:@;
-	MersenneTwister(uint32 iseed);
-	MersenneTwister(const MersenneTwister& mt);
-	virtual ~MersenneTwister() {}
-	uint32 lrand();
-	double drand();
-	double uniform()
-		{@+return drand();@+}
-protected:@;
-	void seed(uint32 iseed);
-	void refresh();
-private:@;
-	@<|MersenneTwister| static inline methods@>;
-};
-
-@ 
-@<|MersenneTwister| static inline methods@>=
-	static uint32 hibit(uint32 u)
-		{return u & 0x80000000UL;}
-	static uint32 lobit(uint32 u)
-		{return u & 0x00000001UL;}
-	static uint32 lobits(uint32 u)
-		{return u & 0x7fffffffUL;}
-	static uint32 mixbits(uint32 u, uint32 v)
-		{return hibit(u) | lobits(v);}
-	static uint32 twist(uint32 m, uint32 s0, uint32 s1)
-		{return m ^ (mixbits(s0,s1)>>1) ^ (-lobit(s1) & 0x9908b0dfUL);}
-
-
-@ 
-@<|MersenneTwister| inline method definitions@>=
-	@<|MersenneTwister| constructor code@>;
-	@<|MersenneTwister| copy constructor code@>;
-	@<|MersenneTwister::lrand| code@>;
-	@<|MersenneTwister::drand| code@>;
-	@<|MersenneTwister::seed| code@>;
-	@<|MersenneTwister::refresh| code@>;
-
-@ 
-@<|MersenneTwister| constructor code@>=
-inline MersenneTwister::MersenneTwister(uint32 iseed)
-{
-	seed(iseed);
-}
-
-@ 
-@<|MersenneTwister| copy constructor code@>=
-inline MersenneTwister::MersenneTwister(const MersenneTwister& mt)
-	: stateptr(mt.stateptr)
-{
-	memcpy(statevec, mt.statevec, sizeof(uint32)*STATE_SIZE);
-}
-
-@ 
-@<|MersenneTwister::lrand| code@>=
-inline MersenneTwister::uint32 MersenneTwister::lrand()
-{
-	if (stateptr >= STATE_SIZE)
-		refresh();
-
-	register uint32 v = statevec[stateptr++];
-	v ^= v >> 11;
-	v ^= (v << 7) & 0x9d2c5680;
-	v ^= (v << 15) & 0xefc60000;
-	return (v ^ (v >> 18));
-}
-
-@ 
-@<|MersenneTwister::drand| code@>=
-inline double MersenneTwister::drand()
-{
-	uint32 a = lrand() >> 5;
-	uint32 b = lrand() >> 6;
-	return (a*67108864.0+b) * (1.0/9007199254740992.0);
-}
-
-@ PRNG of D. Knuth 
-@<|MersenneTwister::seed| code@>=
-inline void MersenneTwister::seed(uint32 iseed)
-{
-	statevec[0] = iseed & 0xffffffffUL;
-    for (int i = 1;  i < STATE_SIZE;  i++) {
-		register uint32 val = statevec[i-1] >> 30;
-		val ^= statevec[i-1];
-		val *= 1812433253ul;
-		val += i;
-		statevec[i] = val & 0xffffffffUL;
-	}
-
-	refresh();
-}
-
-@ 
-@<|MersenneTwister::refresh| code@>=
-inline void MersenneTwister::refresh()
-{
-	register uint32* p = statevec;
-	for (int i = STATE_SIZE-RECUR_OFFSET; i--; ++p)
-		*p = twist(p[RECUR_OFFSET], p[0], p[1]);
-	for (int i = RECUR_OFFSET; --i; ++p)
-		*p = twist(p[RECUR_OFFSET-STATE_SIZE], p[0], p[1]);
-	*p = twist(p[RECUR_OFFSET-STATE_SIZE], p[0], statevec[0]);
-
-	stateptr = 0;
-}
-
-@ End of {\tt mersenne\_twister.h} file.

dynare++/kord/normal_conjugate.cweb

-@q $Id$ @>
-@q Copyright 2007, Ondra Kamenik @>
-
-@ Start of {\tt normal\_conjugate.cpp} file.
-
-@c
-
-#include "normal_conjugate.h"
-#include "kord_exception.h"
-
-@<|NormalConj| diffuse prior constructor@>;
-@<|NormalConj| data update constructor@>;
-@<|NormalConj| copy constructor@>;
-@<|NormalConj::update| one observation code@>;
-@<|NormalConj::update| multiple observations code@>;
-@<|NormalConj::update| with |NormalConj| code@>;
-@<|NormalConj::getVariance| code@>;
-
-@ 
-@<|NormalConj| diffuse prior constructor@>=
-NormalConj::NormalConj(int d)
-	: mu(d), kappa(0), nu(-1), lambda(d,d)
-{
-	mu.zeros();
-	lambda.zeros();
-}
-
-@ 
-@<|NormalConj| data update constructor@>=
-NormalConj::NormalConj(const ConstTwoDMatrix& ydata)
-	: mu(ydata.numRows()), kappa(ydata.numCols()), nu(ydata.numCols()-1),
-	  lambda(ydata.numRows(), ydata.numRows())
-{
-	mu.zeros();
-	for (int i = 0; i < ydata.numCols(); i++)
-		mu.add(1.0/ydata.numCols(), ConstVector(ydata, i));
-
-	lambda.zeros();
-	for (int i = 0; i < ydata.numCols(); i++) {
-		Vector diff(ConstVector(ydata, i));
-		diff.add(-1, mu);
-		lambda.addOuter(diff);
-	}
-}
-
-@ 
-@<|NormalConj| copy constructor@>=
-NormalConj::NormalConj(const NormalConj& nc)
-	: mu(nc.mu), kappa(nc.kappa), nu(nc.nu), lambda(nc.lambda)
-{
-}
-
-@ The method performs the following:
-$$\eqalign{
-  \mu_1 = &\; {\kappa_0\over \kappa_0+1}\mu_0 + {1\over \kappa_0+1}y\cr
-  \kappa_1 = &\; \kappa_0 + 1\cr
-  \nu_1 = &\; \nu_0 + 1\cr
-  \Lambda_1 = &\; \Lambda_0 + {\kappa_0\over\kappa_0+1}(y-\mu_0)(y-\mu_0)^T,
-}$$
-
-@<|NormalConj::update| one observation code@>=
-void NormalConj::update(const ConstVector& y)
-{
-	KORD_RAISE_IF(y.length() != mu.length(),
-				  "Wrong length of a vector in NormalConj::update");
-
-	mu.mult(kappa/(1.0+kappa));
-	mu.add(1.0/(1.0+kappa), y);
-
-	Vector diff(y);
-	diff.add(-1, mu);
-	lambda.addOuter(diff, kappa/(1.0+kappa));
-
-	kappa++;
-	nu++;
-}
-
-@ The method evaluates the formula in the header file.
-
-@<|NormalConj::update| multiple observations code@>=
-void NormalConj::update(const ConstTwoDMatrix& ydata)
-{
-	NormalConj nc(ydata);
-	update(nc);
-}
-
-
-@ 
-@<|NormalConj::update| with |NormalConj| code@>=
-void NormalConj::update(const NormalConj& nc)
-{
-	double wold = ((double)kappa)/(kappa+nc.kappa);
-	double wnew = 1-wold;
-
-	mu.mult(wold);
-	mu.add(wnew, nc.mu);
-
-	Vector diff(nc.mu);
-	diff.add(-1, mu);
-	lambda.add(1.0, nc.lambda);
-	lambda.addOuter(diff);
-
-	kappa = kappa + nc.kappa;
-	nu = nu + nc.kappa;
-}
-
-
-@ This returns ${1\over \nu-d-1}\Lambda$, which is the mean of the
-variance in the posterior distribution. If the number of degrees of
-freedom is less than $d$, then NaNs are returned.
-
-@<|NormalConj::getVariance| code@>=
-void NormalConj::getVariance(TwoDMatrix& v) const
-{
-	if (nu > getDim()+1) {
-		v = (const TwoDMatrix&)lambda;
-		v.mult(1.0/(nu-getDim()-1));
-	} else
-		v.nans();
-}
-
-
-@ End of {\tt normal\_conjugate.cpp} file.

dynare++/kord/normal_conjugate.hweb

-@q $Id$ @>
-@q Copyright 2007, Ondra Kamenik @>
-
-@*2 Conjugate family for normal distribution. Start of {\tt
-normal\_conjugate.h} file.
-
-The main purpose here is to implement a class representing conjugate
-distributions for mean and variance of the normal distribution. The
-class has two main methods: the first one is to update itself with
-respect to one observation, the second one is to update itself with
-respect to anothe object of the class. In the both methods, the
-previous state of the class corresponds to the prior distribution, and
-the final state corresponds to the posterior distribution.
-
-The algrebra can be found in Gelman, Carlin, Stern, Rubin (p.87). It
-goes as follows: Prior conjugate distribution takes the following form:
-$$\eqalign{
-  \Sigma \sim& {\rm InvWishart}_{\nu_0}(\Lambda_0^{-1}) \cr
-  \mu\vert\Sigma \sim& N(\mu_0,\Sigma/\kappa_0)
-}$$
-If the observations are $y_1\ldots y_n$, then the posterior distribution has the
-same form with the following parameters:
-$$\eqalign{
-  \mu_n = &\; {\kappa_0\over \kappa_0+n}\mu_0 + {n\over \kappa_0+n}\bar y\cr
-  \kappa_n = &\; \kappa_0 + n\cr
-  \nu_n = &\; \nu_0 + n\cr
-  \Lambda_n = &\; \Lambda_0 + S + {\kappa_0 n\over\kappa_0+n}(\bar y-\mu_0)(\bar y-\mu_0)^T,
-}$$
-where
-$$\eqalign{
-  \bar y = &\; {1\over n}\sum_{i=1}^ny_i\cr
-  S = &\; \sum_{i=1}^n(y_i-\bar y)(y_i-\bar y)^T
-}$$
-
-@s NormalConj int
-
-@c
-#ifndef NORMAL_CONJUGATE_H
-#define NORMAL_CONJUGATE_H
-
-#include "twod_matrix.h"
-
-@<|NormalConj| class declaration@>;
-
-#endif
-
-@ The class is described by the four parameters: $\mu$, $\kappa$, $\nu$ and $\Lambda$.
-
-@<|NormalConj| class declaration@>=
-class NormalConj {
-protected:@;
-	Vector mu;
-	int kappa;
-	int nu;
-	TwoDMatrix lambda;
-public:@;
-	@<|NormalConj| constructors@>;
-	virtual ~NormalConj() @+{}
-	void update(const ConstVector& y);
-	void update(const ConstTwoDMatrix& ydata);
-	void update(const NormalConj& nc);
-	int getDim() const
-	{@+ return mu.length();@+}
-	const Vector& getMean() const
-	{@+ return mu;@+}
-	void getVariance(TwoDMatrix& v) const;
-};
-
-@ We provide the following constructors: The first constructs diffuse
-(Jeffrey's) prior. It sets $\kappa$, and $\Lambda$ to zeros, $nu$ to
-$-1$ and also the mean $\mu$ to zero (it should not be
-referenced). The second constructs the posterior using the diffuse
-prior and the observed data (columnwise). The third is a copy
-constructor.
-
-@<|NormalConj| constructors@>=
-	NormalConj(int d);
-	NormalConj(const ConstTwoDMatrix& ydata);
-	NormalConj(const NormalConj& nc);
-
-
-@ End of {\tt normal\_conjugate.h} file.

dynare++/kord/random.cweb

-@q $Id: random.cweb 1491 2007-12-19 14:36:53Z kamenik $ @>
-@q Copyright 2007, Ondra Kamenik @>
-
-@ Start of {\tt random.cpp} file.
-@c
-
-#include "random.h"
-
-#include <cstdlib>
-#include <limits>
-#include <cmath>
-
-@<|RandomGenerator::int_uniform| code@>;
-@<|RandomGenerator::normal| code@>;
-SystemRandomGenerator system_random_generator;
-@<|SystemRandomGenerator::uniform| code@>;
-@<|SystemRandomGenerator::initSeed| code@>;
-
-@ 
-@<|RandomGenerator::int_uniform| code@>=
-int RandomGenerator::int_uniform()
-{
-	double s = std::numeric_limits<int>::max()*uniform();
-	return (int)s;
-}
-
-@ This implements Marsaglia Polar Method. 
-@<|RandomGenerator::normal| code@>=
-double RandomGenerator::normal()
-{
-	double x1, x2;
-	double w;
-	do {
-		x1 = 2*uniform()-1;
-		x2 = 2*uniform()-1;
-		w = x1*x1 + x2*x2;
-	} while (w >= 1.0 || w < 1.0e-30);
-	return x1*std::sqrt((-2.0*std::log(w))/w);
-}
-
-@ 
-@<|SystemRandomGenerator::uniform| code@>=
-double SystemRandomGenerator::uniform()
-{
-#if !defined(__MINGW32__)
-	return drand48();
-#else
-	return ((double)rand())/RAND_MAX;
-#endif
-}
-
-@ 
-@<|SystemRandomGenerator::initSeed| code@>=
-void SystemRandomGenerator::initSeed(int seed)
-{
-#if !defined(__MINGW32__)
-	srand48(seed);
-#else
-	srand(seed);
-#endif
-}
-
-@ End of {\tt random.cpp} file.

dynare++/kord/random.hweb

-@q $Id: random.hweb 2335 2009-01-14 10:35:21Z kamenik $ @>
-@q Copyright 2007, Ondra Kamenik @>
-
-@*2 Random number generation. Start of {\tt random.h} file.
-
-@s RandomGenerator int
-@s SystemRandomGenerator int
-
-@c
-#ifndef RANDOM_H
-#define RANDOM_H
-
-@<|RandomGenerator| class declaration@>;
-@<|SystemRandomGenerator| class declaration@>;
-extern SystemRandomGenerator system_random_generator;
-
-#endif
-
-@ This is a general interface to an object able to generate random
-numbers. Subclass needs to implement |uniform| method, other is, by
-default, implemented here.
-@<|RandomGenerator| class declaration@>=
-class RandomGenerator {
-public:@;
-	virtual double uniform() = 0;
-	int int_uniform();
-	double normal();
-};
-
-@ This implements |RandomGenerator| interface with system |drand| or
-|rand|. It is not thread aware.
-@<|SystemRandomGenerator| class declaration@>=
-class SystemRandomGenerator : public RandomGenerator {
-public:@;
-	double uniform();
-	void initSeed(int seed);
-};
-
-@ End of {\tt random.h} file.

dynare++/kord/tests.cpp

-/* $Id: tests.cpp 148 2005-04-19 15:12:26Z kamenik $ */
-/* Copyright 2004, Ondra Kamenik */
-
-#include <cstdlib>
-#include "korder.h"
-#include "SylvException.h"
-
-struct Rand {
-	static void init(int n1, int n2, int n3, int n4, int n5);
-	static double get(double m);
-	static int get(int m);
-	static bool discrete(double prob); // answers true with given probability
-};
-
-void Rand::init(int n1, int n2, int n3, int n4, int n5)
-{
-	long int seed = n1;
-	seed = 256*seed+n2;
-	seed = 256*seed+n3;
-	seed = 256*seed+n4;
-	seed = 256*seed+n5;
-	srand48(seed);
-}
-
-double Rand::get(double m)
-{
-	return 2*m*(drand48()-0.5);
-}
-
-int Rand::get(int m)
-{
-	return (int)(Rand::get(0.9999*m));
-}
-
-bool Rand::discrete(double prob)
-{
-	return drand48() < prob;
-}
-
-struct SparseGenerator {
-	static FSSparseTensor* makeTensor(int dim, int nv, int r,
-									  double fill, double m);
-	static void fillContainer(TensorContainer<FSSparseTensor>& c,
-							  int maxdim, int nv, int r, double m);
-};
-
-FSSparseTensor* SparseGenerator::makeTensor(int dim, int nv, int r,
-											double fill, double m)
-{
-	FSSparseTensor* res = new FSSparseTensor(dim, nv, r);
-	FFSTensor dummy(0, nv, dim);
-	for (Tensor::index fi = dummy.begin(); fi != dummy.end(); ++fi) {
-		for (int i = 0; i < r; i++) {
-			if (Rand::discrete(fill)) {
-				double x = Rand::get(m);
-				res->insert(fi.getCoor(), i, x);
-			}
-		}
-	}
-	return res;
-}
-
-void SparseGenerator::fillContainer(TensorContainer<FSSparseTensor>& c,
-									int maxdim, int nv, int r,
-									double m)
-{
-	Rand::init(maxdim, nv, r, (int)(5*m), 0);
-	double fill = 0.5;
-	for (int d = 1; d <= maxdim; d++) {
-		c.insert(makeTensor(d,nv,r,fill,m));
-		fill *= 0.3;
-	}
-}
-
-const double vdata [] = { // 3x3
-	0.1307870268, 0.1241940078, 0.1356703123,
-	0.1241940078, 0.1986920419, 0.2010160581,
-	0.1356703123, 0.2010160581, 0.2160336975
-};
-
-
-const double gy_data [] = { // 8x4
-	0.3985178619, -0.5688233582, 0.9572900437, -0.6606847776, 0.1453004017,
-	0.3025310675, -0.8627437750, -0.6903410191, 0.4751910580, -0.7270018589,
-	-0.0939612498, -0.1463831989, 0.6742110220, 0.6046671043, 0.5215893126,
-	-1.0412969986, -0.3524898417, -1.0986703430, 0.8006531522, 0.8879776376,
-	-0.1037608317, -0.5587378073, -0.1010366945, 0.9462411248, -0.2439199881,
-	1.3420621236, -0.7820285935, 0.3205293447, 0.3606124791, 0.2975422208,
-	-0.5452861965, 1.6320340279 
-};
-
-const double gu_data [] = { // just some numbers, no structure
-	1.8415286914, -0.2638743845, 1.7690713274, 0.9668585956, 0.2303143646, 
-	-0.2229624279, -0.4381991822, 1.0082401405, -0.3186555860, -0.0624691529, 
-	-0.5189085756, 1.4269672156, 0.1163282969, 1.4020183445, -0.0952660426, 
-	0.2099097124, 0.6912400502, -0.5180935114, 0.5288316624, 0.2188053448, 
-	0.5715516767, 0.7813893410, -0.6385073106, 0.8335131513, 0.3605202168, 
-	-1.1167944865, -1.2263750934, 0.6113636081, 0.6964915482, -0.6451217688, 
-	0.4062810500, -2.0552251116, -1.6383406284, 0.0198915095, 0.0111014458, 
-	-1.2421792262, -1.0724161722, -0.4276904972, 0.1801494950, -2.0716473264
-};
-
-const double vdata2 [] = { // 10x10 positive definite
-	0.79666, -0.15536, 0.05667, -0.21026, 0.20262, 0.28505, 0.60341, -0.09703, 0.32363, 0.13299, 
-	-0.15536, 0.64380, -0.01131, 0.00980, 0.03755, 0.43791, 0.21784, -0.31755, -0.55911, -0.29655, 
-	0.05667, -0.01131, 0.56165, -0.34357, -0.40584, 0.20990, 0.28348, 0.20398, -0.19856, 0.35820, 
-	-0.21026, 0.00980, -0.34357, 0.56147, 0.10972, -0.34146, -0.49906, -0.19685, 0.21088, -0.31560, 
-	0.20262, 0.03755, -0.40584, 0.10972, 0.72278, 0.02155, 0.04089, -0.19696, 0.03446, -0.12919, 
-	0.28505, 0.43791, 0.20990, -0.34146, 0.02155, 0.75867, 0.77699, -0.31125, -0.55141, -0.02155, 
-	0.60341, 0.21784, 0.28348, -0.49906, 0.04089, 0.77699, 1.34553, -0.18613, -0.25811, -0.19016, 
-	-0.09703, -0.31755, 0.20398, -0.19685, -0.19696, -0.31125, -0.18613, 0.59470, 0.08386, 0.41750, 
-	0.32363, -0.55911, -0.19856, 0.21088, 0.03446, -0.55141, -0.25811, 0.08386, 0.98917, -0.12992, 
-	0.13299, -0.29655, 0.35820, -0.31560, -0.12919, -0.02155, -0.19016, 0.41750, -0.12992, 0.89608
-};
-
-const double gy_data2 [] = { // 600 items make gy 30x20, whose gy(6:25,:) has spectrum within unit 
-	0.39414, -0.29766, 0.08948, -0.19204, -0.00750, 0.21159, 0.05494, 0.06225, 0.01771, 0.21913, 
-	-0.01373, 0.20086, -0.06086, -0.10955, 0.14424, -0.08390, 0.03948, -0.14713, 0.11674, 0.05091, 
-	0.24039, 0.28307, -0.11835, 0.13030, 0.11682, -0.27444, -0.19311, -0.16654, 0.12867, 0.25116, 
-	-0.19781, 0.45242, -0.15862, 0.24428, -0.11966, 0.11483, -0.32279, 0.29727, 0.20934, -0.18190, 
-	-0.15080, -0.09477, -0.30551, -0.02672, -0.26919, 0.11165, -0.06390, 0.03449, -0.26622, 0.22197, 
-	0.45141, -0.41683, 0.09760, 0.31094, -0.01652, 0.05809, -0.04514, -0.05645, 0.00554, 0.47980, 
-	0.11726, 0.42459, -0.13136, -0.30902, -0.14648, 0.11455, 0.02947, -0.03835, -0.04044, 0.03559, 
-	-0.26575, -0.01783, 0.31243, -0.14412, -0.13218, -0.05080, 0.18576, 0.13840, -0.05560, 0.35530, 
-	-0.25573, -0.11560, 0.15187, -0.18431, 0.08193, -0.32278, 0.17560, -0.05529, -0.10020, -0.23088, 
-	-0.20979, -0.49245, 0.09915, -0.16909, -0.03443, 0.19497, 0.18473, 0.25662, 0.29605, -0.20531, 
-	-0.39244, -0.43369, 0.05588, 0.24823, -0.14236, -0.08311, 0.16371, -0.19975, 0.30605, -0.17087, 
-	-0.01270, 0.00123, -0.22426, -0.13810, 0.05079, 0.06971, 0.01922, -0.09952, -0.23177, -0.41962, 
-	-0.41991, 0.41430, -0.04247, -0.13706, -0.12048, -0.28906, -0.22813, -0.25057, -0.18579, -0.20642, 
-	-0.47976, 0.25490, -0.05138, -0.30794, 0.31651, 0.02034, 0.12954, -0.20110, 0.13336, -0.40775, 
-	-0.30195, -0.13704, 0.12396, 0.28152, 0.02986, 0.27669, 0.24623, 0.08635, -0.11956, -0.02949, 
-	0.37401, 0.20838, 0.24801, -0.26872, 0.11195, 0.00315, -0.19069, 0.12839, -0.23036, -0.48228, 
-	0.08434, -0.39872, -0.28896, -0.28754, 0.24668, 0.23285, 0.25437, 0.10456, -0.14124, 0.20483, 
-	-0.19117, -0.33836, -0.24875, 0.08207, -0.03930, 0.20364, 0.15384, -0.15270, 0.24372, -0.11199, 
-	-0.46591, 0.30319, 0.05745, 0.09084, 0.06058, 0.31884, 0.05071, -0.28899, -0.30793, -0.03566, 
-	0.02286, 0.28178, 0.00736, -0.31378, -0.18144, -0.22346, -0.27239, 0.31043, -0.26228, 0.22181, 
-	-0.15096, -0.36953, -0.06032, 0.21496, 0.29545, -0.13112, 0.16420, -0.07573, -0.43111, -0.43057, 
-	0.26716, -0.31209, -0.05866, -0.29101, -0.27437, -0.18727, 0.28732, -0.19014, 0.08837, 0.30405, 
-	0.06103, -0.35612, 0.00173, 0.25134, -0.08987, -0.22766, -0.03254, -0.18662, -0.08491, 0.49401, 
-	-0.12145, -0.02961, -0.03668, -0.30043, -0.08555, 0.01701, -0.12544, 0.10969, -0.48202, 0.07245, 
-	0.20673, 0.11408, 0.04343, -0.01815, -0.31594, -0.23632, -0.06258, -0.27474, 0.12180, 0.16613, 
-	-0.37931, 0.30219, 0.15765, 0.25489, 0.17529, -0.17020, -0.30060, 0.22058, -0.02450, -0.42143, 
-	0.49642, 0.46899, -0.28552, -0.22549, -0.01333, 0.21567, 0.22251, 0.21639, -0.19194, -0.19140, 
-	-0.24106, 0.10952, -0.11019, 0.29763, -0.02039, -0.25748, 0.23169, 0.01357, 0.09802, -0.19022, 
-	0.37604, -0.40777, 0.18131, -0.10258, 0.29573, -0.31773, 0.09069, -0.02198, -0.26594, 0.48302, 
-	-0.10041, 0.20210, -0.05609, -0.01169, -0.17339, 0.17862, -0.22502, 0.29009, -0.45160, 0.19771, 
-	0.27634, 0.31695, -0.09993, 0.17167, 0.12394, 0.28088, -0.12502, -0.16967, -0.06296, -0.17036, 
-	0.27320, 0.01595, 0.16955, 0.30146, -0.15173, -0.29807, 0.08178, -0.06811, 0.21655, 0.26348, 
-	0.06316, 0.45661, -0.29756, -0.05742, -0.14715, -0.03037, -0.16656, -0.08768, 0.38078, 0.40679, 
-	-0.32779, -0.09106, 0.16107, -0.07301, 0.07700, -0.22694, -0.15692, -0.02548, 0.38749, -0.12203, 
-	-0.02980, -0.22067, 0.00680, -0.23058, -0.29112, 0.23032, -0.16026, 0.23392, -0.09990, 0.03628, 
-	-0.42592, -0.33474, -0.09499, -0.17442, -0.20110, 0.24618, -0.06418, -0.06715, 0.40754, 0.29377, 
-	0.29543, -0.16832, -0.08468, 0.06491, -0.01410, 0.19988, 0.24950, 0.14626, -0.27851, 0.06079, 
-	0.48134, -0.13475, 0.25398, 0.11738, 0.23369, -0.00661, -0.16811, -0.04557, -0.12030, -0.39527, 
-	-0.35760, 0.01840, -0.15941, 0.03290, 0.09988, -0.08307, 0.06644, -0.24637, 0.34112, -0.08026, 
-	0.00951, 0.27656, 0.16247, 0.28217, 0.17198, -0.16389, -0.03835, -0.02675, -0.08032, -0.21045, 
-	-0.38946, 0.23207, 0.10987, -0.31674, -0.28653, -0.27430, -0.29109, -0.00648, 0.38431, -0.38478, 
-	-0.41195, -0.19364, -0.20977, -0.05524, 0.05558, -0.20109, 0.11803, -0.19884, 0.43318, -0.39255, 
-	0.26612, -0.21771, 0.12471, 0.12856, -0.15104, -0.11676, 0.17582, -0.25330, 0.00298, -0.31712, 
-	0.21532, -0.20319, 0.14507, -0.04588, -0.22995, -0.06470, 0.18849, -0.13444, 0.37107, 0.07387, 
-	-0.14008, 0.09896, 0.13727, -0.28417, -0.09461, -0.18703, 0.04080, 0.02343, -0.49988, 0.17993, 
-	0.23189, -0.30581, -0.18334, -0.09667, -0.27699, -0.05998, 0.09118, -0.32453, 0.46251, 0.41500, 
-	-0.45314, -0.00544, 0.08529, 0.29099, -0.00937, -0.31650, 0.26163, 0.14506, 0.37498, -0.16454, 
-	0.35215, 0.31642, -0.09161, -0.31452, -0.04792, -0.04677, -0.19523, 0.27998, 0.05491, 0.44461, 
-	-0.01258, -0.27887, 0.18361, -0.04539, -0.02977, 0.30821, 0.29454, -0.17932, 0.16193, 0.23934, 
-	0.47923, 0.25373, 0.23258, 0.31484, -0.17958, -0.01136, 0.17681, 0.12869, 0.03235, 0.43762, 
-	0.13734, -0.09433, -0.03735, 0.17949, 0.14122, -0.17814, 0.06359, 0.16044, 0.12249, -0.22314, 
-	0.40775, 0.05147, 0.12389, 0.04290, -0.01642, 0.00082, -0.18056, 0.02875, 0.32690, 0.17712, 
-	0.34001, -0.21581, -0.01086, -0.18180, 0.17480, -0.17774, -0.07503, 0.28438, -0.19747, 0.29595, 
-	-0.28002, -0.02073, -0.16522, -0.18234, -0.20565, 0.29620, 0.07502, 0.01429, -0.31418, 0.43693, 
-	-0.12212, 0.11178, -0.28503, 0.04683, 0.00072, 0.05566, 0.18857, 0.26101, -0.38891, -0.21216, 
-	-0.21850, -0.15147, -0.30749, -0.23762, 0.14984, 0.03535, -0.02862, -0.00105, -0.39907, -0.06909, 
-	-0.36094, 0.21717, 0.15930, -0.18924, 0.13741, 0.01039, 0.13613, 0.00659, 0.07676, -0.13711, 
-	0.24285, -0.07564, -0.28349, -0.15658, 0.03135, -0.30909, -0.22534, 0.17363, -0.19376, 0.26038, 
-	0.05546, -0.22607, 0.32420, -0.02552, -0.05400, 0.13388, 0.04643, -0.31535, -0.06181, 0.30237, 
-	-0.04680, -0.29441, 0.12231, 0.03960, -0.01188, 0.01406, 0.25402, 0.03315, 0.25026, -0.10922
-};
-
-const double gu_data2 [] = { // raw data 300 items
-	0.26599, 0.41329, 0.31846, 0.92590, 0.43050, 0.17466, 0.02322, 0.72621, 0.37921, 0.70597, 
-	0.97098, 0.14023, 0.57619, 0.09938, 0.02281, 0.92341, 0.72654, 0.71000, 0.76687, 0.70182, 
-	0.88752, 0.49524, 0.42549, 0.42806, 0.57615, 0.76051, 0.15341, 0.47457, 0.60066, 0.40880, 
-	0.20668, 0.41949, 0.97620, 0.94318, 0.71491, 0.56402, 0.23553, 0.94387, 0.78567, 0.06362, 
-	0.85252, 0.86262, 0.25190, 0.03274, 0.93216, 0.37971, 0.08797, 0.14596, 0.73871, 0.06574, 
-	0.67447, 0.28575, 0.43911, 0.92133, 0.12327, 0.87762, 0.71060, 0.07141, 0.55443, 0.53310, 
-	0.91529, 0.25121, 0.07593, 0.94490, 0.28656, 0.82174, 0.68887, 0.67337, 0.99291, 0.03316, 
-	0.02849, 0.33891, 0.25594, 0.90071, 0.01248, 0.67871, 0.65953, 0.65369, 0.97574, 0.31578, 
-	0.23678, 0.39220, 0.06706, 0.80943, 0.57694, 0.08220, 0.18151, 0.19969, 0.37096, 0.37858, 
-	0.70153, 0.46816, 0.76511, 0.02520, 0.39387, 0.25527, 0.39050, 0.60141, 0.30322, 0.46195, 
-	0.12025, 0.33616, 0.04174, 0.00196, 0.68886, 0.74445, 0.15869, 0.18994, 0.95195, 0.62874, 
-	0.82874, 0.53369, 0.34383, 0.50752, 0.97023, 0.22695, 0.62407, 0.25840, 0.71279, 0.28785, 
-	0.31611, 0.20391, 0.19702, 0.40760, 0.85158, 0.68369, 0.63760, 0.09879, 0.11924, 0.32920, 
-	0.53052, 0.15900, 0.21229, 0.84080, 0.33933, 0.93651, 0.42705, 0.06199, 0.50092, 0.47192, 
-	0.57152, 0.01818, 0.31404, 0.50173, 0.87725, 0.50530, 0.10717, 0.04035, 0.32901, 0.33538, 
-	0.04780, 0.40984, 0.78216, 0.91288, 0.11314, 0.25248, 0.23823, 0.74001, 0.48089, 0.55531, 
-	0.82486, 0.01058, 0.05409, 0.44357, 0.52641, 0.68188, 0.94629, 0.61627, 0.33037, 0.11961, 
-	0.57988, 0.19653, 0.91902, 0.59838, 0.52974, 0.28364, 0.45767, 0.65836, 0.63045, 0.76140, 
-	0.27918, 0.27256, 0.46035, 0.77418, 0.92918, 0.14095, 0.89645, 0.25146, 0.21172, 0.47910, 
-	0.95451, 0.34377, 0.29927, 0.79220, 0.97654, 0.67591, 0.44385, 0.38434, 0.44860, 0.28170, 
-	0.90712, 0.20337, 0.00292, 0.55046, 0.62255, 0.45127, 0.80896, 0.43965, 0.59145, 0.23801, 
-	0.33601, 0.30119, 0.89935, 0.40850, 0.98226, 0.75430, 0.68318, 0.65407, 0.68067, 0.32942, 
-	0.11756, 0.27626, 0.83879, 0.72174, 0.75430, 0.13702, 0.03402, 0.58781, 0.07393, 0.23067, 
-	0.92537, 0.29445, 0.43437, 0.47685, 0.54548, 0.66082, 0.23805, 0.60208, 0.94337, 0.21363, 
-	0.72637, 0.57181, 0.77679, 0.63931, 0.72860, 0.38901, 0.94920, 0.04535, 0.12863, 0.40550, 
-	0.90095, 0.21418, 0.13953, 0.99639, 0.02526, 0.70018, 0.21828, 0.20294, 0.20191, 0.30954, 
-	0.39490, 0.68955, 0.11506, 0.15748, 0.40252, 0.91680, 0.61547, 0.78443, 0.19693, 0.67630, 
-	0.56552, 0.58556, 0.53554, 0.53507, 0.09831, 0.21229, 0.83135, 0.26375, 0.89287, 0.97069, 
-	0.70615, 0.42041, 0.43117, 0.21291, 0.26086, 0.26978, 0.77340, 0.43833, 0.46179, 0.54418, 
-	0.67878, 0.42776, 0.61454, 0.55915, 0.36363, 0.31999, 0.42442, 0.86649, 0.62513, 0.02047
-};
-
-class TestRunnable {
-	char name[100];
-public:
-	int dim; // dimension of the solved problem
-	int nvar; // number of variable of the solved problem
-	TestRunnable(const char* n, int d, int nv)
-		: dim(d), nvar(nv)
-		{strncpy(name, n, 100);}
-	bool test() const;
-	virtual bool run() const =0;
-	const char* getName() const
-		{return name;}
-protected:
-	static double korder_unfold_fold(int maxdim, int unfold_dim,
-									 int nstat, int npred, int nboth, int forw,
-									 const TwoDMatrix& gy, const TwoDMatrix& gu,
-									 const TwoDMatrix& v);
-};
-
-
-bool TestRunnable::test() const
-{
-	printf("Running test <%s>\n",name);
-	clock_t start = clock();
-	bool passed = run();
-	clock_t end = clock();
-	printf("CPU time %8.4g (CPU seconds)..................",
-		   ((double)(end-start))/CLOCKS_PER_SEC);
-	if (passed) {
-		printf("passed\n\n");
-		return passed;
-	} else {
-		printf("FAILED\n\n");
-		return passed;
-	}
-}
-
-double TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
-										int nstat, int npred, int nboth, int nforw,
-										const TwoDMatrix& gy, const TwoDMatrix& gu,
-										const TwoDMatrix& v)
-{
-	TensorContainer<FSSparseTensor> c(1);
-	int ny = nstat+npred+nboth+nforw;
-	int nu = v.nrows();
-	int nz = nboth+nforw+ny+nboth+npred+nu;
-	SparseGenerator::fillContainer(c, maxdim, nz, ny, 5.0);
-	for (int d = 1; d <= maxdim; d++) {
-		printf("\ttensor fill for dim=%d is:   %3.2f %%\n",
-			   d, c.get(Symmetry(d))->getFillFactor()*100.0);
-	}
-	Journal jr("out.txt");
-	KOrder kord(nstat, npred, nboth, nforw, c, gy, gu, v, jr);
-	// perform unfolded steps until unfold_dim
-	double maxerror = 0.0;
-	for (int d = 2; d <= unfold_dim; d++) {
-		clock_t pertime = clock();
-		kord.performStep<KOrder::unfold>(d);
-		pertime = clock()-pertime;
-		printf("\ttime for unfolded step dim=%d: %8.4g\n",
-			   d, ((double)(pertime))/CLOCKS_PER_SEC);
-		clock_t checktime = clock();
-		double err = kord.check<KOrder::unfold>(d);
-		checktime = clock()-checktime;
-		printf("\ttime for step check dim=%d:    %8.4g\n",
-			   d, ((double)(checktime))/CLOCKS_PER_SEC);
-		printf("\tmax error in step dim=%d:      %10.6g\n",
-			   d, err);
-		if (maxerror < err)
-			maxerror = err;
-	}
-	// perform folded steps until maxdim
-	if (unfold_dim < maxdim) {
-		clock_t swtime = clock();
-		kord.switchToFolded();
-		swtime = clock()-swtime;
-		printf("\ttime for switching dim=%d:   %8.4g\n",
-			   unfold_dim, ((double)(swtime))/CLOCKS_PER_SEC);
-		
-		for (int d = unfold_dim+1; d <= maxdim; d++) {
-			clock_t pertime = clock();
-			kord.performStep<KOrder::fold>(d);
-			pertime = clock()-pertime;
-			printf("\ttime for folded step dim=%d: %8.4g\n",
-				   d, ((double)(pertime))/CLOCKS_PER_SEC);
-			clock_t checktime = clock();
-			double err = kord.check<KOrder::fold>(d);
-			checktime = clock()-checktime;
-			printf("\ttime for step check dim=%d:    %8.4g\n",
-				   d, ((double)(checktime))/CLOCKS_PER_SEC);
-			printf("\tmax error in step dim=%d:      %10.6g\n",
-				   d, err);
-			if (maxerror < err)
-				maxerror = err;
-		}
-	}
-	return maxerror;
-}
-
-class UnfoldKOrderSmall : public TestRunnable {
-public:
-	UnfoldKOrderSmall()
-		: TestRunnable("unfold-3 fold-4 korder (stat=2,pred=3,both=1,forw=2,u=3,dim=4)",
-					   4, 18) {}
-
-	bool run() const
-		{
-			TwoDMatrix gy(8, 4, gy_data);
-			TwoDMatrix gu(8, 3, gu_data);
-			TwoDMatrix v(3, 3, vdata);
-			double err = korder_unfold_fold(4, 3, 2, 3, 1, 2,
-											gy, gu, v);
-
-			return err < 0.08;
-		}
-};
-
-// same dimension as Smets & Wouters
-class UnfoldKOrderSW : public TestRunnable {
-public:
-	UnfoldKOrderSW()
-		: TestRunnable("unfold S&W korder (stat=5,pred=12,both=8,forw=5,u=10,dim=4)",
-					   4, 73) {}
-
-	bool run() const
-		{
-			TwoDMatrix gy(30, 20, gy_data2);
-			TwoDMatrix gu(30, 10, gu_data2);
-			TwoDMatrix v(10, 10, vdata2);
-			v.mult(0.001);
-			gu.mult(.01);
-			double err = korder_unfold_fold(4, 4, 5, 12, 8, 5,
-											gy, gu, v);
-
-			return err < 0.08;
-		}
-};
-
-class UnfoldFoldKOrderSW : public TestRunnable {
-public:
-	UnfoldFoldKOrderSW()
-		: TestRunnable("unfold-2 fold-3 S&W korder (stat=5,pred=12,both=8,forw=5,u=10,dim=3)",
-					   4, 73) {}
-
-	bool run() const
-		{
-			TwoDMatrix gy(30, 20, gy_data2);
-			TwoDMatrix gu(30, 10, gu_data2);
-			TwoDMatrix v(10, 10, vdata2);
-			v.mult(0.001);
-			gu.mult(.01);
-			double err = korder_unfold_fold(4, 3, 5, 12, 8, 5,
-											gy, gu, v);
-
-			return err < 0.08;
-		}
-};
-
-int main()
-{
-	TestRunnable* all_tests[50];
-	// fill in vector of all tests
-	int num_tests = 0;
-	all_tests[num_tests++] = new UnfoldKOrderSmall();
-	all_tests[num_tests++] = new UnfoldKOrderSW();
-	all_tests[num_tests++] = new UnfoldFoldKOrderSW();
-
-	// find maximum dimension and maximum nvar
-	int dmax=0;
-	int nvmax = 0;
-	for (int i = 0; i < num_tests; i++) {
-		if (dmax < all_tests[i]->dim)
-			dmax = all_tests[i]->dim;
-		if (nvmax < all_tests[i]->nvar)
-			nvmax = all_tests[i]->nvar;
-	}
-	tls.init(dmax, nvmax); // initialize library
-
-	// launch the tests
-	int success = 0;
-	for (int i = 0; i < num_tests; i++) {
-		try {
-			if (all_tests[i]->test())
-				success++;
-		} catch (const TLException& e) {
-			printf("Caugth TL 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;
-}

dynare++/parser/cc/Makefile.am

-noinst_LIBRARIES = libparser.a
-
-GENERATED_FILES = assign_tab.cc csv_tab.cc formula_tab.cc matrix_tab.cc namelist_tab.cc assign_tab.hh csv_tab.hh formula_tab.hh matrix_tab.hh namelist_tab.hh assign_ll.cc csv_ll.cc formula_ll.cc matrix_ll.cc namelist_ll.cc
-
-libparser_a_SOURCES = \
-	location.h \
-	namelist.h \
-	atom_assignings.cpp \
-	atom_assignings.h \
-	atom_substitutions.cpp \
-	atom_substitutions.h \
-	csv_parser.cpp \
-	csv_parser.h \
-	dynamic_atoms.cpp \
-	dynamic_atoms.h \
-	fine_atoms.cpp \
-	fine_atoms.h \
-	formula_parser.cpp \
-	formula_parser.h \
-	matrix_parser.cpp \
-	matrix_parser.h \
-	parser_exception.cpp \
-	parser_exception.h \
-	static_atoms.cpp \
-	static_atoms.h \
-	static_fine_atoms.cpp \
-	static_fine_atoms.h \
-	tree.cpp \
-	tree.h \
-	$(GENERATED_FILES)
-
-libparser_a_CPPFLAGS = -I../.. $(BOOST_CPPFLAGS)
-
-BUILT_SOURCES = $(GENERATED_FILES)
-
-EXTRA_DIST = assign.y csv.y formula.y matrix.y namelist.y assign.lex csv.lex formula.lex matrix.lex namelist.lex
-
-%_tab.cc %_tab.hh: %.y
-	$(YACC) -d -o$*_tab.cc $<
-
-%_ll.cc: %.lex
-	$(LEX) -i -o$@ $<

dynare++/parser/cc/assign.lex

-%{
-#include "location.h"
-#include "assign_tab.hh"
-
-	extern YYLTYPE asgn_lloc;
-
-#define YY_USER_ACTION SET_LLOC(asgn_);
-%}
-
-%option nounput
-%option noyy_top_state
-%option stack
-%option yylineno
-%option prefix="asgn_"
-%option never-interactive
-%x CMT
-
-%%
-
- /* comments */
-<*>"/*"            {yy_push_state(CMT);}
-<CMT>[^*\n]*
-<CMT>"*"+[^*/\n]*
-<CMT>"*"+"/"       {yy_pop_state();}
-<CMT>[\n]
-"//".*\n
-
- /* spaces */
-[ \t\r\n]          {return BLANK;}
-
- /* names */
-[A-Za-z_][A-Za-z0-9_]* {
-	asgn_lval.string = asgn_text;
-	return NAME;
-}
-
-;                  {return SEMICOLON;}
-=                  {return EQUAL_SIGN;}
-. {
-	asgn_lval.character = asgn_text[0];
-	return CHARACTER;
-}
-
-%%
-
-int asgn_wrap()
-{
-	return 1;
-}
-
-void asgn__destroy_buffer(void* p)
-{
-	asgn__delete_buffer((YY_BUFFER_STATE)p);
-}

dynare++/parser/cc/assign.y

-%{
-/* Copyright (C) 2006-2011, Ondra Kamenik */
-
-#include "location.h"
-#include "atom_assignings.h"
-#include "assign_tab.hh"
-
-#include <stdio.h>
-
-	void asgn_error(const char*);
-	int asgn_lex(void);
-	extern int asgn_lineno;
-	extern ogp::AtomAssignings* aparser;
-
-%}
-
-%union {
-	int integer;
-	char *string;
-	char character;
-}
-
-%token EQUAL_SIGN SEMICOLON CHARACTER BLANK
-%token <string> NAME;
-
-%name-prefix="asgn_"
-
-%locations
-%error-verbose
-
-%%
-
-root : assignments | ;
-
-assignments : assignments BLANK | assignments assignment | assignment | BLANK;
-
-assignment : NAME EQUAL_SIGN material SEMICOLON {
-	aparser->add_assignment(@1.off, $1, @1.ll, @3.off-@1.off, @3.ll + @4.ll);}
-  | NAME space EQUAL_SIGN material SEMICOLON {
-	aparser->add_assignment(@1.off, $1, @1.ll, @4.off-@1.off, @4.ll + @5.ll);}
-  ;
-
-material : material CHARACTER | material NAME | material BLANK | NAME | CHARACTER | BLANK;
-
-space : space BLANK | BLANK;
-
-%%
-
-void asgn_error(const char* mes)
-{
-	aparser->error(mes);
-}

dynare++/parser/cc/atom_assignings.cpp

-// Copyright (C) 2006, Ondra Kamenik
-
-// $Id: atom_assignings.cpp 92 2007-04-19 11:38:21Z ondra $
-
-#include "atom_assignings.h"
-#include "location.h"
-#include "parser_exception.h"
-
-#include "utils/cc/exception.h"
-
-#include <limits>
-#include <iostream>
-
-using namespace ogp;
-
-AtomAssignings::AtomAssignings(const AtomAssignings& aa, ogp::StaticAtoms& a)
-	: atoms(a), expr(aa.expr, atoms), left_names(aa.left_names),
-	  order(aa.order)
-{
-	// fill the lname2expr
-	for (Tvarintmap::const_iterator it = aa.lname2expr.begin();
-		 it != aa.lname2expr.end(); ++it)
-		lname2expr.insert(Tvarintmap::value_type(left_names.query((*it).first), (*it).second));
-}
-
-/** A global symbol for passing info to the AtomAssignings from
- * asgn_parse(). */
-AtomAssignings* aparser;
-
-/** The declaration of functions defined in asgn_ll.cc and asgn_tab.cc
- * generated from assign.lex assign.y */
-void* asgn__scan_buffer(char*, size_t);
-void asgn__destroy_buffer(void*);
-void asgn_parse();
-extern location_type asgn_lloc;
-
-void AtomAssignings::parse(int length, const char* stream)
-{
-	char* buffer = new char[length+2];
-	strncpy(buffer, stream, length);
-	buffer[length] = '\0';
-	buffer[length+1] = '\0';
-	asgn_lloc.off = 0;
-	asgn_lloc.ll = 0;
-	void* p = asgn__scan_buffer(buffer, (unsigned int)length+2);
-	aparser = this;
-	asgn_parse();
-	delete [] buffer;
-	asgn__destroy_buffer(p);
-}
-
-void AtomAssignings::error(const char* mes)
-{
-	throw ParserException(mes, asgn_lloc.off);
-}
-
-void AtomAssignings::add_assignment_to_double(const char* name, double val)
-{
-	// if left hand side is a registered atom, insert it to tree
-	int t;
-	try {
-		if (atoms.check(name))
-			t = expr.add_nulary(name);
-		else
-			t = -1;
-	} catch (const ParserException& e) {
-		t = -1;
-	}
-	// register left hand side in order
-	order.push_back(t);
-
-	// add the double to the tree
-	char tmp[100];
-	sprintf(tmp, "%30.25g", val);
-	try {
-		expr.parse(strlen(tmp), tmp);
-	} catch (const ParserException& e) {
-		// should never happen
-		throw ParserException(string("Error parsing double ")+tmp+": "+e.message(), 0);
-	}
-
-	// register name of the left hand side and put to lname2expr
-	const char* ss = left_names.insert(name);
-	lname2expr.insert(Tvarintmap::value_type(ss, order.size()-1));
-}
-
-void AtomAssignings::add_assignment(int asgn_off, const char* str, int name_len,
-									int right_off, int right_len)
-{
-	// the order of doing things here is important: since the
-	// FormulaParser requires that all references from the i-th tree
-	// refere to trees with index lass than i, so to capture also a
-	// nulary term for the left hand side, it must be inserted to the
-	// expression tree before the expression is parsed.
-
-	// find the name in the atoms, make copy of name to be able to put
-	// '\0' at the end
-	char* buf = new char[name_len+1];
-	strncpy(buf, str, name_len);
-	buf[name_len] = '\0';
-	// if left hand side is a registered atom, insert it to tree
-	int t;
-	try {
-		t = atoms.check(buf);
-		if (t == -1)
-			t = expr.add_nulary(buf);
-	} catch (const ParserException& e) {
-		atoms.register_name(buf);
-		t = expr.add_nulary(buf);
-	}
-	// register left hand side in order
-	order.push_back(t);
-
-	// parse expression on the right
-	try {
-		expr.parse(right_len, str+right_off);
-	} catch (const ParserException& e) {
-		throw ParserException(e, asgn_off+right_off);
-	}
-
-	// register name of the left hand side and put to lname2expr
-	const char* ss = left_names.insert(buf);
-	if (lname2expr.find(ss) != lname2expr.end()) {
-		// Prevent the occurrence of #415
-		std::cerr << "Changing the value of " << ss << " is not supported. Aborting." << std::endl;
-		exit(EXIT_FAILURE);
-	}
-	lname2expr[ss] = order.size()-1;
-
-	// delete name
-	delete [] buf;
-}
-
-void AtomAssignings::apply_subst(const AtomSubstitutions::Toldnamemap& mm)
-{
-	// go through all old variables and see what are their derived new
-	// variables
-	for (AtomSubstitutions::Toldnamemap::const_iterator it = mm.begin();
-		 it != mm.end(); ++it) {
-		const char* oldname = (*it).first;
-		const AtomSubstitutions::Tshiftnameset& sset = (*it).second;
-		if (! sset.empty()) {
-			int told = atoms.index(oldname);
-			if (told < 0 && ! atoms.get_name_storage().query(oldname))
-				atoms.register_name(oldname);
-			if (told == -1)
-				told = expr.add_nulary(oldname);
-			// at least one substitution here, so make an expression
-			expr.add_formula(told);
-			// say that this expression is not assigned to any atom
-			order.push_back(-1);
-			// now go through all new names derived from the old name and
-			// reference to the newly added formula
-			for (AtomSubstitutions::Tshiftnameset::const_iterator itt = sset.begin();
-				 itt != sset.end(); ++itt) {
-				const char* newname = (*itt).first;
-				const char* nn = left_names.insert(newname);
-				lname2expr.insert(Tvarintmap::value_type(nn, expr.nformulas()-1));
-			}
-		}
-	}
-}
-
-void AtomAssignings::print() const
-{
-	printf("Atom Assignings\nExpressions:\n");
-	expr.print();
-	printf("Left names:\n");
-	for (Tvarintmap::const_iterator it = lname2expr.begin();
-		 it != lname2expr.end(); ++it)
-		printf("%s ==> %d (t=%d)\n", (*it).first, expr.formula((*it).second), order[(*it).second]);
-}
-
-void AtomAsgnEvaluator::setValues(EvalTree& et) const
-{
-	// set values of constants
-	aa.atoms.setValues(et);
-
-	// set values of variables to NaN or to user set values
-	double nan = std::numeric_limits<double>::quiet_NaN();
-	for (int i = 0; i < aa.atoms.nvar(); i++) {
-		const char* ss = aa.atoms.name(i);
-		int t = aa.atoms.index(ss);
-		if (t >= 0) {
-			Tusrvalmap::const_iterator it = user_values.find(t);
-			if (it == user_values.end())
-				et.set_nulary(t, nan);
-			else
-				et.set_nulary(t, (*it).second);
-		}
-	}
-}
-
-void AtomAsgnEvaluator::set_user_value(const char* name, double val)
-{
-	int t = aa.atoms.index(name);
-	if (t >= 0) {
-		Tusrvalmap::iterator it = user_values.find(t);
-		if (it == user_values.end())
-			user_values.insert(Tusrvalmap::value_type(t, val));
-		else
-			(*it).second = val;
-	}
-}
-
-void AtomAsgnEvaluator::load(int i, double res)
-{
-	// set the value
-	operator[](i) = res;
-	// if i-th expression is atom, set its value to this EvalTree
-	int t = aa.order[i];
-	if (t >= 0)
-		etree.set_nulary(t, res);
-}
-
-double AtomAsgnEvaluator::get_value(const char* name) const
-{
-	AtomAssignings::Tvarintmap::const_iterator it = aa.lname2expr.find(name);
-	if (it == aa.lname2expr.end())
-		return std::numeric_limits<double>::quiet_NaN();
-	else
-		return operator[]((*it).second);
-}
-
-

dynare++/parser/cc/atom_assignings.h

-// Copyright (C) 2006, Ondra Kamenik
-
-// $Id: atom_assignings.h 149 2007-04-30 02:11:46Z okamenik $
-
-#ifndef OGP_ATOM_ASSIGNINGS_H
-#define OGP_ATOM_ASSIGNINGS_H
-
-#include "static_atoms.h"
-#include "formula_parser.h"
-#include "atom_substitutions.h"
-
-#include <vector>
-#include <map>
-
-namespace ogp {
-
-	class AtomAsgnEvaluator;
-
-	/** This class represents atom assignments used in parameters
-	 * settings and initval initialization. It maintains atoms of the
-	 * all expressions on the right hand side, the parsed formulas of
-	 * the right hand sides, and the information about the left hand
-	 * sides. See documentation to the order member below. */
-	class AtomAssignings {
-		friend class AtomAsgnEvaluator;
-	protected:
-		typedef std::map<const char*, int, ltstr> Tvarintmap;
-		/** All atoms which should be sufficient for formulas at the
-		 * right hand sides. The atoms should be filled with names
-		 * (preregistered). This is a responsibility of the caller. */
-		StaticAtoms& atoms;
-		/** The formulas of right hand sides. */
-		FormulaParser expr;
-		/** Name storage of the names from left hand sides. */
-		NameStorage left_names;
-		/** Information on left hand sides. This maps a name to the
-		 * index of its assigned expression in expr. More than one
-		 * name may reference to the same expression. */
-		Tvarintmap lname2expr;
-		/** Information on left hand sides. If order[i] >= 0, then it
-		 * says that i-th expression in expr is assigned to atom with
-		 * order[i] tree index. */
-		std::vector<int> order;
-	public:
-		/** Construct the object using the provided static atoms. */
-		AtomAssignings(StaticAtoms& a) : atoms(a), expr(atoms)
-			{}
-		/** Make a copy with provided reference to (posibly different)
-		 * static atoms. */
-		AtomAssignings(const AtomAssignings& aa, StaticAtoms& a);
-		virtual ~AtomAssignings()
-			{}
-		/** Parse the assignments from the given string. */
-		void parse(int length, const char* stream);
-		/** Process a syntax error from bison. */
-		void error(const char* mes);
-		/** Add an assignment of the given name to the given
-		 * double. Can be called by a user, anytime. */
-		void add_assignment_to_double(const char* name, double val);
-		/** Add an assignment. Called from assign.y. */
-		void add_assignment(int asgn_off, const char* str, int name_len,
-							int right_off, int right_len);
-		/** This applies old2new map (possibly from atom
-		 * substitutions) to this object. It registers new variables
-		 * in the atoms, and adds the expressions to expr, and left
-		 * names to lname2expr. The information about dynamical part
-		 * of substitutions is ignored, since we are now in the static
-		 * world. */
-		void apply_subst(const AtomSubstitutions::Toldnamemap& mm);
-		/** Debug print. */
-		void print() const;
-	};
-
-	/** This class basically evaluates the atom assignments
-	 * AtomAssignings, so it inherits from ogp::FormulaEvaluator. It
-	 * is also a storage for the results of the evaluation stored as a
-	 * vector, so the class inherits from std::vector<double> and
-	 * ogp::FormulaEvalLoader. As the expressions for atoms are
-	 * evaluated, the results are values for atoms which will be
-	 * used in subsequent evaluations. For this reason, the class
-	 * inherits also from AtomValues. */
-	class AtomAsgnEvaluator : public FormulaEvalLoader,
-							  public AtomValues,
-							  protected FormulaEvaluator,
-							  public std::vector<double> {
-	protected:
-		typedef std::map<int, double> Tusrvalmap;
-		Tusrvalmap user_values;
-		const AtomAssignings& aa;
-	public:
-		AtomAsgnEvaluator(const AtomAssignings& a)
-			: FormulaEvaluator(a.expr),
-			  std::vector<double>(a.expr.nformulas()), aa(a) {}
-		virtual ~AtomAsgnEvaluator() {}
-		/** This sets all initial values to NaNs, all constants and
-		 * all values set by user by call set_value. This is called by
-		 * FormulaEvaluator::eval() method, which is called by eval()
-		 * method passing this argument as AtomValues. So the
-		 * ogp::EvalTree will be always this->etree. */
-		void setValues(EvalTree& et) const;
-		/** User setting of the values. For example in initval,
-		 * parameters are known and should be set to their values. In
-		 * constrast endogenous variables are set initially to NaNs by
-		 * AtomValues::setValues. */
-		void set_user_value(const char* name, double val);
-		/** This sets the result of i-th expression in aa to res, and
-		 * also checks whether the i-th expression is an atom. If so,
-		 * it sets the value of the atom in ogp::EvalTree
-		 * this->etree. */
-		void load(int i, double res);
-		/** After the user values have been set, the assignments can
-		 * be evaluated. For this purpose we have eval() method. The
-		 * result is that this object as std::vector<double> will
-		 * contain the values. It is ordered given by formulas in
-		 * expr. */
-		void eval()
-			{FormulaEvaluator::eval(*this, *this);}
-		/** This returns a value for a given name. If the name is not
-		 * found among atoms, or there is no assignment for the atom,
-		 * NaN is returned. */
-		double get_value(const char* name) const;
-	};
-
-};
-
-#endif
-
-// Local Variables:
-// mode:C++
-// End:

dynare++/parser/cc/atom_substitutions.cpp

-// Copyright (C) 2006, Ondra Kamenik
-
-// $Id: atom_substitutions.cpp 42 2007-01-22 21:53:24Z ondra $
-
-#include "atom_substitutions.h"
-#include "utils/cc/exception.h"
-
-using namespace ogp;
-
-AtomSubstitutions::AtomSubstitutions(const AtomSubstitutions& as, const FineAtoms& oa,
-									 FineAtoms& na)
-	: old_atoms(oa), new_atoms(na)
-{
-	const NameStorage& ns = na.get_name_storage();
-
-	// fill new2old
-	for (Tshiftmap::const_iterator it = as.new2old.begin();
-		 it != as.new2old.end(); ++it)
-		new2old.insert(Tshiftmap::value_type(ns.query((*it).first),
-											 Tshiftname(ns.query((*it).second.first),
-														(*it).second.second)));
-	// fill old2new
-	for (Toldnamemap::const_iterator it = as.old2new.begin();
-		 it != as.old2new.end(); ++it) {
-		Tshiftnameset sset;
-		for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-			 itt != (*it).second.end(); ++itt)
-			sset.insert(Tshiftname(ns.query((*itt).first), (*itt).second));
-		old2new.insert(Toldnamemap::value_type(ns.query((*it).first), sset));
-	}
-}
-
-
-void AtomSubstitutions::add_substitution(const char* newname, const char* oldname, int tshift)
-{
-	// make sure the storage is from the new_atoms
-	newname = new_atoms.get_name_storage().query(newname);
-	oldname = new_atoms.get_name_storage().query(oldname);
-	if (newname == NULL || oldname == NULL)
-		throw ogu::Exception(__FILE__,__LINE__,
-							 "Bad newname or oldname in AtomSubstitutions::add_substitution");
-
-	// insert to new2old map
-	new2old.insert(Tshiftmap::value_type(newname, Tshiftname(oldname, tshift)));
-	// insert to old2new map
-	Toldnamemap::iterator it = old2new.find(oldname);
-	if (it != old2new.end())
-		(*it).second.insert(Tshiftname(newname, -tshift));
-	else {
-		Tshiftnameset snset;
-		snset.insert(Tshiftname(newname, -tshift));
-		old2new.insert(Toldnamemap::value_type(oldname, snset));
-	}
-
-	// put to info
-	info.num_substs++;
-}
-
-void AtomSubstitutions::substitutions_finished(VarOrdering::ord_type ot)
-{
-	// create an external ordering of new_atoms from old_atoms
-	const vector<const char*>& oa_ext = old_atoms.get_allvar();
-	vector<const char*> na_ext;
-	for (unsigned int i = 0; i < oa_ext.size(); i++) {
-		const char* oname = oa_ext[i];
-		// add the old name itself
-		na_ext.push_back(oname);
-		// add all new names derived from the old name
-		Toldnamemap::const_iterator it = old2new.find(oname);
-		if (it != old2new.end())
-			for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-				 itt != (*it).second.end(); ++itt)
-				na_ext.push_back((*itt).first);
-	}
-
-	// call parsing finished for the new_atoms
-	new_atoms.parsing_finished(ot, na_ext);
-}
-
-const char* AtomSubstitutions::get_new4old(const char* oldname, int tshift) const
-{
-	Toldnamemap::const_iterator it = old2new.find(oldname);
-	if (it != old2new.end()) {
-		const Tshiftnameset& sset = (*it).second;
-		for (Tshiftnameset::const_iterator itt = sset.begin();
-			 itt != sset.end(); ++itt)
-			if ((*itt).second == - tshift)
-				return (*itt).first;
-	}
-	return NULL;
-}
-
-void AtomSubstitutions::print() const
-{
-	printf("Atom Substitutions:\nOld ==> New:\n");
-	for (Toldnamemap::const_iterator it = old2new.begin(); it != old2new.end(); ++it)
-		for (Tshiftnameset::const_iterator itt = (*it).second.begin();
-			 itt != (*it).second.end(); ++itt)
-			printf("    %s ==> [%s, %d]\n", (*it).first, (*itt).first, (*itt).second);
-
-	printf("Old <== New:\n");
-	for (Tshiftmap::const_iterator it = new2old.begin(); it != new2old.end(); ++it)
-		printf("    [%s, %d] <== %s\n", (*it).second.first, (*it).second.second, (*it).first);
-}
-
-void SAtoms::substituteAllLagsAndLeads(FormulaParser& fp, AtomSubstitutions& as)
-{
-	const char* name;
-
-	int mlead, mlag;
-	endovarspan(mlead, mlag);
-
-	// substitute all endo lagged more than 1
-	while (NULL != (name = findEndoWithLeadInInterval(mlag, -2)))
-		makeAuxVariables(name, -1, -2, mlag, fp, as);
-	// substitute all endo leaded more than 1
-	while (NULL != (name = findEndoWithLeadInInterval(2, mlead)))
-		makeAuxVariables(name, 1, 2, mlead, fp, as);
-
-	exovarspan(mlead, mlag);
-	
-	// substitute all lagged exo
-	while (NULL != (name = findExoWithLeadInInterval(mlag, -1)))
-		makeAuxVariables(name, -1, -1, mlag, fp, as);
-	// substitute all leaded exo
-	while (NULL != (name = findExoWithLeadInInterval(1, mlead)))
-		makeAuxVariables(name, 1, 1, mlead, fp, as);
-
-	// notify that substitution have been finished
-	as.substitutions_finished(order_type);
-}
-
-void SAtoms::substituteAllLagsAndExo1Leads(FormulaParser& fp, AtomSubstitutions& as)
-{
-	const char* name;
-
-	int mlead, mlag;
-	endovarspan(mlead, mlag);
-
-	// substitute all endo lagged more than 1
-	while (NULL != (name = findEndoWithLeadInInterval(mlag, -2)))
-		makeAuxVariables(name, -1, -2, mlag, fp, as);
-
-	exovarspan(mlead, mlag);
-	
-	// substitute all lagged exo
-	while (NULL != (name = findExoWithLeadInInterval(mlag, -1)))
-		makeAuxVariables(name, -1, -1, mlag, fp, as);
-	// substitute all leaded exo by 1
-	while (NULL != (name = findExoWithLeadInInterval(1,1)))
-		makeAuxVariables(name, 1, 1, 1, fp, as);
-
-	// notify that substitution have been finished
-	as.substitutions_finished(order_type);
-}
-
-const char* SAtoms::findNameWithLeadInInterval(const vector<const char*>& names,
-											   int ll1, int ll2) const
-{
-	for (unsigned int i = 0; i < names.size(); i++) {
-		const char* name = names[i];
-		DynamicAtoms::Tvarmap::const_iterator it = vars.find(name);
-		if (it != vars.end()) {
-			const DynamicAtoms::Tlagmap& lmap = (*it).second;
-			for (DynamicAtoms::Tlagmap::const_iterator itt = lmap.begin();
-				 itt != lmap.end(); ++itt)
-				if ((*itt).first >= ll1 && (*itt).first <= ll2)
-					return name;
-		}
-	}
-
-	// nothing found
-	return NULL;
-}
-
-void SAtoms::attemptAuxName(const char* str, int ll, string& out) const
-{
-	char c = (ll >= 0)? ((ll == 0)? 'e' : 'p' ) : 'm';
-	char absll[100];
-	sprintf(absll, "%d", std::abs(ll));
-	int iter = 1;
-	do {
-		out = string(str) + '_';
-		for (int i = 0; i < iter; i++)
-			out += c;
-		if (ll != 0)
-			out += absll;
-		iter++;
-	} while (varnames.query(out.c_str()));
-}
-
-void SAtoms::makeAuxVariables(const char* name, int step, int start, int limit_lead,
-							  FormulaParser& fp, AtomSubstitutions& as)
-{
-	if (! (step == 1 || step == -1))
-		throw ogu::Exception(__FILE__,__LINE__,
-							 "Wrong value of step in SAtoms::makeAuxVariables");
-	if (step*start > step*limit_lead)
-		throw ogu::Exception(__FILE__,__LINE__,
-							 "Wrong value of start in SAtoms::makeAuxVariables");
-
-	// make sure that we do not go further than necessary, this is
-	// that the limit lead is not behind maxlead or minlag
-	int mlead, mlag;
-	varspan(name, mlead, mlag);
-	if (step == -1)
-		limit_lead = std::max(limit_lead, mlag);
-	else
-		limit_lead = std::min(limit_lead, mlead);
-
-	// Comment to comments: name="a"; start=-3; step=-1;
-
-	char tmp[500];
-
-	// recover tree index of a previous atom, i.e. set tprev to a tree
-	// index of atom "a(-2)"
-	int tprev = index(name, start-step);
-	if (tprev == -1) {
-		sprintf(tmp, "%s(%d)", name, start-step);
-		tprev = fp.add_nulary(tmp);
-	}
-
-	int ll = start;
-	do {
-		// either create atom "a_m2(0)" with tree index taux and add
-		// equation "a_m2(0)=a(-2)"
-		// or 
-        // check if "a_m2(0)" has not been already created (with
-        // different step), in this case do not add equation "a_m2(0)
-        // = a(-2)"
-		const char* newname;
-		string newname_str;
-		int taux;
-		if (NULL == (newname=as.get_new4old(name, ll-step))) {
-			attemptAuxName(name, ll-step, newname_str);
-			newname = newname_str.c_str();
-			register_uniq_endo(newname);
-			newname = varnames.query(newname);
-			sprintf(tmp, "%s(0)", newname);
-			taux = fp.add_nulary(tmp);
-			// add to substitutions
-			as.add_substitution(newname, name, ll-step);
-
-			// add equation "a_m2(0) = a(-2)", this is taux = tprev
-			fp.add_formula(fp.add_binary(MINUS, taux, tprev));
-		} else {
-			// example: exogenous EPS and occurrence at both EPS(-1)
-			// EPS(+1)
-            // first call makeAuxVariables("EPS",1,1,...) will make endo EPS_p0 = EPS
-            // second call makeAuxVariables("EPS",-1,-1,...) will use this EPS_p0
-			//             to substitute for EPS(-1)
-			taux = index(newname, 0);
-			if (taux < 0)
-				throw ogu::Exception(__FILE__,__LINE__,
-									 "Couldn't find tree index of previously substituted variable");
-		}
-			
-		// create atom "a_m2(-1)" or turn "a(-3)" if any to "a_m2(-1)"; tree index t
-		int t = index(name, ll);
-		if (t == -1) {
-			// no "a(-3)", make t <-> a_m2(-1)
-			sprintf(tmp, "%s(%d)", newname, step);
-			t = fp.add_nulary(tmp);
-		} else {
-			// turn a(-3) to a_m2(-1)
-			unassign_variable(name, ll, t);
-			assign_variable(newname, step, t);
-		}
-
-		// next iteration starts with tprev <-> "a_m2(-1)" (this will be made equal to "a_m3(0)")
-		tprev = t;
-		
-		ll += step;
-	} while (step*ll <= step*limit_lead);
-}