Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4565 deletions
dynare++/extern/R/dynareR.cpp deleted 100644 → 0
View file @ 6ca649c4
// $Id: dynareR.cpp 862 2006-08-04 17:34:56Z tamas $
// Copyright 2006, Tamas K Papp
#include "dynare3.h" // Dynare class
#include "approximation.h" // Approximation class
// exceptions
#include "dynare_exception.h"
#include "parser/cc/parser_exception.h"
#include "utils/cc/exception.h"
#include "SylvException.h"
#include "tl_exception.h"
#include "kord_exception.h"
#include <algorithm>
#include <string.h>
#ifdef DEBUG
#include <stdio.h>
#endif
#include <R_ext/Memory.h>
/** This file containt the C glue functions for an R interface to
* Dynare++. Although written in standard C (except for the use of
* R_alloc), the indexing, calling and memory management conventions
* of the functions in this file were tailored for R.
*
* It is not recommended that you use this interface for anything else
* but R.
*/
/** Error codes: these error codes correspond to possible
* exceptions. */
#define DYNARER_SYLVEXCEPTION 1
#define DYNARER_DYNAREEXCEPTION 2
#define DYNARER_OGUEXCEPTION 3
#define DYNARER_TLEXCEPTION 4
#define DYNARER_KORDEXCEPTION 5
#define DYNARER_NAMESMATCHINGERROR 6
/** Copies the message into a buffer. The buffer is allocated and
* managed by R, ie it will be garbage collected after the .C call
* returns and the contents are duplicated.
*/
char *passmessage(const char *errormessage) {
long l = strlen(errormessage);
char *em = R_alloc(l, 1);
return strcpy(em, errormessage);
}
/** This function puts the mapping between the newtotal items after
* nl[offset] and the items in orig into the buffer perm, which has to
* be at least as long as newtotal. The function uses R indexing,
* that is to say, the first index is 1.
*/
int matchnames(const char **orig, int origtotal,
const NameList &nl, int offset, int newtotal,
int *perm) {
#ifdef DEBUG
printf("matching names (R indexing):\n");
#endif
for (int i=0; i < newtotal; i++) {
int j;
for (j=0; j < origtotal; j++)
if (strcmp(nl.getName(offset+i), *(orig+j))==0) {
*(perm+i) = j+1;
#ifdef DEBUG
printf("%d -> %d\n",i+1,j+1);
#endif
break;
}
if (j==origtotal)
return 1;
}
return 0;
}
/** dynareR is the interface function. The user provides:
* - a list of endogenous and exogenous variables, a list of
* parameters (and the length of each list)
* - the model equations (modeleq, pointer to a 0-terminated string)
* - the order of expansion (ord)
* - journal file name (jnlfile, can be "/dev/null" for no journal)
* - values for the parametes (parval)
* - variance-covariance matrix (vcov, stacked by columns, R does
* this)
* - initial values for finding the steady state (initval)
* - and the number of steps for the approximation algorithm
* (numsteps)
*
* If successful, the interface will write the results to these
* buffers:
* - tensorbuffer for the steady state and the flattened tensors
* - num_state for the number of endogenous variables that ended up in
* the state
* - mappings to variable names (ordering_state, ordering_endo,
* ordering_exo), indices start from 1
* - the deterministic steady state (newinitval)
*
* If dynare throws an exception, the interface tries to catch it and
* return an error code (error), and error message (errormessage), and
* if applicable, information on the stability of the model
* (kordcode). errormessage is allocated into R's memory, and will be
* collected after duplication.
*/
extern "C" {
void dynareR(const char** endo, const int* num_endo,
const char** exo, const int* num_exo,
const char** par, const int* num_par,
const char** equations, const int* ord, const char* jnlfile,
const double *parval, const double *vcov,
const double *initval,
const int *num_steps,
double* tensorbuffer,
int *num_state, int *ordering_state,
int *ordering_endo, int *ordering_exo,
double *newinitval,
int* error, char **errormessage, int *kordcode) {
// construct the model here
try {
#ifdef DEBUG // will print only first var names etc.
printf("eq: %s\nendo: %d %s\nexo: %d %s\npar: %d %s\nord: %d\n",
*equations,*num_endo,*endo,*num_exo,*exo,*num_par,*par,*ord);
#endif
// create journal
Journal journal(jnlfile);
// create Dynare object
Dynare dynare(endo, *num_endo, exo, *num_exo,
par, *num_par, *equations, strlen(*equations),
*ord, journal);
// set Vcov and parameter values
copy(parval,parval+(*num_par),dynare.getParams().base());
#ifdef DEBUG
printf("parameter values (%d):\n",dynare.getParams().length());
dynare.getParams().print();
#endif
copy(vcov,vcov+(*num_exo)*(*num_exo),dynare.getVcov().base());
#ifdef DEBUG
printf("vcov matrix:\n");
dynare.getVcov().print();
#endif
// set initial values
Vector iv(initval,*num_endo);
#ifdef DEBUG
printf("initial values:\n");
iv.print();
#endif
dynare.setInitOuter(iv);
// construct approximation
tls.init(dynare.order(),
dynare.nstat()+2*dynare.npred()+3*dynare.nboth()+
2*dynare.nforw()+dynare.nexog());
Approximation approximation(dynare,journal,*num_steps);
approximation.walkStochSteady();
// write the steady state into the buffer
int ny = dynare.ny();
const Vector ss(dynare.getSteady());
// ss = ConstVector(approximation.getSS(), 0); // FIXME allow
// // for nonzero
int s = dynare.getStateNames().getNum();
int sm = s;
tensorbuffer = copy(ss.base(),ss.base()+ny,tensorbuffer);
// write the tensors into buffer
const UnfoldDecisionRule& udr =
approximation.getUnfoldDecisionRule();
for (int i=1; i <= *ord; i++) {
const UFSTensor* t = udr.get(Symmetry(i));
#ifdef DEBUG
printf("tensor %d:\n", i);
t->print();
#endif
tensorbuffer = copy(t->base(), t->base()+ny*sm, tensorbuffer);
sm *= s;
}
// save number of endogenous states
*num_state = s-(*num_exo);
// ordering
#ifdef DEBUG
printf("all endo names:\n");
dynare.getAllEndoNames().print();
printf("all state names:\n");
dynare.getStateNames().print();
#endif
if (matchnames(endo, *num_endo, dynare.getAllEndoNames(),
0, *num_endo, ordering_endo) ||
matchnames(endo, *num_endo, dynare.getStateNames(),
0, *num_state, ordering_state) ||
matchnames(exo, *num_exo, dynare.getStateNames(),
*num_state, *num_exo, ordering_exo)) {
*error = DYNARER_NAMESMATCHINGERROR;
*errormessage = "There was a problem when matching names. This is weird and should not happen.";
return;
}
// return new init values (first column of SS matrix)
ConstVector newinit((const GeneralMatrix&) approximation.getSS(), 0);
#ifdef DEBUG
printf("new initial values:\n");
newinit.print();
#endif
copy(newinit.base(),newinit.base()+(*num_endo),newinitval);
} catch (const SylvException &e) {
*error = DYNARER_SYLVEXCEPTION;
char errorbuffer[501];
e.printMessage(errorbuffer, 500);
*errormessage = passmessage(errorbuffer);
#ifdef DEBUG
printf("Caught Sylv exception: ");
e.printMessage();
#endif
return;
} catch (const DynareException &e) {
*error = DYNARER_DYNAREEXCEPTION;
*errormessage = passmessage(e.message());
#ifdef DEBUG
printf("Caught Dynare exception: %s\n", e.message());
#endif
return;
} catch (const ogu::Exception &e) {
*error = DYNARER_OGUEXCEPTION;
*errormessage = passmessage(e.message());
#ifdef DEBUG
printf("Caught ogu::Exception: ");
e.print();
#endif
return;
} catch (const TLException &e) {
*error = DYNARER_TLEXCEPTION;
*errormessage = passmessage(e.getmessage());
#ifdef DEBUG
printf("Caugth TL exception: ");
e.print();
#endif
return;
} catch (const KordException &e) {
*error = DYNARER_KORDEXCEPTION;
*errormessage = passmessage(e.getmessage());
*kordcode = e.code(); // Kord error code
#ifdef DEBUG
printf("Caugth Kord exception: ");
e.print();
#endif
return;
}
*error = 0;
return;}
}
dynare++/extern/R/dynareR.r deleted 100644 → 0
View file @ 6ca649c4
## $Id: dynareR.r 862 2006-08-04 17:34:56Z tamas $
## Copyright 2006, Tamas K Papp
dyn.load("dynareR.so") # FIXME: make it platform-independent
## FIXME hide auxiliary functions in a namespace
dynareR.indextensor <- function(ord, nume, nums) {
nume*((nums^ord-1)/(nums-1))
}
dynareR.extracttensor <- function(tensor, ord, nume, nums) {
aperm(array(tensor[dynareR.indextensor(ord,nume,nums)+(1:(nume*nums^ord))],
c(nume,rep(nums,ord))),(ord+1):1)
}
dynareR.errormessages <- c("Sylvester exception",
"Dynare exception",
"OGU exception",
"Tensor library exception",
"K-order expansion library exception",
"Error matching names")
calldynare <- function(modeleq, endo, exo, parameters, expandorder,
parval, vcovmatrix, initval=rep(1,length(endo)),
numsteps=0, jnlfile="/dev/null") {
## check type of parameters
local({
is.charvector <- function(cv) { is.character(cv) && is.vector(cv) }
stopifnot(is.charvector(modeleq) && is.charvector(endo) &&
is.charvector(exo) && is.charvector(parameters) &&
is.charvector(jnlfile))
})
stopifnot(is.numeric(expandorder) && is.vector(expandorder) &&
(length(expandorder) == 1) && (expandorder >= 0))
stopifnot(length(jnlfile) == 1)
local({ # variable names
checkvarname <- function(v) {
stopifnot(length(grep("[^a-zA-Z].*",v)) == 0) # look for strange chars
}
checkvarname(endo)
checkvarname(exo)
checkvarname(parameters)
})
stopifnot(is.vector(parval) && is.numeric(parval))
stopifnot(is.vector(initval) && is.numeric(initval))
stopifnot(is.matrix(vcovmatrix) && is.numeric(vcovmatrix))
stopifnot(is.numeric(numsteps) && is.vector(numsteps) &&
(length(numsteps)==1))
## append semicolons to model equations if necessary
modeleq <- sapply(modeleq, function(s) {
if (length(grep("^.*; *$",s))==1)
s
else
sprintf("%s;",s)
})
## then concatenate into a single string
modeleq <- paste(modeleq, collapse=" ")
## call dynareR
nume <- length(endo)
maxs <- length(endo)+length(exo)
dr <- .C("dynareR",
endo,as.integer(nume),
exo,as.integer(length(exo)),
parameters,as.integer(length(parameters)),
modeleq,as.integer(expandorder),jnlfile,
as.double(parval),as.double(vcovmatrix),
as.double(initval),
as.integer(numsteps),
tensorbuffer=double(dynareR.indextensor(expandorder+1,nume,maxs)),
numstate=integer(1), orderstate=integer(maxs),
orderendo=integer(nume),
orderexo=integer(length(exo)),
newinitval=double(nume),
error=integer(1),
errormessage=character(1),
kordcode=integer(1))
## check for errors
kordcode <- 0
if (dr$error == 0) {
if (dr$error == 5) {
list(kordcode=dr$kordcode - 251) # magic dynare++ constant
} else {
## return result
with(dr, {
nums <- numstate+length(exo)
list(ss=dynareR.extracttensor(dr$tensorbuffer,0,nume,nums), # ss
rule=sapply(1:expandorder,function (o) { # decision rule
dynareR.extracttensor(dr$tensorbuffer,o,nume,nums)
}),
orderstate=orderstate[1:numstate], # state ordering
orderendo=orderendo, # endog. ordering
orderexo=orderexo, # exog. ordering
newinitval=newinitval, # new init values
kordcode=0)
})
}
} else {
stop(sprintf("%s (\"%s\")",dynareR.errormessages[dr$error],
dr$errormessage))
}
}
dynare++/extern/R/dynareR.tex deleted 100644 → 0
View file @ 6ca649c4
%% $Id: dynareR.tex 863 2006-08-04 17:35:21Z tamas $
%% Copyright Tamas K Papp, 2006
%% should compile with any reasonable TeX distribution, I am using tetex
\documentclass[12pt,a4paper]{article}
\usepackage{amsmath}
\usepackage{amsfonts}
%\usepackage[letterpaper,vmargin=1.7in]{geometry}
%\usepackage[letterpaper,left=2cm,right=8cm,bottom=3cm,top=3cm,marginparwidth=4cm]{geometry}
%\usepackage{natbib}
\usepackage{graphicx}
\usepackage{url}
\usepackage{natbib}
\usepackage{color}
\usepackage{paralist} % compactitem
\DeclareMathOperator{\Var}{Var}
\DeclareMathOperator{\Cov}{Cov}
\DeclareMathOperator{\argmin}{argmin}
\DeclareMathOperator{\argmax}{argmax}
\DeclareMathSymbol{\ueps}{\mathord}{letters}{"0F} % ugly epsilon
\renewcommand{\epsilon}{\varepsilon}
\newcommand{\aseq}{\overset{as}=} % almost surely equals
\usepackage{fancyhdr}
\pagestyle{fancy}
\lhead{Tam\'as K Papp} \chead{} \rhead{DynareR}
\cfoot{\thepage}
\renewcommand\floatpagefraction{.9}
\renewcommand\topfraction{.9}
\renewcommand\bottomfraction{.9}
\renewcommand\textfraction{.1}
\usepackage{listings}
\lstset{
language=R,
extendedchars=true,
basicstyle=\footnotesize,
stringstyle=\ttfamily,
commentstyle=\slshape,
% numbers=left,
% stepnumber=5,
% numbersep=6pt,
% numberstyle=\footnotesize,
breaklines=true,
frame=single,
columns=fullflexible,
}
\begin{document}
\title{DynareR}
\author{Tam\'as K Papp (\url{tpapp@princeton.edu})}
\date{\today}
\maketitle
DynareR is an R interface for Ondra Kamen\'ik's Dynare++ program. The
interface is still under development, and the functions might change.
However, I thought that some documentation would help to get users
started.
The purpose of DynareR is to return the transition rule (the
steady state and a list of tensors) for a given model. DynareR
does not simulate, and currently does no checking of the
approximation. Primarily, the interface is to be intended to be used
in Bayesian estimation of DSGE models (via MCMC).
Before you read on, make sure that
\begin{compactitem}
\item you understand what Dynare++ is and how it works,
\item you have compiled Dynare++ and DynareR (see \verb!README! in
\verb!extern/R!), and placed \verb!dynareR.so! and
\verb!dynareR.r! in your load path for R.
\end{compactitem}
The function that performs all the work is called
\lstinline{calldynare}. Its is defined like this:
\begin{lstlisting}
calldynare <- function(modeleq, endo, exo, parameters, expandorder,
parval, vcovmatrix, initval=rep(1,length(endo)),
numsteps=0, jnlfile="/dev/null") {
...
}
\end{lstlisting}
\lstinline{modeleq} is a character vector for the model equations, and
it may have a length longer than one. First, \lstinline{calldynare}
checks if each string in the vector has a terminating semicolon (may
be followed by whitespace), if it doesn't, then it appends one. Then
it concatenates all equations into a single string. Thus, the
following versions of \lstinline{modeleq} give equivalent results:
\begin{lstlisting}
modeleq1 <- c("(c/c(1))^gamma*beta*(alpha*exp(a(1))*k^(alpha-1)+1-delta)=1",
"a=rho*a(-1)+eps",
"k+c=exp(a)*k(-1)^alpha+(1-delta)*k(-1)")
modeleq2 <- c("(c/c(1))^gamma*beta*(alpha*exp(a(1))*k^(alpha-1)+1-delta)=1;",
"a=rho*a(-1)+eps ; ",
"k+c=exp(a)*k(-1)^alpha+(1-delta)*k(-1) \t;\t ")
modeleq3 <- paste(modeleq1, collapse=" ")
\end{lstlisting}
The next three arguments name the endo- and exogenous variables and
the parameters. The names should be character vectors, for example,
\begin{lstlisting}
parameters <- c("beta","gamma","rho","alpha","delta")
varendo <- c("k","c","a")
varexo <- "eps"
\end{lstlisting}
\lstinline{calldynare} also needs the order of the approximation
\lstinline{expandorder} (a nonnegative integer), the parameter values
\lstinline{parval} (should be the same length as
\lstinline{parameters}), a variance-covariance matrix \lstinline{vcov}
(dimensions should match the length of \lstinline{exo}) and initial
values for finding the deterministic steady state
(\lstinline{initval}). If you don't provide initial values,
\lstinline{calldynare} will use a sequence of $1$s, on the assumption
that most variables in economics are positive --- you should always
try to provide a reasonable initial guess for the nonlinear solver if
possible (if you are doing MCMC, chances are that you only have to do
it once, see \lstinline{newinitval} below).
You can also provide the number of steps for calculating the
stochastic steady state (\lstinline{numsteps}, the default is zero,
see the dynare++ tutorial for more information) and the name of the
journal file \lstinline{jnlfile}. If you don't provide a journal
file, the default is \verb!/dev/null!.
Below, you see an example of using dynareR.
\lstinputlisting{test.r}
\lstinline{calldynare} returns the results in a list, variables below
refer to elements of this list. First, you should always check
\lstinline{kordcode}, which tells whether dynare++ could calculate an
approximation. It can have the following values:
\begin{description}
\item[0] the calculation was successful
\item[1] the system is not stable (Blanchard-Kahn)
\item[2] failed to calculate fixed point (infinite values)
\item[3] failed to calculate fixed point (NaN values)
\end{description}
If \lstinline{kordcode} is nonzero, then the list has only this
element.
If \lstinline{kordcode} equals zero, then the list has the following
elements:
\begin{description}
\item[ss] the steady state (ordered by \lstinline{orderendo}), which
is a vector
\item[rule] the transition rule (ordered by \lstinline{orderendo},
\lstinline{orderstate} and \lstinline{orderexo}), a list of arrays
\item[newinitval] the deterministic steady state, you can use this to
initialize the nonlinear solver for a nearby point in the parameter
space (ordered by \lstinline{orderendo})
\item[orderstate] the index of endogenous variables that ended up in
the state
\item[orderendo] the ordering of endogenous variables
\item[orderexo] the ordering of exogenous variables
\item[kordcode] discussed above
\end{description}
An example will illustrate the ordering. To continue the example above,
\begin{lstlisting}
> dd$orderstate
[1] 1 3
> dd$orderendo
[1] 1 3 2
> dd$orderexo
[1] 1
> dd$rule[[1]]
[,1] [,2] [,3]
[1,] 0.9669374 0.0 0.02071077
[2,] 2.4230073 0.9 0.45309125
[3,] 2.6922303 1.0 0.50343473
\end{lstlisting}
Recall that the original ordering of endogenous variables was
\lstinline{k, c, a}. The vectors and matrices of the result are
ordered as \lstinline{varendo[dd$orderendo]}, that is, as
\lstinline{k, a, c}. This is the ordering for the steady state and
the first dimension of the tensors in \lstinline{rule}. The other
dimensions are ordered as
\lstinline{c(varendo[dd$orderstate],varexo[dd$orderexo])}, that is to
say, as \lstinline{k, a, eps}. Use these orderings when calculating
with the tensors and the steady state. Also, remember that the $i$th
tensor is already divided by $i!$.
\lstinline{calldynare} also handles exceptions from dynare. All
exceptions (except KordException, which sets \lstinline{kordcode})
generate an error in the R interface. Normally, when solving a
well-formed model (no typos in the equations, etc), users should not
encounter these exceptions. Having a journal file is useful for
debugging. If you are making long calculations, it is reasonable to
catch errors with \lstinline{try} so that they won't abort the
calculation.
% \bibliographystyle{apalike}
% \bibliography{/home/tpapp/doc/general.bib}
\end{document}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
dynare++/extern/R/test.cpp deleted 100644 → 0
View file @ 6ca649c4
#include "dynareR.cpp"
int main(void) {
const char *parameters[] = {"beta","gamma","rho","alpha","delta"};
const char *varendo[] = {"k","c","a"};
const char *varexo[] = {"eps"};
const int numpar = 5;
const int numendo = 3;
const int numexo = 1;
const int ord = 2;
const int numsteps = 0;
const double parval[] = {.99,2,.9,.3,.025};
const double vcov[] = {0.001};
const double initval[] = {0.066, 0.43, 0.01};
int e;
double tensorbuffer[100];
int num_state;
int ordering_state[] = {0,0,0};
int ordering_endo[] = {0,0,0};
int ordering_exo[] = {0};
double newinitval[] = {0,0,0};
const char *modeleq[] = {"(c/c(1))^gamma*beta*(alpha*exp(a(1))*k^(alpha-1)+1-delta)=1; a=rho*a(-1)+eps; k+c=exp(a)*k(-1)^alpha+(1-delta)*k(-1);"};
dynareR(varendo, &numendo, varexo, &numexo, parameters, &numpar, modeleq,
&ord, "journal", parval, vcov, initval,
&numsteps, tensorbuffer,
&num_state, ordering_state, ordering_endo, ordering_exo,
newinitval,&e);
printf("error code: %d\n", e);
}
dynare++/extern/R/test.r deleted 100644 → 0
View file @ 6ca649c4
source("dynareR.r")
parameters <- c("beta","gamma","rho","alpha","delta")
varendo <- c("k","c","a")
varexo <- "eps"
parval <- c(.99,2,.9,.3,.025)
vcovmatrix <- matrix(1,1,1)
initval <- c(0.066, 0.43, 0.01)
modeleq <- c("(c/c(1))^gamma*beta*(alpha*exp(a(1))*k^(alpha-1)+1-delta)=1",
"a=rho*a(-1)+eps",
"k+c=exp(a)*k(-1)^alpha+(1-delta)*k(-1)")
dd <- calldynare(modeleq,varendo,varexo,parameters,2,parval,vcovmatrix,initval)
dynare++/extern/matlab/dynare_simul.cpp deleted 100644 → 0
View file @ 6ca649c4
// Copyright (C) 2005-2011, Ondra Kamenik
// This is the mexFunction providing interface to
// DecisionRule<>::simulate(). It takes the following input
// parameters:
// order the order of approximation, needs order+1 derivatives
// nstat
// npred
// nboth
// nforw
// nexog
// ystart starting value (full vector of endogenous)
// shocks matrix of shocks (nexog x number of period)
// vcov covariance matrix of shocks (nexog x nexog)
// seed integer seed
// ysteady full vector of decision rule's steady
// ... order+1 matrices of derivatives
// output:
// res simulated results
#include "dynmex.h"
#include "mex.h"
#include "decision_rule.h"
#include "fs_tensor.h"
#include "SylvException.h"
extern "C" {
void mexFunction(int nlhs, mxArray* plhs[],
int nhrs, const mxArray* prhs[])
{
if (nhrs < 12 || nlhs != 2)
DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have at least 12 input parameters and exactly 2 output arguments.\n");
int order = (int)mxGetScalar(prhs[0]);
if (nhrs != 12 + order)
DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have exactly 11+order input parameters.\n");
int nstat = (int)mxGetScalar(prhs[1]);
int npred = (int)mxGetScalar(prhs[2]);
int nboth = (int)mxGetScalar(prhs[3]);
int nforw = (int)mxGetScalar(prhs[4]);
int nexog = (int)mxGetScalar(prhs[5]);
const mxArray* const ystart = prhs[6];
const mxArray* const shocks = prhs[7];
const mxArray* const vcov = prhs[8];
int seed = (int)mxGetScalar(prhs[9]);
const mxArray* const ysteady = prhs[10];
const mwSize* const ystart_dim = mxGetDimensions(ystart);
const mwSize* const shocks_dim = mxGetDimensions(shocks);
const mwSize* const vcov_dim = mxGetDimensions(vcov);
const mwSize* const ysteady_dim = mxGetDimensions(ysteady);
int ny = nstat + npred + nboth + nforw;
if (ny != (int) ystart_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of rows.\n");
if (1 != ystart_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of cols.\n");
int nper = shocks_dim[1];
if (nexog != (int) shocks_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("shocks has a wrong number of rows.\n");
if (nexog != (int) vcov_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of rows.\n");
if (nexog != (int) vcov_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of cols.\n");
if (ny != (int) ysteady_dim[0])
DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of rows.\n");
if (1 != ysteady_dim[1])
DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of cols.\n");
mxArray* res = mxCreateDoubleMatrix(ny, nper, mxREAL);
try {
// initialize tensor library
tls.init(order, npred+nboth+nexog);
// form the polynomial
UTensorPolynomial pol(ny, npred+nboth+nexog);
for (int dim = 0; dim <= order; dim++) {
const mxArray* gk = prhs[11+dim];
const mwSize* const gk_dim = mxGetDimensions(gk);
FFSTensor ft(ny, npred+nboth+nexog, dim);
if (ft.ncols() != (int) gk_dim[1]) {
char buf[1000];
sprintf(buf, "Wrong number of columns for folded tensor: got %d but I want %d\n",
(int) gk_dim[1], ft.ncols());
DYN_MEX_FUNC_ERR_MSG_TXT(buf);
}
if (ft.nrows() != (int) gk_dim[0]) {
char buf[1000];
sprintf(buf, "Wrong number of rows for folded tensor: got %d but I want %d\n",
(int) gk_dim[0], ft.nrows());
DYN_MEX_FUNC_ERR_MSG_TXT(buf);
}
ft.zeros();
ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), mxGetPr(gk));
ft.add(1.0, gk_mat);
UFSTensor* ut = new UFSTensor(ft);
pol.insert(ut);
}
// form the decision rule
UnfoldDecisionRule
dr(pol, PartitionY(nstat, npred, nboth, nforw),
nexog, ConstVector(mxGetPr(ysteady), ny));
// form the shock realization
TwoDMatrix shocks_mat(nexog, nper, (const double*)mxGetPr(shocks));
TwoDMatrix vcov_mat(nexog, nexog, (const double*)mxGetPr(vcov));
GenShockRealization sr(vcov_mat, shocks_mat, seed);
// simulate and copy the results
Vector ystart_vec((const double*)mxGetPr(ystart), ny);
TwoDMatrix* res_mat =
dr.simulate(DecisionRule::horner, nper,
ystart_vec, sr);
TwoDMatrix res_tmp_mat(ny, nper, mxGetPr(res));
res_tmp_mat = (const TwoDMatrix&)(*res_mat);
delete res_mat;
plhs[1] = res;
} catch (const KordException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caugth Kord exception.");
} catch (const TLException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caugth TL exception.");
} catch (SylvException& e) {
DYN_MEX_FUNC_ERR_MSG_TXT("Caught Sylv exception.");
}
plhs[0] = mxCreateDoubleScalar(0);
}
};
dynare++/extern/matlab/dynare_simul.m deleted 100644 → 0
View file @ 6ca649c4
%
% SYNOPSIS
%
% r = dynare_simul(name, shocks)
% r = dynare_simul(name, prefix, shocks)
% r = dynare_simul(name, shocks, start)
% r = dynare_simul(name, prefix, shocks, start)
%
% name name of MAT-file produced by dynare++
% prefix prefix of variables in the MAT-file
% shocks matrix of shocks
% start zero period value
%
% Note that this file requires the dynare_simul_ DLL to be in the path.
% This DLL is distributed with Dynare, under the mex/matlab or mex/octave
% subdirectory.
%
% SEMANTICS
%
% The command reads a decision rule from the MAT-file having the given
% prefix. Then it starts simulating the decision rule with zero time value
% equal to the given start. It uses the given shocks for the simulation. If
% the start is not given, the state about which the decision rule is
% centralized is taken (called fix point, or stochastic steady state, take
% your pick).
%
% prefix Use the prefix with which you called dynare++, the default
% prefix in dynare++ is 'dyn'.
% shocks Number of rows must be a number of exogenous shocks,
% number of columns gives the number of simulated
% periods. NaNs and Infs in the matrix are substitued by
% draws from the normal distribution using the covariance
% matrix given in the model file.
% start Vector of endogenous variables in the ordering given by
% <prefix>_vars.
%
% Seed for random generator is derived from calling rand(1,1). Therefore,
% seeding can be controlled with rand('state') and rand('state',some_seed).
%
% EXAMPLES
%
% All examples suppose that the prefix is 'dyn' and that your_model.mat
% has been loaded into Matlab.
%
% 1. response to permanent negative shock to the third exo var EPS3 for
% 100 periods
%
% shocks = zeros(4,100); % 4 exogenous variables in the model
% shocks(dyn_i_EPS3,:) = -0.1; % the permanent shock to EPS3
% r = dynare_simul('your_model.mat',shocks);
%
% 2. one stochastic simulation for 100 periods
%
% shocks = zeros(4,100)./0; % put NaNs everywhere
% r = dynare_simul('your_model.mat',shocks);
%
% 3. one stochastic simulation starting at 75% undercapitalized economy
%
% shocks = zeros(4,100)./0; % put NaNs everywhere
% ystart = dyn_ss; % get copy of DR fix point
% ystart(dyn_i_K) = 0.75*dyn_ss(dyn_i_K); % scale down the capital
% r = dynare_simul('your_model.mat',shocks,ystart);
%
%
% SEE ALSO
%
% "DSGE Models with Dynare++. A Tutorial.", Ondra Kamenik, 2005
% Copyright (C) 2005-2011, Ondra Kamenik
function r = dynare_simul(varargin)
if exist('dynare_simul_') ~= 3
error('Can''t find dynare_simul_ DLL in the path. The simplest way to add it is to run Dynare once in this session.')
end
% get the file name and load data
fname = varargin{1};
eval(['load ' fname]);
% set prefix, shocks, ystart
if ischar(varargin{2})
prefix = varargin{2};
if length(varargin) == 3
shocks = varargin{3};
ystart = NaN;
elseif length(varargin) == 4
shocks = varargin{3};
ystart = varargin{4};
else
error('Wrong number of parameters.');
end
else
prefix = 'dyn';
if length(varargin) == 2
shocks = varargin{2};
ystart = NaN;
elseif length(varargin) == 3
shocks = varargin{2};
ystart = varargin{3};
else
error('Wrong number of parameters.');
end
end
% load all needed variables but prefix_g_*
if (exist([prefix '_nstat']))
nstat = eval([prefix '_nstat']);
else
error(['Could not find variable ' prefix '_nstat in workspace']);
end
if (exist([prefix '_npred']))
npred = eval([prefix '_npred']);
else
error(['Could not find variable ' prefix '_npred in workspace']);
end
if (exist([prefix '_nboth']))
nboth = eval([prefix '_nboth']);
else
error(['Could not find variable ' prefix '_nboth in workspace']);
end
if (exist([prefix '_nforw']))
nforw = eval([prefix '_nforw']);
else
error(['Could not find variable ' prefix '_nforw in workspace']);
end
if (exist([prefix '_ss']))
ss = eval([prefix '_ss']);
else
error(['Could not find variable ' prefix '_ss in workspace']);
end
if (exist([prefix '_vcov_exo']))
vcov_exo = eval([prefix '_vcov_exo']);
else
error(['Could not find variable ' prefix '_vcov_exo in workspace']);
end
nexog = size(vcov_exo,1);
if isnan(ystart)
ystart = ss;
end
% newer version of dynare++ doesn't return prefix_g_0, we make it here if
% it does not exist in workspace
g_zero = [prefix '_g_0'];
if (~ exist(g_zero))
eval([ g_zero '= zeros(nstat+npred+nboth+nforw,1);']);
end
% make derstr a string of comma seperated existing prefix_g_*
derstr = [',' g_zero];
order = 1;
cont = 1;
while cont == 1
g_ord = [prefix '_g_' num2str(order)];
if (exist(g_ord))
derstr = [derstr ',' g_ord];
order = order + 1;
else
cont = 0;
end
end
% set seed
seed = ceil(10000*rand(1,1));
% call dynare_simul_
command = ['[err,r]=dynare_simul_(' num2str(order-1) ',nstat,npred,nboth,nforw,' ...
'nexog,ystart,shocks,vcov_exo,seed,ss' derstr ');'];
eval(command);
if err
error('Simulation failed')
end
dynare++/integ/Makefile.am deleted 100644 → 0
View file @ 6ca649c4
SUBDIRS = cc src testing
dynare++/integ/cc/Makefile.am deleted 100644 → 0
View file @ 6ca649c4
CWEBSRC = \
product.cweb \
quadrature.cweb \
quasi_mcarlo.cweb \
smolyak.cweb \
vector_function.cweb \
product.hweb \
quadrature.hweb \
quasi_mcarlo.hweb \
smolyak.hweb \
vector_function.hweb
GENERATED_FILES = \
product.cpp \
quadrature.cpp \
quasi_mcarlo.cpp \
smolyak.cpp \
vector_function.cpp \
product.h \
quadrature.h \
quasi_mcarlo.h \
smolyak.h \
vector_function.h
noinst_LIBRARIES = libinteg.a
libinteg_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES) precalc_quadrature.dat
libinteg_a_CPPFLAGS = -I../../sylv/cc -I../../tl/cc -I$(top_srcdir)/mex/sources
libinteg_a_CXXFLAGS = $(PTHREAD_CFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
EXTRA_DIST = main.web dummy.ch
%.cpp: %.cweb dummy.ch
$(CTANGLE) -bhp $< dummy.ch $@
%.h: %.hweb dummy.ch
$(CTANGLE) -bhp $< dummy.ch $@
if HAVE_CWEAVE
if HAVE_PDFTEX
if HAVE_EPLAIN
pdf-local: integ.pdf
integ.pdf: main.web $(CWEBSRC)
$(CWEAVE) -bhp main.web
$(PDFTEX) main
mv main.pdf integ.pdf
endif
endif
endif
CLEANFILES = integ.pdf main.idx main.log main.scn main.tex main.toc
dynare++/integ/cc/dummy.ch deleted 100644 → 0
View file @ 6ca649c4
dynare++/integ/cc/main.web deleted 100644 → 0
View file @ 6ca649c4
@q $Id: main.web 2333 2009-01-14 10:32:55Z kamenik $ @>
@q Copyright 2005, 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 Numerical Integration Module}}
@i ../../c++lib.w
@s Vector int
@s ConstVector int
@s IntSequence int
@s GeneralMatrix int
@s THREAD int
@s THREAD_GROUP int
@s SYNCHRO int
\titletrue
\null\vfill
\centerline{\titlefont Numerical Integration Module}
\vfill\vfill
Copyright \copyright\ 2005 by Ondra Kamenik
\penalty-10000
@i vector_function.hweb
@i vector_function.cweb
@i quadrature.hweb
@i quadrature.cweb
@i product.hweb
@i product.cweb
@i smolyak.hweb
@i smolyak.cweb
@i quasi_mcarlo.hweb
@i quasi_mcarlo.cweb
dynare++/integ/cc/precalc_quadrature.dat deleted 100644 → 0
View file @ 6ca649c4
// $Id: precalc_quadrature.dat 431 2005-08-16 15:41:01Z kamenik $
// Copyright 2005, Ondra Kamenik
// The file contains one dimensional quadrature points and weights for
// a few quadratures. The format of data is clear. There is a class
// OneDPrecalcQuadrature which implements an interface OneDQuadrature
// using the data of this format.
// Gauss-Hermite quadrature; prefix gh
// number of levels
static const int gh_num_levels = 26;
// number of points in each level
static const int gh_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
30, 32, 40, 50, 60, 64
};
// weights, starting with the first level
static const double gh_weights[] = {
// weights 1 = sqrt(pi)
1.77245385090551588191942755656782537698745727539062,
// weights 2
0.886226925452758013649083741671e+00,
0.886226925452758013649083741671e+00,
// weights 3
0.295408975150919337883027913890e+00,
0.118163590060367735153211165556e+01,
0.295408975150919337883027913890e+00,
// weights 4
0.813128354472451771430345571899e-01,
0.804914090005512836506049184481e+00,
0.804914090005512836506049184481e+00,
0.813128354472451771430345571899e-01,
// weights 5
0.199532420590459132077434585942e-01,
0.393619323152241159828495620852e+00,
0.945308720482941881225689324449e+00,
0.393619323152241159828495620852e+00,
0.199532420590459132077434585942e-01,
// weights 6
0.453000990550884564085747256463e-02,
0.157067320322856643916311563508e+00,
0.724629595224392524091914705598e+00,
0.724629595224392524091914705598e+00,
0.157067320322856643916311563508e+00,
0.453000990550884564085747256463e-02,
// weights 7
0.971781245099519154149424255939e-03,
0.545155828191270305921785688417e-01,
0.425607252610127800520317466666e+00,
0.810264617556807326764876563813e+00,
0.425607252610127800520317466666e+00,
0.545155828191270305921785688417e-01,
0.971781245099519154149424255939e-03,
// weights 8
0.199604072211367619206090452544e-03,
0.170779830074134754562030564364e-01,
0.207802325814891879543258620286e+00,
0.661147012558241291030415974496e+00,
0.661147012558241291030415974496e+00,
0.207802325814891879543258620286e+00,
0.170779830074134754562030564364e-01,
0.199604072211367619206090452544e-03,
// weights 9
0.396069772632643819045862946425e-04,
0.494362427553694721722456597763e-02,
0.884745273943765732879751147476e-01,
0.432651559002555750199812112956e+00,
0.720235215606050957124334723389e+00,
0.432651559002555750199812112956e+00,
0.884745273943765732879751147476e-01,
0.494362427553694721722456597763e-02,
0.396069772632643819045862946425e-04,
// weights 10
0.764043285523262062915936785960e-05,
0.134364574678123269220156558585e-02,
0.338743944554810631361647312776e-01,
0.240138611082314686416523295006e+00,
0.610862633735325798783564990433e+00,
0.610862633735325798783564990433e+00,
0.240138611082314686416523295006e+00,
0.338743944554810631361647312776e-01,
0.134364574678123269220156558585e-02,
0.764043285523262062915936785960e-05,
// weights 11
0.143956039371425822033088366032e-05,
0.346819466323345510643413772940e-03,
0.119113954449115324503874202916e-01,
0.117227875167708503381788649308e+00,
0.429359752356125028446073598601e+00,
0.654759286914591779203940657627e+00,
0.429359752356125028446073598601e+00,
0.117227875167708503381788649308e+00,
0.119113954449115324503874202916e-01,
0.346819466323345510643413772940e-03,
0.143956039371425822033088366032e-05,
// weights 12
0.265855168435630160602311400877e-06,
0.857368704358785865456906323153e-04,
0.390539058462906185999438432620e-02,
0.516079856158839299918734423606e-01,
0.260492310264161129233396139765e+00,
0.570135236262479578347113482275e+00,
0.570135236262479578347113482275e+00,
0.260492310264161129233396139765e+00,
0.516079856158839299918734423606e-01,
0.390539058462906185999438432620e-02,
0.857368704358785865456906323153e-04,
0.265855168435630160602311400877e-06,
// weights 13
0.482573185007313108834997332342e-07,
0.204303604027070731248669432937e-04,
0.120745999271938594730924899224e-02,
0.208627752961699392166033805050e-01,
0.140323320687023437762792268873e+00,
0.421616296898543221746893558568e+00,
0.604393187921161642342099068579e+00,
0.421616296898543221746893558568e+00,
0.140323320687023437762792268873e+00,
0.208627752961699392166033805050e-01,
0.120745999271938594730924899224e-02,
0.204303604027070731248669432937e-04,
0.482573185007313108834997332342e-07,
// weights 14
0.862859116812515794532041783429e-08,
0.471648435501891674887688950105e-05,
0.355092613551923610483661076691e-03,
0.785005472645794431048644334608e-02,
0.685055342234652055387163312367e-01,
0.273105609064246603352569187026e+00,
0.536405909712090149794921296776e+00,
0.536405909712090149794921296776e+00,
0.273105609064246603352569187026e+00,
0.685055342234652055387163312367e-01,
0.785005472645794431048644334608e-02,
0.355092613551923610483661076691e-03,
0.471648435501891674887688950105e-05,
0.862859116812515794532041783429e-08,
// weights 15
0.152247580425351702016062666965e-08,
0.105911554771106663577520791055e-05,
0.100004441232499868127296736177e-03,
0.277806884291277589607887049229e-02,
0.307800338725460822286814158758e-01,
0.158488915795935746883839384960e+00,
0.412028687498898627025891079568e+00,
0.564100308726417532852625797340e+00,
0.412028687498898627025891079568e+00,
0.158488915795935746883839384960e+00,
0.307800338725460822286814158758e-01,
0.277806884291277589607887049229e-02,
0.100004441232499868127296736177e-03,
0.105911554771106663577520791055e-05,
0.152247580425351702016062666965e-08,
// weights 16
0.265480747401118224470926366050e-09,
0.232098084486521065338749423185e-06,
0.271186009253788151201891432244e-04,
0.932284008624180529914277305537e-03,
0.128803115355099736834642999312e-01,
0.838100413989858294154207349001e-01,
0.280647458528533675369463335380e+00,
0.507929479016613741913517341791e+00,
0.507929479016613741913517341791e+00,
0.280647458528533675369463335380e+00,
0.838100413989858294154207349001e-01,
0.128803115355099736834642999312e-01,
0.932284008624180529914277305537e-03,
0.271186009253788151201891432244e-04,
0.232098084486521065338749423185e-06,
0.265480747401118224470926366050e-09,
// weights 17
0.458057893079863330580889281222e-10,
0.497707898163079405227863353715e-07,
0.711228914002130958353327376218e-05,
0.298643286697753041151336643059e-03,
0.506734995762753791170069495879e-02,
0.409200341495762798094994877854e-01,
0.172648297670097079217645196219e+00,
0.401826469470411956577635085257e+00,
0.530917937624863560331883103379e+00,
0.401826469470411956577635085257e+00,
0.172648297670097079217645196219e+00,
0.409200341495762798094994877854e-01,
0.506734995762753791170069495879e-02,
0.298643286697753041151336643059e-03,
0.711228914002130958353327376218e-05,
0.497707898163079405227863353715e-07,
0.458057893079863330580889281222e-10,
// weights 18
0.782819977211589102925147471012e-11,
0.104672057957920824443559608435e-07,
0.181065448109343040959702385911e-05,
0.918112686792940352914675407371e-04,
0.188852263026841789438175325426e-02,
0.186400423875446519219315221973e-01,
0.973017476413154293308537234155e-01,
0.284807285669979578595606820713e+00,
0.483495694725455552876410522141e+00,
0.483495694725455552876410522141e+00,
0.284807285669979578595606820713e+00,
0.973017476413154293308537234155e-01,
0.186400423875446519219315221973e-01,
0.188852263026841789438175325426e-02,
0.918112686792940352914675407371e-04,
0.181065448109343040959702385911e-05,
0.104672057957920824443559608435e-07,
0.782819977211589102925147471012e-11,
// weights 19
0.132629709449851575185289154385e-11,
0.216305100986355475019693077221e-08,
0.448824314722312295179447915594e-06,
0.272091977631616257711941025214e-04,
0.670877521407181106194696282100e-03,
0.798886677772299020922211491861e-02,
0.508103869090520673569908110358e-01,
0.183632701306997074156148485766e+00,
0.391608988613030244504042313621e+00,
0.502974888276186530840731361096e+00,
0.391608988613030244504042313621e+00,
0.183632701306997074156148485766e+00,
0.508103869090520673569908110358e-01,
0.798886677772299020922211491861e-02,
0.670877521407181106194696282100e-03,
0.272091977631616257711941025214e-04,
0.448824314722312295179447915594e-06,
0.216305100986355475019693077221e-08,
0.132629709449851575185289154385e-11,
// weights 20
0.222939364553415129252250061603e-12,
0.439934099227318055362885145547e-09,
0.108606937076928169399952456345e-06,
0.780255647853206369414599199965e-05,
0.228338636016353967257145917963e-03,
0.324377334223786183218324713235e-02,
0.248105208874636108821649525589e-01,
0.109017206020023320013755033535e+00,
0.286675505362834129719659706228e+00,
0.462243669600610089650328639861e+00,
0.462243669600610089650328639861e+00,
0.286675505362834129719659706228e+00,
0.109017206020023320013755033535e+00,
0.248105208874636108821649525589e-01,
0.324377334223786183218324713235e-02,
0.228338636016353967257145917963e-03,
0.780255647853206369414599199965e-05,
0.108606937076928169399952456345e-06,
0.439934099227318055362885145547e-09,
0.222939364553415129252250061603e-12,
// weights 30
0.290825470013122622941102747365e-20,
0.281033360275090370876277491534e-16,
0.287860708054870606219239791142e-13,
0.810618629746304420399344796173e-11,
0.917858042437852820850075742492e-09,
0.510852245077594627738963204403e-07,
0.157909488732471028834638794022e-05,
0.293872522892298764150118423412e-04,
0.348310124318685523420995323183e-03,
0.273792247306765846298942568953e-02,
0.147038297048266835152773557787e-01,
0.551441768702342511680754948183e-01,
0.146735847540890099751693643152e+00,
0.280130930839212667413493211293e+00,
0.386394889541813862555601849165e+00,
0.386394889541813862555601849165e+00,
0.280130930839212667413493211293e+00,
0.146735847540890099751693643152e+00,
0.551441768702342511680754948183e-01,
0.147038297048266835152773557787e-01,
0.273792247306765846298942568953e-02,
0.348310124318685523420995323183e-03,
0.293872522892298764150118423412e-04,
0.157909488732471028834638794022e-05,
0.510852245077594627738963204403e-07,
0.917858042437852820850075742492e-09,
0.810618629746304420399344796173e-11,
0.287860708054870606219239791142e-13,
0.281033360275090370876277491534e-16,
0.290825470013122622941102747365e-20,
// weights 32
0.731067642736e-22,
0.923173653649e-18,
0.119734401709e-14,
0.421501021125e-12,
0.593329146300e-10,
0.409883216476e-08,
0.157416779254e-06,
0.365058512955e-05,
0.541658406172e-04,
0.536268365526e-03,
0.365489032664e-02,
0.175534288315e-01,
0.604581309557e-01,
0.151269734076e+00,
0.277458142302e+00,
0.375238352592e+00,
0.375238352592e+00,
0.277458142302e+00,
0.151269734076e+00,
0.604581309557e-01,
0.175534288315e-01,
0.365489032664e-02,
0.536268365526e-03,
0.541658406172e-04,
0.365058512955e-05,
0.157416779254e-06,
0.409883216476e-08,
0.593329146300e-10,
0.421501021125e-12,
0.119734401709e-14,
0.923173653649e-18,
0.731067642736e-22,
// weights 40
0.259104371384e-28,
0.854405696375e-24,
0.256759336540e-20,
0.198918101211e-17,
0.600835878947e-15,
0.880570764518e-13,
0.715652805267e-11,
0.352562079135e-09,
0.112123608322e-07,
0.241114416359e-06,
0.363157615067e-05,
0.393693398108e-04,
0.313853594540e-03,
0.187149682959e-02,
0.846088800823e-02,
0.293125655361e-01,
0.784746058652e-01,
0.163378732713e+00,
0.265728251876e+00,
0.338643277425e+00,
0.338643277425e+00,
0.265728251876e+00,
0.163378732713e+00,
0.784746058652e-01,
0.293125655361e-01,
0.846088800823e-02,
0.187149682959e-02,
0.313853594540e-03,
0.393693398108e-04,
0.363157615067e-05,
0.241114416359e-06,
0.112123608322e-07,
0.352562079135e-09,
0.715652805267e-11,
0.880570764518e-13,
0.600835878947e-15,
0.198918101211e-17,
0.256759336540e-20,
0.854405696375e-24,
0.259104371384e-28,
// weights 50
0.183379404857e-36,
0.167380166790e-31,
0.121524412340e-27,
0.213765830835e-24,
0.141709359957e-21,
0.447098436530e-19,
0.774238295702e-17,
0.809426189344e-15,
0.546594403180e-13,
0.250665552389e-11,
0.811187736448e-10,
0.190904054379e-08,
0.334679340401e-07,
0.445702996680e-06,
0.458168270794e-05,
0.368401905377e-04,
0.234269892109e-03,
0.118901178175e-02,
0.485326382616e-02,
0.160319410684e-01,
0.430791591566e-01,
0.945489354768e-01,
0.170032455676e+00,
0.251130856331e+00,
0.305085129203e+00,
0.305085129203e+00,
0.251130856331e+00,
0.170032455676e+00,
0.945489354768e-01,
0.430791591566e-01,
0.160319410684e-01,
0.485326382616e-02,
0.118901178175e-02,
0.234269892109e-03,
0.368401905377e-04,
0.458168270794e-05,
0.445702996680e-06,
0.334679340401e-07,
0.190904054379e-08,
0.811187736448e-10,
0.250665552389e-11,
0.546594403180e-13,
0.809426189344e-15,
0.774238295702e-17,
0.447098436530e-19,
0.141709359957e-21,
0.213765830835e-24,
0.121524412340e-27,
0.167380166790e-31,
0.183379404857e-36,
// weights 60
0.110958724796e-44,
0.243974758810e-39,
0.377162672698e-35,
0.133255961176e-31,
0.171557314767e-28,
0.102940599693e-25,
0.334575695574e-23,
0.651256725748e-21,
0.815364047300e-19,
0.692324790956e-17,
0.415244410968e-15,
0.181662457614e-13,
0.594843051597e-12,
0.148895734905e-10,
0.289935901280e-09,
0.445682277521e-08,
0.547555461926e-07,
0.543351613419e-06,
0.439428693625e-05,
0.291874190415e-04,
0.160277334681e-03,
0.731773556963e-03,
0.279132482894e-02,
0.893217836028e-02,
0.240612727660e-01,
0.547189709320e-01,
0.105298763697e+00,
0.171776156918e+00,
0.237868904958e+00,
0.279853117522e+00,
0.279853117522e+00,
0.237868904958e+00,
0.171776156918e+00,
0.105298763697e+00,
0.547189709320e-01,
0.240612727660e-01,
0.893217836028e-02,
0.279132482894e-02,
0.731773556963e-03,
0.160277334681e-03,
0.291874190415e-04,
0.439428693625e-05,
0.543351613419e-06,
0.547555461926e-07,
0.445682277521e-08,
0.289935901280e-09,
0.148895734905e-10,
0.594843051597e-12,
0.181662457614e-13,
0.415244410968e-15,
0.692324790956e-17,
0.815364047300e-19,
0.651256725748e-21,
0.334575695574e-23,
0.102940599693e-25,
0.171557314767e-28,
0.133255961176e-31,
0.377162672698e-35,
0.243974758810e-39,
0.110958724796e-44,
// weights 64
0.553570653584e-48,
0.167974799010e-42,
0.342113801099e-38,
0.155739062462e-34,
0.254966089910e-31,
0.192910359546e-28,
0.786179778889e-26,
0.191170688329e-23,
0.298286278427e-21,
0.315225456649e-19,
0.235188471067e-17,
0.128009339117e-15,
0.521862372645e-14,
0.162834073070e-12,
0.395917776693e-11,
0.761521725012e-10,
0.117361674232e-08,
0.146512531647e-07,
0.149553293672e-06,
0.125834025103e-05,
0.878849923082e-05,
0.512592913577e-04,
0.250983698512e-03,
0.103632909950e-02,
0.362258697852e-02,
0.107560405098e-01,
0.272031289536e-01,
0.587399819634e-01,
0.108498349306e+00,
0.171685842349e+00,
0.232994786062e+00,
0.271377424940e+00,
0.271377424940e+00,
0.232994786062e+00,
0.171685842349e+00,
0.108498349306e+00,
0.587399819634e-01,
0.272031289536e-01,
0.107560405098e-01,
0.362258697852e-02,
0.103632909950e-02,
0.250983698512e-03,
0.512592913577e-04,
0.878849923082e-05,
0.125834025103e-05,
0.149553293672e-06,
0.146512531647e-07,
0.117361674232e-08,
0.761521725012e-10,
0.395917776693e-11,
0.162834073070e-12,
0.521862372645e-14,
0.128009339117e-15,
0.235188471067e-17,
0.315225456649e-19,
0.298286278427e-21,
0.191170688329e-23,
0.786179778889e-26,
0.192910359546e-28,
0.254966089910e-31,
0.155739062462e-34,
0.342113801099e-38,
0.167974799010e-42,
0.553570653584e-48
};
// points, starting with the first level
static const double gh_points[] = {
// points 1
0.0,
// points 2
-0.707106781186547524400844362105e+00,
0.707106781186547524400844362105e+00,
// points 3
-0.122474487139158904909864203735e+01,
0.0e+00,
0.122474487139158904909864203735e+01,
// points 4
-0.165068012388578455588334111112e+01,
-0.524647623275290317884060253835e+00,
0.524647623275290317884060253835e+00,
0.165068012388578455588334111112e+01,
// points 5
-0.202018287045608563292872408814e+01,
-0.958572464613818507112770593893e+00,
0.0e+00,
0.958572464613818507112770593893e+00,
0.202018287045608563292872408814e+01,
// points 6
-0.235060497367449222283392198706e+01,
-0.133584907401369694971489528297e+01,
-0.436077411927616508679215948251e+00,
0.436077411927616508679215948251e+00,
0.133584907401369694971489528297e+01,
0.235060497367449222283392198706e+01,
// points 7
-0.265196135683523349244708200652e+01,
-0.167355162876747144503180139830e+01,
-0.816287882858964663038710959027e+00,
0.0e+00,
0.816287882858964663038710959027e+00,
0.167355162876747144503180139830e+01,
0.265196135683523349244708200652e+01,
// points 8
-0.293063742025724401922350270524e+01,
-0.198165675669584292585463063977e+01,
-0.115719371244678019472076577906e+01,
-0.381186990207322116854718885584e+00,
0.381186990207322116854718885584e+00,
0.115719371244678019472076577906e+01,
0.198165675669584292585463063977e+01,
0.293063742025724401922350270524e+01,
// points 9
-0.319099320178152760723004779538e+01,
-0.226658058453184311180209693284e+01,
-0.146855328921666793166701573925e+01,
-0.723551018752837573322639864579e+00,
0.0e+00,
0.723551018752837573322639864579e+00,
0.146855328921666793166701573925e+01,
0.226658058453184311180209693284e+01,
0.319099320178152760723004779538e+01,
// points 10
-0.343615911883773760332672549432e+01,
-0.253273167423278979640896079775e+01,
-0.175668364929988177345140122011e+01,
-0.103661082978951365417749191676e+01,
-0.342901327223704608789165025557e+00,
0.342901327223704608789165025557e+00,
0.103661082978951365417749191676e+01,
0.175668364929988177345140122011e+01,
0.253273167423278979640896079775e+01,
0.343615911883773760332672549432e+01,
// points 11
-0.366847084655958251845837146485e+01,
-0.278329009978165177083671870152e+01,
-0.202594801582575533516591283121e+01,
-0.132655708449493285594973473558e+01,
-0.656809566882099765024611575383e+00,
0.0e+00,
0.656809566882099765024611575383e+00,
0.132655708449493285594973473558e+01,
0.202594801582575533516591283121e+01,
0.278329009978165177083671870152e+01,
0.366847084655958251845837146485e+01,
// points 12
-0.388972489786978191927164274724e+01,
-0.302063702512088977171067937518e+01,
-0.227950708050105990018772856942e+01,
-0.159768263515260479670966277090e+01,
-0.947788391240163743704578131060e+00,
-0.314240376254359111276611634095e+00,
0.314240376254359111276611634095e+00,
0.947788391240163743704578131060e+00,
0.159768263515260479670966277090e+01,
0.227950708050105990018772856942e+01,
0.302063702512088977171067937518e+01,
0.388972489786978191927164274724e+01,
// points 13
-0.410133759617863964117891508007e+01,
-0.324660897837240998812205115236e+01,
-0.251973568567823788343040913628e+01,
-0.185310765160151214200350644316e+01,
-0.122005503659074842622205526637e+01,
-0.605763879171060113080537108602e+00,
0.0e+00,
0.605763879171060113080537108602e+00,
0.122005503659074842622205526637e+01,
0.185310765160151214200350644316e+01,
0.251973568567823788343040913628e+01,
0.324660897837240998812205115236e+01,
0.410133759617863964117891508007e+01,
// points 14
-0.430444857047363181262129810037e+01,
-0.346265693360227055020891736115e+01,
-0.274847072498540256862499852415e+01,
-0.209518325850771681573497272630e+01,
-0.147668273114114087058350654421e+01,
-0.878713787329399416114679311861e+00,
-0.291745510672562078446113075799e+00,
0.291745510672562078446113075799e+00,
0.878713787329399416114679311861e+00,
0.147668273114114087058350654421e+01,
0.209518325850771681573497272630e+01,
0.274847072498540256862499852415e+01,
0.346265693360227055020891736115e+01,
0.430444857047363181262129810037e+01,
// points 15
-0.449999070730939155366438053053e+01,
-0.366995037340445253472922383312e+01,
-0.296716692790560324848896036355e+01,
-0.232573248617385774545404479449e+01,
-0.171999257518648893241583152515e+01,
-0.113611558521092066631913490556e+01,
-0.565069583255575748526020337198e+00,
0.0e+00,
0.565069583255575748526020337198e+00,
0.113611558521092066631913490556e+01,
0.171999257518648893241583152515e+01,
0.232573248617385774545404479449e+01,
0.296716692790560324848896036355e+01,
0.366995037340445253472922383312e+01,
0.449999070730939155366438053053e+01,
// points 16
-0.468873893930581836468849864875e+01,
-0.386944790486012269871942409801e+01,
-0.317699916197995602681399455926e+01,
-0.254620215784748136215932870545e+01,
-0.195178799091625397743465541496e+01,
-0.138025853919888079637208966969e+01,
-0.822951449144655892582454496734e+00,
-0.273481046138152452158280401965e+00,
0.273481046138152452158280401965e+00,
0.822951449144655892582454496734e+00,
0.138025853919888079637208966969e+01,
0.195178799091625397743465541496e+01,
0.254620215784748136215932870545e+01,
0.317699916197995602681399455926e+01,
0.386944790486012269871942409801e+01,
0.468873893930581836468849864875e+01,
// points 17
-0.487134519367440308834927655662e+01,
-0.406194667587547430689245559698e+01,
-0.337893209114149408338327069289e+01,
-0.275776291570388873092640349574e+01,
-0.217350282666662081927537907149e+01,
-0.161292431422123133311288254454e+01,
-0.106764872574345055363045773799e+01,
-0.531633001342654731349086553718e+00,
0.0e+00,
0.531633001342654731349086553718e+00,
0.106764872574345055363045773799e+01,
0.161292431422123133311288254454e+01,
0.217350282666662081927537907149e+01,
0.275776291570388873092640349574e+01,
0.337893209114149408338327069289e+01,
0.406194667587547430689245559698e+01,
0.487134519367440308834927655662e+01,
// points 18
-0.504836400887446676837203757885e+01,
-0.424811787356812646302342016090e+01,
-0.357376906848626607950067599377e+01,
-0.296137750553160684477863254906e+01,
-0.238629908916668600026459301424e+01,
-0.183553160426162889225383944409e+01,
-0.130092085838961736566626555439e+01,
-0.776682919267411661316659462284e+00,
-0.258267750519096759258116098711e+00,
0.258267750519096759258116098711e+00,
0.776682919267411661316659462284e+00,
0.130092085838961736566626555439e+01,
0.183553160426162889225383944409e+01,
0.238629908916668600026459301424e+01,
0.296137750553160684477863254906e+01,
0.357376906848626607950067599377e+01,
0.424811787356812646302342016090e+01,
0.504836400887446676837203757885e+01,
// points 19
-0.522027169053748216460967142500e+01,
-0.442853280660377943723498532226e+01,
-0.376218735196402009751489394104e+01,
-0.315784881834760228184318034120e+01,
-0.259113378979454256492128084112e+01,
-0.204923170985061937575050838669e+01,
-0.152417061939353303183354859367e+01,
-0.101036838713431135136859873726e+01,
-0.503520163423888209373811765050e+00,
0.0e+00,
0.503520163423888209373811765050e+00,
0.101036838713431135136859873726e+01,
0.152417061939353303183354859367e+01,
0.204923170985061937575050838669e+01,
0.259113378979454256492128084112e+01,
0.315784881834760228184318034120e+01,
0.376218735196402009751489394104e+01,
0.442853280660377943723498532226e+01,
0.522027169053748216460967142500e+01,
// points 20
-0.538748089001123286201690041068e+01,
-0.460368244955074427307767524898e+01,
-0.394476404011562521037562880052e+01,
-0.334785456738321632691492452300e+01,
-0.278880605842813048052503375640e+01,
-0.225497400208927552308233334473e+01,
-0.173853771211658620678086566214e+01,
-0.123407621539532300788581834696e+01,
-0.737473728545394358705605144252e+00,
-0.245340708300901249903836530634e+00,
0.245340708300901249903836530634e+00,
0.737473728545394358705605144252e+00,
0.123407621539532300788581834696e+01,
0.173853771211658620678086566214e+01,
0.225497400208927552308233334473e+01,
0.278880605842813048052503375640e+01,
0.334785456738321632691492452300e+01,
0.394476404011562521037562880052e+01,
0.460368244955074427307767524898e+01,
0.538748089001123286201690041068e+01,
// points 30
-6.86334529352989158106110835756e+00,
-6.13827922012393462039499237854e+00,
-5.53314715156749572511833355558e+00,
-4.98891896858994394448649710633e+00,
-4.48305535709251834188703761971e+00,
-4.00390860386122881522787601332e+00,
-3.54444387315534988692540090217e+00,
-3.09997052958644174868873332237e+00,
-2.66713212453561720057110646422e+00,
-2.24339146776150407247297999483e+00,
-1.82674114360368803883588048351e+00,
-1.41552780019818851194072510555e+00,
-1.00833827104672346180498960870e+00,
-0.603921058625552307778155678757e+00,
-0.201128576548871485545763013244e+00,
0.201128576548871485545763013244e+00,
0.603921058625552307778155678757e+00,
1.00833827104672346180498960870e+00,
1.41552780019818851194072510555e+00,
1.82674114360368803883588048351e+00,
2.24339146776150407247297999483e+00,
2.66713212453561720057110646422e+00,
3.09997052958644174868873332237e+00,
3.54444387315534988692540090217e+00,
4.00390860386122881522787601332e+00,
4.48305535709251834188703761971e+00,
4.98891896858994394448649710633e+00,
5.53314715156749572511833355558e+00,
6.13827922012393462039499237854e+00,
6.86334529352989158106110835756e+00,
// points 32
-7.12581390983e+00,
-6.40949814927e+00,
-5.81222594952e+00,
-5.27555098652e+00,
-4.77716450350e+00,
-4.30554795335e+00,
-3.85375548547e+00,
-3.41716749282e+00,
-2.99249082500e+00,
-2.57724953773e+00,
-2.16949918361e+00,
-1.76765410946e+00,
-1.37037641095e+00,
-0.976500463590e+00,
-0.584978765436e+00,
-0.194840741569e+00,
0.194840741569e+00,
0.584978765436e+00,
0.976500463590e+00,
1.37037641095e+00,
1.76765410946e+00,
2.16949918361e+00,
2.57724953773e+00,
2.99249082500e+00,
3.41716749282e+00,
3.85375548547e+00,
4.30554795335e+00,
4.77716450350e+00,
5.27555098652e+00,
5.81222594952e+00,
6.40949814927e+00,
7.12581390983e+00,
// points 40
-8.09876113925e+00,
-7.41158253149e+00,
-6.84023730525e+00,
-6.32825535122e+00,
-5.85409505603e+00,
-5.40665424797e+00,
-4.97926097855e+00,
-4.56750207284e+00,
-4.16825706683e+00,
-3.77920675344e+00,
-3.39855826586e+00,
-3.02487988390e+00,
-2.65699599844e+00,
-2.29391714188e+00,
-1.93479147228e+00,
-1.57886989493e+00,
-1.22548010905e+00,
-0.874006612357e+00,
-0.523874713832e+00,
-0.174537214598e+00,
0.174537214598e+00,
0.523874713832e+00,
0.874006612357e+00,
1.22548010905e+00,
1.57886989493e+00,
1.93479147228e+00,
2.29391714188e+00,
2.65699599844e+00,
3.02487988390e+00,
3.39855826586e+00,
3.77920675344e+00,
4.16825706683e+00,
4.56750207284e+00,
4.97926097855e+00,
5.40665424797e+00,
5.85409505603e+00,
6.32825535122e+00,
6.84023730525e+00,
7.41158253149e+00,
8.09876113925e+00,
// points 50
-9.18240695813e+00,
-8.52277103092e+00,
-7.97562236821e+00,
-7.48640942986e+00,
-7.03432350977e+00,
-6.60864797386e+00,
-6.20295251927e+00,
-5.81299467542e+00,
-5.43578608722e+00,
-5.06911758492e+00,
-4.71129366617e+00,
-4.36097316045e+00,
-4.01706817286e+00,
-3.67867706252e+00,
-3.34503831394e+00,
-3.01549776957e+00,
-2.68948470227e+00,
-2.36649390430e+00,
-2.04607196869e+00,
-1.72780654752e+00,
-1.41131775490e+00,
-1.09625112896e+00,
-0.782271729555e+00,
-0.469059056678e+00,
-0.156302546889e+00,
0.156302546889e+00,
0.469059056678e+00,
0.782271729555e+00,
1.09625112896e+00,
1.41131775490e+00,
1.72780654752e+00,
2.04607196869e+00,
2.36649390430e+00,
2.68948470227e+00,
3.01549776957e+00,
3.34503831394e+00,
3.67867706252e+00,
4.01706817286e+00,
4.36097316045e+00,
4.71129366617e+00,
5.06911758492e+00,
5.43578608722e+00,
5.81299467542e+00,
6.20295251927e+00,
6.60864797386e+00,
7.03432350977e+00,
7.48640942986e+00,
7.97562236821e+00,
8.52277103092e+00,
9.18240695813e+00,
// points 60
-10.1591092462e+00,
-9.52090367701e+00,
-8.99239800140e+00,
-8.52056928412e+00,
-8.08518865425e+00,
-7.67583993750e+00,
-7.28627659440e+00,
-6.91238153219e+00,
-6.55125916706e+00,
-6.20077355799e+00,
-5.85929019639e+00,
-5.52552108614e+00,
-5.19842653458e+00,
-4.87715007747e+00,
-4.56097375794e+00,
-4.24928643596e+00,
-3.94156073393e+00,
-3.63733587617e+00,
-3.33620465355e+00,
-3.03780333823e+00,
-2.74180374807e+00,
-2.44790690231e+00,
-2.15583787123e+00,
-1.86534153123e+00,
-1.57617901198e+00,
-1.28812467487e+00,
-1.00096349956e+00,
-0.714488781673e+00,
-0.428500064221e+00,
-0.142801238703e+00,
0.142801238703e+00,
0.428500064221e+00,
0.714488781673e+00,
1.00096349956e+00,
1.28812467487e+00,
1.57617901198e+00,
1.86534153123e+00,
2.15583787123e+00,
2.44790690231e+00,
2.74180374807e+00,
3.03780333823e+00,
3.33620465355e+00,
3.63733587617e+00,
3.94156073393e+00,
4.24928643596e+00,
4.56097375794e+00,
4.87715007747e+00,
5.19842653458e+00,
5.52552108614e+00,
5.85929019639e+00,
6.20077355799e+00,
6.55125916706e+00,
6.91238153219e+00,
7.28627659440e+00,
7.67583993750e+00,
8.08518865425e+00,
8.52056928412e+00,
8.99239800140e+00,
9.52090367701e+00,
10.1591092462e+00,
// points 64
-10.5261231680e+00,
-9.89528758683e+00,
-9.37315954965e+00,
-8.90724909996e+00,
-8.47752908338e+00,
-8.07368728501e+00,
-7.68954016404e+00,
-7.32101303278e+00,
-6.96524112055e+00,
-6.62011226264e+00,
-6.28401122877e+00,
-5.95566632680e+00,
-5.63405216435e+00,
-5.31832522463e+00,
-5.00777960220e+00,
-4.70181564741e+00,
-4.39991716823e+00,
-4.10163447457e+00,
-3.80657151395e+00,
-3.51437593574e+00,
-3.22473129199e+00,
-2.93735082300e+00,
-2.65197243543e+00,
-2.36835458863e+00,
-2.08627287988e+00,
-1.80551717147e+00,
-1.52588914021e+00,
-1.24720015694e+00,
-0.969269423071e+00,
-0.691922305810e+00,
-0.414988824121e+00,
-0.138302244987e+00,
0.138302244987e+00,
0.414988824121e+00,
0.691922305810e+00,
0.969269423071e+00,
1.24720015694e+00,
1.52588914021e+00,
1.80551717147e+00,
2.08627287988e+00,
2.36835458863e+00,
2.65197243543e+00,
2.93735082300e+00,
3.22473129199e+00,
3.51437593574e+00,
3.80657151395e+00,
4.10163447457e+00,
4.39991716823e+00,
4.70181564741e+00,
5.00777960220e+00,
5.31832522463e+00,
5.63405216435e+00,
5.95566632680e+00,
6.28401122877e+00,
6.62011226264e+00,
6.96524112055e+00,
7.32101303278e+00,
7.68954016404e+00,
8.07368728501e+00,
8.47752908338e+00,
8.90724909996e+00,
9.37315954965e+00,
9.89528758683e+00,
10.5261231680e+00
};
// Gauss-Legendre quadrature; prefix gl
// number of levels
static const int gl_num_levels = 22;
// number of points in each level
static const int gl_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
32, 64
};
// weights, starting with the first level
static const double gl_weights[] = {
// weight 1
2.0e+00,
// weights 2
1.0e+00,
1.0e+00,
// weights 3
0.555555555555555555555555555555e+00,
0.888888888888888888888888888888e+00,
0.555555555555555555555555555555e+00,
// weights 4
0.347854845137453857373063949222e+00,
0.652145154862546142626936050778e+00,
0.652145154862546142626936050778e+00,
0.347854845137453857373063949222e+00,
// weights 5
0.236926885056189087514264040720e+00,
0.478628670499366468041291514836e+00,
0.568888888888888888888888888889e+00,
0.478628670499366468041291514836e+00,
0.236926885056189087514264040720e+00,
// weights 6
0.171324492379170345040296142173e+00,
0.360761573048138607569833513838e+00,
0.467913934572691047389870343990e+00,
0.467913934572691047389870343990e+00,
0.360761573048138607569833513838e+00,
0.171324492379170345040296142173e+00,
// weights 7
0.129484966168869693270611432679e+00,
0.279705391489276667901467771424e+00,
0.381830050505118944950369775489e+00,
0.417959183673469387755102040816e+00,
0.381830050505118944950369775489e+00,
0.279705391489276667901467771424e+00,
0.129484966168869693270611432679e+00,
// weights 8
0.101228536290376259152531354310e+00,
0.222381034453374470544355994426e+00,
0.313706645877887287337962201987e+00,
0.362683783378361982965150449277e+00,
0.362683783378361982965150449277e+00,
0.313706645877887287337962201987e+00,
0.222381034453374470544355994426e+00,
0.101228536290376259152531354310e+00,
// weights 9
0.812743883615744119718921581105e-01,
0.180648160694857404058472031243e+00,
0.260610696402935462318742869419e+00,
0.312347077040002840068630406584e+00,
0.330239355001259763164525069287e+00,
0.312347077040002840068630406584e+00,
0.260610696402935462318742869419e+00,
0.180648160694857404058472031243e+00,
0.812743883615744119718921581105e-01,
// weights 10
0.666713443086881375935688098933e-01,
0.149451349150580593145776339658e+00,
0.219086362515982043995534934228e+00,
0.269266719309996355091226921569e+00,
0.295524224714752870173892994651e+00,
0.295524224714752870173892994651e+00,
0.269266719309996355091226921569e+00,
0.219086362515982043995534934228e+00,
0.149451349150580593145776339658e+00,
0.666713443086881375935688098933e-01,
// weights 11
0.556685671161736664827537204425e-01,
0.125580369464904624634694299224e+00,
0.186290210927734251426097641432e+00,
0.233193764591990479918523704843e+00,
0.262804544510246662180688869891e+00,
0.272925086777900630714483528336e+00,
0.262804544510246662180688869891e+00,
0.233193764591990479918523704843e+00,
0.186290210927734251426097641432e+00,
0.125580369464904624634694299224e+00,
0.556685671161736664827537204425e-01,
// weights 12
0.471753363865118271946159614850e-01,
0.106939325995318430960254718194e+00,
0.160078328543346226334652529543e+00,
0.203167426723065921749064455810e+00,
0.233492536538354808760849898925e+00,
0.249147045813402785000562436043e+00,
0.249147045813402785000562436043e+00,
0.233492536538354808760849898925e+00,
0.203167426723065921749064455810e+00,
0.160078328543346226334652529543e+00,
0.106939325995318430960254718194e+00,
0.471753363865118271946159614850e-01,
// weights 13
0.404840047653158795200215922010e-01,
0.921214998377284479144217759538e-01,
0.138873510219787238463601776869e+00,
0.178145980761945738280046691996e+00,
0.207816047536888502312523219306e+00,
0.226283180262897238412090186040e+00,
0.232551553230873910194589515269e+00,
0.226283180262897238412090186040e+00,
0.207816047536888502312523219306e+00,
0.178145980761945738280046691996e+00,
0.138873510219787238463601776869e+00,
0.921214998377284479144217759538e-01,
0.404840047653158795200215922010e-01,
// weights 14
0.351194603317518630318328761382e-01,
0.801580871597602098056332770629e-01,
0.121518570687903184689414809072e+00,
0.157203167158193534569601938624e+00,
0.185538397477937813741716590125e+00,
0.205198463721295603965924065661e+00,
0.215263853463157790195876443316e+00,
0.215263853463157790195876443316e+00,
0.205198463721295603965924065661e+00,
0.185538397477937813741716590125e+00,
0.157203167158193534569601938624e+00,
0.121518570687903184689414809072e+00,
0.801580871597602098056332770629e-01,
0.351194603317518630318328761382e-01,
// weights 15
0.307532419961172683546283935772e-01,
0.703660474881081247092674164507e-01,
0.107159220467171935011869546686e+00,
0.139570677926154314447804794511e+00,
0.166269205816993933553200860481e+00,
0.186161000015562211026800561866e+00,
0.198431485327111576456118326444e+00,
0.202578241925561272880620199968e+00,
0.198431485327111576456118326444e+00,
0.186161000015562211026800561866e+00,
0.166269205816993933553200860481e+00,
0.139570677926154314447804794511e+00,
0.107159220467171935011869546686e+00,
0.703660474881081247092674164507e-01,
0.307532419961172683546283935772e-01,
// weights 16
0.271524594117540948517805724560e-01,
0.622535239386478928628438369944e-01,
0.951585116824927848099251076022e-01,
0.124628971255533872052476282192e+00,
0.149595988816576732081501730547e+00,
0.169156519395002538189312079030e+00,
0.182603415044923588866763667969e+00,
0.189450610455068496285396723208e+00,
0.189450610455068496285396723208e+00,
0.182603415044923588866763667969e+00,
0.169156519395002538189312079030e+00,
0.149595988816576732081501730547e+00,
0.124628971255533872052476282192e+00,
0.951585116824927848099251076022e-01,
0.622535239386478928628438369944e-01,
0.271524594117540948517805724560e-01,
// weights 17
0.241483028685479319601100262876e-01,
0.554595293739872011294401653582e-01,
0.850361483171791808835353701911e-01,
0.111883847193403971094788385626e+00,
0.135136368468525473286319981702e+00,
0.154045761076810288081431594802e+00,
0.168004102156450044509970663788e+00,
0.176562705366992646325270990113e+00,
0.179446470356206525458265644262e+00,
0.176562705366992646325270990113e+00,
0.168004102156450044509970663788e+00,
0.154045761076810288081431594802e+00,
0.135136368468525473286319981702e+00,
0.111883847193403971094788385626e+00,
0.850361483171791808835353701911e-01,
0.554595293739872011294401653582e-01,
0.241483028685479319601100262876e-01,
// weights 18
0.216160135264833103133427102665e-01,
0.497145488949697964533349462026e-01,
0.764257302548890565291296776166e-01,
0.100942044106287165562813984925e+00,
0.122555206711478460184519126800e+00,
0.140642914670650651204731303752e+00,
0.154684675126265244925418003836e+00,
0.164276483745832722986053776466e+00,
0.169142382963143591840656470135e+00,
0.169142382963143591840656470135e+00,
0.164276483745832722986053776466e+00,
0.154684675126265244925418003836e+00,
0.140642914670650651204731303752e+00,
0.122555206711478460184519126800e+00,
0.100942044106287165562813984925e+00,
0.764257302548890565291296776166e-01,
0.497145488949697964533349462026e-01,
0.216160135264833103133427102665e-01,
// weights 19
0.194617882297264770363120414644e-01,
0.448142267656996003328381574020e-01,
0.690445427376412265807082580060e-01,
0.914900216224499994644620941238e-01,
0.111566645547333994716023901682e+00,
0.128753962539336227675515784857e+00,
0.142606702173606611775746109442e+00,
0.152766042065859666778855400898e+00,
0.158968843393954347649956439465e+00,
0.161054449848783695979163625321e+00,
0.158968843393954347649956439465e+00,
0.152766042065859666778855400898e+00,
0.142606702173606611775746109442e+00,
0.128753962539336227675515784857e+00,
0.111566645547333994716023901682e+00,
0.914900216224499994644620941238e-01,
0.690445427376412265807082580060e-01,
0.448142267656996003328381574020e-01,
0.194617882297264770363120414644e-01,
// weights 20
0.176140071391521183118619623519e-01,
0.406014298003869413310399522749e-01,
0.626720483341090635695065351870e-01,
0.832767415767047487247581432220e-01,
0.101930119817240435036750135480e+00,
0.118194531961518417312377377711e+00,
0.131688638449176626898494499748e+00,
0.142096109318382051329298325067e+00,
0.149172986472603746787828737002e+00,
0.152753387130725850698084331955e+00,
0.152753387130725850698084331955e+00,
0.149172986472603746787828737002e+00,
0.142096109318382051329298325067e+00,
0.131688638449176626898494499748e+00,
0.118194531961518417312377377711e+00,
0.101930119817240435036750135480e+00,
0.832767415767047487247581432220e-01,
0.626720483341090635695065351870e-01,
0.406014298003869413310399522749e-01,
0.176140071391521183118619623519e-01,
// weights 32
0.701861000947009660040706373885e-02,
0.162743947309056706051705622064e-01,
0.253920653092620594557525897892e-01,
0.342738629130214331026877322524e-01,
0.428358980222266806568786466061e-01,
0.509980592623761761961632446895e-01,
0.586840934785355471452836373002e-01,
0.658222227763618468376500637069e-01,
0.723457941088485062253993564785e-01,
0.781938957870703064717409188283e-01,
0.833119242269467552221990746043e-01,
0.876520930044038111427714627518e-01,
0.911738786957638847128685771116e-01,
0.938443990808045656391802376681e-01,
0.956387200792748594190820022041e-01,
0.965400885147278005667648300636e-01,
0.965400885147278005667648300636e-01,
0.956387200792748594190820022041e-01,
0.938443990808045656391802376681e-01,
0.911738786957638847128685771116e-01,
0.876520930044038111427714627518e-01,
0.833119242269467552221990746043e-01,
0.781938957870703064717409188283e-01,
0.723457941088485062253993564785e-01,
0.658222227763618468376500637069e-01,
0.586840934785355471452836373002e-01,
0.509980592623761761961632446895e-01,
0.428358980222266806568786466061e-01,
0.342738629130214331026877322524e-01,
0.253920653092620594557525897892e-01,
0.162743947309056706051705622064e-01,
0.701861000947009660040706373885e-02,
// weights 64
0.178328072169643294729607914497e-02,
0.414703326056246763528753572855e-02,
0.650445796897836285611736039998e-02,
0.884675982636394772303091465973e-02,
0.111681394601311288185904930192e-01,
0.134630478967186425980607666860e-01,
0.157260304760247193219659952975e-01,
0.179517157756973430850453020011e-01,
0.201348231535302093723403167285e-01,
0.222701738083832541592983303842e-01,
0.243527025687108733381775504091e-01,
0.263774697150546586716917926252e-01,
0.283396726142594832275113052002e-01,
0.302346570724024788679740598195e-01,
0.320579283548515535854675043479e-01,
0.338051618371416093915654821107e-01,
0.354722132568823838106931467152e-01,
0.370551285402400460404151018096e-01,
0.385501531786156291289624969468e-01,
0.399537411327203413866569261283e-01,
0.412625632426235286101562974736e-01,
0.424735151236535890073397679088e-01,
0.435837245293234533768278609737e-01,
0.445905581637565630601347100309e-01,
0.454916279274181444797709969713e-01,
0.462847965813144172959532492323e-01,
0.469681828162100173253262857546e-01,
0.475401657148303086622822069442e-01,
0.479993885964583077281261798713e-01,
0.483447622348029571697695271580e-01,
0.485754674415034269347990667840e-01,
0.486909570091397203833653907347e-01,
0.486909570091397203833653907347e-01,
0.485754674415034269347990667840e-01,
0.483447622348029571697695271580e-01,
0.479993885964583077281261798713e-01,
0.475401657148303086622822069442e-01,
0.469681828162100173253262857546e-01,
0.462847965813144172959532492323e-01,
0.454916279274181444797709969713e-01,
0.445905581637565630601347100309e-01,
0.435837245293234533768278609737e-01,
0.424735151236535890073397679088e-01,
0.412625632426235286101562974736e-01,
0.399537411327203413866569261283e-01,
0.385501531786156291289624969468e-01,
0.370551285402400460404151018096e-01,
0.354722132568823838106931467152e-01,
0.338051618371416093915654821107e-01,
0.320579283548515535854675043479e-01,
0.302346570724024788679740598195e-01,
0.283396726142594832275113052002e-01,
0.263774697150546586716917926252e-01,
0.243527025687108733381775504091e-01,
0.222701738083832541592983303842e-01,
0.201348231535302093723403167285e-01,
0.179517157756973430850453020011e-01,
0.157260304760247193219659952975e-01,
0.134630478967186425980607666860e-01,
0.111681394601311288185904930192e-01,
0.884675982636394772303091465973e-02,
0.650445796897836285611736039998e-02,
0.414703326056246763528753572855e-02,
0.178328072169643294729607914497e-02
};
// points, starting with the first level
static const double gl_points[] = {
// points 1
0.0e+00,
// points 2
-0.577350269189625764509148780502e+00,
0.577350269189625764509148780502e+00,
// points 3
-0.774596669241483377035853079956e+00,
0.0e+00,
0.774596669241483377035853079956e+00,
// points 4
-0.861136311594052575223946488893e+00,
-0.339981043584856264802665759103e+00,
0.339981043584856264802665759103e+00,
0.861136311594052575223946488893e+00,
// points 5
-0.906179845938663992797626878299e+00,
-0.538469310105683091036314420700e+00,
0.0e+00,
0.538469310105683091036314420700e+00,
0.906179845938663992797626878299e+00,
// points 6
-0.932469514203152027812301554494e+00,
-0.661209386466264513661399595020e+00,
-0.238619186083196908630501721681e+00,
0.238619186083196908630501721681e+00,
0.661209386466264513661399595020e+00,
0.932469514203152027812301554494e+00,
// points 7
-0.949107912342758524526189684048e+00,
-0.741531185599394439863864773281e+00,
-0.405845151377397166906606412077e+00,
0.0e+00,
0.405845151377397166906606412077e+00,
0.741531185599394439863864773281e+00,
0.949107912342758524526189684048e+00,
// points 8
-0.960289856497536231683560868569e+00,
-0.796666477413626739591553936476e+00,
-0.525532409916328985817739049189e+00,
-0.183434642495649804939476142360e+00,
0.183434642495649804939476142360e+00,
0.525532409916328985817739049189e+00,
0.796666477413626739591553936476e+00,
0.960289856497536231683560868569e+00,
// points 9
-0.968160239507626089835576202904e+00,
-0.836031107326635794299429788070e+00,
-0.613371432700590397308702039341e+00,
-0.324253423403808929038538014643e+00,
0.0e+00,
0.324253423403808929038538014643e+00,
0.613371432700590397308702039341e+00,
0.836031107326635794299429788070e+00,
0.968160239507626089835576202904e+00,
// points 10
-0.973906528517171720077964012084e+00,
-0.865063366688984510732096688423e+00,
-0.679409568299024406234327365115e+00,
-0.433395394129247190799265943166e+00,
-0.148874338981631210884826001130e+00,
0.148874338981631210884826001130e+00,
0.433395394129247190799265943166e+00,
0.679409568299024406234327365115e+00,
0.865063366688984510732096688423e+00,
0.973906528517171720077964012084e+00,
// points 11
-0.978228658146056992803938001123e+00,
-0.887062599768095299075157769304e+00,
-0.730152005574049324093416252031e+00,
-0.519096129206811815925725669459e+00,
-0.269543155952344972331531985401e+00,
0.0e+00,
0.269543155952344972331531985401e+00,
0.519096129206811815925725669459e+00,
0.730152005574049324093416252031e+00,
0.887062599768095299075157769304e+00,
0.978228658146056992803938001123e+00,
// points 12
-0.981560634246719250690549090149e+00,
-0.904117256370474856678465866119e+00,
-0.769902674194304687036893833213e+00,
-0.587317954286617447296702418941e+00,
-0.367831498998180193752691536644e+00,
-0.125233408511468915472441369464e+00,
0.125233408511468915472441369464e+00,
0.367831498998180193752691536644e+00,
0.587317954286617447296702418941e+00,
0.769902674194304687036893833213e+00,
0.904117256370474856678465866119e+00,
0.981560634246719250690549090149e+00,
// points 13
-0.984183054718588149472829448807e+00,
-0.917598399222977965206547836501e+00,
-0.801578090733309912794206489583e+00,
-0.642349339440340220643984606996e+00,
-0.448492751036446852877912852128e+00,
-0.230458315955134794065528121098e+00,
0.0e+00,
0.230458315955134794065528121098e+00,
0.448492751036446852877912852128e+00,
0.642349339440340220643984606996e+00,
0.801578090733309912794206489583e+00,
0.917598399222977965206547836501e+00,
0.984183054718588149472829448807e+00,
// points 14
-0.986283808696812338841597266704e+00,
-0.928434883663573517336391139378e+00,
-0.827201315069764993189794742650e+00,
-0.687292904811685470148019803019e+00,
-0.515248636358154091965290718551e+00,
-0.319112368927889760435671824168e+00,
-0.108054948707343662066244650220e+00,
0.108054948707343662066244650220e+00,
0.319112368927889760435671824168e+00,
0.515248636358154091965290718551e+00,
0.687292904811685470148019803019e+00,
0.827201315069764993189794742650e+00,
0.928434883663573517336391139378e+00,
0.986283808696812338841597266704e+00,
// points 15
-0.987992518020485428489565718587e+00,
-0.937273392400705904307758947710e+00,
-0.848206583410427216200648320774e+00,
-0.724417731360170047416186054614e+00,
-0.570972172608538847537226737254e+00,
-0.394151347077563369897207370981e+00,
-0.201194093997434522300628303395e+00,
0.0e+00,
0.201194093997434522300628303395e+00,
0.394151347077563369897207370981e+00,
0.570972172608538847537226737254e+00,
0.724417731360170047416186054614e+00,
0.848206583410427216200648320774e+00,
0.937273392400705904307758947710e+00,
0.987992518020485428489565718587e+00,
// points 16
-0.989400934991649932596154173450e+00,
-0.944575023073232576077988415535e+00,
-0.865631202387831743880467897712e+00,
-0.755404408355003033895101194847e+00,
-0.617876244402643748446671764049e+00,
-0.458016777657227386342419442984e+00,
-0.281603550779258913230460501460e+00,
-0.950125098376374401853193354250e-01,
0.950125098376374401853193354250e-01,
0.281603550779258913230460501460e+00,
0.458016777657227386342419442984e+00,
0.617876244402643748446671764049e+00,
0.755404408355003033895101194847e+00,
0.865631202387831743880467897712e+00,
0.944575023073232576077988415535e+00,
0.989400934991649932596154173450e+00,
// points 17
-0.990575475314417335675434019941e+00,
-0.950675521768767761222716957896e+00,
-0.880239153726985902122955694488e+00,
-0.781514003896801406925230055520e+00,
-0.657671159216690765850302216643e+00,
-0.512690537086476967886246568630e+00,
-0.351231763453876315297185517095e+00,
-0.178484181495847855850677493654e+00,
0.0e+00,
0.178484181495847855850677493654e+00,
0.351231763453876315297185517095e+00,
0.512690537086476967886246568630e+00,
0.657671159216690765850302216643e+00,
0.781514003896801406925230055520e+00,
0.880239153726985902122955694488e+00,
0.950675521768767761222716957896e+00,
0.990575475314417335675434019941e+00,
// points 18
-0.991565168420930946730016004706e+00,
-0.955823949571397755181195892930e+00,
-0.892602466497555739206060591127e+00,
-0.803704958972523115682417455015e+00,
-0.691687043060353207874891081289e+00,
-0.559770831073947534607871548525e+00,
-0.411751161462842646035931793833e+00,
-0.251886225691505509588972854878e+00,
-0.847750130417353012422618529358e-01,
0.847750130417353012422618529358e-01,
0.251886225691505509588972854878e+00,
0.411751161462842646035931793833e+00,
0.559770831073947534607871548525e+00,
0.691687043060353207874891081289e+00,
0.803704958972523115682417455015e+00,
0.892602466497555739206060591127e+00,
0.955823949571397755181195892930e+00,
0.991565168420930946730016004706e+00,
// points 19
-0.992406843843584403189017670253e+00,
-0.960208152134830030852778840688e+00,
-0.903155903614817901642660928532e+00,
-0.822714656537142824978922486713e+00,
-0.720966177335229378617095860824e+00,
-0.600545304661681023469638164946e+00,
-0.464570741375960945717267148104e+00,
-0.316564099963629831990117328850e+00,
-0.160358645640225375868096115741e+00,
0.0e+00,
0.160358645640225375868096115741e+00,
0.316564099963629831990117328850e+00,
0.464570741375960945717267148104e+00,
0.600545304661681023469638164946e+00,
0.720966177335229378617095860824e+00,
0.822714656537142824978922486713e+00,
0.903155903614817901642660928532e+00,
0.960208152134830030852778840688e+00,
0.992406843843584403189017670253e+00,
// points 20
-0.993128599185094924786122388471e+00,
-0.963971927277913791267666131197e+00,
-0.912234428251325905867752441203e+00,
-0.839116971822218823394529061702e+00,
-0.746331906460150792614305070356e+00,
-0.636053680726515025452836696226e+00,
-0.510867001950827098004364050955e+00,
-0.373706088715419560672548177025e+00,
-0.227785851141645078080496195369e+00,
-0.765265211334973337546404093988e-01,
0.765265211334973337546404093988e-01,
0.227785851141645078080496195369e+00,
0.373706088715419560672548177025e+00,
0.510867001950827098004364050955e+00,
0.636053680726515025452836696226e+00,
0.746331906460150792614305070356e+00,
0.839116971822218823394529061702e+00,
0.912234428251325905867752441203e+00,
0.963971927277913791267666131197e+00,
0.993128599185094924786122388471e+00,
// points 32
-0.997263861849481563544981128665e+00,
-0.985611511545268335400175044631e+00,
-0.964762255587506430773811928118e+00,
-0.934906075937739689170919134835e+00,
-0.896321155766052123965307243719e+00,
-0.849367613732569970133693004968e+00,
-0.794483795967942406963097298970e+00,
-0.732182118740289680387426665091e+00,
-0.663044266930215200975115168663e+00,
-0.587715757240762329040745476402e+00,
-0.506899908932229390023747474378e+00,
-0.421351276130635345364119436172e+00,
-0.331868602282127649779916805730e+00,
-0.239287362252137074544603209166e+00,
-0.144471961582796493485186373599e+00,
-0.483076656877383162348125704405e-01,
0.483076656877383162348125704405e-01,
0.144471961582796493485186373599e+00,
0.239287362252137074544603209166e+00,
0.331868602282127649779916805730e+00,
0.421351276130635345364119436172e+00,
0.506899908932229390023747474378e+00,
0.587715757240762329040745476402e+00,
0.663044266930215200975115168663e+00,
0.732182118740289680387426665091e+00,
0.794483795967942406963097298970e+00,
0.849367613732569970133693004968e+00,
0.896321155766052123965307243719e+00,
0.934906075937739689170919134835e+00,
0.964762255587506430773811928118e+00,
0.985611511545268335400175044631e+00,
0.997263861849481563544981128665e+00,
// points 64
-0.999305041735772139456905624346e+00,
-0.996340116771955279346924500676e+00,
-0.991013371476744320739382383443e+00,
-0.983336253884625956931299302157e+00,
-0.973326827789910963741853507352e+00,
-0.961008799652053718918614121897e+00,
-0.946411374858402816062481491347e+00,
-0.929569172131939575821490154559e+00,
-0.910522137078502805756380668008e+00,
-0.889315445995114105853404038273e+00,
-0.865999398154092819760783385070e+00,
-0.840629296252580362751691544696e+00,
-0.813265315122797559741923338086e+00,
-0.783972358943341407610220525214e+00,
-0.752819907260531896611863774886e+00,
-0.719881850171610826848940217832e+00,
-0.685236313054233242563558371031e+00,
-0.648965471254657339857761231993e+00,
-0.611155355172393250248852971019e+00,
-0.571895646202634034283878116659e+00,
-0.531279464019894545658013903544e+00,
-0.489403145707052957478526307022e+00,
-0.446366017253464087984947714759e+00,
-0.402270157963991603695766771260e+00,
-0.357220158337668115950442615046e+00,
-0.311322871990210956157512698560e+00,
-0.264687162208767416373964172510e+00,
-0.217423643740007084149648748989e+00,
-0.169644420423992818037313629748e+00,
-0.121462819296120554470376463492e+00,
-0.729931217877990394495429419403e-01,
-0.243502926634244325089558428537e-01,
0.243502926634244325089558428537e-01,
0.729931217877990394495429419403e-01,
0.121462819296120554470376463492e+00,
0.169644420423992818037313629748e+00,
0.217423643740007084149648748989e+00,
0.264687162208767416373964172510e+00,
0.311322871990210956157512698560e+00,
0.357220158337668115950442615046e+00,
0.402270157963991603695766771260e+00,
0.446366017253464087984947714759e+00,
0.489403145707052957478526307022e+00,
0.531279464019894545658013903544e+00,
0.571895646202634034283878116659e+00,
0.611155355172393250248852971019e+00,
0.648965471254657339857761231993e+00,
0.685236313054233242563558371031e+00,
0.719881850171610826848940217832e+00,
0.752819907260531896611863774886e+00,
0.783972358943341407610220525214e+00,
0.813265315122797559741923338086e+00,
0.840629296252580362751691544696e+00,
0.865999398154092819760783385070e+00,
0.889315445995114105853404038273e+00,
0.910522137078502805756380668008e+00,
0.929569172131939575821490154559e+00,
0.946411374858402816062481491347e+00,
0.961008799652053718918614121897e+00,
0.973326827789910963741853507352e+00,
0.983336253884625956931299302157e+00,
0.991013371476744320739382383443e+00,
0.996340116771955279346924500676e+00,
0.999305041735772139456905624346e+00
};
// this is the positive half of normal inverse cum distribution
// function starting at 0.5 and ending at 0.998, with step 0.002
static const int normal_icdf_num = 250;
static const double normal_icdf_end = 0.998;
static const double normal_icdf_step = 0.002;
static const double normal_icdf_data[] = {
0, 5.013277548926632e-03, 1.002668110027482e-02,
1.504033667863573e-02, 2.005437035295075e-02, 2.506890825871118e-02,
3.008407662018906e-02, 3.510000177270896e-02, 4.011681018496811e-02,
4.513462848142118e-02, 5.015358346473358e-02, 5.517380213831685e-02,
6.019541172895673e-02, 6.521853970954372e-02, 7.024331382191684e-02,
7.526986209982979e-02, 8.029831289205518e-02, 8.532879488562921e-02,
9.036143712925872e-02, 9.539636905689193e-02, 1.004337205114700e-01,
1.054736217688682e-01, 1.105162035620419e-01, 1.155615971053833e-01,
1.206099341193073e-01, 1.256613468550742e-01, 1.307159681198632e-01,
1.357739313021116e-01, 1.408353703971274e-01, 1.459004200329941e-01,
1.509692154967774e-01, 1.560418927610502e-01, 1.611185885107454e-01,
1.661994401703590e-01, 1.712845859315068e-01, 1.763741647808615e-01,
1.814683165284770e-01, 1.865671818365194e-01, 1.916709022484199e-01,
1.967796202184666e-01, 2.018934791418509e-01, 2.070126233851871e-01,
2.121371983175242e-01, 2.172673503418634e-01, 2.224032269272064e-01,
2.275449766411493e-01, 2.326927491830447e-01, 2.378466954177492e-01,
2.430069674099821e-01, 2.481737184593126e-01, 2.533471031357997e-01,
2.585272773163098e-01, 2.637143982215299e-01, 2.689086244537098e-01,
2.741101160351471e-01, 2.793190344474543e-01, 2.845355426716215e-01,
2.897598052289143e-01, 2.949919882226262e-01, 3.002322593807220e-01,
3.054807880993972e-01, 3.107377454875922e-01, 3.160033044124830e-01,
3.212776395459965e-01, 3.265609274123727e-01, 3.318533464368166e-01,
3.371550769952773e-01, 3.424663014653906e-01, 3.477872042786273e-01,
3.531179719736894e-01, 3.584587932511938e-01, 3.638098590296960e-01,
3.691713625030897e-01, 3.745434991994428e-01, 3.799264670413076e-01,
3.853204664075677e-01, 3.907257001968699e-01, 3.961423738926983e-01,
4.015706956301487e-01, 4.070108762644656e-01, 4.124631294414047e-01,
4.179276716694820e-01, 4.234047223941831e-01, 4.288945040742017e-01,
4.343972422597815e-01, 4.399131656732339e-01, 4.454425062917200e-01,
4.509854994323708e-01, 4.565423838398405e-01, 4.621134017763774e-01,
4.676987991145082e-01, 4.732988254324370e-01, 4.789137341122557e-01,
4.845437824410792e-01, 4.901892317152095e-01, 4.958503473474533e-01,
5.015273989777081e-01, 5.072206605869456e-01, 5.129304106147284e-01,
5.186569320803909e-01, 5.244005127080407e-01, 5.301614450555191e-01,
5.359400266474903e-01, 5.417365601128169e-01, 5.475513533264015e-01,
5.533847195556728e-01, 5.592369776119069e-01, 5.651084520065839e-01,
5.709994731129874e-01, 5.769103773332714e-01, 5.828415072712163e-01,
5.887932119109195e-01, 5.947658468016782e-01, 6.007597742493188e-01,
6.067753635142652e-01, 6.128129910166273e-01, 6.188730405486286e-01,
6.249559034946875e-01, 6.310619790594989e-01, 6.371916745044747e-01,
6.433454053929173e-01, 6.495235958443252e-01, 6.557266787982537e-01,
6.619550962881621e-01, 6.682092997257233e-01, 6.744897501960819e-01,
6.807969187645747e-01, 6.871312867954694e-01, 6.934933462832894e-01,
6.998836001973414e-01, 7.063025628400875e-01, 7.127507602200432e-01,
7.192287304399239e-01, 7.257370241008051e-01, 7.322762047230997e-01,
7.388468491852137e-01, 7.454495481807891e-01, 7.520849066954916e-01,
7.587535445043710e-01, 7.654560966908778e-01, 7.721932141886847e-01,
7.789655643475453e-01, 7.857738315244843e-01, 7.926187177017122e-01,
7.995009431327367e-01, 8.064212470182405e-01, 8.133803882134047e-01,
8.203791459684610e-01, 8.274183207043821e-01, 8.344987348257406e-01,
8.416212335729145e-01, 8.487866859159668e-01, 8.559959854926823e-01,
8.632500515934207e-01, 8.705498301956541e-01, 8.778962950512290e-01,
8.852904488296417e-01, 8.927333243208563e-01, 9.002259857014339e-01,
9.077695298680560e-01, 9.153650878428145e-01, 9.230138262549803e-01,
9.307169489043392e-01, 9.384756984115684e-01, 9.462913579615760e-01,
9.541652531461944e-01, 9.620987539131418e-01, 9.700932766287370e-01,
9.781502862624715e-01, 9.862712987022384e-01, 9.944578832097529e-01,
1.002711665026549e+00, 1.011034328141817e+00, 1.019427618234370e+00,
1.027893345802143e+00, 1.036433389493790e+00, 1.045049699658389e+00,
1.053744302130666e+00, 1.062519302270867e+00, 1.071376889280213e+00,
1.080319340814956e+00, 1.089349027924277e+00, 1.098468420339863e+00,
1.107680092147800e+00, 1.116986727876610e+00, 1.126391129038801e+00,
1.135896221167312e+00, 1.145505061392697e+00, 1.155220846611952e+00,
1.165046922305602e+00, 1.174986792066090e+00, 1.185044127907810e+00,
1.195222781437427e+00, 1.205526795972518e+00, 1.215960419707319e+00,
1.226528120036610e+00, 1.237234599162827e+00, 1.248084811127547e+00,
1.259083980427072e+00, 1.270237622393149e+00, 1.281551565544601e+00,
1.293031976144243e+00, 1.304685385228790e+00, 1.316518718418261e+00,
1.328539328856810e+00, 1.340755033690217e+00, 1.353174154548003e+00,
1.365805562572272e+00, 1.378658728623277e+00, 1.391743779396326e+00,
1.405071560309632e+00, 1.418653706172739e+00, 1.432502720825812e+00,
1.446632067158978e+00, 1.461056269186906e+00, 1.475791028179170e+00,
1.490853355246661e+00, 1.506261723278244e+00, 1.522036241735856e+00,
1.538198858584064e+00, 1.554773594596853e+00, 1.571786816509860e+00,
1.589267557051392e+00, 1.607247891900218e+00, 1.625763386233235e+00,
1.644853626951473e+00, 1.664562861202721e+00, 1.684940767871913e+00,
1.706043396888962e+00, 1.727934322388419e+00, 1.750686071252170e+00,
1.774381910344958e+00, 1.799118106837967e+00, 1.825006821146403e+00,
1.852179858769047e+00, 1.880793608151250e+00, 1.911035647549119e+00,
1.943133751105067e+00, 1.977368428181947e+00, 2.014090812018140e+00,
2.053748910631823e+00, 2.096927429164343e+00, 2.144410620911840e+00,
2.197286376641053e+00, 2.257129244486226e+00, 2.326347874040842e+00,
2.408915545815460e+00, 2.512144327930459e+00, 2.652069807902199e+00,
2.878161739095476e+00
};
// Local Variables:
// mode:C++
// End:
dynare++/integ/cc/product.cweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: product.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@ This is {\tt product.cpp} file.
@c
#include "product.h"
#include "symmetry.h"
@<|prodpit| empty constructor@>;
@<|prodpit| regular constructor@>;
@<|prodpit| copy constructor@>;
@<|prodpit| destructor@>;
@<|prodpit::operator==| code@>;
@<|prodpit::operator=| code@>;
@<|prodpit::operator++| code@>;
@<|prodpit::setPointAndWeight| code@>;
@<|prodpit::print| code@>;
@<|ProductQuadrature| constructor@>;
@<|ProductQuadrature::begin| code@>;
@<|ProductQuadrature::designLevelForEvals| code@>;
@
@<|prodpit| empty constructor@>=
prodpit::prodpit()
: prodq(NULL), level(0), npoints(0), jseq(NULL),
end_flag(true), sig(NULL), p(NULL)
{
}
@ This constructs a product iterator corresponding to index $(j0,0\ldots,0)$.
@<|prodpit| regular constructor@>=
prodpit::prodpit(const ProductQuadrature& q, int j0, int l)
: prodq(&q), level(l), npoints(q.uquad.numPoints(l)), jseq(new IntSequence(q.dimen(), 0)),
end_flag(false), sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
{
if (j0 < npoints) {
(*jseq)[0] = j0;
setPointAndWeight();
} else {
end_flag = true;
}
}
@ Copy constructor, clear.
@<|prodpit| copy constructor@>=
prodpit::prodpit(const prodpit& ppit)
: prodq(ppit.prodq), level(ppit.level), npoints(ppit.npoints),
end_flag(ppit.end_flag), w(ppit.w)
{
if (ppit.jseq)
jseq = new IntSequence(*(ppit.jseq));
else
jseq = NULL;
if (ppit.sig)
sig = new ParameterSignal(*(ppit.sig));
else
sig = NULL;
if (ppit.p)
p = new Vector(*(ppit.p));
else
p = NULL;
}
@
@<|prodpit| destructor@>=
prodpit::~prodpit()
{
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
}
@
@<|prodpit::operator==| code@>=
bool prodpit::operator==(const prodpit& ppit) const
{
bool ret = true;
ret = ret & prodq == ppit.prodq;
ret = ret & end_flag == ppit.end_flag;
ret = ret & ((jseq==NULL && ppit.jseq==NULL) ||
(jseq!=NULL && ppit.jseq!=NULL && *jseq == *(ppit.jseq)));
return ret;
}
@
@<|prodpit::operator=| code@>=
const prodpit& prodpit::operator=(const prodpit& ppit)
{
prodq = ppit.prodq;
end_flag = ppit.end_flag;
w = ppit.w;
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
if (ppit.jseq)
jseq = new IntSequence(*(ppit.jseq));
else
jseq = NULL;
if (ppit.sig)
sig = new ParameterSignal(*(ppit.sig));
else
sig = NULL;
if (ppit.p)
p = new Vector(*(ppit.p));
else
p = NULL;
return *this;
}
@
@<|prodpit::operator++| code@>=
prodpit& prodpit::operator++()
{
// todo: throw if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or |end_flag==true|
int i = prodq->dimen()-1;
(*jseq)[i]++;
while (i >= 0 && (*jseq)[i] == npoints) {
(*jseq)[i] = 0;
i--;
if (i >= 0)
(*jseq)[i]++;
}
sig->signalAfter(std::max(i,0));
if (i == -1)
end_flag = true;
if (! end_flag)
setPointAndWeight();
return *this;
}
@ This calculates the weight and sets point coordinates from the indices.
@<|prodpit::setPointAndWeight| code@>=
void prodpit::setPointAndWeight()
{
// todo: raise if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or
// |p==NULL| or |end_flag==true|
w = 1.0;
for (int i = 0; i < prodq->dimen(); i++) {
(*p)[i] = (prodq->uquad).point(level, (*jseq)[i]);
w* = (prodq->uquad).weight(level, (*jseq)[i]);
}
}
@ Debug print.
@<|prodpit::print| code@>=
void prodpit::print() const
{
printf("j=[");
for (int i = 0; i < prodq->dimen(); i++)
printf("%2d ", (*jseq)[i]);
printf("] %+4.3f*(",w);
for (int i = 0; i < prodq->dimen()-1; i++)
printf("%+4.3f ", (*p)[i]);
printf("%+4.3f)\n",(*p)[prodq->dimen()-1]);
}
@
@<|ProductQuadrature| constructor@>=
ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature& uq)
: QuadratureImpl<prodpit>(d), uquad(uq)
{
// todo: check |d>=1|
}
@ This calls |prodpit| constructor to return an iterator which points
approximatelly at |ti|-th portion out of |tn| portions. First we find
out how many points are in the level, and then construct an interator
$(j0,0,\ldots,0)$ where $j0=$|ti*npoints/tn|.
@<|ProductQuadrature::begin| code@>=
prodpit ProductQuadrature::begin(int ti, int tn, int l) const
{
// todo: raise is |l<dimen()|
// todo: check |l<=uquad.numLevels()|
int npoints = uquad.numPoints(l);
return prodpit(*this, ti*npoints/tn, l);
}
@ This just starts at the first level and goes to a higher level as
long as a number of evaluations (which is $n_k^d$ for $k$ being the
level) is less than the given number of evaluations.
@<|ProductQuadrature::designLevelForEvals| code@>=
void ProductQuadrature::designLevelForEvals(int max_evals, int& lev, int& evals) const
{
int last_evals;
evals = 1;
lev = 1;
do {
lev++;
last_evals = evals;
evals = numEvals(lev);
} while (lev < uquad.numLevels()-2 && evals < max_evals);
lev--;
evals = last_evals;
}
@ End of {\tt product.cpp} file
dynare++/integ/cc/product.hweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: product.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@*2 Product quadrature. This is {\tt product.h} file
This file defines a product multidimensional quadrature. If $Q_k$
denotes the one dimensional quadrature, then the product quadrature
$Q$ of $k$ level and dimension $d$ takes the form
$$Qf=\sum_{i_1=1}^{n_k}\ldots\sum_{i_d=1}^{n^k}w_{i_1}\cdot\ldots\cdot w_{i_d}
f(x_{i_1},\ldots,x_{i_d})$$
which can be written in terms of the one dimensional quadrature $Q_k$ as
$$Qf=(Q_k\otimes\ldots\otimes Q_k)f$$
Here we define the product quadrature iterator |prodpit| and plug it
into |QuadratureImpl| to obtains |ProductQuadrature|.
@s prodpit int
@s ProductQuadrature int
@c
#ifndef PRODUCT_H
#define PRODUCT_H
#include "int_sequence.h"
#include "vector_function.h"
#include "quadrature.h"
@<|prodpit| class declaration@>;
@<|ProductQuadrature| class declaration@>;
#endif
@ This defines a product point iterator. We have to maintain the
following: a pointer to product quadrature in order to know the
dimension and the underlying one dimensional quadrature, then level,
number of points in the level, integer sequence of indices, signal,
the coordinates of the point and the weight.
The point indices, signal, and point coordinates are implmented as
pointers in order to allow for empty constructor.
The constructor |prodpit(const ProductQuadrature& q, int j0, int l)|
constructs an iterator pointing to $(j0,0,\ldots,0)$, which is used by
|begin| dictated by |QuadratureImpl|.
@<|prodpit| class declaration@>=
class ProductQuadrature;
class prodpit {
protected:@;
const ProductQuadrature* prodq;
int level;
int npoints;
IntSequence* jseq;
bool end_flag;
ParameterSignal* sig;
Vector* p;
double w;
public:@;
prodpit();
prodpit(const ProductQuadrature& q, int j0, int l);
prodpit(const prodpit& ppit);
~prodpit();
bool operator==(const prodpit& ppit) const;
bool operator!=(const prodpit& ppit) const
{@+ return ! operator==(ppit);@+}
const prodpit& operator=(const prodpit& spit);
prodpit& operator++();
const ParameterSignal& signal() const
{@+ return *sig;@+}
const Vector& point() const
{@+ return *p;@+}
double weight() const
{@+ return w;@+}
void print() const;
protected:@;
void setPointAndWeight();
};
@ The product quadrature is just |QuadratureImpl| with the product
iterator plugged in. The object is constructed by just giving the
underlying one dimensional quadrature, and the dimension. The only
extra method is |designLevelForEvals| which for the given maximum
number of evaluations (and dimension and underlying quadrature from
the object) returns a maximum level yeilding number of evaluations
less than the given number.
@<|ProductQuadrature| class declaration@>=
class ProductQuadrature : public QuadratureImpl<prodpit> {
friend class prodpit;
const OneDQuadrature& uquad;
public:@;
ProductQuadrature(int d, const OneDQuadrature& uq);
virtual ~ProductQuadrature()@+ {}
int numEvals(int l) const
{
int res = 1;
for (int i = 0; i < dimen(); i++)
res *= uquad.numPoints(l);
return res;
}
void designLevelForEvals(int max_eval, int& lev, int& evals) const;
protected:@;
prodpit begin(int ti, int tn, int level) const;
};
@ End of {\tt product.h} file
dynare++/integ/cc/quadrature.cweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: quadrature.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@ This is {\tt quadrature.cpp} file.
@c
#include "quadrature.h"
#include "precalc_quadrature.dat"
#include <cmath>
@<|OneDPrecalcQuadrature::calcOffsets| code@>;
@<|GaussHermite| constructor code@>;
@<|GaussLegendre| constructor code@>;
@<|NormalICDF| get code@>;
@
@<|OneDPrecalcQuadrature::calcOffsets| code@>=
void OneDPrecalcQuadrature::calcOffsets()
{
offsets[0] = 0;
for (int i = 1; i < num_levels; i++)
offsets[i] = offsets[i-1] + num_points[i-1];
}
@
@<|GaussHermite| constructor code@>=
GaussHermite::GaussHermite()
: OneDPrecalcQuadrature(gh_num_levels, gh_num_points, gh_weights, gh_points)@+ {}
@
@<|GaussLegendre| constructor code@>=
GaussLegendre::GaussLegendre()
: OneDPrecalcQuadrature(gl_num_levels, gl_num_points, gl_weights, gl_points)@+ {}
@ Here we transform a draw from univariate $\langle 0,1\rangle$ to the
draw from Gaussina $N(0,1)$. This is done by a table lookup, the table
is given by |normal_icdf_step|, |normal_icfd_data|, |normal_icdf_num|,
and a number |normal_icdf_end|. In order to avoid wrong tails for lookups close
to zero or one, we rescale input |x| by $(1-2*(1-end))=2*end-1$.
@<|NormalICDF| get code@>=
double NormalICDF::get(double x)
{
double xx = (2*normal_icdf_end-1)*std::abs(x-0.5);
int i = (int)floor(xx/normal_icdf_step);
double xx1 = normal_icdf_step*i;
double yy1 = normal_icdf_data[i];
double y;
if (i < normal_icdf_num-1) {
double yy2 = normal_icdf_data[i+1];
y = yy1 + (yy2-yy1)*(xx-xx1)/normal_icdf_step;
} else { // this should never happen
y = yy1;
}
if (x > 0.5)
return y;
else
return -y;
}
@ End of {\tt quadrature.cpp} file
dynare++/integ/cc/quadrature.hweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: quadrature.hweb 2269 2008-11-23 14:33:22Z michel $ @>
@q Copyright 2005, Ondra Kamenik @>
@*2 Quadrature. This is {\tt quadrature.h} file
This file defines an interface for one dimensional (non-nested) quadrature
|OneDQuadrature|, and a parent for all multi-dimensional
quadratures. This parent class |Quadrature| presents a general concept of
quadrature, this is
$$\int f(x){\rm d}x \approx\sum_{i=1}^N w_ix_i$$
The class |Quadrature| just declares this concept. The concept is
implemented by class |QuadratureImpl| which paralelizes the
summation. All implementations therefore wishing to use the parallel
implementation should inherit from |QuadratureImpl| and integration is
done.
The integration concept relies on a point iterator, which goes through
all $x_i$ and $w_i$ for $i=1,\ldots,N$. All the iterators must be able
to go through only a portion of the set $i=1,\ldots,N$. This enables
us to implement paralelism, for two threads for example, one iterator
goes from the beginning to the (approximately) half, and the other
goes from the half to the end.
Besides this concept of the general quadrature, this file defines also
one dimensional quadrature, which is basically a scheme of points and
weights for different levels. The class |OneDQuadrature| is a parent
of all such objects, the classes |GaussHermite| and |GaussLegendre|
are specific implementations for Gauss--Hermite and Gauss--Legendre
quadratures resp.
@s OneDQuadrature int
@s Quadrature int
@s IntegrationWorker int
@s QuadratureImpl int
@s OneDPrecalcQuadrature int
@s GaussHermite int
@s GaussLegendre int
@s NormalICDF int
@s _Tpit int
@c
#ifndef QUADRATURE_H
#define QUADRATURE_H
#include <cstdlib>
#include "vector_function.h"
#include "int_sequence.h"
#include "sthread.h"
@<|OneDQuadrature| class declaration@>;
@<|Quadrature| class declaration@>;
@<|IntegrationWorker| class declaration@>;
@<|QuadratureImpl| class declaration@>;
@<|OneDPrecalcQuadrature| class declaration@>;
@<|GaussHermite| class declaration@>;
@<|GaussLegendre| class declaration@>;
@<|NormalICDF| class declaration@>;
#endif
@ This pure virtual class represents a concept of one-dimensional
(non-nested) quadrature. So, one dimensional quadrature must return
number of levels, number of points in a given level, and then a point
and a weight in a given level and given order.
@<|OneDQuadrature| class declaration@>=
class OneDQuadrature {
public:@;
virtual ~OneDQuadrature()@+ {}
virtual int numLevels() const =0;
virtual int numPoints(int level) const =0;
virtual double point(int level, int i) const =0;
virtual double weight(int lelel, int i) const =0;
};
@ This is a general concept of multidimensional quadrature. at this
general level, we maintain only a dimension, and declare virtual
functions for integration. The function take two forms; first takes a
constant |VectorFunction| as an argument, creates locally
|VectorFunctionSet| and do calculation, second one takes as an
argument |VectorFunctionSet|.
Part of the interface is a method returning a number of evaluations
for a specific level. Note two things: this assumes that the number of
evaluations is known apriori and thus it is not applicable for
adaptive quadratures, second for Monte Carlo type of quadrature, the
level is a number of evaluations.
@<|Quadrature| class declaration@>=
class Quadrature {
protected:@;
int dim;
public:@;
Quadrature(int d) : dim(d)@+ {}
virtual ~Quadrature()@+ {}
int dimen() const
{@+ return dim;@+}
virtual void integrate(const VectorFunction& func, int level,
int tn, Vector& out) const =0;
virtual void integrate(VectorFunctionSet& fs, int level, Vector& out) const =0;
virtual int numEvals(int level) const =0;
};
@ This is just an integration worker, which works over a given
|QuadratureImpl|. It also needs the function, level, a specification
of the subgroup of points, and output vector.
See |@<|QuadratureImpl| class declaration@>| for details.
@<|IntegrationWorker| class declaration@>=
template <typename _Tpit>
class QuadratureImpl;
template <typename _Tpit>
class IntegrationWorker : public THREAD {
const QuadratureImpl<_Tpit>& quad;
VectorFunction& func;
int level;
int ti;
int tn;
Vector& outvec;
public:@;
IntegrationWorker(const QuadratureImpl<_Tpit>& q, VectorFunction& f, int l,
int tii, int tnn, Vector& out)
: quad(q), func(f), level(l), ti(tii), tn(tnn), outvec(out) @+{}
@<|IntegrationWorker::operator()()| code@>;
};
@ This integrates the given portion of the integral. We obtain first
and last iterators for the portion (|beg| and |end|). Then we iterate
through the portion. and finally we add the intermediate result to the
result |outvec|.
This method just everything up as it is coming. This might be imply
large numerical errors, perhaps in future I will implement something
smarter.
@<|IntegrationWorker::operator()()| code@>=
void operator()() {
_Tpit beg = quad.begin(ti, tn, level);
_Tpit end = quad.begin(ti+1, tn, level);
Vector tmpall(outvec.length());
tmpall.zeros();
Vector tmp(outvec.length());
// note that since beg came from begin, it has empty signal
// and first evaluation gets no signal
for (_Tpit run = beg; run != end; ++run) {
func.eval(run.point(), run.signal(), tmp);
tmpall.add(run.weight(), tmp);
}
{
SYNCHRO@, syn(&outvec, "IntegrationWorker");
outvec.add(1.0, tmpall);
}
}
@ This is the class which implements the integration. The class is
templated by the iterator type. We declare a method |begin| returning
an iterator to the beginnning of the |ti|-th portion out of total |tn|
portions for a given level.
In addition, we define a method which saves all the points to a given
file. Only for debugging purposes.
@<|QuadratureImpl| class declaration@>=
template <typename _Tpit>
class QuadratureImpl : public Quadrature {
friend class IntegrationWorker<_Tpit>;
public:@;
QuadratureImpl(int d) : Quadrature(d)@+ {}
@<|QuadratureImpl::integrate| code@>;
void integrate(const VectorFunction& func,
int level, int tn, Vector& out) const {
VectorFunctionSet fs(func, tn);
integrate(fs, level, out);
}
@<|Quadrature::savePoints| code@>;
_Tpit start(int level) const
{@+ return begin(0,1,level);@+}
_Tpit end(int level) const
{@+ return begin(1,1,level);@+}
protected:@;
virtual _Tpit begin(int ti, int tn, int level) const =0;
};
@ Just fill a thread group with workes and run it.
@<|QuadratureImpl::integrate| code@>=
void integrate(VectorFunctionSet& fs, int level, Vector& out) const {
// todo: out.length()==func.outdim()
// todo: dim == func.indim()
out.zeros();
THREAD_GROUP@, gr;
for (int ti = 0; ti < fs.getNum(); ti++) {
gr.insert(new IntegrationWorker<_Tpit>(*this, fs.getFunc(ti),
level, ti, fs.getNum(), out));
}
gr.run();
}
@ Just for debugging.
@<|Quadrature::savePoints| code@>=
void savePoints(const char* fname, int level) const
{
FILE* fd;
if (NULL==(fd = fopen(fname,"w"))) {
// todo: raise
fprintf(stderr, "Cannot open file %s for writing.\n", fname);
exit(1);
}
_Tpit beg = begin(0,1,level);
_Tpit end = begin(1,1,level);
for (_Tpit run = beg; run != end; ++run) {
fprintf(fd, "%16.12g", run.weight());
for (int i = 0; i < dimen(); i++)
fprintf(fd, "\t%16.12g", run.point()[i]);
fprintf(fd, "\n");
}
fclose(fd);
}
@ This is only an interface to a precalculated data in file {\tt
precalc\_quadrature.dat} which is basically C coded static data. It
implements |OneDQuadrature|. The data file is supposed to define the
following data: number of levels, array of number of points at each
level, an array of weights and array of points. The both latter array
store data level by level. An offset for a specific level is stored in
|offsets| integer sequence.
The implementing subclasses just fill the necessary data from the
file, the rest is calculated here.
@<|OneDPrecalcQuadrature| class declaration@>=
class OneDPrecalcQuadrature : public OneDQuadrature {
int num_levels;
const int* num_points;
const double* weights;
const double* points;
IntSequence offsets;
public:@;
OneDPrecalcQuadrature(int nlevels, const int* npoints,
const double* wts, const double* pts)
: num_levels(nlevels), num_points(npoints),
weights(wts), points(pts), offsets(num_levels)
{@+ calcOffsets();@+}
virtual ~OneDPrecalcQuadrature()@+ {}
int numLevels() const
{@+ return num_levels;@+}
int numPoints(int level) const
{@+ return num_points[level-1];@+}
double point(int level, int i) const
{@+ return points[offsets[level-1]+i];@+}
double weight(int level, int i) const
{@+ return weights[offsets[level-1]+i];@+}
protected:@;
void calcOffsets();
};
@ Just precalculated Gauss--Hermite quadrature. This quadrature integrates exactly integrals
$$\int_{-\infty}^{\infty} x^ke^{-x^2}{\rm d}x$$
for level $k$.
Note that if pluging this one-dimensional quadrature to product or
Smolyak rule in order to integrate a function $f$ through normally
distributed inputs, one has to wrap $f$ to
|GaussConverterFunction| and apply the product or Smolyak rule to the
new function.
Check {\tt precalc\_quadrature.dat} for available levels.
@<|GaussHermite| class declaration@>=
class GaussHermite : public OneDPrecalcQuadrature {
public:@;
GaussHermite();
};
@ Just precalculated Gauss--Legendre quadrature. This quadrature integrates exactly integrals
$$\int_0^1x^k{\rm d}x$$
for level $k$.
Check {\tt precalc\_quadrature.dat} for available levels.
@<|GaussLegendre| class declaration@>=
class GaussLegendre : public OneDPrecalcQuadrature {
public:@;
GaussLegendre();
};
@ This is just an inverse cummulative density function of normal
distribution. Its only method |get| returns for a given number
$x\in(0,1)$ a number $y$ such that $P(z<y)=x$, where the probability
is taken over normal distribution $N(0,1)$.
Currently, the implementation is done by a table lookup which implies
that the tails had to be chopped off. This further implies that Monte
Carlo quadratures using this transformation to draw points from
multidimensional $N(0,I)$ fail to integrate with satisfactory
precision polynomial functions, for which the tails matter.
@<|NormalICDF| class declaration@>=
class NormalICDF {
public:@;
static double get(double x);
};
@ End of {\tt quadrature.h} file
dynare++/integ/cc/quasi_mcarlo.cweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: quasi_mcarlo.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@ This is {\tt quasi\_mcarlo.cpp} file.
@c
#include "quasi_mcarlo.h"
#include <cmath>
@<|RadicalInverse| constructor code@>;
@<|RadicalInverse::eval| code@>;
@<|RadicalInverse::increase| code@>;
@<|RadicalInverse::print| code@>;
@<|HaltonSequence| static data@>;
@<|HaltonSequence| constructor code@>;
@<|HaltonSequence::operator=| code@>;
@<|HaltonSequence::increase| code@>;
@<|HaltonSequence::eval| code@>;
@<|HaltonSequence::print| code@>;
@<|qmcpit| empty constructor code@>;
@<|qmcpit| regular constructor code@>;
@<|qmcpit| copy constructor code@>;
@<|qmcpit| destructor@>;
@<|qmcpit::operator==| code@>;
@<|qmcpit::operator=| code@>;
@<|qmcpit::operator++| code@>;
@<|qmcpit::weight| code@>;
@<|qmcnpit| empty constructor code@>;
@<|qmcnpit| regular constructor code@>;
@<|qmcnpit| copy constructor code@>;
@<|qmcnpit| destructor@>;
@<|qmcnpit::operator=| code@>;
@<|qmcnpit::operator++| code@>;
@<|WarnockPerScheme::permute| code@>;
@<|ReversePerScheme::permute| code@>;
@ Here in the constructor, we have to calculate a maximum length of
|coeff| array for a given |base| and given maximum |maxn|. After
allocation, we calculate the coefficients.
@<|RadicalInverse| constructor code@>=
RadicalInverse::RadicalInverse(int n, int b, int mxn)
: num(n), base(b), maxn(mxn),
coeff((int)(floor(log((double)maxn)/log((double)b))+2), 0)
{
int nr = num;
j = -1;
do {
j++;
coeff[j] = nr % base;
nr = nr / base;
} while (nr > 0);
}
@ This evaluates the radical inverse. If there was no permutation, we have to calculate
$$
{c_0\over b}+{c_1\over b^2}+\ldots+{c_j\over b^{j+1}}
$$
which is evaluated as
$$
\left(\ldots\left(\left({c_j\over b}\cdot{1\over b}+{c_{j-1}\over b}\right)
\cdot{1\over b}+{c_{j-2}\over b}\right)
\ldots\right)\cdot{1\over b}+{c_0\over b}
$$
Now with permutation $\pi$, we have
$$
\left(\ldots\left(\left({\pi(c_j)\over b}\cdot{1\over b}+
{\pi(c_{j-1})\over b}\right)\cdot{1\over b}+
{\pi(c_{j-2})\over b}\right)
\ldots\right)\cdot{1\over b}+{\pi(c_0)\over b}
$$
@<|RadicalInverse::eval| code@>=
double RadicalInverse::eval(const PermutationScheme& p) const
{
double res = 0;
for (int i = j; i >= 0; i--) {
int cper = p.permute(i, base, coeff[i]);
res = (cper + res)/base;
}
return res;
}
@ We just add 1 to the lowest coefficient and check for overflow with respect to the base.
@<|RadicalInverse::increase| code@>=
void RadicalInverse::increase()
{
// todo: raise if |num+1 > maxn|
num++;
int i = 0;
coeff[i]++;
while (coeff[i] == base) {
coeff[i] = 0;
coeff[++i]++;
}
if (i > j)
j = i;
}
@ Debug print.
@<|RadicalInverse::print| code@>=
void RadicalInverse::print() const
{
printf("n=%d b=%d c=", num, base);
coeff.print();
}
@ Here we have the first 170 primes. This means that we are not able
to integrate dimensions greater than 170.
@<|HaltonSequence| static data@>=
int HaltonSequence::num_primes = 170;
int HaltonSequence::primes[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
};
@ This takes first |dim| primes and constructs |dim| radical inverses
and calls |eval|.
@<|HaltonSequence| constructor code@>=
HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme& p)
: num(n), maxn(mxn), per(p), pt(dim)
{
// todo: raise if |dim > num_primes|
// todo: raise if |n > mxn|
for (int i = 0; i < dim; i++)
ri.push_back(RadicalInverse(num, primes[i], maxn));
eval();
}
@
@<|HaltonSequence::operator=| code@>=
const HaltonSequence& HaltonSequence::operator=(const HaltonSequence& hs)
{
num = hs.num;
maxn = hs.maxn;
ri.clear();
for (unsigned int i = 0; i < hs.ri.size(); i++)
ri.push_back(RadicalInverse(hs.ri[i]));
pt = hs.pt;
return *this;
}
@ This calls |RadicalInverse::increase| for all radical inverses and
calls |eval|.
@<|HaltonSequence::increase| code@>=
void HaltonSequence::increase()
{
for (unsigned int i = 0; i < ri.size(); i++)
ri[i].increase();
num++;
if (num <= maxn)
eval();
}
@ This sets point |pt| to radical inverse evaluations in each dimension.
@<|HaltonSequence::eval| code@>=
void HaltonSequence::eval()
{
for (unsigned int i = 0; i < ri.size(); i++)
pt[i] = ri[i].eval(per);
}
@ Debug print.
@<|HaltonSequence::print| code@>=
void HaltonSequence::print() const
{
for (unsigned int i = 0; i < ri.size(); i++)
ri[i].print();
printf("point=[ ");
for (unsigned int i = 0; i < ri.size(); i++)
printf("%7.6f ", pt[i]);
printf("]\n");
}
@
@<|qmcpit| empty constructor code@>=
qmcpit::qmcpit()
: spec(NULL), halton(NULL), sig(NULL)@+ {}
@
@<|qmcpit| regular constructor code@>=
qmcpit::qmcpit(const QMCSpecification& s, int n)
: spec(&s), halton(new HaltonSequence(n, s.level(), s.dimen(), s.getPerScheme())),
sig(new ParameterSignal(s.dimen()))
{
}
@
@<|qmcpit| copy constructor code@>=
qmcpit::qmcpit(const qmcpit& qpit)
: spec(qpit.spec), halton(NULL), sig(NULL)
{
if (qpit.halton)
halton = new HaltonSequence(*(qpit.halton));
if (qpit.sig)
sig = new ParameterSignal(qpit.spec->dimen());
}
@
@<|qmcpit| destructor@>=
qmcpit::~qmcpit()
{
if (halton)
delete halton;
if (sig)
delete sig;
}
@
@<|qmcpit::operator==| code@>=
bool qmcpit::operator==(const qmcpit& qpit) const
{
return (spec == qpit.spec) &&
((halton == NULL && qpit.halton == NULL) ||
(halton != NULL && qpit.halton != NULL && halton->getNum() == qpit.halton->getNum()));
}
@
@<|qmcpit::operator=| code@>=
const qmcpit& qmcpit::operator=(const qmcpit& qpit)
{
spec = qpit.spec;
if (halton)
delete halton;
if (qpit.halton)
halton = new HaltonSequence(*(qpit.halton));
else
halton = NULL;
return *this;
}
@
@<|qmcpit::operator++| code@>=
qmcpit& qmcpit::operator++()
{
// todo: raise if |halton == null || qmcq == NULL|
halton->increase();
return *this;
}
@
@<|qmcpit::weight| code@>=
double qmcpit::weight() const
{
return 1.0/spec->level();
}
@
@<|qmcnpit| empty constructor code@>=
qmcnpit::qmcnpit()
: qmcpit(), pnt(NULL)@+ {}
@
@<|qmcnpit| regular constructor code@>=
qmcnpit::qmcnpit(const QMCSpecification& s, int n)
: qmcpit(s, n), pnt(new Vector(s.dimen()))
{
}
@
@<|qmcnpit| copy constructor code@>=
qmcnpit::qmcnpit(const qmcnpit& qpit)
: qmcpit(qpit), pnt(NULL)
{
if (qpit.pnt)
pnt = new Vector(*(qpit.pnt));
}
@
@<|qmcnpit| destructor@>=
qmcnpit::~qmcnpit()
{
if (pnt)
delete pnt;
}
@
@<|qmcnpit::operator=| code@>=
const qmcnpit& qmcnpit::operator=(const qmcnpit& qpit)
{
qmcpit::operator=(qpit);
if (pnt)
delete pnt;
if (qpit.pnt)
pnt = new Vector(*(qpit.pnt));
else
pnt = NULL;
return *this;
}
@ Here we inccrease a point in Halton sequence ant then store images
of the points in |NormalICDF| function.
@<|qmcnpit::operator++| code@>=
qmcnpit& qmcnpit::operator++()
{
qmcpit::operator++();
for (int i = 0; i < halton->point().length(); i++)
(*pnt)[i] = NormalICDF::get(halton->point()[i]);
return *this;
}
@ Clear from code.
@<|WarnockPerScheme::permute| code@>=
int WarnockPerScheme::permute(int i, int base, int c) const
{
return (c+i) % base;
}
@ Clear from code.
@<|ReversePerScheme::permute| code@>=
int ReversePerScheme::permute(int i, int base, int c) const
{
return (base-c) % base;
}
@ End of {\tt quasi\_mcarlo.cpp} file.
\ No newline at end of file
dynare++/integ/cc/quasi_mcarlo.hweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: quasi_mcarlo.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@*2 Quasi Monte Carlo quadrature. This is {\tt quasi\_mcarlo.h} file.
This defines quasi Monte Carlo quadratures for cube and for a function
multiplied by normal density. The quadrature for a cube is named
|QMCarloCubeQuadrature| and integrates:
$$\int_{\langle 0,1\rangle^n}f(x){\rm d}x$$
The quadrature for a function of normally distributed parameters is
named |QMCarloNormalQuadrature| and integrates:
$${1\over\sqrt{(2\pi)^n}}\int_{(-\infty,\infty)^n}f(x)e^{-{1\over 2}x^Tx}{\rm d}x$$
For a cube we define |qmcpit| as iterator of |QMCarloCubeQuadrature|,
and for the normal density multiplied function we define |qmcnpit| as
iterator of |QMCarloNormalQuadrature|.
The quasi Monte Carlo method generates low discrepancy points with
equal weights. The one dimensional low discrepancy sequences are
generated by |RadicalInverse| class, the sequences are combined for
higher dimensions by |HaltonSequence| class. The Halton sequence can
use a permutation scheme; |PermutattionScheme| is an abstract class
for all permutaton schemes. We have three implementations:
|WarnockPerScheme|, |ReversePerScheme|, and |IdentityPerScheme|.
@s PermutationScheme int
@s RadicalInverse int
@s HaltonSequence int
@s QMCSpecification int
@s qmcpit int
@s QMCarloCubeQuadrature int
@s qmcnpit int
@s QMCarloNormalQuadrature int
@s WarnockPerScheme int
@s ReversePerScheme int
@s IdentityPerScheme int
@c
#ifndef QUASI_MCARLO_H
#define QUASI_MCARLO_H
#include "int_sequence.h"
#include "quadrature.h"
#include "Vector.h"
#include <vector>
@<|PermutationScheme| class declaration@>;
@<|RadicalInverse| class declaration@>;
@<|HaltonSequence| class declaration@>;
@<|QMCSpecification| class declaration@>;
@<|qmcpit| class declaration@>;
@<|QMCarloCubeQuadrature| class declaration@>;
@<|qmcnpit| class declaration@>;
@<|QMCarloNormalQuadrature| class declaration@>;
@<|WarnockPerScheme| class declaration@>;
@<|ReversePerScheme| class declaration@>;
@<|IdentityPerScheme| class declaration@>;
#endif
@ This abstract class declares |permute| method which permutes
coefficient |c| having index of |i| fro the base |base| and returns
the permuted coefficient which must be in $0,\ldots,base-1$.
@<|PermutationScheme| class declaration@>=
class PermutationScheme {
public:@;
PermutationScheme()@+ {}
virtual ~PermutationScheme()@+ {}
virtual int permute(int i, int base, int c) const =0;
};
@ This class represents an integer number |num| as
$c_0+c_1b+c_2b^2+\ldots+c_jb^j$, where $b$ is |base| and
$c_0,\ldots,c_j$ is stored in |coeff|. The size of |IntSequence| coeff
does not grow with growing |num|, but is fixed from the very beginning
and is set according to supplied maximum |maxn|.
The basic method is |eval| which evaluates the |RadicalInverse| with a
given permutation scheme and returns the point, and |increase| which
increases |num| and recalculates the coefficients.
@<|RadicalInverse| class declaration@>=
class RadicalInverse {
int num;
int base;
int maxn;
int j;
IntSequence coeff;
public:@;
RadicalInverse(int n, int b, int mxn);
RadicalInverse(const RadicalInverse& ri)
: num(ri.num), base(ri.base), maxn(ri.maxn), j(ri.j), coeff(ri.coeff)@+ {}
const RadicalInverse& operator=(const RadicalInverse& radi)
{
num = radi.num; base = radi.base; maxn = radi.maxn;
j = radi.j; coeff = radi.coeff;
return *this;
}
double eval(const PermutationScheme& p) const;
void increase();
void print() const;
};
@ This is a vector of |RadicalInverse|s, each |RadicalInverse| has a
different prime as its base. The static members |primes| and
|num_primes| define a precalculated array of primes. The |increase|
method of the class increases indices in all |RadicalInverse|s and
sets point |pt| to contain the points in each dimension.
@<|HaltonSequence| class declaration@>=
class HaltonSequence {
private:@;
static int primes[];
static int num_primes;
protected:@;
int num;
int maxn;
vector<RadicalInverse> ri;
const PermutationScheme& per;
Vector pt;
public:@;
HaltonSequence(int n, int mxn, int dim, const PermutationScheme& p);
HaltonSequence(const HaltonSequence& hs)
: num(hs.num), maxn(hs.maxn), ri(hs.ri), per(hs.per), pt(hs.pt)@+ {}
const HaltonSequence& operator=(const HaltonSequence& hs);
void increase();
const Vector& point() const
{@+ return pt;@+}
const int getNum() const
{@+ return num;@+}
void print() const;
protected:@;
void eval();
};
@ This is a specification of quasi Monte Carlo quadrature. It consists
of dimension |dim|, number of points (or level) |lev|, and the
permutation scheme. This class is common to all quasi Monte Carlo
classes.
@<|QMCSpecification| class declaration@>=
class QMCSpecification {
protected:@;
int dim;
int lev;
const PermutationScheme& per_scheme;
public:@;
QMCSpecification(int d, int l, const PermutationScheme& p)
: dim(d), lev(l), per_scheme(p)@+ {}
virtual ~QMCSpecification() {}
int dimen() const
{@+ return dim;@+}
int level() const
{@+ return lev;@+}
const PermutationScheme& getPerScheme() const
{@+ return per_scheme;@+}
};
@ This is an iterator for quasi Monte Carlo over a cube
|QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of
the same dimension as given by the specification. An iterator can be
constructed from a given number |n|, or by a copy constructor. For
technical reasons, there is also an empty constructor; for that
reason, every member is a pointer.
@<|qmcpit| class declaration@>=
class qmcpit {
protected:@;
const QMCSpecification* spec;
HaltonSequence* halton;
ParameterSignal* sig;
public:@;
qmcpit();
qmcpit(const QMCSpecification& s, int n);
qmcpit(const qmcpit& qpit);
~qmcpit();
bool operator==(const qmcpit& qpit) const;
bool operator!=(const qmcpit& qpit) const
{@+ return ! operator==(qpit);@+}
const qmcpit& operator=(const qmcpit& qpit);
qmcpit& operator++();
const ParameterSignal& signal() const
{@+ return *sig;@+}
const Vector& point() const
{@+ return halton->point();@+}
double weight() const;
void print() const
{@+ halton->print();@+}
};
@ This is an easy declaration of quasi Monte Carlo quadrature for a
cube. Everything important has been done in its iterator |qmcpit|, so
we only inherit from general |Quadrature| and reimplement |begin| and
|numEvals|.
@<|QMCarloCubeQuadrature| class declaration@>=
class QMCarloCubeQuadrature : public QuadratureImpl<qmcpit>, public QMCSpecification {
public:@;
QMCarloCubeQuadrature(int d, int l, const PermutationScheme& p)
: QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)@+ {}
virtual ~QMCarloCubeQuadrature()@+ {}
int numEvals(int l) const
{@+ return l;@+}
protected:@;
qmcpit begin(int ti, int tn, int lev) const
{@+ return qmcpit(*this, ti*level()/tn + 1);@+}
};
@ This is an iterator for |QMCarloNormalQuadrature|. It is equivalent
to an iterator for quasi Monte Carlo cube quadrature but a point. The
point is obtained from a point of |QMCarloCubeQuadrature| by a
transformation $\langle
0,1\rangle\rightarrow\langle-\infty,\infty\rangle$ aplied to all
dimensions. The transformation yields a normal distribution from a
uniform distribution on $\langle0,1\rangle$. It is in fact
|NormalICDF|.
@<|qmcnpit| class declaration@>=
class qmcnpit : public qmcpit {
protected:@;
Vector* pnt;
public:@;
qmcnpit();
qmcnpit(const QMCSpecification& spec, int n);
qmcnpit(const qmcnpit& qpit);
~qmcnpit();
bool operator==(const qmcnpit& qpit) const
{@+ return qmcpit::operator==(qpit);@+}
bool operator!=(const qmcnpit& qpit) const
{@+ return ! operator==(qpit);@+}
const qmcnpit& operator=(const qmcnpit& qpit);
qmcnpit& operator++();
const ParameterSignal& signal() const
{@+ return *sig;@+}
const Vector& point() const
{@+ return *pnt;@+}
void print() const
{@+ halton->print();pnt->print();@+}
};
@ This is an easy declaration of quasi Monte Carlo quadrature for a
cube. Everything important has been done in its iterator |qmcnpit|, so
we only inherit from general |Quadrature| and reimplement |begin| and
|numEvals|.
@<|QMCarloNormalQuadrature| class declaration@>=
class QMCarloNormalQuadrature : public QuadratureImpl<qmcnpit>, public QMCSpecification {
public:@;
QMCarloNormalQuadrature(int d, int l, const PermutationScheme& p)
: QuadratureImpl<qmcnpit>(d), QMCSpecification(d, l, p)@+ {}
virtual ~QMCarloNormalQuadrature()@+ {}
int numEvals(int l) const
{@+ return l;@+}
protected:@;
qmcnpit begin(int ti, int tn, int lev) const
{@+ return qmcnpit(*this, ti*level()/tn + 1);@+}
};
@ Declares Warnock permutation scheme.
@<|WarnockPerScheme| class declaration@>=
class WarnockPerScheme : public PermutationScheme {
public:@;
int permute(int i, int base, int c) const;
};
@ Declares reverse permutation scheme.
@<|ReversePerScheme| class declaration@>=
class ReversePerScheme : public PermutationScheme {
public:@;
int permute(int i, int base, int c) const;
};
@ Declares no permutation (identity) scheme.
@<|IdentityPerScheme| class declaration@>=
class IdentityPerScheme : public PermutationScheme {
public:@;
int permute(int i, int base, int c) const
{@+ return c;@+}
};
@ End of {\tt quasi\_mcarlo.h} file
dynare++/integ/cc/smolyak.cweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: smolyak.cweb 1208 2007-03-19 21:33:12Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@ This is {\tt smolyak.cpp} file.
@c
#include "smolyak.h"
#include "symmetry.h"
@<|smolpit| empty constructor@>;
@<|smolpit| regular constructor@>;
@<|smolpit| copy constructor@>;
@<|smolpit| destructor@>;
@<|smolpit::operator==| code@>;
@<|smolpit::operator=| code@>;
@<|smolpit::operator++| code@>;
@<|smolpit::setPointAndWeight| code@>;
@<|smolpit::print| code@>;
@<|SmolyakQuadrature| constructor@>;
@<|SmolyakQuadrature::numEvals| code@>;
@<|SmolyakQuadrature::begin| code@>;
@<|SmolyakQuadrature::calcNumEvaluations| code@>;
@<|SmolyakQuadrature::designLevelForEvals| code@>;
@
@<|smolpit| empty constructor@>=
smolpit::smolpit()
: smolq(NULL), isummand(0), jseq(NULL), sig(NULL), p(NULL)
{
}
@ This constructs a beginning of |isum| summand in |smolq|. We must be
careful here, since |isum| can be past-the-end, so no reference to
vectors in |smolq| by |isum| must be done in this case.
@<|smolpit| regular constructor@>=
smolpit::smolpit(const SmolyakQuadrature& q, unsigned int isum)
: smolq(&q), isummand(isum), jseq(new IntSequence(q.dimen(), 0)),
sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
{
if (isummand < q.numSummands()) {
setPointAndWeight();
}
}
@
@<|smolpit| copy constructor@>=
smolpit::smolpit(const smolpit& spit)
: smolq(spit.smolq), isummand(spit.isummand), w(spit.w)
{
if (spit.jseq)
jseq = new IntSequence(*(spit.jseq));
else
jseq = NULL;
if (spit.sig)
sig = new ParameterSignal(*(spit.sig));
else
sig = NULL;
if (spit.p)
p = new Vector(*(spit.p));
else
p = NULL;
}
@
@<|smolpit| destructor@>=
smolpit::~smolpit()
{
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
}
@
@<|smolpit::operator==| code@>=
bool smolpit::operator==(const smolpit& spit) const
{
bool ret = true;
ret = ret & smolq == spit.smolq;
ret = ret & isummand == spit.isummand;
ret = ret & ((jseq==NULL && spit.jseq==NULL) ||
(jseq!=NULL && spit.jseq!=NULL && *jseq == *(spit.jseq)));
return ret;
}
@
@<|smolpit::operator=| code@>=
const smolpit& smolpit::operator=(const smolpit& spit)
{
smolq = spit.smolq;
isummand = spit.isummand;
w = spit.w;
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
if (spit.jseq)
jseq = new IntSequence(*(spit.jseq));
else
jseq = NULL;
if (spit.sig)
sig = new ParameterSignal(*(spit.sig));
else
sig = NULL;
if (spit.p)
p = new Vector(*(spit.p));
else
p = NULL;
return *this;
}
@ We first try to increase index within the current summand. If we are
at maximum, we go to a subsequent summand. Note that in this case all
indices in |jseq| will be zero, so no change is needed.
@<|smolpit::operator++| code@>=
smolpit& smolpit::operator++()
{
// todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL|
const IntSequence& levpts = smolq->levpoints[isummand];
int i = smolq->dimen()-1;
(*jseq)[i]++;
while (i >= 0 && (*jseq)[i] == levpts[i]) {
(*jseq)[i] = 0;
i--;
if (i >= 0)
(*jseq)[i]++;
}
sig->signalAfter(std::max(i,0));
if (i < 0)
isummand++;
if (isummand < smolq->numSummands())
setPointAndWeight();
return *this;
}
@ Here we set the point coordinates according to |jseq| and
|isummand|. Also the weight is set here.
@<|smolpit::setPointAndWeight| code@>=
void smolpit::setPointAndWeight()
{
// todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or
// |p==NULL| or |isummand>=smolq->numSummands()|
int l = smolq->level;
int d = smolq->dimen();
int sumk = (smolq->levels[isummand]).sum();
int m1exp = l + d - sumk - 1;
w = (2*(m1exp/2)==m1exp)? 1.0 : -1.0;
w *= smolq->psc.noverk(d-1, sumk-l);
for (int i = 0; i < d; i++) {
int ki = (smolq->levels[isummand])[i];
(*p)[i] = (smolq->uquad).point(ki, (*jseq)[i]);
w *= (smolq->uquad).weight(ki, (*jseq)[i]);
}
}
@ Debug print.
@<|smolpit::print| code@>=
void smolpit::print() const
{
printf("isum=%-3d: [", isummand);
for (int i = 0; i < smolq->dimen(); i++)
printf("%2d ", (smolq->levels[isummand])[i]);
printf("] j=[");
for (int i = 0; i < smolq->dimen(); i++)
printf("%2d ", (*jseq)[i]);
printf("] %+4.3f*(",w);
for (int i = 0; i < smolq->dimen()-1; i++)
printf("%+4.3f ", (*p)[i]);
printf("%+4.3f)\n",(*p)[smolq->dimen()-1]);
}
@ Here is the constructor of |SmolyakQuadrature|. We have to setup
|levels|, |levpoints| and |cumevals|. We have to go through all
$d$-dimensional sequences $k$, such that $l\leq \vert k\vert\leq
l+d-1$ and all $k_i$ are positive integers. This is equivalent to
going through all $k$ such that $l-d\leq\vert k\vert\leq l-1$ and all
$k_i$ are non-negative integers. This is equivalent to going through
$d+1$ dimensional sequences $(k,x)$ such that $\vert(k,x)\vert =l-1$
and $x=0,\ldots,d-1$. The resulting sequence of positive integers is
obtained by adding $1$ to all $k_i$.
@<|SmolyakQuadrature| constructor@>=
SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature& uq)
: QuadratureImpl<smolpit>(d), level(l), uquad(uq), psc(d-1,d-1)
{
// todo: check |l>1|, |l>=d|
// todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
int cum = 0;
SymmetrySet ss(l-1, d+1);
for (symiterator si(ss); !si.isEnd(); ++si) {
if ((*si)[d] <= d-1) {
IntSequence lev((const IntSequence&)*si, 0, d);
lev.add(1);
levels.push_back(lev);
IntSequence levpts(d);
for (int i = 0; i < d; i++)
levpts[i] = uquad.numPoints(lev[i]);
levpoints.push_back(levpts);
cum += levpts.mult();
cumevals.push_back(cum);
}
}
}
@ Here we return a number of evalutions of the quadrature for the
given level. If the given level is the current one, we simply return
the maximum cumulative number of evaluations. Otherwise we call costly
|calcNumEvaluations| method.
@<|SmolyakQuadrature::numEvals| code@>=
int SmolyakQuadrature::numEvals(int l) const
{
if (l != level)
return calcNumEvaluations(l);
else
return cumevals[numSummands()-1];
}
@ This divides all the evaluations to |tn| approximately equal groups,
and returns the beginning of the specified group |ti|. The granularity
of divisions are summands as listed by |levels|.
@<|SmolyakQuadrature::begin| code@>=
smolpit SmolyakQuadrature::begin(int ti, int tn, int l) const
{
// todo: raise is |level!=l|
if (ti == tn)
return smolpit(*this, numSummands());
int totevals = cumevals[numSummands()-1];
int evals = (totevals*ti)/tn;
unsigned int isum = 0;
while (isum+1 < numSummands() && cumevals[isum+1] < evals)
isum++;
return smolpit(*this, isum);
}
@ This is the same in a structure as |@<|SmolyakQuadrature| constructor@>|.
We have to go through all summands and calculate
a number of evaluations in each summand.
@<|SmolyakQuadrature::calcNumEvaluations| code@>=
int SmolyakQuadrature::calcNumEvaluations(int lev) const
{
int cum = 0;
SymmetrySet ss(lev-1, dim+1);
for (symiterator si(ss); !si.isEnd(); ++si) {
if ((*si)[dim] <= dim-1) {
IntSequence lev((const IntSequence&)*si, 0, dim);
lev.add(1);
IntSequence levpts(dim);
for (int i = 0; i < dim; i++)
levpts[i] = uquad.numPoints(lev[i]);
cum += levpts.mult();
}
}
return cum;
}
@ This returns a maximum level such that the number of evaluations is
less than the given number.
@<|SmolyakQuadrature::designLevelForEvals| code@>=
void SmolyakQuadrature::designLevelForEvals(int max_evals, int& lev, int& evals) const
{
int last_evals;
evals = 1;
lev = 1;
do {
lev++;
last_evals = evals;
evals = calcNumEvaluations(lev);
} while (lev < uquad.numLevels() && evals <= max_evals);
lev--;
evals = last_evals;
}
@ End of {\tt smolyak.cpp} file
dynare++/integ/cc/smolyak.hweb deleted 100644 → 0
View file @ 6ca649c4
@q $Id: smolyak.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@*2 Smolyak quadrature. This is {\tt smolyak.h} file
This file defines Smolyak (sparse grid) multidimensional quadrature
for non-nested underlying one dimensional quadrature. Let $Q^1_l$ denote
the one dimensional quadrature of $l$ level. Let $n_l$ denote a
number of points in the $l$ level. Than the Smolyak quadrature can be
defined as
$$Q^df=\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert k\vert-1}\left(\matrix{d-1\cr
\vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f,$$
where $d$ is the dimension, $k$ is $d$-dimensional sequence of
integers, and $\vert k\vert$ denotes a sum of the sequence.
Here we define |smolpit| as Smolyak iterator and |SmolyakQuadrature|.
@s smolpit int
@s SmolyakQuadrature int
@s PascalTriangle int
@s SymmetrySet int
@s symiterator int
@c
#ifndef SMOLYAK_H
#define SMOLYAK_H
#include "int_sequence.h"
#include "tl_static.h"
#include "vector_function.h"
#include "quadrature.h"
@<|smolpit| class declaration@>;
@<|SmolyakQuadrature| class declaration@>;
#endif
@ Here we define the Smolyak point iterator. The Smolyak formula can
be broken to a sum of product quadratures with various combinations of
levels. The iterator follows this pattern. It maintains an index to a
summand and then a point coordinates within the summand (product
quadrature). The array of summands to which the |isummand| points is
maintained by the |SmolyakQuadrature| class to which the object knows
the pointer |smolq|.
We provide a constructor which points to the beginning of the given
summand. This constructor is used in |SmolyakQuadrature::begin| method
which approximately divideds all the iterators to subsets of equal
size.
@<|smolpit| class declaration@>=
class SmolyakQuadrature;
class smolpit {
protected:@;
const SmolyakQuadrature* smolq;
unsigned int isummand;
IntSequence* jseq;
ParameterSignal* sig;
Vector* p;
double w;
public:@;
smolpit();
smolpit(const SmolyakQuadrature& q, unsigned int isum);
smolpit(const smolpit& spit);
~smolpit();
bool operator==(const smolpit& spit) const;
bool operator!=(const smolpit& spit) const
{@+ return ! operator==(spit);@+}
const smolpit& operator=(const smolpit& spit);
smolpit& operator++();
const ParameterSignal& signal() const
{@+ return *sig;@+}
const Vector& point() const
{@+ return *p;@+}
double weight() const
{@+ return w;@+}
void print() const;
protected:@;
void setPointAndWeight();
};
@ Here we define the class |SmolyakQuadrature|. It maintains an array
of summands of the Smolyak quadrature formula:
$$\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert
k\vert-1}\left(\matrix{d-1\cr
\vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f$$
Each summand is fully specified by sequence $k$. The summands are here
represented (besides $k$) also by sequence of number of points in each
level selected by $k$, and also by a cummulative number of
evaluations. The latter two are added only for conveniency.
The summands in the code are given by |levels|, which is a vector of
$k$ sequences, further by |levpoints| which is a vector of sequences
of nuber of points in each level, and by |cumevals| which is the
cumulative number of points, this is $\sum_k\prod_{i=1}^dn_{k_i}$,
where the sum is done through all $k$ before the current.
The |levels| and |levpoints| vectors are used by |smolpit|.
@<|SmolyakQuadrature| class declaration@>=
class SmolyakQuadrature : public QuadratureImpl<smolpit> {
friend class smolpit;
int level;
const OneDQuadrature& uquad;
vector<IntSequence> levels;
vector<IntSequence> levpoints;
vector<int> cumevals;
PascalTriangle psc;
public:@;
SmolyakQuadrature(int d, int l, const OneDQuadrature& uq);
virtual ~SmolyakQuadrature()@+ {}
virtual int numEvals(int level) const;
void designLevelForEvals(int max_eval, int& lev, int& evals) const;
protected:@;
smolpit begin(int ti, int tn, int level) const;
unsigned int numSummands() const
{@+ return levels.size();@+}
private:@;
int calcNumEvaluations(int level) const;
};
@ End of {\tt smolyak.h} file