diff --git a/README.md b/README.md
index 36ea89126..4c59c93c0 100644
--- a/README.md
+++ b/README.md
@@ -78,9 +78,9 @@ A number of tools and libraries are needed in order to recompile everything. You
- [Flex](https://github.com/westes/flex), version 2.5.4 or later (only if you get the source through Git)
- [Autoconf](https://www.gnu.org/software/autoconf/), version 2.62 or later (only if you get the source through Git)
- [Automake](https://www.gnu.org/software/automake/), version 1.11.2 or later (only if you get the source through Git)
-- [MAT File I/O library](https://sourceforge.net/projects/matio/), version 1.5 or later (if you want to compile Markov-Switching code, the estimation DLL, and the k-order DLL)
+- [MAT File I/O library](https://sourceforge.net/projects/matio/), version 1.5 or later (if you want to compile Markov-Switching SBVAR code)
- [SLICOT](http://www.slicot.org) (if you want to compile the Kalman steady state DLL)
-- [GSL library](https://www.gnu.org/software/gsl/) (if you want to compile Markov-Switching code)
+- [GSL library](https://www.gnu.org/software/gsl/) (if you want to compile Markov-Switching SBVAR code)
- A decent LaTeX distribution (if you want to compile PDF documentation),
ideally with Beamer
- For building the reference manual:
diff --git a/license.txt b/license.txt
index 34140b637..06f5bd8c5 100644
--- a/license.txt
+++ b/license.txt
@@ -302,9 +302,8 @@ Copyright: 2010 Yannick Kalantzis
License: GPL-3+
Files: mex/sources/gensylv/gensylv.cc
- mex/sources/libkorder/kord/* mex/sources/libkorder/integ/*
- mex/sources/libkorder/sylv/* mex/sources/libkorder/tl/*
- mex/sources/libkorder/utils/*
+ mex/sources/libkorder/kord/* mex/sources/libkorder/sylv/*
+ mex/sources/libkorder/tl/* mex/sources/libkorder/utils/*
Copyright: 2004-2011 Ondra Kamenik
2019-2023 Dynare Team
License: GPL-3+
diff --git a/mex/build/libkorder.am b/mex/build/libkorder.am
index 954ce0226..6c7d573fc 100644
--- a/mex/build/libkorder.am
+++ b/mex/build/libkorder.am
@@ -2,7 +2,7 @@ noinst_LIBRARIES = libkorder.a
TOPDIR = $(top_srcdir)/../../sources/libkorder
-libkorder_a_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR) -I$(TOPDIR)/tl -I$(TOPDIR)/sylv -I$(TOPDIR)/integ -I$(TOPDIR)/kord -I$(TOPDIR)/utils $(CPPFLAGS_MATIO)
+libkorder_a_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR) -I$(TOPDIR)/tl -I$(TOPDIR)/sylv -I$(TOPDIR)/kord -I$(TOPDIR)/utils $(CPPFLAGS_MATIO)
# TODO: fix the codebase so that the following line is no longer needed
libkorder_a_CXXFLAGS = $(AM_CXXFLAGS) -Wno-unused-parameter
@@ -13,12 +13,9 @@ KORD_SRCS = \
dynamic_model.cc \
faa_di_bruno.cc \
first_order.cc \
- global_check.cc \
korder.cc \
korder_stoch.cc \
- journal.cc \
- normal_conjugate.cc \
- seed_generator.cc
+ journal.cc
SYLV_SRCS = \
BlockDiagonal.cc \
@@ -61,18 +58,10 @@ TL_SRCS = \
tl_static.cc \
twod_matrix.cc
-INTEG_SRCS = \
- quadrature.cc \
- quasi_mcarlo.cc \
- product.cc \
- smolyak.cc \
- vector_function.cc
-
UTILS_SRCS = \
pascal_triangle.cc \
int_power.cc \
- sthread.cc \
- nlsolve.cc
+ sthread.cc
TOPDIR_SRCS = \
k_ord_dynare.cc \
@@ -83,7 +72,6 @@ nodist_libkorder_a_SOURCES = \
$(KORD_SRCS) \
$(TL_SRCS) \
$(SYLV_SRCS) \
- $(INTEG_SRCS) \
$(UTILS_SRCS) \
$(TOPDIR_SRCS)
@@ -98,7 +86,5 @@ $(TL_SRCS): %.cc: $(TOPDIR)/tl/%.cc
$(LN_S) -f $< $@
$(SYLV_SRCS): %.cc: $(TOPDIR)/sylv/%.cc
$(LN_S) -f $< $@
-$(INTEG_SRCS): %.cc: $(TOPDIR)/integ/%.cc
- $(LN_S) -f $< $@
$(UTILS_SRCS): %.cc: $(TOPDIR)/utils/%.cc
$(LN_S) -f $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index d822bc07b..78134214b 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -1,11 +1,6 @@
ACLOCAL_AMFLAGS = -I ../../../m4
-SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update
-
-# libkorder must come first
-if ENABLE_MEX_GENSYLV_KORDER
-SUBDIRS += libkorder gensylv k_order_perturbation k_order_welfare
-endif
+SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update libkorder gensylv k_order_perturbation k_order_welfare
if ENABLE_MEX_MS_SBVAR
SUBDIRS += ms_sbvar
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index db678163c..54cf3c717 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -64,9 +64,6 @@ AC_SUBST([LIBADD_DLOPEN])
AC_ARG_ENABLE([mex-ms-sbvar], AS_HELP_STRING([--disable-mex-ms-sbvar], [disable compilation of the MS-SBVAR MEX]), [], [enable_mex_ms_sbvar=yes])
AM_CONDITIONAL([ENABLE_MEX_MS_SBVAR], [test "$enable_mex_ms_sbvar" = yes])
-AC_ARG_ENABLE([mex-gensylv-korder], AS_HELP_STRING([--disable-mex-gensylv-korder], [disable compilation of gensylv and k-order MEX]), [], [enable_mex_gensylv_korder=yes])
-AM_CONDITIONAL([ENABLE_MEX_GENSYLV_KORDER], [test "$enable_mex_gensylv_korder" = yes])
-
AC_ARG_ENABLE([mex-kalman-steady-state], AS_HELP_STRING([--disable-mex-kalman-steady-state], [disable compilation of the kalman_steady_state MEX]), [], [enable_mex_kalman_steady_state=yes])
AM_CONDITIONAL([ENABLE_MEX_KALMAN_STEADY_STATE], [test "$enable_mex_kalman_steady_state" = yes])
@@ -76,12 +73,6 @@ if test "$enable_mex_ms_sbvar" = yes; then
test "$has_gsl" != yes && AC_MSG_ERROR([GSL cannot be found. If you want to skip the compilation of the MS-SBVAR MEX, pass the --disable-mex-ms-sbvar flag.])
fi
-# Check for libmatio, needed by MEX files using Dynare++ code
-if test "$enable_mex_gensylv_korder" = yes; then
- AX_MATIO
- test "$has_matio" != yes && AC_MSG_ERROR([libmatio cannot be found. If you want to skip the compilation of gensylv and k-order MEX, pass the --disable-mex-gensylv-korder flag.])
-fi
-
# Check for libslicot, needed by kalman_steady_state
if test "$enable_mex_kalman_steady_state" = yes; then
# FCFLAGS must be temporarily modified, because otherwise -fno-underscoring is not
@@ -98,7 +89,6 @@ fi
case ${host_os} in
*mingw32*)
GSL_LIBS="-Wl,-Bstatic $GSL_LIBS -Wl,-Bdynamic"
- LIBADD_MATIO="-Wl,-Bstatic $LIBADD_MATIO -Wl,-Bdynamic"
LIBADD_SLICOT="-Wl,-Bstatic $LIBADD_SLICOT -Wl,-Bdynamic"
;;
esac
@@ -115,12 +105,6 @@ AC_SUBST([M2HTML])
AM_CONDITIONAL([HAVE_M2HTML], [test "x$M2HTML" != "x"])
# Construct final output message
-if test "$enable_mex_gensylv_korder" = yes; then
- BUILD_GENSYLV_KORDER_MEX_MATLAB="yes"
-else
- BUILD_GENSYLV_KORDER_MEX_MATLAB="no"
-fi
-
if test "$enable_mex_kalman_steady_state" = yes; then
BUILD_KALMAN_STEADY_STATE_MATLAB="yes"
else
@@ -139,7 +123,6 @@ Dynare is now configured for building the following components...
Binaries (with "make"):
MEX files for MATLAB (except those listed below): yes
- Gensylv and k-order MEX files for MATLAB: $BUILD_GENSYLV_KORDER_MEX_MATLAB
MS-SBVAR MEX files for MATLAB: $BUILD_MS_SBVAR_MEX_MATLAB
Kalman Steady State MEX file for MATLAB: $BUILD_KALMAN_STEADY_STATE_MATLAB
M2HTML documentation: $BUILD_M2HTML
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 7c47da953..8e93bdd0a 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -1,11 +1,6 @@
ACLOCAL_AMFLAGS = -I ../../../m4
-SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update
-
-# libkorder must come first
-if ENABLE_MEX_GENSYLV_KORDER
-SUBDIRS += libkorder gensylv k_order_perturbation k_order_welfare
-endif
+SUBDIRS = mjdgges kronecker bytecode sobol perfect_foresight_problem num_procs block_trust_region disclyap_fast libkordersim local_state_space_iterations folded_to_unfolded_dr k_order_simul k_order_mean cycle_reduction logarithmic_reduction riccati_update libkorder gensylv k_order_perturbation k_order_welfare
if ENABLE_MEX_MS_SBVAR
SUBDIRS += ms_sbvar
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 2b0f46f48..56cbbaf4d 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -66,9 +66,6 @@ AC_SUBST([LIBADD_DLOPEN])
AC_ARG_ENABLE([mex-ms-sbvar], AS_HELP_STRING([--disable-mex-ms-sbvar], [disable compilation of the MS-SBVAR MEX]), [], [enable_mex_ms_sbvar=yes])
AM_CONDITIONAL([ENABLE_MEX_MS_SBVAR], [test "$enable_mex_ms_sbvar" = yes])
-AC_ARG_ENABLE([mex-gensylv-korder], AS_HELP_STRING([--disable-mex-gensylv-korder], [disable compilation of gensylv and k-order MEX]), [], [enable_mex_gensylv_korder=yes])
-AM_CONDITIONAL([ENABLE_MEX_GENSYLV_KORDER], [test "$enable_mex_gensylv_korder" = yes])
-
AC_ARG_ENABLE([mex-kalman-steady-state], AS_HELP_STRING([--disable-mex-kalman-steady-state], [disable compilation of the kalman_steady_state MEX]), [], [enable_mex_kalman_steady_state=yes])
AM_CONDITIONAL([ENABLE_MEX_KALMAN_STEADY_STATE], [test "$enable_mex_kalman_steady_state" = yes])
@@ -78,10 +75,10 @@ if test "$enable_mex_ms_sbvar" = yes; then
test "$has_gsl" != yes && AC_MSG_ERROR([GSL cannot be found. If you want to skip the compilation of the MS-SBVAR MEX, pass the --disable-mex-ms-sbvar flag.])
fi
-# Check for libmatio, needed by MEX files using Dynare++ code, and by ms-sbvar (the latter only under Octave, as an alternative to MATLAB's libmat)
-if test "$enable_mex_gensylv_korder" = yes -o "$enable_mex_ms_sbvar" = yes; then
+# Check for libmatio, needed by ms-sbvar (the latter only under Octave, as an alternative to MATLAB's libmat)
+if test "$enable_mex_ms_sbvar" = yes; then
AX_MATIO
- test "$has_matio" != yes && AC_MSG_ERROR([libmatio cannot be found. If you want to skip the compilation of MS-SBVAR, gensylv and k-order MEX, pass the --disable-mex-gensylv-korder and --disable-mex-ms-sbvar flags.])
+ test "$has_matio" != yes && AC_MSG_ERROR([libmatio cannot be found. If you want to skip the compilation of MS-SBVAR, pass the --disable-mex-ms-sbvar flags.])
fi
# Check for libslicot, needed by kalman_steady_state
@@ -127,12 +124,6 @@ esac
AM_CONDITIONAL([LINK_OCTAVE_LIBS], [test ${link_octave_libs} = yes])
# Construct final output message
-if test "$enable_mex_gensylv_korder" = yes; then
- BUILD_GENSYLV_KORDER_MEX_OCTAVE="yes"
-else
- BUILD_GENSYLV_KORDER_MEX_OCTAVE="no"
-fi
-
if test "$enable_mex_kalman_steady_state" = yes; then
BUILD_KALMAN_STEADY_STATE_OCTAVE="yes"
else
@@ -151,7 +142,6 @@ Dynare is now configured for building the following components...
Binaries (with "make"):
MEX files for Octave (except those listed below): yes
- Gensylv and k-order MEX for Octave: $BUILD_GENSYLV_KORDER_MEX_OCTAVE
MS-SBVAR MEX files for Octave: $BUILD_MS_SBVAR_MEX_OCTAVE
Kalman Steady State MEX file for Octave: $BUILD_KALMAN_STEADY_STATE_OCTAVE
diff --git a/mex/sources/libkorder/integ/precalc_quadrature.hh b/mex/sources/libkorder/integ/precalc_quadrature.hh
deleted file mode 100644
index fca4d80e3..000000000
--- a/mex/sources/libkorder/integ/precalc_quadrature.hh
+++ /dev/null
@@ -1,1748 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// 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 = √π
- 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
- };
diff --git a/mex/sources/libkorder/integ/product.cc b/mex/sources/libkorder/integ/product.cc
deleted file mode 100644
index ad9f7296d..000000000
--- a/mex/sources/libkorder/integ/product.cc
+++ /dev/null
@@ -1,143 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "product.hh"
-#include "symmetry.hh"
-
-#include
-#include
-
-/* This constructs a product iterator corresponding to index (j0,0,…,0). */
-
-prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
- : prodq(q), level(l), npoints(q.uquad.numPoints(l)),
- jseq(q.dimen(), 0),
- end_flag(false),
- sig{q.dimen()},
- p{q.dimen()}
-{
- if (j0 < npoints)
- {
- jseq[0] = j0;
- setPointAndWeight();
- }
- else
- end_flag = true;
-}
-
-bool
-prodpit::operator==(const prodpit &ppit) const
-{
- return &prodq == &ppit.prodq && end_flag == ppit.end_flag && jseq == ppit.jseq;
-}
-
-prodpit &
-prodpit::operator++()
-{
- 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. */
-
-void
-prodpit::setPointAndWeight()
-{
- 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. */
-
-void
-prodpit::print() const
-{
- auto ff = std::cout.flags();
- std::cout << "j=[";
- for (int i = 0; i < prodq.dimen(); i++)
- std::cout << std::setw(2) << jseq[i];
- std::cout << std::showpos << std::fixed << std::setprecision(3)
- << "] " << std::setw(4) << w << "*(";
- for (int i = 0; i < prodq.dimen()-1; i++)
- std::cout << std::setw(4) << p[i] << ' ';
- std::cout << std::setw(4) << p[prodq.dimen()-1] << ')' << std::endl;
- std::cout.flags(ff);
-}
-
-ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq)
- : QuadratureImpl(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,…,0) where j0=ti·npoints/tn. */
-
-prodpit
-ProductQuadrature::begin(int ti, int tn, int l) const
-{
- // TODO: raise if l.
- */
-
-// Product quadrature.
-
-/* This file defines a product multidimensional quadrature. If Qₖ$ denotes the
- one dimensional quadrature, then the product quadrature Q of k level and
- dimension d takes the form:
-
- nₖ nₖ
- Qf = ∑ … ∑ w_i₁·…·w_{i_d} f(x_i₁,…,x_{i_d})
- i₁=1 i_d=1
-
- which can be written in terms of the one dimensional quadrature Qₖ as
-
- Qf=(Qₖ⊗…⊗Qₖ)f
-
- Here we define the product quadrature iterator prodpit and plug it into
- QuadratureImpl to obtains ProductQuadrature. */
-
-#ifndef PRODUCT_H
-#define PRODUCT_H
-
-#include "int_sequence.hh"
-#include "vector_function.hh"
-#include "quadrature.hh"
-
-/* 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,…,0), which is used by begin()
- dictated by QuadratureImpl. */
-
-class ProductQuadrature;
-
-class prodpit
-{
-protected:
- const ProductQuadrature &prodq;
- int level{0};
- int npoints{0};
- IntSequence jseq;
- bool end_flag{true};
- ParameterSignal sig;
- Vector p;
- double w;
-public:
- prodpit() = default;
- prodpit(const ProductQuadrature &q, int j0, int l);
- prodpit(const prodpit &ppit) = default;
- ~prodpit() = default;
- bool operator==(const prodpit &ppit) const;
- bool
- operator!=(const prodpit &ppit) const
- {
- return !operator==(ppit);
- }
- prodpit &operator=(const prodpit &spit) = delete;
- 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. */
-
-class ProductQuadrature : public QuadratureImpl
-{
- friend class prodpit;
- const OneDQuadrature &uquad;
-public:
- ProductQuadrature(int d, const OneDQuadrature &uq);
- ~ProductQuadrature() override = default;
- int
- numEvals(int l) const override
- {
- 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 override;
-};
-
-#endif
diff --git a/mex/sources/libkorder/integ/quadrature.cc b/mex/sources/libkorder/integ/quadrature.cc
deleted file mode 100644
index 9053240e4..000000000
--- a/mex/sources/libkorder/integ/quadrature.cc
+++ /dev/null
@@ -1,42 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "quadrature.hh"
-#include "precalc_quadrature.hh"
-
-#include
-
-void
-OneDPrecalcQuadrature::calcOffsets()
-{
- offsets[0] = 0;
- for (int i = 1; i < num_levels; i++)
- offsets[i] = offsets[i-1] + num_points[i-1];
-}
-
-GaussHermite::GaussHermite()
- : OneDPrecalcQuadrature(gh_num_levels, gh_num_points, gh_weights, gh_points)
-{
-}
-
-GaussLegendre::GaussLegendre()
- : OneDPrecalcQuadrature(gl_num_levels, gl_num_points, gl_weights, gl_points)
-{
-}
diff --git a/mex/sources/libkorder/integ/quadrature.hh b/mex/sources/libkorder/integ/quadrature.hh
deleted file mode 100644
index 5f1bff9fc..000000000
--- a/mex/sources/libkorder/integ/quadrature.hh
+++ /dev/null
@@ -1,329 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Quadrature.
-
-/* 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
-
- N
- ∫ f(x)dx ≃ ∑ wᵢxᵢ
- ⁱ⁼¹
-
- 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ᵢ and wᵢ for i=1,…,N. All the iterators must be able to go through only a
- portion of the set i=1,…,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. */
-
-#ifndef QUADRATURE_H
-#define QUADRATURE_H
-
-#include
-#include
-#include
-
-#include "vector_function.hh"
-#include "int_sequence.hh"
-#include "sthread.hh"
-
-/* 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. */
-
-class OneDQuadrature
-{
-public:
- virtual ~OneDQuadrature() = default;
- 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. */
-
-class Quadrature
-{
-protected:
- int dim;
-public:
- Quadrature(int d) : dim(d)
- {
- }
- virtual ~Quadrature() = default;
- 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. */
-
-template
-class QuadratureImpl;
-
-template
-class IntegrationWorker : public sthread::detach_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)
- {
- }
-
- /* 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 something smarter should be implemented. */
-
- void
- operator()(std::mutex &mut) override
- {
- _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);
- }
-
- {
- std::unique_lock lk{mut};
- 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. */
-
-template
-class QuadratureImpl : public Quadrature
-{
- friend class IntegrationWorker<_Tpit>;
-public:
- QuadratureImpl(int d) : Quadrature(d)
- {
- }
-
- /* Just fill a thread group with workes and run it. */
- void
- integrate(VectorFunctionSet &fs, int level, Vector &out) const override
- {
- // TODO: out.length()==func.outdim()
- // TODO: dim == func.indim()
- out.zeros();
- sthread::detach_thread_group gr;
- for (int ti = 0; ti < fs.getNum(); ti++)
- gr.insert(std::make_unique>(*this, fs.getFunc(ti),
- level, ti, fs.getNum(), out));
- gr.run();
- }
- void
- integrate(const VectorFunction &func,
- int level, int tn, Vector &out) const override
- {
- VectorFunctionSet fs(func, tn);
- integrate(fs, level, out);
- }
-
- /* Just for debugging. */
- void
- savePoints(const std::string &fname, int level) const
- {
- std::ofstream fd{fname, std::ios::out | std::ios::trunc};
- if (fd.fail())
- {
- // TODO: raise
- std::cerr << "Cannot open file " << fname << " for writing." << std::endl;
- exit(EXIT_FAILURE);
- }
- _Tpit beg = begin(0, 1, level);
- _Tpit end = begin(1, 1, level);
- fd << std::setprecision(12);
- for (_Tpit run = beg; run != end; ++run)
- {
- fd << std::setw(16) << run.weight();
- for (int i = 0; i < dimen(); i++)
- fd << '\t' << std::setw(16) << run.point()[i];
- fd << std::endl;
- }
- fd.close();
- }
-
- _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;
-};
-
-/* This is only an interface to a precalculated data in file
- precalc_quadrature.hh 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. */
-
-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();
- }
- ~OneDPrecalcQuadrature() override = default;
- int
- numLevels() const override
- {
- return num_levels;
- }
- int
- numPoints(int level) const override
- {
- return num_points[level-1];
- }
- double
- point(int level, int i) const override
- {
- return points[offsets[level-1]+i];
- }
- double
- weight(int level, int i) const override
- {
- return weights[offsets[level-1]+i];
- }
-protected:
- void calcOffsets();
-};
-
-/* Just precalculated Gauss-Hermite quadrature. This quadrature integrates
- exactly integrals
-
- +∞
- ∫ xᵏe^{−x²}dx
- −∞
-
- 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 precalc_quadrature.hh for available levels. */
-
-class GaussHermite : public OneDPrecalcQuadrature
-{
-public:
- GaussHermite();
-};
-
-/* Just precalculated Gauss-Legendre quadrature. This quadrature integrates
- exactly integrals
-
- ₁
- ∫ xᵏdx
- ⁰
-
- for level k.
-
- Check precalc_quadrature.hh for available levels. */
-
-class GaussLegendre : public OneDPrecalcQuadrature
-{
-public:
- GaussLegendre();
-};
-
-#endif
diff --git a/mex/sources/libkorder/integ/quasi_mcarlo.cc b/mex/sources/libkorder/integ/quasi_mcarlo.cc
deleted file mode 100644
index 4cdde2ce9..000000000
--- a/mex/sources/libkorder/integ/quasi_mcarlo.cc
+++ /dev/null
@@ -1,214 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "quasi_mcarlo.hh"
-
-#include
-#include
-#include
-#include
-
-/* 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::RadicalInverse(int n, int b, int mxn)
- : num(n), base(b), maxn(mxn),
- coeff(static_cast(floor(log(static_cast(maxn))/log(static_cast(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₀ c₁ cⱼ
- ── + ── + … + ────
- b b² bʲ⁺¹
-
- which is evaluated as:
-
- ⎛ ⎛⎛cⱼ 1 cⱼ₋₁⎞ 1 cⱼ₋₂⎞ ⎞ 1 c₀
- ⎢…⎢⎢──·─ + ────⎥·─ + ────⎥…⎥·─ + ──
- ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
-
-
- Now with permutation π, we have:
-
- ⎛ ⎛⎛π(cⱼ) 1 π(cⱼ₋₁)⎞ 1 π(cⱼ₋₂)⎞ ⎞ 1 π(c₀)
- ⎢…⎢⎢─────·─ + ───────⎥·─ + ───────⎥…⎥·─ + ─────
- ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
-*/
-
-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. */
-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. */
-void
-RadicalInverse::print() const
-{
- std::cout << "n=" << num << " b=" << base << " c=";
- coeff.print();
-}
-
-/* Here we have the first 170 primes. This means that we are not able to
- integrate dimensions greater than 170. */
-
-std::array 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::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.emplace_back(num, primes[i], maxn);
- eval();
-}
-
-/* This calls RadicalInverse::increase() for all radical inverses and calls
- eval(). */
-
-void
-HaltonSequence::increase()
-{
- for (auto &i : ri)
- i.increase();
- num++;
- if (num <= maxn)
- eval();
-}
-
-/* This sets point ‘pt’ to radical inverse evaluations in each dimension. */
-void
-HaltonSequence::eval()
-{
- for (unsigned int i = 0; i < ri.size(); i++)
- pt[i] = ri[i].eval(per);
-}
-
-/* Debug print. */
-void
-HaltonSequence::print() const
-{
- auto ff = std::cout.flags();
- for (const auto &i : ri)
- i.print();
- std::cout << "point=[ "
- << std::fixed << std::setprecision(6);
- for (unsigned int i = 0; i < ri.size(); i++)
- std::cout << std::setw(7) << pt[i] << ' ';
- std::cout << ']' << std::endl;
- std::cout.flags(ff);
-}
-
-qmcpit::qmcpit(const QMCSpecification &s, int n)
- : spec(s), halton{n, s.level(), s.dimen(), s.getPerScheme()},
- sig{s.dimen()}
-{
-}
-
-bool
-qmcpit::operator==(const qmcpit &qpit) const
-{
- return &spec == &qpit.spec && halton.getNum() == qpit.halton.getNum();
-}
-
-qmcpit &
-qmcpit::operator++()
-{
- halton.increase();
- return *this;
-}
-
-double
-qmcpit::weight() const
-{
- return 1.0/spec.level();
-}
-
-int
-WarnockPerScheme::permute(int i, int base, int c) const
-{
- return (c+i) % base;
-}
-
-int
-ReversePerScheme::permute(int i, int base, int c) const
-{
- return (base-c) % base;
-}
diff --git a/mex/sources/libkorder/integ/quasi_mcarlo.hh b/mex/sources/libkorder/integ/quasi_mcarlo.hh
deleted file mode 100644
index 7a50e8ffe..000000000
--- a/mex/sources/libkorder/integ/quasi_mcarlo.hh
+++ /dev/null
@@ -1,256 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019-2022 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Quasi Monte Carlo quadrature.
-
-/* 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:
-
- ∫ f(x)dx
- [0,1]ⁿ
-
- The quadrature for a function of normally distributed parameters is named
- QMCarloNormalQuadrature and integrates:
-
- 1
- ──────── ∫ f(x)e^{−½xᵀx}dx
- √{(2π)ⁿ} [−∞,+∞]ⁿ
-
- 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. */
-
-#ifndef QUASI_MCARLO_H
-#define QUASI_MCARLO_H
-
-#include "int_sequence.hh"
-#include "quadrature.hh"
-
-#include "Vector.hh"
-
-#include
-
-/* 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,…,base−1. */
-
-class PermutationScheme
-{
-public:
- PermutationScheme() = default;
- virtual ~PermutationScheme() = default;
- virtual int permute(int i, int base, int c) const = 0;
-};
-
-/* This class represents an integer number ‘num’ as c₀+c₁b+c₂b²+…+cⱼbʲ, where b
- is ‘base’ and c₀,…,cⱼ 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. */
-
-class RadicalInverse
-{
- int num;
- int base;
- int maxn;
- int j;
- IntSequence coeff;
-public:
- RadicalInverse(int n, int b, int mxn);
- RadicalInverse(const RadicalInverse &ri) = default;
- RadicalInverse &operator=(const RadicalInverse &radi) = default;
- double eval(const PermutationScheme &p) const;
- void increase();
- void print() const;
-};
-
-/* This is a vector of RadicalInverses, 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. */
-
-class HaltonSequence
-{
-private:
- static std::array primes;
-protected:
- int num;
- int maxn;
- std::vector ri;
- const PermutationScheme &per;
- Vector pt;
-public:
- HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p);
- HaltonSequence(const HaltonSequence &hs) = default;
- HaltonSequence &operator=(const HaltonSequence &hs) = delete;
- void increase();
- const Vector &
- point() const
- {
- return pt;
- }
- 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. */
-
-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() = default;
- 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. */
-
-class qmcpit
-{
-protected:
- const QMCSpecification &spec;
- HaltonSequence halton;
- ParameterSignal sig;
-public:
- qmcpit(const QMCSpecification &s, int n);
- qmcpit(const qmcpit &qpit) = default;
- virtual ~qmcpit() = default;
- bool operator==(const qmcpit &qpit) const;
- bool
- operator!=(const qmcpit &qpit) const
- {
- return !operator==(qpit);
- }
- qmcpit &operator=(const qmcpit &qpit) = delete;
- 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(). */
-
-class QMCarloCubeQuadrature : public QuadratureImpl, public QMCSpecification
-{
-public:
- QMCarloCubeQuadrature(int d, int l, const PermutationScheme &p)
- : QuadratureImpl(d), QMCSpecification(d, l, p)
- {
- }
- ~QMCarloCubeQuadrature() override = default;
- int
- numEvals(int l) const override
- {
- return l;
- }
-protected:
- qmcpit
- begin(int ti, int tn, int lev) const override
- {
- return qmcpit(*this, ti*level()/tn + 1);
- }
-};
-
-/* Declares Warnock permutation scheme. */
-class WarnockPerScheme : public PermutationScheme
-{
-public:
- int permute(int i, int base, int c) const override;
-};
-
-/* Declares reverse permutation scheme. */
-class ReversePerScheme : public PermutationScheme
-{
-public:
- int permute(int i, int base, int c) const override;
-};
-
-/* Declares no permutation (identity) scheme. */
-class IdentityPerScheme : public PermutationScheme
-{
-public:
- int
- permute(int i, int base, int c) const override
- {
- return c;
- }
-};
-
-#endif
diff --git a/mex/sources/libkorder/integ/smolyak.cc b/mex/sources/libkorder/integ/smolyak.cc
deleted file mode 100644
index df123785c..000000000
--- a/mex/sources/libkorder/integ/smolyak.cc
+++ /dev/null
@@ -1,220 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "smolyak.hh"
-#include "symmetry.hh"
-
-#include
-#include
-
-/* 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::smolpit(const SmolyakQuadrature &q, unsigned int isum)
- : smolq(q), isummand(isum),
- jseq(q.dimen(), 0),
- sig{q.dimen()},
- p{q.dimen()}
-{
- if (isummand < q.numSummands())
- setPointAndWeight();
-}
-
-bool
-smolpit::operator==(const smolpit &spit) const
-{
- return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq;
-}
-
-/* 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 &
-smolpit::operator++()
-{
- 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. */
-
-void
-smolpit::setPointAndWeight()
-{
- // todo: raise if 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 *= PascalTriangle::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. */
-void
-smolpit::print() const
-{
- auto ff = std::cout.flags();
- std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": [";
- for (int i = 0; i < smolq.dimen(); i++)
- std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' ';
- std::cout << "] j=[";
- for (int i = 0; i < smolq.dimen(); i++)
- std::cout << std::setw(2) << jseq[i] << ' ';
- std::cout << std::showpos << std::fixed << std::setprecision(3)
- << "] " << std::setw(4) << w << "*(";
- for (int i = 0; i < smolq.dimen()-1; i++)
- std::cout << std::setw(4) << p[i] << ' ';
- std::cout << std::setw(4) << p[smolq.dimen()-1] << ')' << std::endl;
- std::cout.flags(ff);
-}
-
-/* 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≤|k|≤l+d−1 and all kᵢ are positive integers. This is
- equivalent to going through all k such that l−d≤|k|≤l−1 and all kᵢ are
- non-negative integers. This is equivalent to going through d+1 dimensional
- sequences (k,x) such that |(k,x)|=l−1 and x=0,…,d−1. The resulting sequence
- of positive integers is obtained by adding 1 to all kᵢ. */
-
-SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
- : QuadratureImpl(d), level(l), uquad(uq)
-{
- // TODO: check l>1, l≥d
- // TODO: check l≥uquad.miLevel(), l≤uquad.maxLevel()
- int cum = 0;
- for (const auto &si : SymmetrySet(l-1, d+1))
- {
- if (si[d] <= d-1)
- {
- IntSequence lev(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. */
-
-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’. */
-
-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. */
-
-int
-SmolyakQuadrature::calcNumEvaluations(int lev) const
-{
- int cum = 0;
- for (const auto &si : SymmetrySet(lev-1, dim+1))
- {
- if (si[dim] <= dim-1)
- {
- IntSequence lev(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. */
-
-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;
-}
diff --git a/mex/sources/libkorder/integ/smolyak.hh b/mex/sources/libkorder/integ/smolyak.hh
deleted file mode 100644
index 561433593..000000000
--- a/mex/sources/libkorder/integ/smolyak.hh
+++ /dev/null
@@ -1,149 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Smolyak quadrature.
-
-/* This file defines Smolyak (sparse grid) multidimensional quadrature for
- non-nested underlying one dimensional quadrature. Let Q¹ₗ denote the one
- dimensional quadrature of l level. Let nₗ denote a number of points in the l
- level. Than the Smolyak quadrature can be defined as
-
- ⎛ d−1 ⎞
- Qᵈf = ∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
- l≤|k|≤l+d−1 ⎝|k|−l⎠
-
- where d is the dimension, k is d-dimensional sequence of integers, and |k|
- denotes the sum of the sequence.
-
- Here we define smolpit as Smolyak iterator and SmolyakQuadrature. */
-
-#ifndef SMOLYAK_H
-#define SMOLYAK_H
-
-#include "int_sequence.hh"
-#include "tl_static.hh"
-#include "vector_function.hh"
-#include "quadrature.hh"
-#include "pascal_triangle.hh"
-
-/* 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. */
-
-class SmolyakQuadrature;
-
-class smolpit
-{
-protected:
- const SmolyakQuadrature &smolq;
- unsigned int isummand{0};
- IntSequence jseq;
- ParameterSignal sig;
- Vector p;
- double w;
-public:
- smolpit(const SmolyakQuadrature &q, unsigned int isum);
- smolpit(const smolpit &spit) = default;
- ~smolpit() = default;
- bool operator==(const smolpit &spit) const;
- bool
- operator!=(const smolpit &spit) const
- {
- return !operator==(spit);
- }
- smolpit &operator=(const smolpit &spit) = delete;
- 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:
-
- ⎛ d−1 ⎞
- ∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
- l≤|k|≤l+d−1 ⎝|k|−l⎠
-
- 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:
-
- d
- ∑ ∏ n_kᵢ
- ᵏ ⁱ⁼¹
-
- where the sum is done through all k before the current.
-
- The ‘levels’ and ‘levpoints’ vectors are used by smolpit. */
-
-class SmolyakQuadrature : public QuadratureImpl
-{
- friend class smolpit;
- int level;
- const OneDQuadrature &uquad;
- std::vector levels;
- std::vector levpoints;
- std::vector cumevals;
-public:
- SmolyakQuadrature(int d, int l, const OneDQuadrature &uq);
- ~SmolyakQuadrature() override = default;
- int numEvals(int level) const override;
- void designLevelForEvals(int max_eval, int &lev, int &evals) const;
-protected:
- smolpit begin(int ti, int tn, int level) const override;
- unsigned int
- numSummands() const
- {
- return levels.size();
- }
-private:
- int calcNumEvaluations(int level) const;
-};
-
-#endif
diff --git a/mex/sources/libkorder/integ/tests/Makefile.am b/mex/sources/libkorder/integ/tests/Makefile.am
deleted file mode 100644
index 2b7dcf48b..000000000
--- a/mex/sources/libkorder/integ/tests/Makefile.am
+++ /dev/null
@@ -1,9 +0,0 @@
-check_PROGRAMS = tests
-
-tests_SOURCES = tests.cc
-tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I../../utils/cc -I$(top_srcdir)/mex/sources
-tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
-tests_LDADD = ../../sylv/cc/libsylv.a ../cc/libinteg.a ../../tl/cc/libtl.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_MATIO)
-
-check-local:
- ./tests
diff --git a/mex/sources/libkorder/integ/tests/tests.cc b/mex/sources/libkorder/integ/tests/tests.cc
deleted file mode 100644
index a3f778885..000000000
--- a/mex/sources/libkorder/integ/tests/tests.cc
+++ /dev/null
@@ -1,550 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "GeneralMatrix.hh"
-#include
-#include "SylvException.hh"
-
-#include "rfs_tensor.hh"
-#include "normal_moments.hh"
-
-#include "vector_function.hh"
-#include "quadrature.hh"
-#include "smolyak.hh"
-#include "product.hh"
-#include "quasi_mcarlo.hh"
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-/* Evaluates unfolded (Dx)ᵏ power, where x is a vector, D is a Cholesky factor
- (lower triangular) */
-class MomentFunction : public VectorFunction
-{
- GeneralMatrix D;
- int k;
-public:
- MomentFunction(const GeneralMatrix &inD, int kk)
- : VectorFunction(inD.nrows(), UFSTensor::calcMaxOffset(inD.nrows(), kk)),
- D(inD), k(kk)
- {
- }
- MomentFunction(const MomentFunction &func) = default;
-
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- if (point.length() != indim() || out.length() != outdim())
- {
- std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl;
- std::exit(EXIT_FAILURE);
- }
- Vector y(point);
- y.zeros();
- D.multaVec(y, point);
- URSingleTensor ypow(y, k);
- out.zeros();
- out.add(1.0, ypow.getData());
-}
-
-class TensorPower : public VectorFunction
-{
- int k;
-public:
- TensorPower(int nvar, int kk)
- : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
- {
- }
- TensorPower(const TensorPower &func) = default;
-
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- if (point.length() != indim() || out.length() != outdim())
- {
- std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl;
- std::exit(EXIT_FAILURE);
- }
- URSingleTensor ypow(point, k);
- out.zeros();
- out.add(1.0, ypow.getData());
-}
-
-/* Evaluates (1+1/d)ᵈ(x₁·…·x_d)^(1/d), its integral over [0,1]ᵈ
- is 1.0, and its variation grows exponentially */
-class Function1 : public VectorFunction
-{
- int dim;
-public:
- Function1(int d)
- : VectorFunction(d, 1), dim(d)
- {
- }
- Function1(const Function1 &f)
- : VectorFunction(f.indim(), f.outdim()), dim(f.dim)
- {
- }
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- if (point.length() != dim || out.length() != 1)
- {
- std::cerr << "Wrong length of vectors in Function1::eval" << std::endl;
- std::exit(EXIT_FAILURE);
- }
- double r = 1;
- for (int i = 0; i < dim; i++)
- r *= point[i];
- r = pow(r, 1.0/dim);
- r *= pow(1.0 + 1.0/dim, static_cast(dim));
- out[0] = r;
-}
-
-// Evaluates Function1 but with transformation xᵢ=0.5(yᵢ+1)
-// This makes the new function integrate over [−1,1]ᵈ to 1.0
-class Function1Trans : public Function1
-{
-public:
- Function1Trans(int d)
- : Function1(d)
- {
- }
- Function1Trans(const Function1Trans &func) = default;
-
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- Vector p(point.length());
- for (int i = 0; i < p.length(); i++)
- p[i] = 0.5*(point[i]+1);
- Function1::eval(p, sig, out);
- out.mult(pow(0.5, indim()));
-}
-
-/* WallTimer class. Constructor saves the wall time, destructor cancels the
- current time from the saved, and prints the message with time information */
-class WallTimer
-{
- std::string mes;
- std::chrono::time_point start;
- bool new_line;
-public:
- WallTimer(std::string m, bool nl = true)
- : mes{m}, start{std::chrono::high_resolution_clock::now()}, new_line{nl}
- {
- }
- ~WallTimer()
- {
- auto end = std::chrono::high_resolution_clock::now();
- std::chrono::duration duration = end - start;
- std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count();
- if (new_line)
- std::cout << std::endl;
- }
-};
-
-/****************************************************/
-/* declaration of TestRunnable class */
-/****************************************************/
-class TestRunnable
-{
-public:
- const std::string name;
- int dim; // dimension of the solved problem
- int nvar; // number of variables of the solved problem
- TestRunnable(std::string name_arg, int d, int nv)
- : name{move(name_arg)}, dim(d), nvar(nv)
- {
- }
- virtual ~TestRunnable() = default;
- bool test() const;
- virtual bool run() const = 0;
-protected:
- static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level);
- static bool product_normal_moments(const GeneralMatrix &m, int imom, int level);
- static bool qmc_normal_moments(const GeneralMatrix &m, int imom, int level);
- static bool smolyak_product_cube(const VectorFunction &func, const Vector &res,
- double tol, int level);
- static bool qmc_cube(const VectorFunction &func, double res, double tol, int level);
-};
-
-bool
-TestRunnable::test() const
-{
- std::cout << "Running test <" << name << ">" << std::endl;
- bool passed;
- {
- WallTimer tim("Wall clock time ", false);
- passed = run();
- }
- if (passed)
- {
- std::cout << "............................ passed" << std::endl << std::endl;
- return passed;
- }
- else
- {
- std::cout << "............................ FAILED" << std::endl << std::endl;
- return passed;
- }
-}
-
-/****************************************************/
-/* definition of TestRunnable static methods */
-/****************************************************/
-bool
-TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
-{
- // First make m·mᵀ and then Cholesky factor
- GeneralMatrix msq(m * transpose(m));
-
- // Make vector function
- int dim = m.nrows();
- TensorPower tp(dim, imom);
- GaussConverterFunction func(tp, msq);
-
- // Smolyak quadrature
- Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
- {
- WallTimer tim("\tSmolyak quadrature time: ");
- GaussHermite gs;
- SmolyakQuadrature quad(dim, level, gs);
- quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, smol_out);
- std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
- }
-
- // Check against theoretical moments
- UNormalMoments moments(imom, msq);
- smol_out.add(-1.0, moments.get(Symmetry{imom}).getData());
- std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
- return smol_out.getMax() < 1.e-7;
-}
-
-bool
-TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
-{
- // First make m·mᵀ and then Cholesky factor
- GeneralMatrix msq(m * transpose(m));
-
- // Make vector function
- int dim = m.nrows();
- TensorPower tp(dim, imom);
- GaussConverterFunction func(tp, msq);
-
- // Product quadrature
- Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
- {
- WallTimer tim("\tProduct quadrature time: ");
- GaussHermite gs;
- ProductQuadrature quad(dim, gs);
- quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, prod_out);
- std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
- }
-
- // Check against theoretical moments
- UNormalMoments moments(imom, msq);
- prod_out.add(-1.0, moments.get(Symmetry{imom}).getData());
- std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
- return prod_out.getMax() < 1.e-7;
-}
-
-bool
-TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res,
- double tol, int level)
-{
- if (res.length() != func.outdim())
- {
- std::cerr << "Incompatible dimensions of check value and function." << std::endl;
- std::exit(EXIT_FAILURE);
- }
-
- GaussLegendre glq;
- Vector out(func.outdim());
- double smol_error;
- double prod_error;
- {
- WallTimer tim("\tSmolyak quadrature time: ");
- SmolyakQuadrature quad(func.indim(), level, glq);
- quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
- out.add(-1.0, res);
- smol_error = out.getMax();
- std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
- std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_error << std::endl;
- }
- {
- WallTimer tim("\tProduct quadrature time: ");
- ProductQuadrature quad(func.indim(), glq);
- quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
- out.add(-1.0, res);
- prod_error = out.getMax();
- std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
- std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_error << std::endl;
- }
-
- return smol_error < tol && prod_error < tol;
-}
-
-bool
-TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int level)
-{
- Vector r(1);
- double error1;
- {
- WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time: ");
- WarnockPerScheme wps;
- QMCarloCubeQuadrature qmc(func.indim(), level, wps);
- qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
- error1 = std::max(res - r[0], r[0] - res);
- std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl;
- }
- double error2;
- {
- WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time: ");
- ReversePerScheme rps;
- QMCarloCubeQuadrature qmc(func.indim(), level, rps);
- qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
- error2 = std::max(res - r[0], r[0] - res);
- std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl;
- }
- double error3;
- {
- WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time: ");
- IdentityPerScheme ips;
- QMCarloCubeQuadrature qmc(func.indim(), level, ips);
- qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
- error3 = std::max(res - r[0], r[0] - res);
- std::cout << "\tQuasi-Monte Carlo (no scrambling) error: " << std::setw(16) << std::setprecision(12) << error3 << std::endl;
- }
-
- return error1 < tol && error2 < tol && error3 < tol;
-}
-
-/****************************************************/
-/* definition of TestRunnable subclasses */
-/****************************************************/
-class SmolyakNormalMom1 : public TestRunnable
-{
-public:
- SmolyakNormalMom1()
- : TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2)
- {
- }
-
- bool
- run() const override
- {
- GeneralMatrix m(2, 2);
- m.zeros();
- m.get(0, 0) = 1;
- m.get(1, 1) = 1;
- return smolyak_normal_moments(m, 4, 4);
- }
-};
-
-class SmolyakNormalMom2 : public TestRunnable
-{
-public:
- SmolyakNormalMom2()
- : TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3)
- {
- }
-
- bool
- run() const override
- {
- GeneralMatrix m(3, 3);
- m.zeros();
- m.get(0, 0) = 1;
- m.get(0, 2) = 0.5;
- m.get(1, 1) = 1;
- m.get(1, 0) = 0.5;
- m.get(2, 2) = 2;
- m.get(2, 1) = 4;
- return smolyak_normal_moments(m, 8, 8);
- }
-};
-
-class ProductNormalMom1 : public TestRunnable
-{
-public:
- ProductNormalMom1()
- : TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2)
- {
- }
-
- bool
- run() const override
- {
- GeneralMatrix m(2, 2);
- m.zeros();
- m.get(0, 0) = 1;
- m.get(1, 1) = 1;
- return product_normal_moments(m, 4, 4);
- }
-};
-
-class ProductNormalMom2 : public TestRunnable
-{
-public:
- ProductNormalMom2()
- : TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3)
- {
- }
-
- bool
- run() const override
- {
- GeneralMatrix m(3, 3);
- m.zeros();
- m.get(0, 0) = 1;
- m.get(0, 2) = 0.5;
- m.get(1, 1) = 1;
- m.get(1, 0) = 0.5;
- m.get(2, 2) = 2;
- m.get(2, 1) = 4;
- return product_normal_moments(m, 8, 8);
- }
-};
-
-// Note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
-class F1GaussLegendre : public TestRunnable
-{
-public:
- F1GaussLegendre()
- : TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1)
- {
- }
-
- bool
- run() const override
- {
- Function1Trans f1(6);
- Vector res(1);
- res[0] = 1.0;
- return smolyak_product_cube(f1, res, 1e-2, 13);
- }
-};
-
-class F1QuasiMCarlo : public TestRunnable
-{
-public:
- F1QuasiMCarlo()
- : TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1)
- {
- }
-
- bool
- run() const override
- {
- Function1 f1(6);
- return qmc_cube(f1, 1.0, 1.e-4, 1000000);
- }
-};
-
-int
-main()
-{
- std::vector> all_tests;
- // Fill in vector of all tests
- all_tests.push_back(std::make_unique());
- all_tests.push_back(std::make_unique());
- all_tests.push_back(std::make_unique());
- all_tests.push_back(std::make_unique());
- all_tests.push_back(std::make_unique());
- all_tests.push_back(std::make_unique());
-
- // Find maximum dimension and maximum nvar
- int dmax = 0;
- int nvmax = 0;
- for (const auto &test : all_tests)
- {
- dmax = std::max(dmax, test->dim);
- nvmax = std::max(nvmax, test->nvar);
- }
- TLStatic::init(dmax, nvmax); // initialize library
-
- // Launch the tests
- int success = 0;
- for (const auto &test : all_tests)
- {
- try
- {
- if (test->test())
- success++;
- }
- catch (const TLException &e)
- {
- std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
- e.print();
- }
- catch (SylvException &e)
- {
- std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
- e.printMessage();
- }
- }
-
- int nfailed = all_tests.size() - success;
- std::cout << "There were " << nfailed << " tests that failed out of "
- << all_tests.size() << " tests run." << std::endl;
-
- if (nfailed)
- return EXIT_FAILURE;
- else
- return EXIT_SUCCESS;
-}
diff --git a/mex/sources/libkorder/integ/vector_function.cc b/mex/sources/libkorder/integ/vector_function.cc
deleted file mode 100644
index b639a6a42..000000000
--- a/mex/sources/libkorder/integ/vector_function.cc
+++ /dev/null
@@ -1,147 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "vector_function.hh"
-
-#include
-
-#include
-#include
-
-/* Just an easy constructor of sequence of booleans defaulting to change
- everywhere. */
-
-ParameterSignal::ParameterSignal(int n)
- : data(n, true)
-{
-}
-
-/* This sets ‘false’ (no change) before a given parameter, and ‘true’ (change)
- after the given parameter (including). */
-
-void
-ParameterSignal::signalAfter(int l)
-{
- for (size_t i = 0; i < std::min(static_cast(l), data.size()); i++)
- data[i] = false;
- for (size_t i = l; i < data.size(); i++)
- data[i] = true;
-}
-
-/* This constructs a function set hardcopying also the first. */
-VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
- : funcs(n)
-{
- for (int i = 0; i < n; i++)
- {
- func_copies.push_back(f.clone());
- funcs[i] = func_copies.back().get();
- }
-}
-
-/* This constructs a function set with shallow copy in the first and hard
- copies in others. */
-VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
- : funcs(n)
-{
- if (n > 0)
- funcs[0] = &f;
- for (int i = 1; i < n; i++)
- {
- func_copies.push_back(f.clone());
- funcs[i] = func_copies.back().get();
- }
-}
-
-/* Here we construct the object from the given function f and given
- variance-covariance matrix Σ=vcov. The matrix A is calculated as lower
- triangular and yields Σ=AAᵀ. */
-
-GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
- : VectorFunction(f), func(&f), A(vcov.nrows(), vcov.nrows()),
- multiplier(calcMultiplier())
-{
- // TODO: raise if A.nrows() ≠ indim()
- calcCholeskyFactor(vcov);
-}
-
-GaussConverterFunction::GaussConverterFunction(std::unique_ptr f, const GeneralMatrix &vcov)
- : VectorFunction(*f), func_storage{move(f)}, func{func_storage.get()}, A(vcov.nrows(), vcov.nrows()),
- multiplier(calcMultiplier())
-{
- // TODO: raise if A.nrows() ≠ indim()
- calcCholeskyFactor(vcov);
-}
-
-GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
- : VectorFunction(f), func_storage{f.func->clone()}, func{func_storage.get()}, A(f.A),
- multiplier(f.multiplier)
-{
-}
-
-/* Here we evaluate the function
-
- g(y) = 1/√(πⁿ) f(√2·Ay).
-
- Since the matrix A is lower triangular, the change signal for the function f
- will look like (0,…,0,1,…,1) where the first 1 is in the same position as
- the first change in the given signal ‘sig’ of the input y=point. */
-
-void
-GaussConverterFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- ParameterSignal s(sig);
- int i = 0;
- while (i < indim() && !sig[i])
- i++;
- s.signalAfter(i);
-
- Vector x(indim());
- x.zeros();
- A.multaVec(x, point);
- x.mult(sqrt(2.0));
-
- func->eval(x, s, out);
-
- out.mult(multiplier);
-}
-
-/* This returns 1/√(πⁿ). */
-
-double
-GaussConverterFunction::calcMultiplier() const
-{
- return sqrt(pow(M_PI, -1*indim()));
-}
-
-void
-GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
-{
- A = vcov;
-
- lapack_int rows = A.nrows(), lda = A.getLD();
- for (int i = 0; i < rows; i++)
- for (int j = i+1; j < rows; j++)
- A.get(i, j) = 0.0;
-
- lapack_int info;
- dpotrf("L", &rows, A.base(), &lda, &info);
- // TODO: raise if info≠1
-}
diff --git a/mex/sources/libkorder/integ/vector_function.hh b/mex/sources/libkorder/integ/vector_function.hh
deleted file mode 100644
index 8beee94dc..000000000
--- a/mex/sources/libkorder/integ/vector_function.hh
+++ /dev/null
@@ -1,186 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Vector function.
-
-/* This file defines interface for functions taking a vector as an input and
- returning a vector (with a different size) as an output. We are also
- introducing a parameter signalling; it is a boolean vector which tracks
- parameters which were changed from the previous call. The VectorFunction
- implementation can exploit this information and evaluate the function more
- efficiently. The information can be completely ignored.
-
- From the signalling reason, and from other reasons, the function evaluation
- is not const. */
-
-#ifndef VECTOR_FUNCTION_H
-#define VECTOR_FUNCTION_H
-
-#include "Vector.hh"
-#include "GeneralMatrix.hh"
-
-#include
-#include
-
-/* This is a simple class representing a vector of booleans. The items night be
- retrieved or changed, or can be set ‘true’ after some point. This is useful
- when we multiply the vector with lower triangular matrix.
-
- ‘true’ means that a parameter was changed. */
-
-class ParameterSignal
-{
-protected:
- std::vector data;
-public:
- ParameterSignal(int n);
- ParameterSignal(const ParameterSignal &sig) = default;
- ~ParameterSignal() = default;
- void signalAfter(int l);
- bool
- operator[](int i) const
- {
- return data[i];
- }
- std::vector::reference
- operator[](int i)
- {
- return data[i];
- }
-};
-
-/* This is the abstract class for vector function. At this level of abstraction
- we only need to know size of input vector and a size of output vector.
-
- The important thing here is a clone method, we will need to make hard copies
- of vector functions since the evaluations are not const. The hardcopies
- apply for parallelization. */
-
-class VectorFunction
-{
-protected:
- int in_dim;
- int out_dim;
-public:
- VectorFunction(int idim, int odim)
- : in_dim(idim), out_dim(odim)
- {
- }
- VectorFunction(const VectorFunction &func) = default;
- virtual ~VectorFunction() = default;
- virtual std::unique_ptr clone() const = 0;
- virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0;
- int
- indim() const
- {
- return in_dim;
- }
- int
- outdim() const
- {
- return out_dim;
- }
-};
-
-/* This makes ‘n’ copies of VectorFunction. The first constructor make exactly
- ‘n’ new copies, the second constructor copies only the pointer to the first
- and others are hard (real) copies.
-
- The class is useful for making a given number of copies at once, and this
- set can be reused many times if we need mupliple copis of the function (for
- example for paralelizing the code). */
-
-class VectorFunctionSet
-{
-private:
- // Stores the hard copies made by the class
- std::vector> func_copies;
-protected:
- std::vector funcs;
-public:
- VectorFunctionSet(const VectorFunction &f, int n);
- VectorFunctionSet(VectorFunction &f, int n);
- ~VectorFunctionSet() = default;
- VectorFunction &
- getFunc(int i)
- {
- return *(funcs[i]);
- }
- int
- getNum() const
- {
- return funcs.size();
- }
-};
-
-/* This class wraps another VectorFunction to allow integration of a function
- through normally distributed inputs. Namely, if one wants to integrate
-
- 1
- ─────────── ∫ f(x)e^{−½xᵀ|Σ|⁻¹x}dx
- √{(2π)ⁿ|Σ|}
-
- then if we write Σ=AAᵀ and x=√2·Ay, we get integral
-
- 1 1
- ─────────── ∫ f(√2·Ay)e^{−½yᵀy} √(2ⁿ)|A|dy = ───── ∫ f(√2·Ay)e^{−½yᵀy}dy
- √{(2π)ⁿ|Σ|} √(πⁿ)
-
- which means that a given function f we have to wrap to yield a function
-
- g(y) = 1/√(πⁿ) f(√2·Ay).
-
- This is exactly what this class is doing. This transformation is useful
- since the Gauss-Hermite points and weights are defined for weighting
- function e^{−y²}, so this transformation allows using Gauss-Hermite
- quadratures seemlessly in a context of integration through normally
- distributed inputs.
-
- The class maintains a pointer to the function f. When the object is
- constructed by the first constructor, the f is assumed to be owned by the
- caller. If the object of this class is copied, then f is copied and hence
- stored in a std::unique_ptr. The second constructor takes a smart pointer to
- the function and in that case the class takes ownership of f. */
-
-class GaussConverterFunction : public VectorFunction
-{
-private:
- std::unique_ptr func_storage;
-protected:
- VectorFunction *func;
- GeneralMatrix A;
- double multiplier;
-public:
- GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov);
- GaussConverterFunction(std::unique_ptr f, const GeneralMatrix &vcov);
- GaussConverterFunction(const GaussConverterFunction &f);
- ~GaussConverterFunction() override = default;
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-private:
- double calcMultiplier() const;
- void calcCholeskyFactor(const GeneralMatrix &vcov);
-};
-
-#endif
diff --git a/mex/sources/libkorder/kord/decision_rule.cc b/mex/sources/libkorder/kord/decision_rule.cc
index 9330a73aa..b433a7265 100644
--- a/mex/sources/libkorder/kord/decision_rule.cc
+++ b/mex/sources/libkorder/kord/decision_rule.cc
@@ -1,6 +1,6 @@
/*
* Copyright © 2004 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@@ -21,14 +21,9 @@
#include "kord_exception.hh"
#include "decision_rule.hh"
#include "dynamic_model.hh"
-#include "seed_generator.hh"
#include "SymSchurDecomp.hh"
-#include
-
-#include
-#include
#include
// FoldDecisionRule conversion from UnfoldDecisionRule
@@ -48,539 +43,3 @@ UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
for (const auto &it : fdr)
insert(std::make_unique::Ttensym>(*(it.second)));
}
-
-/* This runs simulations with an output to journal file. Note that we
- report how many simulations had to be thrown out due to Nan or Inf. */
-
-void
-SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov, Journal &journal)
-{
- JournalRecordPair paa(journal);
- paa << "Performing " << num_sim << " stochastic simulations for "
- << num_per << " periods burning " << num_burn << " initial periods" << endrec;
- simulate(num_sim, dr, start, vcov);
- int thrown = num_sim - data.size();
- if (thrown > 0)
- {
- JournalRecord rec(journal);
- rec << "I had to throw " << thrown << " simulations away due to Nan or Inf" << endrec;
- }
-}
-
-/* This runs a given number of simulations by creating
- SimulationWorker for each simulation and inserting them to the
- thread group. */
-
-void
-SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov)
-{
- std::vector rsrs;
- rsrs.reserve(num_sim);
-
- sthread::detach_thread_group gr;
- for (int i = 0; i < num_sim; i++)
- {
- RandomShockRealization sr(vcov, seed_generator::get_new_seed());
- rsrs.push_back(sr);
- gr.insert(std::make_unique(*this, dr, DecisionRule::emethod::horner,
- num_per+num_burn, start, rsrs.back()));
- }
- gr.run();
-}
-
-/* This adds the data with the realized shocks. It takes only periods
- which are not to be burnt. If the data is not finite, the both data
- and shocks are thrown away. */
-
-bool
-SimResults::addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st)
-{
- KORD_RAISE_IF(d.nrows() != num_y,
- "Incompatible number of rows for SimResults::addDataSets");
- KORD_RAISE_IF(d.ncols() != num_per+num_burn,
- "Incompatible number of cols for SimResults::addDataSets");
- bool ret = false;
- if (d.isFinite())
- {
- data.emplace_back(d, num_burn, num_per);
- shocks.emplace_back(ConstTwoDMatrix(sr.getShocks(), num_burn, num_per));
- if (num_burn == 0)
- start.emplace_back(st);
- else
- start.emplace_back(d.getCol(num_burn-1));
- ret = true;
- }
-
- return ret;
-}
-
-void
-SimResults::writeMat(const std::string &base, const std::string &lname) const
-{
- std::string matfile_name = base + ".mat";
- mat_t *matfd = Mat_Create(matfile_name.c_str(), nullptr);
- if (matfd)
- {
- writeMat(matfd, lname);
- Mat_Close(matfd);
- }
-}
-
-/* This save the results as matrices with given prefix and with index
- appended. If there is only one matrix, the index is not appended. */
-
-void
-SimResults::writeMat(mat_t *fd, const std::string &lname) const
-{
- for (int i = 0; i < getNumSets(); i++)
- {
- std::string tmp = lname + "_data";
- if (getNumSets() > 1)
- tmp += std::to_string(i+1);
- data[i].writeMat(fd, tmp);
- }
-}
-
-void
-SimResultsStats::simulate(int num_sim, const DecisionRule &dr,
- const Vector &start,
- const TwoDMatrix &vcov, Journal &journal)
-{
- SimResults::simulate(num_sim, dr, start, vcov, journal);
- {
- JournalRecordPair paa(journal);
- paa << "Calculating means from the simulations." << endrec;
- calcMean();
- }
- {
- JournalRecordPair paa(journal);
- paa << "Calculating covariances from the simulations." << endrec;
- calcVcov();
- }
-}
-
-/* Here we do not save the data itself, we save only mean and vcov. */
-void
-SimResultsStats::writeMat(mat_t *fd, const std::string &lname) const
-{
- ConstTwoDMatrix(num_y, 1, mean).writeMat(fd, lname + "_mean");;
- vcov.writeMat(fd, lname + "_vcov");
-}
-
-void
-SimResultsStats::calcMean()
-{
- mean.zeros();
- if (data.size()*num_per > 0)
- {
- double mult = 1.0/data.size()/num_per;
- for (const auto &i : data)
- {
- for (int j = 0; j < num_per; j++)
- {
- ConstVector col{i.getCol(j)};
- mean.add(mult, col);
- }
- }
- }
-}
-
-void
-SimResultsStats::calcVcov()
-{
- if (data.size()*num_per > 1)
- {
- vcov.zeros();
- double mult = 1.0/(data.size()*num_per - 1);
- for (const auto &d : data)
- for (int j = 0; j < num_per; j++)
- for (int m = 0; m < num_y; m++)
- for (int n = m; n < num_y; n++)
- {
- double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
- vcov.get(m, n) += mult*s;
- if (m != n)
- vcov.get(n, m) += mult*s;
- }
- }
- else
- vcov.infs();
-}
-
-void
-SimResultsDynamicStats::simulate(int num_sim, const DecisionRule &dr,
- const Vector &start,
- const TwoDMatrix &vcov, Journal &journal)
-{
- SimResults::simulate(num_sim, dr, start, vcov, journal);
- {
- JournalRecordPair paa(journal);
- paa << "Calculating means of the conditional simulations." << endrec;
- calcMean();
- }
- {
- JournalRecordPair paa(journal);
- paa << "Calculating variances of the conditional simulations." << endrec;
- calcVariance();
- }
-}
-
-void
-SimResultsDynamicStats::writeMat(mat_t *fd, const std::string &lname) const
-{
- mean.writeMat(fd, lname + "_cond_mean");
- variance.writeMat(fd, lname + "_cond_variance");
-}
-
-void
-SimResultsDynamicStats::calcMean()
-{
- mean.zeros();
- if (data.size() > 0)
- {
- double mult = 1.0/data.size();
- for (int j = 0; j < num_per; j++)
- {
- Vector meanj{mean.getCol(j)};
- for (const auto &i : data)
- {
- ConstVector col{i.getCol(j)};
- meanj.add(mult, col);
- }
- }
- }
-}
-
-void
-SimResultsDynamicStats::calcVariance()
-{
- if (data.size() > 1)
- {
- variance.zeros();
- double mult = 1.0/(data.size()-1);
- for (int j = 0; j < num_per; j++)
- {
- ConstVector meanj{mean.getCol(j)};
- Vector varj{variance.getCol(j)};
- for (const auto &i : data)
- {
- Vector col{i.getCol(j)};
- col.add(-1.0, meanj);
- for (int k = 0; k < col.length(); k++)
- col[k] = col[k]*col[k];
- varj.add(mult, col);
- }
- }
- }
- else
- variance.infs();
-}
-
-void
-SimResultsIRF::simulate(const DecisionRule &dr, Journal &journal)
-{
- JournalRecordPair paa(journal);
- paa << "Performing " << control.getNumSets() << " IRF simulations for "
- << num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
- simulate(dr);
- int thrown = control.getNumSets() - data.size();
- if (thrown > 0)
- {
- JournalRecord rec(journal);
- rec << "I had to throw " << thrown
- << " simulations away due to Nan or Inf" << endrec;
- }
- calcMeans();
- calcVariances();
-}
-
-void
-SimResultsIRF::simulate(const DecisionRule &dr)
-{
- sthread::detach_thread_group gr;
- for (int idata = 0; idata < control.getNumSets(); idata++)
- gr.insert(std::make_unique(*this, dr, DecisionRule::emethod::horner,
- num_per, idata, ishock, imp));
- gr.run();
-}
-
-void
-SimResultsIRF::calcMeans()
-{
- means.zeros();
- if (data.size() > 0)
- {
- for (const auto &i : data)
- means.add(1.0, i);
- means.mult(1.0/data.size());
- }
-}
-
-void
-SimResultsIRF::calcVariances()
-{
- if (data.size() > 1)
- {
- variances.zeros();
- for (const auto &i : data)
- {
- TwoDMatrix d(i);
- d.add(-1.0, means);
- for (int j = 0; j < d.nrows(); j++)
- for (int k = 0; k < d.ncols(); k++)
- variances.get(j, k) += d.get(j, k)*d.get(j, k);
- d.mult(1.0/(data.size()-1));
- }
- }
- else
- variances.infs();
-}
-
-void
-SimResultsIRF::writeMat(mat_t *fd, const std::string &lname) const
-{
- means.writeMat(fd, lname + "_mean");
- variances.writeMat(fd, lname + "_var");
-}
-
-void
-RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &v, Journal &journal)
-{
- JournalRecordPair paa(journal);
- paa << "Performing " << num_sim << " real-time stochastic simulations for "
- << num_per << " periods" << endrec;
- simulate(num_sim, dr, start, v);
- mean = nc.getMean();
- mean.add(1.0, dr.getSteady());
- nc.getVariance(vcov);
- if (thrown_periods > 0)
- {
- JournalRecord rec(journal);
- rec << "I had to throw " << thrown_periods << " periods away due to Nan or Inf" << endrec;
- JournalRecord rec1(journal);
- rec1 << "This affected " << incomplete_simulations << " out of "
- << num_sim << " simulations" << endrec;
- }
-}
-
-void
-RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov)
-{
- std::vector rsrs;
- rsrs.reserve(num_sim);
-
- sthread::detach_thread_group gr;
- for (int i = 0; i < num_sim; i++)
- {
- RandomShockRealization sr(vcov, seed_generator::get_new_seed());
- rsrs.push_back(sr);
- gr.insert(std::make_unique(*this, dr, DecisionRule::emethod::horner,
- num_per, start, rsrs.back()));
- }
- gr.run();
-}
-
-void
-RTSimResultsStats::writeMat(mat_t *fd, const std::string &lname)
-{
- ConstTwoDMatrix(nc.getDim(), 1, mean).writeMat(fd, lname + "_rt_mean");
- vcov.writeMat(fd, lname + "_rt_vcov");
-}
-
-IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
- const SimResults &control, std::vector ili,
- Journal &journal)
- : model(mod), irf_list_ind(std::move(ili))
-{
- int num_per = control.getNumPer();
- JournalRecordPair pa(journal);
- pa << "Calculating IRFs against control for " << static_cast(irf_list_ind.size()) << " shocks and for "
- << num_per << " periods" << endrec;
- const TwoDMatrix &vcov = mod.getVcov();
- for (int ishock : irf_list_ind)
- {
- double stderror = sqrt(vcov.get(ishock, ishock));
- irf_res.emplace_back(control, model.numeq(), num_per,
- ishock, stderror);
- irf_res.emplace_back(control, model.numeq(), num_per,
- ishock, -stderror);
- }
-
- for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
- {
- irf_res[2*ii].simulate(dr, journal);
- irf_res[2*ii+1].simulate(dr, journal);
- }
-}
-
-void
-IRFResults::writeMat(mat_t *fd, const std::string &prefix) const
-{
- for (unsigned int i = 0; i < irf_list_ind.size(); i++)
- {
- int ishock = irf_list_ind[i];
- auto shockname = model.getExogNames().getName(ishock);
- irf_res[2*i].writeMat(fd, prefix + "_irfp_" + shockname);
- irf_res[2*i+1].writeMat(fd, prefix + "_irfm_" + shockname);
- }
-}
-
-void
-SimulationWorker::operator()(std::mutex &mut)
-{
- ExplicitShockRealization esr(sr, np);
- TwoDMatrix m{dr.simulate(em, np, st, esr)};
- {
- std::unique_lock lk{mut};
- res.addDataSet(m, esr, st);
- }
-}
-
-/* Here we create a new instance of ExplicitShockRealization of the
- corresponding control, add the impulse, and simulate. */
-
-void
-SimulationIRFWorker::operator()(std::mutex &mut)
-{
- ExplicitShockRealization esr(res.control.getShocks(idata));
- esr.addToShock(ishock, 0, imp);
- TwoDMatrix m{dr.simulate(em, np, res.control.getStart(idata), esr)};
- m.add(-1.0, res.control.getData(idata));
- {
- std::unique_lock lk{mut};
- res.addDataSet(m, esr, res.control.getStart(idata));
- }
-}
-
-void
-RTSimulationWorker::operator()(std::mutex &mut)
-{
- NormalConj nc(res.nc.getDim());
- const PartitionY &ypart = dr.getYPart();
- int nu = dr.nexog();
- const Vector &ysteady = dr.getSteady();
-
- // initialize vectors and subvectors for simulation
- Vector dyu(ypart.nys()+nu);
- ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
- ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
- Vector dy(dyu, 0, ypart.nys());
- Vector u(dyu, ypart.nys(), nu);
- Vector y(nc.getDim());
- ConstVector ypred(y, ypart.nstat, ypart.nys());
-
- // simulate the first real-time period
- int ip = 0;
- dy = ystart_pred;
- dy.add(-1.0, ysteady_pred);
- sr.get(ip, u);
- dr.eval(em, y, dyu);
- if (ip >= res.num_burn)
- nc.update(y);
-
- // simulate other real-time periods
- while (y.isFinite() && ip < res.num_burn + res.num_per)
- {
- ip++;
- dy = ypred;
- sr.get(ip, u);
- dr.eval(em, y, dyu);
- if (ip >= res.num_burn)
- nc.update(y);
- }
- {
- std::unique_lock lk{mut};
- res.nc.update(nc);
- if (res.num_per-ip > 0)
- {
- res.incomplete_simulations++;
- res.thrown_periods += res.num_per-ip;
- }
- }
-}
-
-/* This calculates factorization FFᵀ=V in the Cholesky way. It does
- not work for semidefinite matrices. */
-
-void
-RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
-{
- factor = v;
- lapack_int rows = factor.nrows(), lda = factor.getLD();
- for (int i = 0; i < rows; i++)
- for (int j = i+1; j < rows; j++)
- factor.get(i, j) = 0.0;
- lapack_int info;
-
- dpotrf("L", &rows, factor.base(), &lda, &info);
- KORD_RAISE_IF(info != 0,
- "Info!=0 in RandomShockRealization::choleskyFactor");
-}
-
-/* This calculates FFᵀ=V factorization by symmetric Schur
- decomposition. It works for semidefinite matrices. */
-
-void
-RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
-{
- SymSchurDecomp(v).getFactor(factor);
-}
-
-void
-RandomShockRealization::get(int n, Vector &out)
-{
- KORD_RAISE_IF(out.length() != numShocks(),
- "Wrong length of out vector in RandomShockRealization::get");
- Vector d(out.length());
- for (int i = 0; i < d.length(); i++)
- d[i] = dis(mtwister);
- out.zeros();
- factor.multaVec(out, ConstVector(d));
-}
-
-ExplicitShockRealization::ExplicitShockRealization(ShockRealization &sr,
- int num_per)
- : shocks(sr.numShocks(), num_per)
-{
- for (int j = 0; j < num_per; j++)
- {
- Vector jcol{shocks.getCol(j)};
- sr.get(j, jcol);
- }
-}
-
-void
-ExplicitShockRealization::get(int n, Vector &out)
-{
- KORD_RAISE_IF(out.length() != numShocks(),
- "Wrong length of out vector in ExplicitShockRealization::get");
- int i = n % shocks.ncols();
- ConstVector icol{shocks.getCol(i)};
- out = icol;
-}
-
-void
-ExplicitShockRealization::addToShock(int ishock, int iper, double val)
-{
- KORD_RAISE_IF(ishock < 0 || ishock > numShocks(),
- "Wrong index of shock in ExplicitShockRealization::addToShock");
- int j = iper % shocks.ncols();
- shocks.get(ishock, j) += val;
-}
-
-void
-GenShockRealization::get(int n, Vector &out)
-{
- KORD_RAISE_IF(out.length() != numShocks(),
- "Wrong length of out vector in GenShockRealization::get");
- ExplicitShockRealization::get(n, out);
- Vector r(numShocks());
- RandomShockRealization::get(n, r);
- for (int j = 0; j < numShocks(); j++)
- if (!std::isfinite(out[j]))
- out[j] = r[j];
-}
diff --git a/mex/sources/libkorder/kord/decision_rule.hh b/mex/sources/libkorder/kord/decision_rule.hh
index 2736ab4cc..1e07e1d09 100644
--- a/mex/sources/libkorder/kord/decision_rule.hh
+++ b/mex/sources/libkorder/kord/decision_rule.hh
@@ -1,6 +1,6 @@
/*
* Copyright © 2004 Ondra Kamenik
- * Copyright © 2019-2021 Dynare Team
+ * Copyright © 2019-2023 Dynare Team
*
* This file is part of Dynare.
*
@@ -35,27 +35,13 @@
#ifndef DECISION_RULE_H
#define DECISION_RULE_H
-#include
-
#include "kord_exception.hh"
#include "korder.hh"
-#include "normal_conjugate.hh"
#include
#include
#include
-/* This is a general interface to a shock realizations. The interface has only
- one method returning the shock realizations at the given time. This method
- is not constant, since it may change a state of the object. */
-class ShockRealization
-{
-public:
- virtual ~ShockRealization() = default;
- virtual void get(int n, Vector &out) = 0;
- virtual int numShocks() const = 0;
-};
-
/* This class is an abstract interface to decision rule. Its main purpose is to
define a common interface for simulation of a decision rule. We need only a
simulate, evaluate, centralized clone and output method. */
@@ -65,10 +51,6 @@ public:
enum class emethod { horner, trad };
virtual ~DecisionRule() = default;
- // simulates the rule for a given realization of the shocks
- virtual TwoDMatrix simulate(emethod em, int np, const ConstVector &ystart,
- ShockRealization &sr) const = 0;
-
/* primitive evaluation (it takes a vector of state variables (predetermined,
both and shocks) and returns the next period variables. Both input and
output are in deviations from the rule's steady. */
@@ -79,9 +61,6 @@ public:
virtual void evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const = 0;
- // writes the decision rule to the MAT file
- virtual void writeMat(mat_t *fd, const std::string &prefix) const = 0;
-
/* returns a new copy of the decision rule, which is centralized about
provided fix-point */
virtual std::unique_ptr centralizedClone(const Vector &fixpoint) const = 0;
@@ -167,12 +146,9 @@ public:
{
return ysteady;
}
- TwoDMatrix simulate(emethod em, int np, const ConstVector &ystart,
- ShockRealization &sr) const override;
void evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const override;
std::unique_ptr centralizedClone(const Vector &fixpoint) const override;
- void writeMat(mat_t *fd, const std::string &prefix) const override;
int
nexog() const override
@@ -352,79 +328,6 @@ DecisionRuleImpl::centralize(const DecisionRuleImpl &dr)
}
}
-/* Here we evaluate repeatedly the polynomial storing results in the created
- matrix. For exogenous shocks, we use ShockRealization class, for
- predetermined variables, we use ‘ystart’ as the first state. The ‘ystart’
- vector is required to be all state variables ypart.ny(), although only the
- predetermined part of ‘ystart’ is used.
-
- We simulate in terms of Δy, this is, at the beginning the ‘ysteady’ is
- canceled from ‘ystart’, we simulate, and at the end ‘ysteady’ is added to
- all columns of the result. */
-
-template
-TwoDMatrix
-DecisionRuleImpl::simulate(emethod em, int np, const ConstVector &ystart,
- ShockRealization &sr) const
-{
- KORD_RAISE_IF(ysteady.length() != ystart.length(),
- "Start and steady lengths differ in DecisionRuleImpl::simulate");
- TwoDMatrix res(ypart.ny(), np);
-
- // initialize vectors and subvectors for simulation
- /* Here allocate the stack vector (Δy*,u), define the subvectors ‘dy’, and
- ‘u’, then we pickup predetermined parts of ‘ystart’ and ‘ysteady’. */
- Vector dyu(ypart.nys()+nu);
- ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
- ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
- Vector dy(dyu, 0, ypart.nys());
- Vector u(dyu, ypart.nys(), nu);
-
- // perform the first step of simulation
- /* We cancel ‘ysteady’ from ‘ystart’, get realization to ‘u’, and evaluate
- the polynomial. */
- dy = ystart_pred;
- dy.add(-1.0, ysteady_pred);
- sr.get(0, u);
- Vector out{res.getCol(0)};
- eval(em, out, dyu);
-
- // perform all other steps of simulations
- /* Also clear. If the result at some period is not finite, we pad the rest of
- the matrix with zeros. */
- int i = 1;
- while (i < np)
- {
- ConstVector ym{res.getCol(i-1)};
- ConstVector dym(ym, ypart.nstat, ypart.nys());
- dy = dym;
- sr.get(i, u);
- Vector out{res.getCol(i)};
- eval(em, out, dyu);
- if (!out.isFinite())
- {
- if (i+1 < np)
- {
- TwoDMatrix rest(res, i+1, np-i-1);
- rest.zeros();
- }
- break;
- }
- i++;
- }
-
- // add the steady state to columns of ‘res’
- /* Even clearer. We add the steady state to the numbers computed above and
- leave the padded columns to zero. */
- for (int j = 0; j < i; j++)
- {
- Vector col{res.getCol(j)};
- col.add(1.0, ysteady);
- }
-
- return res;
-}
-
/* This is one period evaluation of the decision rule. The simulation is a
sequence of repeated one period evaluations with a difference, that the
steady state (fix point) is cancelled and added once. Hence we have two
@@ -473,18 +376,6 @@ DecisionRuleImpl::eval(emethod em, Vector &out, const ConstVector &v) const
_Tpol::evalTrad(out, v);
}
-/* Write the decision rule and steady state to the MAT file. */
-
-template
-void
-DecisionRuleImpl::writeMat(mat_t *fd, const std::string &prefix) const
-{
- ctraits::Tpol::writeMat(fd, prefix);
- TwoDMatrix dum(ysteady.length(), 1);
- dum.getData() = ysteady;
- ConstTwoDMatrix(dum).writeMat(fd, prefix + "_ss");
-}
-
/* This is exactly the same as DecisionRuleImpl. The only
difference is that we have a conversion from UnfoldDecisionRule, which is
exactly DecisionRuleImpl. */
@@ -786,354 +677,4 @@ DRFixPoint::calcFixPoint(emethod em, Vector &out)
return converged;
}
-/* This is a basically a number of matrices of the same dimensions, which can
- be obtained as simulation results from a given decision rule and shock
- realizations. We also store the realizations of shocks and the starting
- point of each simulation. */
-
-class ExplicitShockRealization;
-class SimResults
-{
-protected:
- int num_y;
- int num_per;
- int num_burn;
- std::vector data;
- std::vector shocks;
- std::vector start;
-public:
- SimResults(int ny, int nper, int nburn = 0)
- : num_y(ny), num_per(nper), num_burn(nburn)
- {
- }
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov, Journal &journal);
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov);
- int
- getNumPer() const
- {
- return num_per;
- }
- int
- getNumBurn() const
- {
- return num_burn;
- }
- int
- getNumSets() const
- {
- return static_cast(data.size());
- }
- const TwoDMatrix &
- getData(int i) const
- {
- return data[i];
- }
- const ExplicitShockRealization &
- getShocks(int i) const
- {
- return shocks[i];
- }
- const ConstVector &
- getStart(int i) const
- {
- return start[i];
- }
-
- bool addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st);
- void writeMat(const std::string &base, const std::string &lname) const;
- void writeMat(mat_t *fd, const std::string &lname) const;
-};
-
-/* This does the same as SimResults plus it calculates means and covariances of
- the simulated data. */
-
-class SimResultsStats : public SimResults
-{
-protected:
- Vector mean;
- TwoDMatrix vcov;
-public:
- SimResultsStats(int ny, int nper, int nburn = 0)
- : SimResults(ny, nper, nburn), mean(ny), vcov(ny, ny)
- {
- }
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov, Journal &journal);
- void writeMat(mat_t *fd, const std::string &lname) const;
-protected:
- void calcMean();
- void calcVcov();
-};
-
-/* This does the similar thing as SimResultsStats but the statistics are not
- calculated over all periods but only within each period. Then we do not
- calculate covariances with periods but only variances. */
-
-class SimResultsDynamicStats : public SimResults
-{
-protected:
- TwoDMatrix mean;
- TwoDMatrix variance;
-public:
- SimResultsDynamicStats(int ny, int nper, int nburn = 0)
- : SimResults(ny, nper, nburn), mean(ny, nper), variance(ny, nper)
- {
- }
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov, Journal &journal);
- void writeMat(mat_t *fd, const std::string &lname) const;
-protected:
- void calcMean();
- void calcVariance();
-};
-
-/* This goes through control simulation results, and for each control it adds a
- given impulse to a given shock and runs a simulation. The control simulation
- is then cancelled and the result is stored. After that these results are
- averaged with variances calculated.
-
- The means and the variances are then written to the MAT file. */
-
-class SimulationIRFWorker;
-class SimResultsIRF : public SimResults
-{
- friend class SimulationIRFWorker;
-protected:
- const SimResults &control;
- int ishock;
- double imp;
- TwoDMatrix means;
- TwoDMatrix variances;
-public:
- SimResultsIRF(const SimResults &cntl, int ny, int nper, int i, double impulse)
- : SimResults(ny, nper, 0), control(cntl),
- ishock(i), imp(impulse),
- means(ny, nper), variances(ny, nper)
- {
- }
- void simulate(const DecisionRule &dr, Journal &journal);
- void simulate(const DecisionRule &dr);
- void writeMat(mat_t *fd, const std::string &lname) const;
-protected:
- void calcMeans();
- void calcVariances();
-};
-
-/* This simulates and gathers all statistics from the real time simulations. In
- the simulate() method, it runs RTSimulationWorker’s which accummulate
- information from their own estimates. The estimation is done by means of
- NormalConj class, which is a conjugate family of densities for normal
- distibutions. */
-
-class RTSimulationWorker;
-class RTSimResultsStats
-{
- friend class RTSimulationWorker;
-protected:
- Vector mean;
- TwoDMatrix vcov;
- int num_per;
- int num_burn;
- NormalConj nc;
- int incomplete_simulations;
- int thrown_periods;
-public:
- RTSimResultsStats(int ny, int nper, int nburn = 0)
- : mean(ny), vcov(ny, ny),
- num_per(nper), num_burn(nburn), nc(ny),
- incomplete_simulations(0), thrown_periods(0)
- {
- }
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov, Journal &journal);
- void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
- const TwoDMatrix &vcov);
- void writeMat(mat_t *fd, const std::string &lname);
-};
-
-/* For each shock, this simulates plus and minus impulse. The class maintains a
- vector of simulation results, each gets a particular shock and sign
- (positive/negative). The results of type SimResultsIRF are stored in a
- vector so that even ones are positive, odd ones are negative.
-
- The constructor takes a reference to the control simulations, which must be
- finished before the constructor is called. The control simulations are
- passed to all SimResultsIRF’s.
-
- The constructor also takes the vector of indices of exogenous variables
- (‘ili’) for which the IRFs are generated. The list is kept (as
- ‘irf_list_ind’) for other methods. */
-
-class DynamicModel;
-class IRFResults
-{
- std::vector irf_res;
- const DynamicModel &model;
- std::vector irf_list_ind;
-public:
- IRFResults(const DynamicModel &mod, const DecisionRule &dr,
- const SimResults &control, std::vector ili,
- Journal &journal);
- void writeMat(mat_t *fd, const std::string &prefix) const;
-};
-
-/* This worker simulates the given decision rule and inserts the result to
- SimResults. */
-
-class SimulationWorker : public sthread::detach_thread
-{
-protected:
- SimResults &res;
- const DecisionRule &dr;
- DecisionRule::emethod em;
- int np;
- const Vector &st;
- ShockRealization &sr;
-public:
- SimulationWorker(SimResults &sim_res,
- const DecisionRule &dec_rule,
- DecisionRule::emethod emet, int num_per,
- const Vector &start, ShockRealization &shock_r)
- : res(sim_res), dr(dec_rule), em(emet), np(num_per), st(start), sr(shock_r)
- {
- }
- void operator()(std::mutex &mut) override;
-};
-
-/* This worker simulates a given impulse ‘imp’ to a given shock ‘ishock’ based
- on a given control simulation with index ‘idata’. The control simulations
- are contained in SimResultsIRF which is passed to the constructor. */
-
-class SimulationIRFWorker : public sthread::detach_thread
-{
- SimResultsIRF &res;
- const DecisionRule &dr;
- DecisionRule::emethod em;
- int np;
- int idata;
- int ishock;
- double imp;
-public:
- SimulationIRFWorker(SimResultsIRF &sim_res,
- const DecisionRule &dec_rule,
- DecisionRule::emethod emet, int num_per,
- int id, int ishck, double impulse)
- : res(sim_res), dr(dec_rule), em(emet), np(num_per),
- idata(id), ishock(ishck), imp(impulse)
- {
- }
- void operator()(std::mutex &mut) override;
-};
-
-/* This class does the real time simulation job for RTSimResultsStats. It
- simulates the model period by period. It accummulates the information in
- ‘RTSimResultsStats::nc’. If NaN or Inf is observed, it ends the simulation
- and adds to the ‘thrown_periods’ of RTSimResultsStats. */
-
-class RTSimulationWorker : public sthread::detach_thread
-{
-protected:
- RTSimResultsStats &res;
- const DecisionRule &dr;
- DecisionRule::emethod em;
- int np;
- const Vector &ystart;
- ShockRealization &sr;
-public:
- RTSimulationWorker(RTSimResultsStats &sim_res,
- const DecisionRule &dec_rule,
- DecisionRule::emethod emet, int num_per,
- const Vector &start, ShockRealization &shock_r)
- : res(sim_res), dr(dec_rule), em(emet), np(num_per), ystart(start), sr(shock_r)
- {
- }
- void operator()(std::mutex &mut) override;
-};
-
-/* This class generates draws from Gaussian distribution with zero mean and the
- given variance-covariance matrix. It stores the factor of vcov V matrix,
- yielding FFᵀ = V. */
-
-class RandomShockRealization : virtual public ShockRealization
-{
-protected:
- std::mt19937 mtwister;
- std::normal_distribution<> dis;
- TwoDMatrix factor;
-public:
- RandomShockRealization(const ConstTwoDMatrix &v, decltype(mtwister)::result_type iseed)
- : mtwister(iseed), factor(v.nrows(), v.nrows())
- {
- schurFactor(v);
- }
- void get(int n, Vector &out) override;
- int
- numShocks() const override
- {
- return factor.nrows();
- }
-protected:
- void choleskyFactor(const ConstTwoDMatrix &v);
- void schurFactor(const ConstTwoDMatrix &v);
-};
-
-/* This is just a matrix of finite numbers. It can be constructed from any
- ShockRealization with a given number of periods. */
-
-class ExplicitShockRealization : virtual public ShockRealization
-{
- TwoDMatrix shocks;
-public:
- explicit ExplicitShockRealization(const ConstTwoDMatrix &sh)
- : shocks(sh)
- {
- }
- ExplicitShockRealization(ShockRealization &sr, int num_per);
- void get(int n, Vector &out) override;
- int
- numShocks() const override
- {
- return shocks.nrows();
- }
- const TwoDMatrix &
- getShocks() const
- {
- return shocks;
- }
- void addToShock(int ishock, int iper, double val);
- void
- print() const
- {
- shocks.print();
- }
-};
-
-/* This represents a user given shock realization. The first matrix of the
- constructor is a covariance matrix of shocks, the second matrix is a
- rectangular matrix, where columns correspond to periods, rows to shocks. If
- an element of the matrix is NaN or ±∞, then the random shock is taken
- instead of that element.
-
- In this way it is a generalization of both RandomShockRealization and
- ExplicitShockRealization. */
-
-class GenShockRealization : public RandomShockRealization, public ExplicitShockRealization
-{
-public:
- GenShockRealization(const ConstTwoDMatrix &v, const ConstTwoDMatrix &sh, int seed)
- : RandomShockRealization(v, seed), ExplicitShockRealization(sh)
- {
- KORD_RAISE_IF(sh.nrows() != v.nrows() || v.nrows() != v.ncols(),
- "Wrong dimension of input matrix in GenShockRealization constructor");
- }
- void get(int n, Vector &out) override;
- int
- numShocks() const override
- {
- return RandomShockRealization::numShocks();
- }
-};
-
#endif
diff --git a/mex/sources/libkorder/kord/dynamic_model.cc b/mex/sources/libkorder/kord/dynamic_model.cc
index 99b5d525c..f39214447 100644
--- a/mex/sources/libkorder/kord/dynamic_model.cc
+++ b/mex/sources/libkorder/kord/dynamic_model.cc
@@ -29,44 +29,3 @@ NameList::print() const
for (int i = 0; i < getNum(); i++)
std::cout << getName(i) << '\n';
}
-
-void
-NameList::writeMat(mat_t *fd, const std::string &vname) const
-{
- int maxlen = 0;
- for (int i = 0; i < getNum(); i++)
- maxlen = std::max(maxlen, static_cast(getName(i).size()));
-
- if (maxlen == 0)
- return;
-
- auto m = std::make_unique(getNum()*maxlen);
-
- for (int i = 0; i < getNum(); i++)
- for (int j = 0; j < maxlen; j++)
- if (j < static_cast(getName(i).size()))
- m[j*getNum()+i] = getName(i)[j];
- else
- m[j*getNum()+i] = ' ';
-
- size_t dims[2];
- dims[0] = getNum();
- dims[1] = maxlen;
-
- matvar_t *v = Mat_VarCreate(vname.c_str(), MAT_C_CHAR, MAT_T_UINT8, 2, dims, m.get(), 0);
-
- Mat_VarWrite(fd, v, MAT_COMPRESSION_NONE);
-
- Mat_VarFree(v);
-}
-
-void
-NameList::writeMatIndices(mat_t *fd, const std::string &prefix) const
-{
- TwoDMatrix aux(1, 1);
- for (int i = 0; i < getNum(); i++)
- {
- aux.get(0, 0) = i+1;
- aux.writeMat(fd, prefix + "_i_" + getName(i));
- }
-}
diff --git a/mex/sources/libkorder/kord/dynamic_model.hh b/mex/sources/libkorder/kord/dynamic_model.hh
index 9ceb02cac..2171027d5 100644
--- a/mex/sources/libkorder/kord/dynamic_model.hh
+++ b/mex/sources/libkorder/kord/dynamic_model.hh
@@ -48,8 +48,6 @@ public:
virtual int getNum() const = 0;
virtual const std::string &getName(int i) const = 0;
void print() const;
- void writeMat(mat_t *fd, const std::string &vname) const;
- void writeMatIndices(mat_t *fd, const std::string &prefix) const;
};
/* This is the interface to an information on a generic DSGE model. It is
diff --git a/mex/sources/libkorder/kord/global_check.cc b/mex/sources/libkorder/kord/global_check.cc
deleted file mode 100644
index 0992c8ad1..000000000
--- a/mex/sources/libkorder/kord/global_check.cc
+++ /dev/null
@@ -1,384 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "SymSchurDecomp.hh"
-
-#include "global_check.hh"
-#include "seed_generator.hh"
-
-#include "smolyak.hh"
-#include "product.hh"
-#include "quasi_mcarlo.hh"
-
-#include
-#include
-
-/* Here we just set a reference to the approximation, and create a new
- DynamicModel. */
-
-ResidFunction::ResidFunction(const Approximation &app)
- : VectorFunction(app.getModel().nexog(), app.getModel().numeq()), approx(app),
- model(app.getModel().clone())
-{
-}
-
-ResidFunction::ResidFunction(const ResidFunction &rf)
- : VectorFunction(rf), approx(rf.approx), model(rf.model->clone())
-{
- if (rf.yplus)
- yplus = std::make_unique(*(rf.yplus));
- if (rf.ystar)
- ystar = std::make_unique(*(rf.ystar));
- if (rf.u)
- u = std::make_unique(*(rf.u));
- if (rf.hss)
- hss = std::make_unique(*(rf.hss));
-}
-
-/* This sets y* and u. We have to create ‘ystar’, ‘u’, ‘yplus’ and ‘hss’. */
-
-void
-ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
-{
- ystar = std::make_unique(ys);
- u = std::make_unique(xx);
- yplus = std::make_unique(model->numeq());
- approx.getFoldDecisionRule().evaluate(DecisionRule::emethod::horner,
- *yplus, *ystar, *u);
-
- // make a tensor polynomial of in-place subtensors from decision rule
- /* Note that the non-const polynomial will be used for a construction of
- ‘hss’ and will be used in a const context. So this const cast is safe.
-
- Note, that there is always a folded decision rule in Approximation. */
- const FoldDecisionRule &dr = approx.getFoldDecisionRule();
- FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
- const_cast(dr));
-
- // make ‘ytmp_star’ be a difference of ‘yplus’ from steady
- Vector ytmp_star(ConstVector(*yplus, model->nstat(), model->npred()+model->nboth()));
- ConstVector ysteady_star(dr.getSteady(), model->nstat(),
- model->npred()+model->nboth());
- ytmp_star.add(-1.0, ysteady_star);
-
- // make ‘hss’ and add steady to it
- /* Here is the const context of ‘dr_ss’. */
- hss = std::make_unique(dr_ss, ytmp_star);
- ConstVector ysteady_ss(dr.getSteady(), model->nstat()+model->npred(),
- model->nboth()+model->nforw());
- if (hss->check(Symmetry{0}))
- hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
- else
- {
- auto ten = std::make_unique(hss->nrows(), hss->nvars(), 0);
- ten->getData() = ysteady_ss;
- hss->insert(std::move(ten));
- }
-}
-
-/* Here we evaluate the residual F(y*,u,u′). We have to evaluate ‘hss’ for
- u′=point and then we evaluate the system f. */
-
-void
-ResidFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
- KORD_RAISE_IF(point.length() != hss->nvars(),
- "Wrong dimension of input vector in ResidFunction::eval");
- KORD_RAISE_IF(out.length() != model->numeq(),
- "Wrong dimension of output vector in ResidFunction::eval");
- Vector yss(hss->nrows());
- hss->evalHorner(yss, point);
- model->evaluateSystem(out, *ystar, *yplus, yss, *u);
-}
-
-/* This checks the 𝔼[F(y*,u,u′)] for a given y* and u by integrating with a
- given quadrature. Note that the input ‘ys’ is y* not whole y. */
-
-void
-GlobalChecker::check(const Quadrature &quad, int level,
- const ConstVector &ys, const ConstVector &x, Vector &out)
-{
- for (int ifunc = 0; ifunc < vfs.getNum(); ifunc++)
- dynamic_cast(vfs.getFunc(ifunc)).setYU(ys, x);
- quad.integrate(vfs, level, out);
-}
-
-/* This method is a bulk version of GlobalChecker::check() vector code. It
- decides between Smolyak and product quadrature according to ‘max_evals’
- constraint.
-
- Note that ‘y’ can be either full (all endogenous variables including static
- and forward looking), or just y* (state variables). The method is able to
- recognize it. */
-
-void
-GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
- const ConstTwoDMatrix &x, TwoDMatrix &out)
-{
- JournalRecordPair pa(journal);
- pa << "Checking approximation error for " << y.ncols()
- << " states with at most " << max_evals << " evaluations" << endrec;
-
- // Decide about which type of quadrature
- GaussHermite gh;
-
- SmolyakQuadrature dummy_sq(model.nexog(), 1, gh);
- int smol_evals;
- int smol_level;
- dummy_sq.designLevelForEvals(max_evals, smol_level, smol_evals);
-
- ProductQuadrature dummy_pq(model.nexog(), gh);
- int prod_evals;
- int prod_level;
- dummy_pq.designLevelForEvals(max_evals, prod_level, prod_evals);
-
- bool take_smolyak = (smol_evals < prod_evals) && (smol_level >= prod_level-1);
-
- std::unique_ptr quad;
- int lev;
-
- // Create the quadrature and report the decision
- if (take_smolyak)
- {
- quad = std::make_unique(model.nexog(), smol_level, gh);
- lev = smol_level;
- JournalRecord rec(journal);
- rec << "Selected Smolyak (level,evals)=(" << smol_level << ","
- << smol_evals << ") over product (" << prod_level << ","
- << prod_evals << ")" << endrec;
- }
- else
- {
- quad = std::make_unique(model.nexog(), gh);
- lev = prod_level;
- JournalRecord rec(journal);
- rec << "Selected product (level,evals)=(" << prod_level << ","
- << prod_evals << ") over Smolyak (" << smol_level << ","
- << smol_evals << ")" << endrec;
- }
-
- // Check all columns of ‘y’ and ‘x’
- int first_row = (y.nrows() == model.numeq()) ? model.nstat() : 0;
- ConstTwoDMatrix ysmat(y, first_row, 0, model.npred()+model.nboth(), y.ncols());
- for (int j = 0; j < y.ncols(); j++)
- {
- ConstVector yj{ysmat.getCol(j)};
- ConstVector xj{x.getCol(j)};
- Vector outj{out.getCol(j)};
- check(*quad, lev, yj, xj, outj);
- }
-}
-
-/* This method checks an error of the approximation by evaluating residual
- 𝔼[F(y*,u,u′) | y*,u] for y* equal to the steady state, and changing u. We go
- through all elements of u and vary them from −mult·σ to mult·σ in ‘m’
- steps. */
-
-void
-GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const std::string &prefix,
- int m, double mult, int max_evals)
-{
- JournalRecordPair pa(journal);
- pa << "Calculating errors along shocks +/- "
- << mult << " std errors, granularity " << m << endrec;
-
- // Setup ‘y_mat’ of steady states for checking
- TwoDMatrix y_mat(model.numeq(), 2*m*model.nexog()+1);
- for (int j = 0; j < 2*m*model.nexog()+1; j++)
- {
- Vector yj{y_mat.getCol(j)};
- yj = model.getSteady();
- }
-
- // Setup ‘exo_mat’ for checking
- TwoDMatrix exo_mat(model.nexog(), 2*m*model.nexog()+1);
- exo_mat.zeros();
- for (int ishock = 0; ishock < model.nexog(); ishock++)
- {
- double max_sigma = sqrt(model.getVcov().get(ishock, ishock));
- for (int j = 0; j < 2*m; j++)
- {
- int jmult = (j < m) ? j-m : j-m+1;
- exo_mat.get(ishock, 1+2*m*ishock+j) = mult*jmult*max_sigma/m;
- }
- }
-
- TwoDMatrix errors(model.numeq(), 2*m*model.nexog()+1);
- check(max_evals, y_mat, exo_mat, errors);
-
- // Report errors along shock and save them
- TwoDMatrix res(model.nexog(), 2*m+1);
- JournalRecord rec(journal);
- rec << "Shock value error" << endrec;
- ConstVector err0{errors.getCol(0)};
- for (int ishock = 0; ishock < model.nexog(); ishock++)
- {
- TwoDMatrix err_out(model.numeq(), 2*m+1);
- for (int j = 0; j < 2*m+1; j++)
- {
- int jj;
- Vector error{err_out.getCol(j)};
- if (j != m)
- {
- if (j < m)
- jj = 1 + 2*m*ishock+j;
- else
- jj = 1 + 2*m*ishock+j-1;
- ConstVector coljj{errors.getCol(jj)};
- error = coljj;
- }
- else
- {
- jj = 0;
- error = err0;
- }
- JournalRecord rec1(journal);
- std::string shockname{model.getExogNames().getName(ishock)};
- shockname.resize(8, ' ');
- rec1 << shockname << ' ' << exo_mat.get(ishock, jj)
- << "\t" << error.getMax() << endrec;
- }
- err_out.writeMat(fd, prefix + "_shock_" + model.getExogNames().getName(ishock) + "_errors");
- }
-}
-
-/* This method checks errors on ellipse of endogenous states (predetermined
- variables). The ellipse is shaped according to covariance matrix of
- endogenous variables based on the first order approximation and scaled by
- ‘mult’. The points on the ellipse are chosen as polar images of the low
- discrepancy grid in a cube.
-
- The method works as follows. First we calculate symmetric Schur factor of
- covariance matrix of the states. Second we generate low discrepancy points
- on the unit sphere. Third we transform the sphere with the
- variance-covariance matrix factor and multiplier ‘mult’ and initialize
- matrix of uₜ to zeros. Fourth we run the check() method and save the
- results. */
-
-void
-GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
- int m, double mult, int max_evals)
-{
- JournalRecordPair pa(journal);
- pa << "Calculating errors at " << m
- << " ellipse points scaled by " << mult << endrec;
-
- // Make factor of covariance of variables
- /* Here we set ‘ycovfac’ to the symmetric Schur decomposition factor of a
- submatrix of covariances of all endogenous variables. The submatrix
- corresponds to state variables (predetermined plus both). */
- TwoDMatrix ycov{approx.calcYCov()};
- TwoDMatrix ycovpred(const_cast(ycov), model.nstat(), model.nstat(),
- model.npred()+model.nboth(), model.npred()+model.nboth());
- SymSchurDecomp ssd(ycovpred);
- ssd.correctDefinitness(1.e-05);
- TwoDMatrix ycovfac(ycovpred.nrows(), ycovpred.ncols());
- KORD_RAISE_IF(!ssd.isPositiveSemidefinite(),
- "Covariance matrix of the states not positive \
- semidefinite in GlobalChecker::checkOnEllipseAndSave");
- ssd.getFactor(ycovfac);
-
- // Put low discrepancy sphere points to ‘ymat’
- /* Here we first calculate dimension ‘d’ of the sphere, which is a number of
- state variables minus one. We go through the ‘d’-dimensional cube [0,1]ᵈ
- by QMCarloCubeQuadrature and make a polar transformation to the sphere.
- The polar transformation fⁱ can be written recursively w.r.t. the
- dimension i as:
-
- f⁰() = [1]
-
- ⎡cos(2πxᵢ)·fⁱ⁻¹(x₁,…,xᵢ₋₁)⎤
- fⁱ(x₁,…,xᵢ) = ⎣ sin(2πxᵢ) ⎦
- */
- int d = model.npred()+model.nboth()-1;
- TwoDMatrix ymat(model.npred()+model.nboth(), (d == 0) ? 2 : m);
- if (d == 0)
- {
- ymat.get(0, 0) = 1;
- ymat.get(0, 1) = -1;
- }
- else
- {
- int icol = 0;
- ReversePerScheme ps;
- QMCarloCubeQuadrature qmc(d, m, ps);
- qmcpit beg = qmc.start(m);
- qmcpit end = qmc.end(m);
- for (qmcpit run = beg; run != end; ++run, icol++)
- {
- Vector ycol{ymat.getCol(icol)};
- Vector x(run.point());
- x.mult(2*M_PI);
- ycol[0] = 1;
- for (int i = 0; i < d; i++)
- {
- Vector subsphere(ycol, 0, i+1);
- subsphere.mult(cos(x[i]));
- ycol[i+1] = sin(x[i]);
- }
- }
- }
-
- // Transform sphere ‘ymat’ and prepare ‘umat’ for checking
- /* Here we multiply the sphere points in ‘ymat’ with the Cholesky factor to
- obtain the ellipse, scale the ellipse by the given ‘mult’, and initialize
- matrix of shocks ‘umat’ to zero. */
- TwoDMatrix umat(model.nexog(), ymat.ncols());
- umat.zeros();
- ymat.mult(mult);
- ymat.multLeft(ycovfac);
- ConstVector ys(model.getSteady(), model.nstat(),
- model.npred()+model.nboth());
- for (int icol = 0; icol < ymat.ncols(); icol++)
- {
- Vector ycol{ymat.getCol(icol)};
- ycol.add(1.0, ys);
- }
-
- // Check on ellipse and save
- /* Here we check the points and save the results to MAT-4 file. */
- TwoDMatrix out(model.numeq(), ymat.ncols());
- check(max_evals, ymat, umat, out);
-
- ymat.writeMat(fd, prefix + "_ellipse_points");
- out.writeMat(fd, prefix + "_ellipse_errors");
-}
-
-/* Here we check the errors along a simulation. We simulate, then set ‘x’ to
- zeros, check and save results. */
-
-void
-GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
- int m, int max_evals)
-{
- JournalRecordPair pa(journal);
- pa << "Calculating errors at " << m
- << " simulated points" << endrec;
- RandomShockRealization sr(model.getVcov(), seed_generator::get_new_seed());
- TwoDMatrix y{approx.getFoldDecisionRule().simulate(DecisionRule::emethod::horner,
- m, model.getSteady(), sr)};
- TwoDMatrix x(model.nexog(), m);
- x.zeros();
- TwoDMatrix out(model.numeq(), m);
- check(max_evals, y, x, out);
-
- y.writeMat(fd, prefix + "_simul_points");
- out.writeMat(fd, prefix + "_simul_errors");
-}
diff --git a/mex/sources/libkorder/kord/global_check.hh b/mex/sources/libkorder/kord/global_check.hh
deleted file mode 100644
index e5dae5057..000000000
--- a/mex/sources/libkorder/kord/global_check.hh
+++ /dev/null
@@ -1,171 +0,0 @@
-/*
- * Copyright © 2005 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Global check
-
-/* The purpose of this file is to provide classes for checking error of
- approximation. If yₜ=g(y*ₜ₋₁,u) is an approximate solution, then we check
- for the error of residuals of the system equations. Let
- F(y*,u,u′)=f(g**(g*(y*,u′),u),g(y*,u),y*,u), then we calculate the integral:
-
- 𝔼ₜ[F(y*,u,u′)]
-
- which we want to be zero for all y* and u.
-
- There are a few possibilities how and where the integral is evaluated.
- Currently we offer the following ones:
-
- — Along shocks. The y* is set to the steady state, and u is set to zero but
- one element is going from minus through plus shocks in few steps. The user
- gives the scaling factor, for instance the interval [−3σ,3σ] (where σ is
- one standard error of the shock), and a number of steps. This is repeated
- for each shock (element of the u vector).
-
- — Along simulation. Some random simulation is run, and for each realization
- of y* and u along the path we evaluate the residual.
-
- — On ellipse. Let V=AAᵀ be a covariance matrix of the predetermined
- variables y* based on linear approximation, then we calculate integral for
- points on the ellipse { Ax | ‖x‖₂=1 }. The points are selected by means of
- low discrepancy method and polar transformation. The shock u are zeros.
-
- — Unconditional distribution.
-*/
-
-#ifndef GLOBAL_CHECK_H
-#define GLOBAL_CHECK_H
-
-#include
-
-#include
-
-#include "vector_function.hh"
-#include "quadrature.hh"
-
-#include "dynamic_model.hh"
-#include "journal.hh"
-#include "approximation.hh"
-
-/* This is a class for implementing the VectorFunction interface evaluating the
- residual of equations, this is F(y*,u,u′)=f(g**(g*(y*,u),u′),y*,u) is
- written as a function of u′.
-
- When the object is constructed, one has to specify (y*,u), this is done by
- the setYU() method. The object has basically two states. One is after
- construction and before the call to setYU(). The second is after the call to
- setYU(). We distinguish between the two states, an object in the second
- state contains ‘yplus’, ‘ystar’, ‘u’, and ‘hss’.
-
- The vector ‘yplus’ is g*(y*,u). ‘ystar’ is y*, and polynomial ‘hss’ is
- partially evaluated g**(yplus, u).
-
- The pointer to DynamicModel is important, since the DynamicModel evaluates
- the function f. When copying the object, we have to make also a copy of
- DynamicModel. */
-
-class ResidFunction : public VectorFunction
-{
-protected:
- const Approximation ≈
- std::unique_ptr model;
- std::unique_ptr yplus, ystar, u;
- std::unique_ptr hss;
-public:
- ResidFunction(const Approximation &app);
- ResidFunction(const ResidFunction &rf);
-
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
- void setYU(const ConstVector &ys, const ConstVector &xx);
-};
-
-/* This is a ResidFunction wrapped with GaussConverterFunction. */
-
-class GResidFunction : public GaussConverterFunction
-{
-public:
- GResidFunction(const Approximation &app)
- : GaussConverterFunction(std::make_unique(app), app.getModel().getVcov())
- {
- }
- std::unique_ptr
- clone() const override
- {
- return std::make_unique(*this);
- }
- void
- setYU(const ConstVector &ys, const ConstVector &xx)
- {
- dynamic_cast(func)->setYU(ys, xx);
- }
-};
-
-/* This is a class encapsulating checking algorithms. Its core routine is
- check(), which calculates integral 𝔼[F(y*,u,u′) | y*,u] for given
- realizations of y* and u. The both are given in matrices. The methods
- checking along shocks, on ellipse and anlong a simulation path, just fill
- the matrices and call the core check().
-
- The method checkUnconditionalAndSave() evaluates unconditional 𝔼[F(y,u,u′)].
-
- The object also maintains a set of GResidFunction functions ‘vfs’ in order
- to save (possibly expensive) copying of DynamicModel’s. */
-
-class GlobalChecker
-{
- const Approximation ≈
- const DynamicModel &model;
- Journal &journal;
- GResidFunction rf;
- VectorFunctionSet vfs;
-public:
- GlobalChecker(const Approximation &app, int n, Journal &jr)
- : approx(app), model(approx.getModel()), journal(jr),
- rf(approx), vfs(rf, n)
- {
- }
- void check(int max_evals, const ConstTwoDMatrix &y,
- const ConstTwoDMatrix &x, TwoDMatrix &out);
- void checkAlongShocksAndSave(mat_t *fd, const std::string &prefix,
- int m, double mult, int max_evals);
- void checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
- int m, double mult, int max_evals);
- void checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
- int m, int max_evals);
- void checkUnconditionalAndSave(mat_t *fd, const std::string &prefix,
- int m, int max_evals);
-protected:
- void check(const Quadrature &quad, int level,
- const ConstVector &y, const ConstVector &x, Vector &out);
-};
-
-/* Signalled resid function. Not implemented yet. todo: */
-
-class ResidFunctionSig : public ResidFunction
-{
-public:
- ResidFunctionSig(const Approximation &app, const Vector &ys, const Vector &xx);
-};
-
-#endif
diff --git a/mex/sources/libkorder/kord/normal_conjugate.cc b/mex/sources/libkorder/kord/normal_conjugate.cc
deleted file mode 100644
index 3c5fd8b13..000000000
--- a/mex/sources/libkorder/kord/normal_conjugate.cc
+++ /dev/null
@@ -1,123 +0,0 @@
-/*
- * Copyright © 2007 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "normal_conjugate.hh"
-#include "kord_exception.hh"
-
-// NormalConj diffuse prior constructor
-NormalConj::NormalConj(int d)
- : mu(d), kappa(0), nu(-1), lambda(d, d)
-{
- mu.zeros();
- lambda.zeros();
-}
-
-// NormalConj data update constructor
-NormalConj::NormalConj(const ConstTwoDMatrix &ydata)
- : mu(ydata.nrows()), kappa(ydata.ncols()), nu(ydata.ncols()-1),
- lambda(ydata.nrows(), ydata.nrows())
-{
- mu.zeros();
- for (int i = 0; i < ydata.ncols(); i++)
- mu.add(1.0/ydata.ncols(), ydata.getCol(i));
-
- lambda.zeros();
- for (int i = 0; i < ydata.ncols(); i++)
- {
- Vector diff{ydata.getCol(i)};
- diff.add(-1, mu);
- lambda.addOuter(diff);
- }
-}
-
-// NormalConj::update() one observation code
-/* The method performs the following:
-
- κ₀ 1
- μ₁ = ──── μ₀ + ──── y
- κ₀+1 κ₀+1
-
- κ₁ = κ₀ + 1
-
- ν₁ = ν₀ + 1
-
- κ₀
- Λ₁ = Λ₀ + ──── (y − μ₀)(y − μ₀)ᵀ
- κ₀+1
-*/
-void
-NormalConj::update(const ConstVector &y)
-{
- KORD_RAISE_IF(y.length() != mu.length(),
- "Wrong length of a vector in NormalConj::update");
-
- mu.mult(kappa/(1.0+kappa));
- mu.add(1.0/(1.0+kappa), y);
-
- Vector diff(y);
- diff.add(-1, mu);
- lambda.addOuter(diff, kappa/(1.0+kappa));
-
- kappa++;
- nu++;
-}
-
-// NormalConj::update() multiple observations code
-/* The method evaluates the formula in the header file. */
-void
-NormalConj::update(const ConstTwoDMatrix &ydata)
-{
- NormalConj nc(ydata);
- update(nc);
-}
-
-// NormalConj::update() with NormalConj code
-void
-NormalConj::update(const NormalConj &nc)
-{
- double wold = static_cast(kappa)/(kappa+nc.kappa);
- double wnew = 1-wold;
-
- mu.mult(wold);
- mu.add(wnew, nc.mu);
-
- Vector diff(nc.mu);
- diff.add(-1, mu);
- lambda.add(1.0, nc.lambda);
- lambda.addOuter(diff);
-
- kappa = kappa + nc.kappa;
- nu = nu + nc.kappa;
-}
-
-/* This returns 1/(ν−d−1)·Λ, which is the mean of the variance in the posterior
- distribution. If the number of degrees of freedom is less than d, then NaNs
- are returned. */
-void
-NormalConj::getVariance(TwoDMatrix &v) const
-{
- if (nu > getDim()+1)
- {
- v = const_cast(lambda);
- v.mult(1.0/(nu-getDim()-1));
- }
- else
- v.nans();
-}
diff --git a/mex/sources/libkorder/kord/normal_conjugate.hh b/mex/sources/libkorder/kord/normal_conjugate.hh
deleted file mode 100644
index 2f6268087..000000000
--- a/mex/sources/libkorder/kord/normal_conjugate.hh
+++ /dev/null
@@ -1,105 +0,0 @@
-/*
- * Copyright © 2007 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-// Conjugate family for normal distribution
-
-/* The main purpose here is to implement a class representing conjugate
- distributions for mean and variance of the normal distribution. The class
- has two main methods: the first one is to update itself with respect to one
- observation, the second one is to update itself with respect to anothe
- object of the class. In the both methods, the previous state of the class
- corresponds to the prior distribution, and the final state corresponds to
- the posterior distribution.
-
- The algebra can be found in Gelman, Carlin, Stern, Rubin (p.87). It goes as
- follows. Prior conjugate distribution takes the following form:
-
- Σ ↝ InvWishart_ν₀(Λ₀⁻¹)
- μ|Σ ↝ 𝒩(μ₀,Σ/κ₀)
-
- If the observations are y₁…yₙ, then the posterior distribution has the same
- form with the following parameters:
-
- κ₀ n
- μₙ = ──── μ₀ + ──── ȳ
- κ₀+n κ₀+n
-
- κₙ = κ₀ + n
-
- νₙ = ν₀ + n
-
- κ₀·n
- Λₙ = Λ₀ + S + ──── (ȳ − μ₀)(ȳ − μ₀)ᵀ
- κ₀+n
-
- where
-
- 1 ₙ
- ȳ = ─ ∑ yᵢ
- n ⁱ⁼¹
-
- ₙ
- S = ∑ (yᵢ − ȳ)(yᵢ − ȳ)ᵀ
- ⁱ⁼¹
-*/
-
-#ifndef NORMAL_CONJUGATE_H
-#define NORMAL_CONJUGATE_H
-
-#include "twod_matrix.hh"
-
-/* The class is described by the four parameters: μ, κ, ν and Λ. */
-
-class NormalConj
-{
-protected:
- Vector mu;
- int kappa;
- int nu;
- TwoDMatrix lambda;
-public:
- /* We provide the following constructors: The first constructs diffuse
- (Jeffrey’s) prior. It sets κ and Λ to zeros, ν to −1 and also the mean μ
- to zero (it should not be referenced). The second constructs the posterior
- using the diffuse prior and the observed data (columnwise). The third is a
- copy constructor. */
- NormalConj(int d);
- NormalConj(const ConstTwoDMatrix &ydata);
- NormalConj(const NormalConj &) = default;
- NormalConj(NormalConj &&) = default;
-
- virtual ~NormalConj() = default;
- void update(const ConstVector &y);
- void update(const ConstTwoDMatrix &ydata);
- void update(const NormalConj &nc);
- int
- getDim() const
- {
- return mu.length();
- }
- const Vector &
- getMean() const
- {
- return mu;
- }
- void getVariance(TwoDMatrix &v) const;
-};
-
-#endif
diff --git a/mex/sources/libkorder/kord/seed_generator.cc b/mex/sources/libkorder/kord/seed_generator.cc
deleted file mode 100644
index ccc138c51..000000000
--- a/mex/sources/libkorder/kord/seed_generator.cc
+++ /dev/null
@@ -1,47 +0,0 @@
-/*
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "seed_generator.hh"
-
-#include
-#include
-
-namespace seed_generator
-{
- std::mutex mut;
-
- std::mt19937 rng;
-
- std::uniform_int_distribution seed_generator(std::numeric_limits::min(),
- std::numeric_limits::max());
-
- std::mt19937::result_type
- get_new_seed()
- {
- std::lock_guard lk{mut};
- return seed_generator(rng);
- }
-
- void
- set_meta_seed(std::mt19937::result_type s)
- {
- std::lock_guard lk{mut};
- rng.seed(s);
- }
-};
diff --git a/mex/sources/libkorder/kord/seed_generator.hh b/mex/sources/libkorder/kord/seed_generator.hh
deleted file mode 100644
index 5679ad2a9..000000000
--- a/mex/sources/libkorder/kord/seed_generator.hh
+++ /dev/null
@@ -1,34 +0,0 @@
-/*
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#ifndef RANDOM_H
-#define RANDOM_H
-
-#include
-
-namespace seed_generator
-{
- // Produces seeds that can be used with Mersenne-Twister generators (thread-safe)
- std::mt19937::result_type get_new_seed();
-
- // Sets the seed for the seed generator (!)
- void set_meta_seed(std::mt19937::result_type s);
-};
-
-#endif
diff --git a/mex/sources/libkorder/tl/t_container.hh b/mex/sources/libkorder/tl/t_container.hh
index 884131eb8..8aef209f3 100644
--- a/mex/sources/libkorder/tl/t_container.hh
+++ b/mex/sources/libkorder/tl/t_container.hh
@@ -75,8 +75,6 @@
#include
#include
-#include
-
// ltsym predicate
/* We need a predicate on strict weak ordering of
symmetries. */
@@ -225,21 +223,6 @@ public:
}
}
- /* Output to the MAT file. */
- void
- writeMat(mat_t *fd, const std::string &prefix) const
- {
- for (auto &it : *this)
- {
- std::string lname = prefix + "_g";
- const Symmetry &sym = it.first;
- for (int i = 0; i < sym.num(); i++)
- lname += '_' + std::to_string(sym[i]);
- ConstTwoDMatrix m(*(it.second));
- m.writeMat(fd, lname);
- }
- }
-
/* Output to the Memory Map. */
void
writeMMap(std::map &mm, const std::string &prefix) const
diff --git a/mex/sources/libkorder/tl/twod_matrix.cc b/mex/sources/libkorder/tl/twod_matrix.cc
index 327445eaf..1e308dc08 100644
--- a/mex/sources/libkorder/tl/twod_matrix.cc
+++ b/mex/sources/libkorder/tl/twod_matrix.cc
@@ -51,25 +51,6 @@ ConstTwoDMatrix::ConstTwoDMatrix(int first_row, int num, const ConstTwoDMatrix &
{
}
-void
-ConstTwoDMatrix::writeMat(mat_t *fd, const std::string &vname) const
-{
- size_t dims[2];
- dims[0] = nrows();
- dims[1] = ncols();
- auto data = std::make_unique(nrows()*ncols());
-
- for (int j = 0; j < ncols(); j++)
- for (int i = 0; i < nrows(); i++)
- data[j*nrows()+i] = get(i, j);
-
- matvar_t *v = Mat_VarCreate(vname.c_str(), MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, data.get(), 0);
-
- Mat_VarWrite(fd, v, MAT_COMPRESSION_NONE);
-
- Mat_VarFree(v);
-}
-
TwoDMatrix &
TwoDMatrix::operator=(const ConstTwoDMatrix &m)
{
diff --git a/mex/sources/libkorder/tl/twod_matrix.hh b/mex/sources/libkorder/tl/twod_matrix.hh
index 1099c031f..728a7691c 100644
--- a/mex/sources/libkorder/tl/twod_matrix.hh
+++ b/mex/sources/libkorder/tl/twod_matrix.hh
@@ -34,8 +34,6 @@
#include "GeneralMatrix.hh"
-#include
-
#include
#include
@@ -75,8 +73,6 @@ public:
ConstTwoDMatrix &operator=(const ConstTwoDMatrix &v) = delete;
ConstTwoDMatrix &operator=(ConstTwoDMatrix &&v) = delete;
-
- void writeMat(mat_t *fd, const std::string &vname) const;
};
class TwoDMatrix : public GeneralMatrix
@@ -200,12 +196,6 @@ public:
// Saves the matrix to a text file
void save(const std::string &fname) const;
-
- void
- writeMat(mat_t *fd, const std::string &vname) const
- {
- ConstTwoDMatrix(*this).writeMat(fd, vname);
- }
};
#endif
diff --git a/mex/sources/libkorder/utils/nlsolve.cc b/mex/sources/libkorder/utils/nlsolve.cc
deleted file mode 100644
index 3286652c3..000000000
--- a/mex/sources/libkorder/utils/nlsolve.cc
+++ /dev/null
@@ -1,277 +0,0 @@
-/*
- * Copyright © 2006 Ondra Kamenik
- * Copyright © 2019-2022 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include "nlsolve.hh"
-#include "dynare_exception.hh"
-
-#include
-#include
-
-using namespace ogu;
-
-double GoldenSectionSearch::golden = (3.-std::sqrt(5.))/2;
-
-double
-GoldenSectionSearch::search(OneDFunction &f, double x1, double x2)
-{
- double b;
- if (init_bracket(f, x1, x2, b))
- {
- double fb = f.eval(b);
- double f1 = f.eval(x1);
- f.eval(x2);
- double dx;
- do
- {
- double w = (b-x1)/(x2-x1);
- dx = std::abs((1-2*w)*(x2-x1));
- double x;
- if (b-x1 > x2-b)
- x = b - dx;
- else
- x = b + dx;
- double fx = f.eval(x);
- if (!std::isfinite(fx))
- return x1;
- if (b-x1 > x2-b)
- {
- // x is on the left from b
- if (f1 > fx && fx < fb)
- {
- // pickup bracket [f1,fx,fb]
- x2 = b;
- fb = fx;
- b = x;
- }
- else
- {
- // pickup bracket [fx,fb,fx2]
- f1 = fx;
- x1 = x;
- }
- }
- else
- {
- // x is on the right from b
- if (f1 > fb && fb < fx)
- // pickup bracket [f1,fb,fx]
- x2 = x;
- else
- {
- // pickup bracket [fb,fx,f2]
- f1 = fb;
- x1 = b;
- fb = fx;
- b = x;
- }
- }
- }
- while (dx > tol);
- }
- return b;
-}
-
-bool
-GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double &b)
-{
- double f1 = f.eval(x1);
- if (!std::isfinite(f1))
- throw DynareException(__FILE__, __LINE__,
- "Safer point not finite in GoldenSectionSearch::init_bracket");
-
- int cnt = 0;
- bool bracket_found = false;
- do
- {
- bool finite_found = search_for_finite(f, x1, x2, b);
- if (!finite_found)
- {
- b = x1;
- return false;
- }
- double f2 = f.eval(x2);
- double fb = f.eval(b);
- double bsym = 2*x2 - b;
- double fbsym = f.eval(bsym);
- // now we know that f1, f2, and fb are finite
- if (std::isfinite(fbsym))
- {
- /* we have four numbers f1, fb, f2, fbsym, we test for the following
- combinations to find the bracket: [f1,f2,fbsym], [f1,fb,fbsym] and
- [f1,fb,fbsym] */
- if (f1 > f2 && f2 < fbsym)
- {
- bracket_found = true;
- b = x2;
- x2 = bsym;
- }
- else if (f1 > fb && fb < fbsym)
- {
- bracket_found = true;
- x2 = bsym;
- }
- else if (f1 > fb && fb < f2)
- bracket_found = true;
- else
- {
- double newx2 = b;
- // choose the smallest value in case we end
- if (f1 > fbsym)
- {
- /* the smallest value is on the other end, we do not want to
- continue */
- b = bsym;
- return false;
- }
- else
- b = x1;
- // move x2 to b in case we continue
- x2 = newx2;
- }
- }
- else
- {
- /* we have only three numbers, we test for the bracket, and if not
- found, we set b as potential result and shorten x2 as potential
- init value for next cycle */
- if (f1 > fb && fb < f2)
- bracket_found = true;
- else
- {
- double newx2 = b;
- // choose the smaller value in case we end
- if (f1 > f2)
- b = x2;
- else
- b = x1;
- // move x2 to b in case we continue
- x2 = newx2;
- }
- }
- cnt++;
- }
- while (!bracket_found && cnt < 5);
-
- return bracket_found;
-}
-
-/* This moves x2 toward to x1 until the function at x2 is finite and b as a
- golden section between x1 and x2 yields also finite f. */
-bool
-GoldenSectionSearch::search_for_finite(OneDFunction &f, double x1, double &x2, double &b)
-{
- int cnt = 0;
- bool found = false;
- do
- {
- double f2 = f.eval(x2);
- b = (1-golden)*x1 + golden*x2;
- double fb = f.eval(b);
- found = std::isfinite(f2) && std::isfinite(fb);
- if (!found)
- x2 = b;
- cnt++;
- }
- while (!found && cnt < 5);
-
- return found;
-}
-
-void
-VectorFunction::check_for_eval(const ConstVector &in, Vector &out) const
-{
- if (inDim() != in.length() || outDim() != out.length())
- throw DynareException(__FILE__, __LINE__,
- "Wrong dimensions in VectorFunction::check_for_eval");
-}
-
-double
-NLSolver::eval(double lambda)
-{
- Vector xx(const_cast(x));
- xx.add(1-lambda, xcauchy);
- xx.add(lambda, xnewton);
- Vector ff(func.outDim());
- func.eval(xx, ff);
- return ff.dot(ff);
-}
-
-bool
-NLSolver::solve(Vector &xx, int &iter)
-{
- JournalRecord rec(journal);
- rec << "Iter lambda residual" << endrec;
- JournalRecord rec1(journal);
- rec1 << "───────────────────────────" << endrec;
-
- x = const_cast(xx);
- iter = 0;
- // setup fx
- Vector fx(func.outDim());
- func.eval(x, fx);
- if (!fx.isFinite())
- throw DynareException(__FILE__, __LINE__,
- "Initial guess does not yield finite residual in NLSolver::solve");
- bool converged = fx.getMax() < tol;
- JournalRecord rec2(journal);
- auto format_double = [](double v)
- {
- std::ostringstream buf;
- buf << std::setw(11) << v;
- return buf.str();
- };
- rec2 << iter << " N/A " << format_double(fx.getMax()) << endrec;
- while (!converged && iter < max_iter)
- {
- // setup Jacobian
- jacob.eval(x);
- // calculate cauchy step
- Vector g(func.inDim());
- g.zeros();
- ConstTwoDMatrix(jacob).multaVecTrans(g, fx);
- Vector Jg(func.inDim());
- Jg.zeros();
- ConstTwoDMatrix(jacob).multaVec(Jg, g);
- double m = -g.dot(g)/Jg.dot(Jg);
- xcauchy = const_cast(g);
- xcauchy.mult(m);
- // calculate newton step
- xnewton = const_cast(fx);
- ConstTwoDMatrix(jacob).multInvLeft(xnewton);
- xnewton.mult(-1);
-
- // line search
- double lambda = GoldenSectionSearch::search(*this, 0, 1);
- x.add(1-lambda, xcauchy);
- x.add(lambda, xnewton);
- // evaluate func
- func.eval(x, fx);
- converged = fx.getMax() < tol;
-
- // iter
- iter++;
-
- JournalRecord rec3(journal);
- rec3 << iter << " " << lambda << " " << format_double(fx.getMax()) << endrec;
- }
- xx = const_cast(x);
-
- return converged;
-}
diff --git a/mex/sources/libkorder/utils/nlsolve.hh b/mex/sources/libkorder/utils/nlsolve.hh
deleted file mode 100644
index b1e1bc21f..000000000
--- a/mex/sources/libkorder/utils/nlsolve.hh
+++ /dev/null
@@ -1,113 +0,0 @@
-/*
- * Copyright © 2006 Ondra Kamenik
- * Copyright © 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#ifndef OGU_NLSOLVE_H
-#define OGU_NLSOLVE_H
-
-#include "twod_matrix.hh"
-#include "journal.hh"
-
-#include
-
-namespace ogu
-{
- class OneDFunction
- {
- public:
- virtual ~OneDFunction() = default;
- virtual double eval(double) = 0;
- };
-
- class GoldenSectionSearch
- {
- protected:
- constexpr static double tol = 1.e-4;
- /* This is equal to the golden section ratio. */
- static double golden;
- public:
- static double search(OneDFunction &f, double x1, double x2);
- protected:
- /* This initializes a bracket by moving x2 and b (as a golden section of
- x1,x2) so that f(x1)>f(b) && f(b)