diff options
Diffstat (limited to 'gsl-1.9/monte')
| -rw-r--r-- | gsl-1.9/monte/ChangeLog | 166 | ||||
| -rw-r--r-- | gsl-1.9/monte/Makefile.am | 21 | ||||
| -rw-r--r-- | gsl-1.9/monte/Makefile.in | 547 | ||||
| -rw-r--r-- | gsl-1.9/monte/README | 37 | ||||
| -rw-r--r-- | gsl-1.9/monte/TODO | 25 | ||||
| -rw-r--r-- | gsl-1.9/monte/gsl_monte.h | 55 | ||||
| -rw-r--r-- | gsl-1.9/monte/gsl_monte_miser.h | 83 | ||||
| -rw-r--r-- | gsl-1.9/monte/gsl_monte_plain.h | 65 | ||||
| -rw-r--r-- | gsl-1.9/monte/gsl_monte_vegas.h | 105 | ||||
| -rw-r--r-- | gsl-1.9/monte/miser.c | 681 | ||||
| -rw-r--r-- | gsl-1.9/monte/plain.c | 151 | ||||
| -rw-r--r-- | gsl-1.9/monte/test.c | 373 | ||||
| -rw-r--r-- | gsl-1.9/monte/test_main.c | 70 | ||||
| -rw-r--r-- | gsl-1.9/monte/vegas.c | 864 |
14 files changed, 3243 insertions, 0 deletions
diff --git a/gsl-1.9/monte/ChangeLog b/gsl-1.9/monte/ChangeLog new file mode 100644 index 0000000..2ae49d2 --- /dev/null +++ b/gsl-1.9/monte/ChangeLog @@ -0,0 +1,166 @@ +2004-06-02 Brian Gough <bjg@network-theory.co.uk> + + * test_main.c: handle the case where sd==0 && error[i] !=0 more + carefully + +Mon Apr 23 13:23:47 2001 Brian Gough <bjg@network-theory.co.uk> + + * test_main.c: removed unused status variable + +Sat Jan 6 19:56:49 2001 Brian Gough <bjg@network-theory.co.uk> + + * miser.c (gsl_monte_miser_free): fixed memory leak for s->hits_{r,l} + + * vegas.c (gsl_monte_vegas_free): fixed memory leak for s->x + +Fri Dec 22 21:43:04 2000 Brian Gough <bjg@network-theory.co.uk> + + * miser.c: removed max-min estimation method for subvolumes. We + will use the standard variance method and try to use a sufficient + number of points to estimate the variance reliably. + +Wed Dec 20 21:32:40 2000 Brian Gough <bjg@network-theory.co.uk> + + * vegas.c: tidied up the algorithm, deal with cases of sigma = 0 + explicitly. + +Sat Dec 9 14:19:53 2000 Brian Gough <bjg@network-theory.co.uk> + + * reorganization and clean up + +Thu Nov 16 19:50:27 2000 Brian Gough <bjg@network-theory.co.uk> + + * miser.c: rewrite to fix overflows and make calling conventions + consistent + + * plain.c: rewrite to fix overflows and make calling conventions + consistent + +Thu Oct 26 20:06:36 2000 Brian Gough <bjg@network-theory.co.uk> + + * plain.c (gsl_monte_plain_integrate): integer factor + calls*(calls-1) used in numerical expression can easily overflow, + changed to calls*(calls-1.0). + +Sat Oct 21 20:36:06 2000 Brian Gough <bjg@network-theory.co.uk> + + * miser.c (gsl_monte_miser_integrate): fixed bug where hits_l was + used in place of hits_r + +Tue Sep 19 19:16:37 2000 Brian Gough <bjg@network-theory.co.uk> + + * plain.c (gsl_monte_plain_alloc): initialise check_done to avoid + warning about use of unitialised memory + + * vegas.c (gsl_monte_vegas_alloc): as above + + * miser.c (gsl_monte_miser_alloc): as above + + * plain.c miser.c: removed use of sprintf for error messages + +Wed May 31 19:50:19 2000 Brian Gough <bjg@network-theory.co.uk> + + * miser_test.c (main): increased tolerances to allow tests to pass + with different compilers + +Mon May 15 15:26:22 2000 Brian Gough <bjg@network-theory.co.uk> + + * vegas_test.c (main): added gsl_ieee_env_setup() to allow change + of precision in tests + + * miser_test.c (main): ditto + + * plain_test.c (main): ditto + +Fri Feb 26 14:49:56 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * Makefile.am: removed ..._LDFLAGS = -static since this is + specific to gcc + + + +Wed Nov 18 10:59:56 1998 Brian Gough <bjg@vvv.lanl.gov> + + * use standard headers templates_on.h and templates_off.h instead + of source.h + +Tue Nov 17 16:49:12 1998 Brian Gough <bjg@vvv.lanl.gov> + + * added #include <config.h> to all top-level source files + + * plain.c (gsl_monte_plain_integrate): replaced myMAX by GSL_MAX + + * utils.c: move macros around to avoid double definition + +Fri Aug 14 10:12:06 1998 Brian Gough <bjg@vvv.lanl.gov> + + * Makefile.am: I needed to add utils.h to libgslmonte_a_SOURCES to + get it to work with my automake + +Thu Jul 30 17:31:29 1998 booth <booth@planck.pha.jhu.edu> + + * gsl_monte_miser.h, miser.c: + Turn off the annoying warning in miser unless the user requests it. + +Wed Jul 29 20:24:54 1998 bjg <bjg@vvv.lanl.gov> + + * Makefile.am, Makefile.in: some fixes to pass make distcheck + + * Makefile.am, Makefile.in: + experimenting with new top level makefile.am + +Tue Jul 28 17:05:20 1998 booth <booth@planck.pha.jhu.edu> + + * plain.c, vegas.c, miser.c: + make all the _free functions check their argument. Also, the init functions + now throw EFAULT if given a null pointer. + + * gsl_monte_vegas.h, vegas.c: + vegas1, vegas2 and vegas3 all go away. Get them by setting the "stage" + variabe appropriately. Also, make _free check its argument. + + * vegas.c, gsl_monte_vegas.h: minor cleanup prior to be change. + + * vegas.c: minor changes, commiting by accident. + +Mon Jul 27 22:52:49 1998 bjg <bjg@vvv.lanl.gov> + + * Makefile.in, Makefile.am: + fixed some of the include requirements for make dist + +Mon Jul 27 15:19:54 1998 booth <booth@planck.pha.jhu.edu> + + * vegas_print.h, vegas_test.c, miser_test.c, vegas-print.c, vegas.c, Attic/gsl_vegas.h, Attic/gsl_vegas_print.h, gsl_monte_vegas.h, miser.c, Attic/gsl_miser.h, Makefile.am, Makefile.in, gsl_monte_miser.h: + Renamed public header files to follow convention correctly. + + * vegas.c, vegas_test.c, miser.c, miser_test.c, plain.c, plain_test.c, Attic/gsl_miser.h, Attic/gsl_vegas.h, gsl_monte_plain.h: + All the integration functions now end with _integrate (except for + vegas1, vegas2 and vegas3 which are going away RSN). + +Tue Jul 21 21:54:33 1998 booth <booth@planck.pha.jhu.edu> + + * vegas.c, vegas-print.c, Attic/gsl_vegas_print.h, Attic/gsl_vegas.h: + trivial stuff: eliminate compiler warnings, eliminate some unneeded variables, + change some types to make consistent with plain and miser. + + * gsl_monte_plain.h, plain.c, plain_test.c: + Make "plain" conform to same style as miser and vegas. + +Fri Jul 17 02:23:40 1998 jungman <jungman@nnn.lanl.gov> + + * Makefile.in: we're gonna make it... + +Thu Jul 16 16:23:45 1998 booth <booth@planck.pha.jhu.edu> + + * vegas-print.c, vegas.c, Attic/gsl_vegas_print.h: + Have now completely eliminated all static variables. + + * vegas_test.c, vegas.c, Attic/gsl_vegas.h, Attic/gsl_vegas_print.h, vegas-print.c: + Continuing to remove all static variables from vegas. + +Wed Jul 15 19:10:24 1998 booth <booth@planck.pha.jhu.edu> + + * vegas.c, vegas_test.c, Attic/gsl_vegas.h: + Do the state-structure thing for vegas. Not finished, so far the only + real content is the rng structure. + diff --git a/gsl-1.9/monte/Makefile.am b/gsl-1.9/monte/Makefile.am new file mode 100644 index 0000000..fec4ff4 --- /dev/null +++ b/gsl-1.9/monte/Makefile.am @@ -0,0 +1,21 @@ +noinst_LTLIBRARIES = libgslmonte.la + +libgslmonte_la_SOURCES = miser.c plain.c gsl_monte_vegas.h gsl_monte_miser.h gsl_monte_plain.h gsl_monte.h vegas.c + +pkginclude_HEADERS = gsl_monte.h gsl_monte_vegas.h gsl_monte_miser.h gsl_monte_plain.h + +INCLUDES= -I$(top_builddir) -I$(top_srcdir) + +TESTS = $(check_PROGRAMS) + +check_PROGRAMS = test #demo + +test_SOURCES = test.c +test_LDADD = libgslmonte.la ../rng/libgslrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + +noinst_HEADERS = test_main.c + +#demo_SOURCES= demo.c +#demo_LDADD = libgslmonte.la ../rng/libgslrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la + + diff --git a/gsl-1.9/monte/Makefile.in b/gsl-1.9/monte/Makefile.in new file mode 100644 index 0000000..acbd39e --- /dev/null +++ b/gsl-1.9/monte/Makefile.in @@ -0,0 +1,547 @@ +# Makefile.in generated by automake 1.9.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005 Free Software Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +top_builddir = .. +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +INSTALL = @INSTALL@ +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +check_PROGRAMS = test$(EXEEXT) +subdir = monte +DIST_COMMON = README $(noinst_HEADERS) $(pkginclude_HEADERS) \ + $(srcdir)/Makefile.am $(srcdir)/Makefile.in ChangeLog TODO +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(SHELL) $(top_srcdir)/mkinstalldirs +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +LTLIBRARIES = $(noinst_LTLIBRARIES) +libgslmonte_la_LIBADD = +am_libgslmonte_la_OBJECTS = miser.lo plain.lo vegas.lo +libgslmonte_la_OBJECTS = $(am_libgslmonte_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslmonte.la ../rng/libgslrng.la \ + ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la \ + ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir) +depcomp = +am__depfiles_maybe = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +SOURCES = $(libgslmonte_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslmonte_la_SOURCES) $(test_SOURCES) +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = `echo $$p | sed -e 's|^.*/||'`; +am__installdirs = "$(DESTDIR)$(pkgincludedir)" +pkgincludeHEADERS_INSTALL = $(INSTALL_HEADER) +HEADERS = $(noinst_HEADERS) $(pkginclude_HEADERS) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +ECHO = @ECHO@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GSL_LT_CBLAS_VERSION = @GSL_LT_CBLAS_VERSION@ +GSL_LT_VERSION = @GSL_LT_VERSION@ +HAVE_AIX_IEEE_INTERFACE = @HAVE_AIX_IEEE_INTERFACE@ +HAVE_DARWIN86_IEEE_INTERFACE = @HAVE_DARWIN86_IEEE_INTERFACE@ +HAVE_DARWIN_IEEE_INTERFACE = @HAVE_DARWIN_IEEE_INTERFACE@ +HAVE_EXTENDED_PRECISION_REGISTERS = @HAVE_EXTENDED_PRECISION_REGISTERS@ +HAVE_FREEBSD_IEEE_INTERFACE = @HAVE_FREEBSD_IEEE_INTERFACE@ +HAVE_GNUM68K_IEEE_INTERFACE = @HAVE_GNUM68K_IEEE_INTERFACE@ +HAVE_GNUPPC_IEEE_INTERFACE = @HAVE_GNUPPC_IEEE_INTERFACE@ +HAVE_GNUSPARC_IEEE_INTERFACE = @HAVE_GNUSPARC_IEEE_INTERFACE@ +HAVE_GNUX86_IEEE_INTERFACE = @HAVE_GNUX86_IEEE_INTERFACE@ +HAVE_HPUX11_IEEE_INTERFACE = @HAVE_HPUX11_IEEE_INTERFACE@ +HAVE_HPUX_IEEE_INTERFACE = @HAVE_HPUX_IEEE_INTERFACE@ +HAVE_IEEE_COMPARISONS = @HAVE_IEEE_COMPARISONS@ +HAVE_IEEE_DENORMALS = @HAVE_IEEE_DENORMALS@ +HAVE_INLINE = @HAVE_INLINE@ +HAVE_IRIX_IEEE_INTERFACE = @HAVE_IRIX_IEEE_INTERFACE@ +HAVE_NETBSD_IEEE_INTERFACE = @HAVE_NETBSD_IEEE_INTERFACE@ +HAVE_OPENBSD_IEEE_INTERFACE = @HAVE_OPENBSD_IEEE_INTERFACE@ +HAVE_OS2EMX_IEEE_INTERFACE = @HAVE_OS2EMX_IEEE_INTERFACE@ +HAVE_PRINTF_LONGDOUBLE = @HAVE_PRINTF_LONGDOUBLE@ +HAVE_SOLARIS_IEEE_INTERFACE = @HAVE_SOLARIS_IEEE_INTERFACE@ +HAVE_SUNOS4_IEEE_INTERFACE = @HAVE_SUNOS4_IEEE_INTERFACE@ +HAVE_TRU64_IEEE_INTERFACE = @HAVE_TRU64_IEEE_INTERFACE@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAINT = @MAINT@ +MAINTAINER_MODE_FALSE = @MAINTAINER_MODE_FALSE@ +MAINTAINER_MODE_TRUE = @MAINTAINER_MODE_TRUE@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +RELEASED = @RELEASED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_STRIP = @ac_ct_STRIP@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +datadir = @datadir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +noinst_LTLIBRARIES = libgslmonte.la +libgslmonte_la_SOURCES = miser.c plain.c gsl_monte_vegas.h gsl_monte_miser.h gsl_monte_plain.h gsl_monte.h vegas.c +pkginclude_HEADERS = gsl_monte.h gsl_monte_vegas.h gsl_monte_miser.h gsl_monte_plain.h +INCLUDES = -I$(top_builddir) -I$(top_srcdir) +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c +test_LDADD = libgslmonte.la ../rng/libgslrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +noinst_HEADERS = test_main.c +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \ + && exit 0; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps monte/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps monte/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +clean-noinstLTLIBRARIES: + -test -z "$(noinst_LTLIBRARIES)" || rm -f $(noinst_LTLIBRARIES) + @list='$(noinst_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +libgslmonte.la: $(libgslmonte_la_OBJECTS) $(libgslmonte_la_DEPENDENCIES) + $(LINK) $(libgslmonte_la_LDFLAGS) $(libgslmonte_la_OBJECTS) $(libgslmonte_la_LIBADD) $(LIBS) + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; for p in $$list; do \ + f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f $$p $$f"; \ + rm -f $$p $$f ; \ + done +test$(EXEEXT): $(test_OBJECTS) $(test_DEPENDENCIES) + @rm -f test$(EXEEXT) + $(LINK) $(test_LDFLAGS) $(test_OBJECTS) $(test_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +distclean-libtool: + -rm -f libtool +uninstall-info-am: +install-pkgincludeHEADERS: $(pkginclude_HEADERS) + @$(NORMAL_INSTALL) + test -z "$(pkgincludedir)" || $(mkdir_p) "$(DESTDIR)$(pkgincludedir)" + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \ + f=$(am__strip_dir) \ + echo " $(pkgincludeHEADERS_INSTALL) '$$d$$p' '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + $(pkgincludeHEADERS_INSTALL) "$$d$$p" "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +uninstall-pkgincludeHEADERS: + @$(NORMAL_UNINSTALL) + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + f=$(am__strip_dir) \ + echo " rm -f '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + rm -f "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +check-TESTS: $(TESTS) + @failed=0; all=0; xfail=0; xpass=0; skip=0; \ + srcdir=$(srcdir); export srcdir; \ + list='$(TESTS)'; \ + if test -n "$$list"; then \ + for tst in $$list; do \ + if test -f ./$$tst; then dir=./; \ + elif test -f $$tst; then dir=; \ + else dir="$(srcdir)/"; fi; \ + if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xpass=`expr $$xpass + 1`; \ + failed=`expr $$failed + 1`; \ + echo "XPASS: $$tst"; \ + ;; \ + *) \ + echo "PASS: $$tst"; \ + ;; \ + esac; \ + elif test $$? -ne 77; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xfail=`expr $$xfail + 1`; \ + echo "XFAIL: $$tst"; \ + ;; \ + *) \ + failed=`expr $$failed + 1`; \ + echo "FAIL: $$tst"; \ + ;; \ + esac; \ + else \ + skip=`expr $$skip + 1`; \ + echo "SKIP: $$tst"; \ + fi; \ + done; \ + if test "$$failed" -eq 0; then \ + if test "$$xfail" -eq 0; then \ + banner="All $$all tests passed"; \ + else \ + banner="All $$all tests behaved as expected ($$xfail expected failures)"; \ + fi; \ + else \ + if test "$$xpass" -eq 0; then \ + banner="$$failed of $$all tests failed"; \ + else \ + banner="$$failed of $$all tests did not behave as expected ($$xpass unexpected passes)"; \ + fi; \ + fi; \ + dashes="$$banner"; \ + skipped=""; \ + if test "$$skip" -ne 0; then \ + skipped="($$skip tests were not run)"; \ + test `echo "$$skipped" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$skipped"; \ + fi; \ + report=""; \ + if test "$$failed" -ne 0 && test -n "$(PACKAGE_BUGREPORT)"; then \ + report="Please report to $(PACKAGE_BUGREPORT)"; \ + test `echo "$$report" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$report"; \ + fi; \ + dashes=`echo "$$dashes" | sed s/./=/g`; \ + echo "$$dashes"; \ + echo "$$banner"; \ + test -z "$$skipped" || echo "$$skipped"; \ + test -z "$$report" || echo "$$report"; \ + echo "$$dashes"; \ + test "$$failed" -eq 0; \ + else :; fi + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-TESTS +check: check-am +all-am: Makefile $(LTLIBRARIES) $(HEADERS) +installdirs: + for dir in "$(DESTDIR)$(pkgincludedir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-libtool distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: install-pkgincludeHEADERS + +install-exec-am: + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-info-am uninstall-pkgincludeHEADERS + +.PHONY: CTAGS GTAGS all all-am check check-TESTS check-am clean \ + clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES ctags distclean distclean-compile \ + distclean-generic distclean-libtool distclean-tags distdir dvi \ + dvi-am html html-am info info-am install install-am \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man \ + install-pkgincludeHEADERS install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-info-am \ + uninstall-pkgincludeHEADERS + + +#demo_SOURCES= demo.c +#demo_LDADD = libgslmonte.la ../rng/libgslrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gsl-1.9/monte/README b/gsl-1.9/monte/README new file mode 100644 index 0000000..9658d7b --- /dev/null +++ b/gsl-1.9/monte/README @@ -0,0 +1,37 @@ +G. P. Lepage, "VEGAS..." J. Comp. Phys. 27, 192(1978). +W.H. Press, G.R. Farrar, "Miser..." Computers in Physics, v4 (1990), pp190-195 + +References from PVEGAS by Richard Kreckel + +[1]: G.P. Lepage: A New Algorithm for Adaptive Multidimensional + Integration; Journal of Computational Physics 27, + 192-203, (1978) + +[2]: W. Press, S. Teukolsky, W. Vetterling, B. Flannery: Numerical + Recipes in C, (second edition) Cambridge University Press, + 1992. + +[3]: R.C. Tausworthe: Random numbers generated by linear + recurrence modulo two, Math. Comput. 19, 201-209, (1965) + +[4]: I. Deak: Uniform random number generators for parallel + computers; Parallel Computing, 15, 155-164, (1990) + +[5]: R. Kreckel: Parallelization of adaptive MC integrators, + MZ-TH/97-30, physics/9710028; Comp. Phys. Comm., 106, 258-266, + (1997) (A preliminary version of this article can be retrieved + via anonymous ftp from ftpthep.physik.uni-mainz.de in the directory + /pub/pvegas/.) + +[6]: S. Veseli: Multidimensional integration in a heterogeneous + network environment; FERMILAB-PUB-97/271-T; physics/9710017 + +[7]: R. Kreckel: Addendum: Parallelization of adaptive MC integrators, + (Available via anonymous ftp from ftpthep.physik.uni-mainz.de in + the directory /pub/pvegas/.) + +[8]: MPI-Forum: MPI: A Message-Passing Interface Standard, University + of Tennessee, Knoxville, Tennessee, (1995) + +[9]: T. Ohl: Vegas Revisited: Adaptive Monte Carlo Integration Beyond + Factorization; hep-ph/9806432 diff --git a/gsl-1.9/monte/TODO b/gsl-1.9/monte/TODO new file mode 100644 index 0000000..29d2eb4 --- /dev/null +++ b/gsl-1.9/monte/TODO @@ -0,0 +1,25 @@ +* Add Lattice Method and Quasi-random Method to monte carlo integration + +* Fix the "No-points in left/right half space" error in miser. Random +errors like that are discouraged in a library. The routine should +iterate over the dimensions choosing a point on each side of the +bisection to ensure that the error does not occur. + +* VEGAS could estimate its roundoff error, caused by dividing up into +many boxes and then summing (I know this is ridiculous for any +realistic case, but it shows up on the tests when integrating a +constant!). In particular, there are some failures reported for the +VEGAS constant tests on Windows. + +* VEGAS gives a negative chisq for some cases where there are wildly +inconsistent values. Calculation should be improved so that chisq>=0. + +Also, when there are inconsistent values it might make sense to scale +the error by chisq, as is often done experimentally -- at least that +way the error would be increased. + +* The original VEGAS allowed to user to calculate other weighted +quantities. The appropriate way to implement it in GSL would be to +provide a separate routine which returns sample points and weights so +that the user can caculate any quantity. + diff --git a/gsl-1.9/monte/gsl_monte.h b/gsl-1.9/monte/gsl_monte.h new file mode 100644 index 0000000..dd9c9d8 --- /dev/null +++ b/gsl-1.9/monte/gsl_monte.h @@ -0,0 +1,55 @@ +/* monte/gsl_monte.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Some things common to all the Monte-Carlo implementations */ +/* Author: MJB */ + +#ifndef __GSL_MONTE_H__ +#define __GSL_MONTE_H__ + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +/* Hide the function type in a typedef so that we can use it in all our + integration functions, and make it easy to change things. +*/ + +struct gsl_monte_function_struct { + double (*f)(double * x_array, size_t dim, void * params); + size_t dim; + void * params; +}; + +typedef struct gsl_monte_function_struct gsl_monte_function; + +#define GSL_MONTE_FN_EVAL(F,x) (*((F)->f))(x,(F)->dim,(F)->params) + + +__END_DECLS + +#endif /* __GSL_MONTE_H__ */ diff --git a/gsl-1.9/monte/gsl_monte_miser.h b/gsl-1.9/monte/gsl_monte_miser.h new file mode 100644 index 0000000..d08cb8c --- /dev/null +++ b/gsl-1.9/monte/gsl_monte_miser.h @@ -0,0 +1,83 @@ +/* monte/gsl_monte_miser.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: MJB */ + +#ifndef __GSL_MONTE_MISER_H__ +#define __GSL_MONTE_MISER_H__ + +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte.h> +#include <gsl/gsl_monte_plain.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +typedef struct { + size_t min_calls; + size_t min_calls_per_bisection; + double dither; + double estimate_frac; + double alpha; + size_t dim; + int estimate_style; + int depth; + int verbose; + double * x; + double * xmid; + double * sigma_l; + double * sigma_r; + double * fmax_l; + double * fmax_r; + double * fmin_l; + double * fmin_r; + double * fsum_l; + double * fsum_r; + double * fsum2_l; + double * fsum2_r; + size_t * hits_l; + size_t * hits_r; +} gsl_monte_miser_state; + +int gsl_monte_miser_integrate(gsl_monte_function * f, + const double xl[], const double xh[], + size_t dim, size_t calls, + gsl_rng *r, + gsl_monte_miser_state* state, + double *result, double *abserr); + +gsl_monte_miser_state* gsl_monte_miser_alloc(size_t dim); + +int gsl_monte_miser_init(gsl_monte_miser_state* state); + +void gsl_monte_miser_free(gsl_monte_miser_state* state); + + +__END_DECLS + +#endif /* __GSL_MONTE_MISER_H__ */ diff --git a/gsl-1.9/monte/gsl_monte_plain.h b/gsl-1.9/monte/gsl_monte_plain.h new file mode 100644 index 0000000..c3d2b28 --- /dev/null +++ b/gsl-1.9/monte/gsl_monte_plain.h @@ -0,0 +1,65 @@ +/* monte/gsl_monte_plain.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Plain Monte-Carlo. */ + +/* Author: MJB */ + +#ifndef __GSL_MONTE_PLAIN_H__ +#define __GSL_MONTE_PLAIN_H__ + +#include <stdio.h> +#include <gsl/gsl_monte.h> +#include <gsl/gsl_rng.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +typedef struct { + size_t dim; + double *x; +} gsl_monte_plain_state; + +int +gsl_monte_plain_integrate (const gsl_monte_function * f, + const double xl[], const double xu[], + const size_t dim, + const size_t calls, + gsl_rng * r, + gsl_monte_plain_state * state, + double *result, double *abserr); + +gsl_monte_plain_state* gsl_monte_plain_alloc(size_t dim); + +int gsl_monte_plain_init(gsl_monte_plain_state* state); + +void gsl_monte_plain_free (gsl_monte_plain_state* state); + +__END_DECLS + +#endif /* __GSL_MONTE_PLAIN_H__ */ diff --git a/gsl-1.9/monte/gsl_monte_vegas.h b/gsl-1.9/monte/gsl_monte_vegas.h new file mode 100644 index 0000000..b328307 --- /dev/null +++ b/gsl-1.9/monte/gsl_monte_vegas.h @@ -0,0 +1,105 @@ +/* monte/gsl_monte_vegas.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* header for the gsl "vegas" routines. Mike Booth, May 1998 */ + +#ifndef __GSL_MONTE_VEGAS_H__ +#define __GSL_MONTE_VEGAS_H__ + +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +enum {GSL_VEGAS_MODE_IMPORTANCE = 1, + GSL_VEGAS_MODE_IMPORTANCE_ONLY = 0, + GSL_VEGAS_MODE_STRATIFIED = -1}; + +typedef struct { + /* grid */ + size_t dim; + size_t bins_max; + unsigned int bins; + unsigned int boxes; /* these are both counted along the axes */ + double * xi; + double * xin; + double * delx; + double * weight; + double vol; + + double * x; + int * bin; + int * box; + + /* distribution */ + double * d; + + /* control variables */ + double alpha; + int mode; + int verbose; + unsigned int iterations; + int stage; + + /* scratch variables preserved between calls to vegas1/2/3 */ + double jac; + double wtd_int_sum; + double sum_wgts; + double chi_sum; + double chisq; + + double result; + double sigma; + + unsigned int it_start; + unsigned int it_num; + unsigned int samples; + unsigned int calls_per_box; + + FILE * ostream; + +} gsl_monte_vegas_state; + +int gsl_monte_vegas_integrate(gsl_monte_function * f, + double xl[], double xu[], + size_t dim, size_t calls, + gsl_rng * r, + gsl_monte_vegas_state *state, + double* result, double* abserr); + +gsl_monte_vegas_state* gsl_monte_vegas_alloc(size_t dim); + +int gsl_monte_vegas_init(gsl_monte_vegas_state* state); + +void gsl_monte_vegas_free (gsl_monte_vegas_state* state); + +__END_DECLS + +#endif /* __GSL_MONTE_VEGAS_H__ */ + diff --git a/gsl-1.9/monte/miser.c b/gsl-1.9/monte/miser.c new file mode 100644 index 0000000..92e9215 --- /dev/null +++ b/gsl-1.9/monte/miser.c @@ -0,0 +1,681 @@ +/* monte/miser.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* MISER. Based on the algorithm described in the following article, + + W.H. Press, G.R. Farrar, "Recursive Stratified Sampling for + Multidimensional Monte Carlo Integration", Computers in Physics, + v4 (1990), pp190-195. + +*/ + +/* Author: MJB */ +/* Modified by Brian Gough 12/2000 */ + +#include <config.h> +#include <math.h> +#include <stdlib.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte.h> +#include <gsl/gsl_monte_miser.h> + +static int +estimate_corrmc (gsl_monte_function * f, + const double xl[], const double xu[], + size_t dim, size_t calls, + gsl_rng * r, + gsl_monte_miser_state * state, + double *result, double *abserr, + const double xmid[], double sigma_l[], double sigma_r[]); + + +int +gsl_monte_miser_integrate (gsl_monte_function * f, + const double xl[], const double xu[], + size_t dim, size_t calls, + gsl_rng * r, + gsl_monte_miser_state * state, + double *result, double *abserr) +{ + size_t n, estimate_calls, calls_l, calls_r; + const size_t min_calls = state->min_calls; + size_t i; + size_t i_bisect; + int found_best; + + double res_est = 0, err_est = 0; + double res_r = 0, err_r = 0, res_l = 0, err_l = 0; + double xbi_l, xbi_m, xbi_r, s; + + double vol; + double weight_l, weight_r; + + double *x = state->x; + double *xmid = state->xmid; + double *sigma_l = state->sigma_l, *sigma_r = state->sigma_r; + + if (dim != state->dim) + { + GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL); + } + + for (i = 0; i < dim; i++) + { + if (xu[i] <= xl[i]) + { + GSL_ERROR ("xu must be greater than xl", GSL_EINVAL); + } + + if (xu[i] - xl[i] > GSL_DBL_MAX) + { + GSL_ERROR ("Range of integration is too large, please rescale", + GSL_EINVAL); + } + } + + if (state->alpha < 0) + { + GSL_ERROR ("alpha must be non-negative", GSL_EINVAL); + } + + /* Compute volume */ + + vol = 1; + + for (i = 0; i < dim; i++) + { + vol *= xu[i] - xl[i]; + } + + if (calls < state->min_calls_per_bisection) + { + double m = 0.0, q = 0.0; + + if (calls < 2) + { + GSL_ERROR ("insufficient calls for subvolume", GSL_EFAILED); + } + + for (n = 0; n < calls; n++) + { + /* Choose a random point in the integration region */ + + for (i = 0; i < dim; i++) + { + x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]); + } + + { + double fval = GSL_MONTE_FN_EVAL (f, x); + + /* recurrence for mean and variance */ + + double d = fval - m; + m += d / (n + 1.0); + q += d * d * (n / (n + 1.0)); + } + } + + *result = vol * m; + + *abserr = vol * sqrt (q / (calls * (calls - 1.0))); + + return GSL_SUCCESS; + } + + estimate_calls = GSL_MAX (min_calls, calls * (state->estimate_frac)); + + if (estimate_calls < 4 * dim) + { + GSL_ERROR ("insufficient calls to sample all halfspaces", GSL_ESANITY); + } + + /* Flip coins to bisect the integration region with some fuzz */ + + for (i = 0; i < dim; i++) + { + s = (gsl_rng_uniform (r) - 0.5) >= 0.0 ? state->dither : -state->dither; + state->xmid[i] = (0.5 + s) * xl[i] + (0.5 - s) * xu[i]; + } + + /* The idea is to chose the direction to bisect based on which will + give the smallest total variance. We could (and may do so later) + use MC to compute these variances. But the NR guys simply estimate + the variances by finding the min and max function values + for each half-region for each bisection. */ + + estimate_corrmc (f, xl, xu, dim, estimate_calls, + r, state, &res_est, &err_est, xmid, sigma_l, sigma_r); + + /* We have now used up some calls for the estimation */ + + calls -= estimate_calls; + + /* Now find direction with the smallest total "variance" */ + + { + double best_var = GSL_DBL_MAX; + double beta = 2.0 / (1.0 + state->alpha); + found_best = 0; + i_bisect = 0; + weight_l = weight_r = 1.0; + + for (i = 0; i < dim; i++) + { + if (sigma_l[i] >= 0 && sigma_r[i] >= 0) + { + /* estimates are okay */ + double var = pow (sigma_l[i], beta) + pow (sigma_r[i], beta); + + if (var <= best_var) + { + found_best = 1; + best_var = var; + i_bisect = i; + weight_l = pow (sigma_l[i], beta); + weight_r = pow (sigma_r[i], beta); + } + } + else + { + if (sigma_l[i] < 0) + { + GSL_ERROR ("no points in left-half space!", GSL_ESANITY); + } + if (sigma_r[i] < 0) + { + GSL_ERROR ("no points in right-half space!", GSL_ESANITY); + } + } + } + } + + if (!found_best) + { + /* All estimates were the same, so chose a direction at random */ + + i_bisect = gsl_rng_uniform_int (r, dim); + } + + xbi_l = xl[i_bisect]; + xbi_m = xmid[i_bisect]; + xbi_r = xu[i_bisect]; + + /* Get the actual fractional sizes of the two "halves", and + distribute the remaining calls among them */ + + { + double fraction_l = fabs ((xbi_m - xbi_l) / (xbi_r - xbi_l)); + double fraction_r = 1 - fraction_l; + + double a = fraction_l * weight_l; + double b = fraction_r * weight_r; + + calls_l = min_calls + (calls - 2 * min_calls) * a / (a + b); + calls_r = min_calls + (calls - 2 * min_calls) * b / (a + b); + } + + /* Compute the integral for the left hand side of the bisection */ + + /* Due to the recursive nature of the algorithm we must allocate + some new memory for each recursive call */ + + { + int status; + + double *xu_tmp = (double *) malloc (dim * sizeof (double)); + + if (xu_tmp == 0) + { + GSL_ERROR_VAL ("out of memory for left workspace", GSL_ENOMEM, 0); + } + + for (i = 0; i < dim; i++) + { + xu_tmp[i] = xu[i]; + } + + xu_tmp[i_bisect] = xbi_m; + + status = gsl_monte_miser_integrate (f, xl, xu_tmp, + dim, calls_l, r, state, + &res_l, &err_l); + free (xu_tmp); + + if (status != GSL_SUCCESS) + { + return status; + } + } + + /* Compute the integral for the right hand side of the bisection */ + + { + int status; + + double *xl_tmp = (double *) malloc (dim * sizeof (double)); + + if (xl_tmp == 0) + { + GSL_ERROR_VAL ("out of memory for right workspace", GSL_ENOMEM, 0); + } + + for (i = 0; i < dim; i++) + { + xl_tmp[i] = xl[i]; + } + + xl_tmp[i_bisect] = xbi_m; + + status = gsl_monte_miser_integrate (f, xl_tmp, xu, + dim, calls_r, r, state, + &res_r, &err_r); + free (xl_tmp); + + if (status != GSL_SUCCESS) + { + return status; + } + } + + *result = res_l + res_r; + *abserr = sqrt (err_l * err_l + err_r * err_r); + + return GSL_SUCCESS; +} + +gsl_monte_miser_state * +gsl_monte_miser_alloc (size_t dim) +{ + gsl_monte_miser_state *s = + (gsl_monte_miser_state *) malloc (sizeof (gsl_monte_miser_state)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for miser state struct", + GSL_ENOMEM, 0); + } + + s->x = (double *) malloc (dim * sizeof (double)); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->xmid = (double *) malloc (dim * sizeof (double)); + + if (s->xmid == 0) + { + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for xmid", GSL_ENOMEM, 0); + } + + s->sigma_l = (double *) malloc (dim * sizeof (double)); + + if (s->sigma_l == 0) + { + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for sigma_l", GSL_ENOMEM, 0); + } + + s->sigma_r = (double *) malloc (dim * sizeof (double)); + + if (s->sigma_r == 0) + { + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for sigma_r", GSL_ENOMEM, 0); + } + + s->fmax_l = (double *) malloc (dim * sizeof (double)); + + if (s->fmax_l == 0) + { + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fmax_l", GSL_ENOMEM, 0); + } + + s->fmax_r = (double *) malloc (dim * sizeof (double)); + + if (s->fmax_r == 0) + { + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fmax_r", GSL_ENOMEM, 0); + } + + s->fmin_l = (double *) malloc (dim * sizeof (double)); + + if (s->fmin_l == 0) + { + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fmin_l", GSL_ENOMEM, 0); + } + + s->fmin_r = (double *) malloc (dim * sizeof (double)); + + if (s->fmin_r == 0) + { + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fmin_r", GSL_ENOMEM, 0); + } + + s->fsum_l = (double *) malloc (dim * sizeof (double)); + + if (s->fsum_l == 0) + { + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum_l", GSL_ENOMEM, 0); + } + + s->fsum_r = (double *) malloc (dim * sizeof (double)); + + if (s->fsum_r == 0) + { + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum_r", GSL_ENOMEM, 0); + } + + s->fsum2_l = (double *) malloc (dim * sizeof (double)); + + if (s->fsum2_l == 0) + { + free (s->fsum_r); + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum2_l", GSL_ENOMEM, 0); + } + + s->fsum2_r = (double *) malloc (dim * sizeof (double)); + + if (s->fsum2_r == 0) + { + free (s->fsum2_l); + free (s->fsum_r); + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0); + } + + + s->hits_r = (size_t *) malloc (dim * sizeof (size_t)); + + if (s->hits_r == 0) + { + free (s->fsum2_r); + free (s->fsum2_l); + free (s->fsum_r); + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0); + } + + s->hits_l = (size_t *) malloc (dim * sizeof (size_t)); + + if (s->hits_l == 0) + { + free (s->hits_r); + free (s->fsum2_r); + free (s->fsum2_l); + free (s->fsum_r); + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for fsum2_r", GSL_ENOMEM, 0); + } + + s->dim = dim; + + gsl_monte_miser_init (s); + + return s; +} + +int +gsl_monte_miser_init (gsl_monte_miser_state * s) +{ + /* We use 8 points in each halfspace to estimate the variance. There are + 2*dim halfspaces. A variance estimate requires a minimum of 2 points. */ + s->min_calls = 16 * s->dim; + s->min_calls_per_bisection = 32 * s->min_calls; + s->estimate_frac = 0.1; + s->alpha = 2.0; + s->dither = 0.0; + + return GSL_SUCCESS; +} + +void +gsl_monte_miser_free (gsl_monte_miser_state * s) +{ + free (s->hits_r); + free (s->hits_l); + free (s->fsum2_r); + free (s->fsum2_l); + free (s->fsum_r); + free (s->fsum_l); + free (s->fmin_r); + free (s->fmin_l); + free (s->fmax_r); + free (s->fmax_l); + free (s->sigma_r); + free (s->sigma_l); + free (s->xmid); + free (s->x); + free (s); +} + +static int +estimate_corrmc (gsl_monte_function * f, + const double xl[], const double xu[], + size_t dim, size_t calls, + gsl_rng * r, + gsl_monte_miser_state * state, + double *result, double *abserr, + const double xmid[], double sigma_l[], double sigma_r[]) +{ + size_t i, n; + + double *x = state->x; + double *fsum_l = state->fsum_l; + double *fsum_r = state->fsum_r; + double *fsum2_l = state->fsum2_l; + double *fsum2_r = state->fsum2_r; + size_t *hits_l = state->hits_l; + size_t *hits_r = state->hits_r; + + double m = 0.0, q = 0.0; + double vol = 1.0; + + for (i = 0; i < dim; i++) + { + vol *= xu[i] - xl[i]; + hits_l[i] = hits_r[i] = 0; + fsum_l[i] = fsum_r[i] = 0.0; + fsum2_l[i] = fsum2_r[i] = 0.0; + sigma_l[i] = sigma_r[i] = -1; + } + + for (n = 0; n < calls; n++) + { + double fval; + + unsigned int j = (n/2) % dim; + unsigned int side = (n % 2); + + for (i = 0; i < dim; i++) + { + double z = gsl_rng_uniform_pos (r) ; + + if (i != j) + { + x[i] = xl[i] + z * (xu[i] - xl[i]); + } + else + { + if (side == 0) + { + x[i] = xmid[i] + z * (xu[i] - xmid[i]); + } + else + { + x[i] = xl[i] + z * (xmid[i] - xl[i]); + } + } + } + + fval = GSL_MONTE_FN_EVAL (f, x); + + /* recurrence for mean and variance */ + { + double d = fval - m; + m += d / (n + 1.0); + q += d * d * (n / (n + 1.0)); + } + + /* compute the variances on each side of the bisection */ + for (i = 0; i < dim; i++) + { + if (x[i] <= xmid[i]) + { + fsum_l[i] += fval; + fsum2_l[i] += fval * fval; + hits_l[i]++; + } + else + { + fsum_r[i] += fval; + fsum2_r[i] += fval * fval; + hits_r[i]++; + } + } + } + + for (i = 0; i < dim; i++) + { + double fraction_l = (xmid[i] - xl[i]) / (xu[i] - xl[i]); + + if (hits_l[i] > 0) + { + fsum_l[i] /= hits_l[i]; + sigma_l[i] = sqrt (fsum2_l[i] - fsum_l[i] * fsum_l[i] / hits_l[i]); + sigma_l[i] *= fraction_l * vol / hits_l[i]; + } + + if (hits_r[i] > 0) + { + fsum_r[i] /= hits_r[i]; + sigma_r[i] = sqrt (fsum2_r[i] - fsum_r[i] * fsum_r[i] / hits_r[i]); + sigma_r[i] *= (1 - fraction_l) * vol / hits_r[i]; + } + } + + *result = vol * m; + + if (calls < 2) + { + *abserr = GSL_POSINF; + } + else + { + *abserr = vol * sqrt (q / (calls * (calls - 1.0))); + } + + return GSL_SUCCESS; +} + diff --git a/gsl-1.9/monte/plain.c b/gsl-1.9/monte/plain.c new file mode 100644 index 0000000..cab3ae2 --- /dev/null +++ b/gsl-1.9/monte/plain.c @@ -0,0 +1,151 @@ +/* monte/plain.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Plain Monte-Carlo. */ + +/* Author: MJB */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte_plain.h> + +int +gsl_monte_plain_integrate (const gsl_monte_function * f, + const double xl[], const double xu[], + const size_t dim, + const size_t calls, + gsl_rng * r, + gsl_monte_plain_state * state, + double *result, double *abserr) +{ + double vol, m = 0, q = 0; + double *x = state->x; + size_t n, i; + + if (dim != state->dim) + { + GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL); + } + + for (i = 0; i < dim; i++) + { + if (xu[i] <= xl[i]) + { + GSL_ERROR ("xu must be greater than xl", GSL_EINVAL); + } + + if (xu[i] - xl[i] > GSL_DBL_MAX) + { + GSL_ERROR ("Range of integration is too large, please rescale", + GSL_EINVAL); + } + } + + /* Compute the volume of the region */ + + vol = 1; + + for (i = 0; i < dim; i++) + { + vol *= xu[i] - xl[i]; + } + + for (n = 0; n < calls; n++) + { + /* Choose a random point in the integration region */ + + for (i = 0; i < dim; i++) + { + x[i] = xl[i] + gsl_rng_uniform_pos (r) * (xu[i] - xl[i]); + } + + { + double fval = GSL_MONTE_FN_EVAL (f, x); + + /* recurrence for mean and variance */ + + double d = fval - m; + m += d / (n + 1.0); + q += d * d * (n / (n + 1.0)); + } + } + + *result = vol * m; + + if (calls < 2) + { + *abserr = GSL_POSINF; + } + else + { + *abserr = vol * sqrt (q / (calls * (calls - 1.0))); + } + + return GSL_SUCCESS; +} + +gsl_monte_plain_state * +gsl_monte_plain_alloc (size_t dim) +{ + gsl_monte_plain_state *s = + (gsl_monte_plain_state *) malloc (sizeof (gsl_monte_plain_state)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for state struct", + GSL_ENOMEM, 0); + } + + s->x = (double *) malloc (dim * sizeof (double)); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for working vector", + GSL_ENOMEM, 0); + } + + s->dim = dim; + + return s; +} + +/* Set some default values and whatever */ + +int +gsl_monte_plain_init (gsl_monte_plain_state * s) +{ + size_t i; + + for (i = 0; i < s->dim; i++) + { + s->x[i] = 0.0; + } + + return GSL_SUCCESS; +} + +void +gsl_monte_plain_free (gsl_monte_plain_state * s) +{ + free (s->x); + free (s); +} diff --git a/gsl-1.9/monte/test.c b/gsl-1.9/monte/test.c new file mode 100644 index 0000000..7eb0a7a --- /dev/null +++ b/gsl-1.9/monte/test.c @@ -0,0 +1,373 @@ +/* monte/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdlib.h> +#include <math.h> +#include <stdio.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte_plain.h> +#include <gsl/gsl_monte_miser.h> +#include <gsl/gsl_monte_vegas.h> + +#define CONSTANT +#define PRODUCT +#define GAUSSIAN +#define DBLGAUSSIAN +#define TSUDA + +#define PLAIN +#define MISER +#define VEGAS + +double xl[11] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; +double xu[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }; +double xu2[11] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 }; +double xu3[2] = { GSL_DBL_MAX, GSL_DBL_MAX }; + +double fconst (double x[], size_t d, void *params); +double f0 (double x[], size_t d, void *params); +double f1 (double x[], size_t d, void *params); +double f2 (double x[], size_t d, void *params); +double f3 (double x[], size_t d, void *params); + +void my_error_handler (const char *reason, const char *file, + int line, int err); + +struct problem { + gsl_monte_function * f; + double * xl; + double * xu; + size_t dim; + size_t calls; + double expected_result; + double expected_error; + char * description; +} ; + +gsl_monte_function +make_function (double (*f)(double *, size_t, void *), size_t d, void * p); + +gsl_monte_function +make_function (double (*f)(double *, size_t, void *), size_t d, void * p) +{ + gsl_monte_function f_new; + + f_new.f = f; + f_new.dim = d; + f_new.params = p; + + return f_new; +} + + +void +add (struct problem * problems, int * n, + gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls, + double result, double err, char * description); + +void +add (struct problem * problems, int * n, + gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls, + double result, double err, char * description) +{ + int i = *n; + problems[i].f = f; + problems[i].xl = xl; + problems[i].xu = xu; + problems[i].dim = dim; + problems[i].calls = calls; + problems[i].expected_result = result; + problems[i].expected_error = err; + problems[i].description = description; + (*n)++; +} + +#define TRIALS 10 + +int +main (void) +{ + double result[TRIALS], error[TRIALS]; + double a = 0.1; + double c = (1.0 + sqrt (10.0)) / 9.0; + + gsl_monte_function Fc = make_function(&fconst, 0, 0); + gsl_monte_function F0 = make_function(&f0, 0, &a); + gsl_monte_function F1 = make_function(&f1, 0, &a); + gsl_monte_function F2 = make_function(&f2, 0, &a); + gsl_monte_function F3 = make_function(&f3, 0, &c); + + /* The relationship between the variance of the function itself, the + error on the integral and the number of calls is, + + sigma = sqrt(variance/N) + + where the variance is the <(f - <f>)^2> where <.> denotes the + volume average (integral over the integration region divided by + the volume) */ + + int n = 0; + struct problem * I; + struct problem problems[256]; + +#ifdef CONSTANT + /* variance(Fc) = 0 */ + + add(problems,&n, &Fc, xl, xu, 1, 1000, 1.0, 0.0, "constant, 1d"); + add(problems,&n, &Fc, xl, xu, 2, 1000, 1.0, 0.0, "constant, 2d"); + add(problems,&n, &Fc, xl, xu, 3, 1000, 1.0, 0.0, "constant, 3d"); + add(problems,&n, &Fc, xl, xu, 4, 1000, 1.0, 0.0, "constant, 4d"); + add(problems,&n, &Fc, xl, xu, 5, 1000, 1.0, 0.0, "constant, 5d"); + add(problems,&n, &Fc, xl, xu, 6, 1000, 1.0, 0.0, "constant, 6d"); + add(problems,&n, &Fc, xl, xu, 7, 1000, 1.0, 0.0, "constant, 7d"); + add(problems,&n, &Fc, xl, xu, 8, 1000, 1.0, 0.0, "constant, 8d"); + add(problems,&n, &Fc, xl, xu, 9, 1000, 1.0, 0.0, "constant, 9d"); + add(problems,&n, &Fc, xl, xu, 10, 1000, 1.0, 0.0, "constant, 10d"); +#endif + +#ifdef PRODUCT + /* variance(F0) = (4/3)^d - 1 */ + + add(problems,&n, &F0, xl, xu, 1, 3333, 1.0, 0.01, "product, 1d" ); + add(problems,&n, &F0, xl, xu, 2, 7777, 1.0, 0.01, "product, 2d" ); + add(problems,&n, &F0, xl, xu, 3, 13703, 1.0, 0.01, "product, 3d" ); + add(problems,&n, &F0, xl, xu, 4, 21604, 1.0, 0.01, "product, 4d" ); + add(problems,&n, &F0, xl, xu, 5, 32139, 1.0, 0.01, "product, 5d" ); + add(problems,&n, &F0, xl, xu, 6, 46186, 1.0, 0.01, "product, 6d" ); + add(problems,&n, &F0, xl, xu, 7, 64915, 1.0, 0.01, "product, 7d" ); + add(problems,&n, &F0, xl, xu, 8, 89887, 1.0, 0.01, "product, 8d" ); + add(problems,&n, &F0, xl, xu, 9, 123182, 1.0, 0.01, "product, 9d" ); + add(problems,&n, &F0, xl, xu, 10, 167577, 1.0, 0.01, "product, 10d" ); +#endif + +#ifdef GAUSSIAN + /* variance(F1) = (1/(a sqrt(2 pi)))^d - 1 */ + + add(problems,&n, &F1, xl, xu, 1, 298, 1.0, 0.1, "gaussian, 1d" ); + add(problems,&n, &F1, xl, xu, 2, 1492, 1.0, 0.1, "gaussian, 2d" ); + add(problems,&n, &F1, xl, xu, 3, 6249, 1.0, 0.1, "gaussian, 3d" ); + add(problems,&n, &F1, xl, xu, 4, 25230, 1.0, 0.1, "gaussian, 4d" ); + add(problems,&n, &F1, xl, xu, 5, 100953, 1.0, 0.1, "gaussian, 5d" ); + add(problems,&n, &F1, xl, xu, 6, 44782, 1.0, 0.3, "gaussian, 6d" ); + add(problems,&n, &F1, xl, xu, 7, 178690, 1.0, 0.3, "gaussian, 7d" ); + add(problems,&n, &F1, xl, xu, 8, 712904, 1.0, 0.3, "gaussian, 8d" ); + add(problems,&n, &F1, xl, xu, 9, 2844109, 1.0, 0.3, "gaussian, 9d" ); + add(problems,&n, &F1, xl, xu, 10, 11346390, 1.0, 0.3, "gaussian, 10d" ); +#endif + +#ifdef DBLGAUSSIAN + /* variance(F2) = 0.5 * (1/(a sqrt(2 pi)))^d - 1 */ + + add(problems,&n, &F2, xl, xu, 1, 9947, 1.0, 0.01, "double gaussian, 1d" ); + add(problems,&n, &F2, xl, xu, 2, 695, 1.0, 0.1, "double gaussian, 2d" ); + add(problems,&n, &F2, xl, xu, 3, 3074, 1.0, 0.1, "double gaussian, 3d" ); + add(problems,&n, &F2, xl, xu, 4, 12565, 1.0, 0.1, "double gaussian, 4d" ); + add(problems,&n, &F2, xl, xu, 5, 50426, 1.0, 0.1, "double gaussian, 5d" ); + add(problems,&n, &F2, xl, xu, 6, 201472, 1.0, 0.1, "double gaussian, 6d" ); + add(problems,&n, &F2, xl, xu, 7, 804056, 1.0, 0.1, "double gaussian, 7d" ); + add(problems,&n, &F2, xl, xu, 8, 356446, 1.0, 0.3, "double gaussian, 8d" ); + add(problems,&n, &F2, xl, xu, 9, 1422049, 1.0, 0.3, "double gaussian, 9d" ); + add(problems,&n, &F2, xl, xu, 10, 5673189, 1.0, 0.3, "double gaussian, 10d" ); +#endif + +#ifdef TSUDA + /* variance(F3) = ((c^2 + c + 1/3)/(c(c+1)))^d - 1 */ + + add(problems,&n, &F3, xl, xu, 1, 4928, 1.0, 0.01, "tsuda function, 1d" ); + add(problems,&n, &F3, xl, xu, 2, 12285, 1.0, 0.01, "tsuda function, 2d" ); + add(problems,&n, &F3, xl, xu, 3, 23268, 1.0, 0.01, "tsuda function, 3d" ); + add(problems,&n, &F3, xl, xu, 4, 39664, 1.0, 0.01, "tsuda function, 4d" ); + add(problems,&n, &F3, xl, xu, 5, 64141, 1.0, 0.01, "tsuda function, 5d" ); + add(problems,&n, &F3, xl, xu, 6, 100680, 1.0, 0.01, "tsuda function, 6d" ); + add(problems,&n, &F3, xl, xu, 7, 155227, 1.0, 0.01, "tsuda function, 7d" ); + add(problems,&n, &F3, xl, xu, 8, 236657, 1.0, 0.01, "tsuda function, 8d" ); + add(problems,&n, &F3, xl, xu, 9, 358219, 1.0, 0.01, "tsuda function, 9d" ); + add(problems,&n, &F3, xl, xu, 10, 539690, 1.0, 0.01, "tsuda function, 10d" ); +#endif + + add(problems,&n, 0, 0, 0, 0, 0, 0, 0, 0 ); + + + /* gsl_set_error_handler (&my_error_handler); */ + gsl_ieee_env_setup (); + gsl_rng_env_setup (); + +#ifdef A + printf ("testing allocation/input checks\n"); + + status = gsl_monte_plain_validate (s, xl, xu3, 1, 1); + gsl_test (status != 0, "error if limits too large"); + status = gsl_monte_plain_validate (s, xl, xu, 0, 10); + gsl_test (status != 0, "error if num_dim = 0"); + status = gsl_monte_plain_validate (s, xl, xu, 1, 0); + gsl_test (status != 0, "error if calls = 0"); + status = gsl_monte_plain_validate (s, xu, xl, 1, 10); + gsl_test (status != 0, "error if xu < xl"); +#endif + +#ifdef PLAIN +#define NAME "plain" +#define MONTE_STATE gsl_monte_plain_state +#define MONTE_ALLOC gsl_monte_plain_alloc +#define MONTE_INTEGRATE gsl_monte_plain_integrate +#define MONTE_FREE gsl_monte_plain_free +#define MONTE_SPEEDUP 1 +#define MONTE_ERROR_TEST(err,expected) gsl_test_factor(err,expected, 5.0, NAME ", %s, abserr[%d]", I->description, i) +#include "test_main.c" +#undef NAME +#undef MONTE_STATE +#undef MONTE_ALLOC +#undef MONTE_INTEGRATE +#undef MONTE_FREE +#undef MONTE_ERROR_TEST +#undef MONTE_SPEEDUP +#endif + +#ifdef MISER +#define NAME "miser" +#define MONTE_STATE gsl_monte_miser_state +#define MONTE_ALLOC gsl_monte_miser_alloc +#define MONTE_INTEGRATE gsl_monte_miser_integrate +#define MONTE_FREE gsl_monte_miser_free +#define MONTE_SPEEDUP 2 +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 5.0 * expected, NAME ", %s, abserr[%d] (obs %g vs plain %g)", I->description, i, err, expected) +#include "test_main.c" +#undef NAME +#undef MONTE_STATE +#undef MONTE_ALLOC +#undef MONTE_INTEGRATE +#undef MONTE_FREE +#undef MONTE_ERROR_TEST +#undef MONTE_SPEEDUP +#endif + +#ifdef VEGAS +#define NAME "vegas" +#define MONTE_STATE gsl_monte_vegas_state +#define MONTE_ALLOC gsl_monte_vegas_alloc +#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; if (s->chisq < 0.5 || s->chisq > 2) gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err); } +#define MONTE_FREE gsl_monte_vegas_free +#define MONTE_SPEEDUP 3 +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected) +#include "test_main.c" +#undef NAME +#undef MONTE_STATE +#undef MONTE_ALLOC +#undef MONTE_INTEGRATE +#undef MONTE_FREE +#undef MONTE_ERROR_TEST +#undef MONTE_SPEEDUP +#endif + + exit (gsl_test_summary ()); +} + +/* Simple constant function */ +double +fconst (double x[], size_t num_dim, void *params) +{ + return 1; +} + +/* Simple product function */ +double +f0 (double x[], size_t num_dim, void *params) +{ + double prod = 1.0; + unsigned int i; + + for (i = 0; i < num_dim; ++i) + { + prod *= 2.0 * x[i]; + } + + return prod; +} + +/* Gaussian centered at 1/2. */ + +double +f1 (double x[], size_t num_dim, void *params) +{ + double a = *(double *)params; + double sum = 0.; + + unsigned int i; + for (i = 0; i < num_dim; i++) + { + double dx = x[i] - 0.5; + sum += dx * dx; + } + return (pow (M_2_SQRTPI / (2. * a), (double) num_dim) * + exp (-sum / (a * a))); +} + +/* double gaussian */ +double +f2 (double x[], size_t num_dim, void *params) +{ + double a = *(double *)params; + double sum1 = 0.; + double sum2 = 0.; + + unsigned int i; + for (i = 0; i < num_dim; i++) + { + double dx1 = x[i] - 1. / 3.; + double dx2 = x[i] - 2. / 3.; + sum1 += dx1 * dx1; + sum2 += dx2 * dx2; + } + return 0.5 * pow (M_2_SQRTPI / (2. * a), num_dim) + * (exp (-sum1 / (a * a)) + exp (-sum2 / (a * a))); +} + +/* Tsuda's example */ +double +f3 (double x[], size_t num_dim, void *params) +{ + double c = *(double *)params; + + double prod = 1.; + + unsigned int i; + + for (i = 0; i < num_dim; i++) + { + prod *= c / (c + 1) * pow((c + 1) / (c + x[i]), 2.0); + } + + return prod; +} + + +void +my_error_handler (const char *reason, const char *file, int line, int err) +{ + if (0) + printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err); +} diff --git a/gsl-1.9/monte/test_main.c b/gsl-1.9/monte/test_main.c new file mode 100644 index 0000000..d35f357 --- /dev/null +++ b/gsl-1.9/monte/test_main.c @@ -0,0 +1,70 @@ +for (I = problems ; I->f != 0; I++) +{ + size_t i; + double sum = 0, mean, sumd2 = 0, sd, res, err; + + gsl_rng * r; + + if (I->dim > 3) + { + continue ; + } + + r = gsl_rng_alloc (gsl_rng_default); + + for (i = 0; i < TRIALS ; i++) + { + MONTE_STATE *s = MONTE_ALLOC (I->dim); + + I->f->dim = I->dim; + + MONTE_INTEGRATE (I->f, I->xl, I->xu, + I->dim, I->calls / MONTE_SPEEDUP, r, s, + &res, &err); + + gsl_test_abs (res, I->expected_result, + 5 * GSL_MAX(err, 1024*GSL_DBL_EPSILON), + NAME ", %s, result[%d]", I->description, i); + + MONTE_ERROR_TEST (err, I->expected_error); + + result[i] = res; + error[i] = err; + + MONTE_FREE (s); + } + + for (i = 0; i < TRIALS; i++) + { + sum += result[i]; + } + + mean = sum / TRIALS ; + + for (i = 0; i < TRIALS; i++) + { + sumd2 += pow(result[i] - mean, 2.0); + } + + sd = sqrt(sumd2 / (TRIALS-1.0)) ; + + if (sd < TRIALS * GSL_DBL_EPSILON * fabs (mean)) + { + sd = 0; + } + + for (i = 0; i < TRIALS; i++) + { + if (sd == 0 && fabs(error[i]) < GSL_DBL_EPSILON * fabs(result[i])) + { + error[i] = 0.0; + } + + gsl_test_factor (error[i], sd, 5.0, + NAME ", %s, abserr[%d] vs sd", I->description, i); + } + + + gsl_rng_free (r); +} + diff --git a/gsl-1.9/monte/vegas.c b/gsl-1.9/monte/vegas.c new file mode 100644 index 0000000..1a3eb21 --- /dev/null +++ b/gsl-1.9/monte/vegas.c @@ -0,0 +1,864 @@ +/* monte/vegas.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth + * + * This program 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 2 of the License, or (at + * your option) any later version. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: MJB */ +/* Modified by: Brian Gough, 12/2000 */ + +/* This is an implementation of the adaptive Monte-Carlo algorithm + of G. P. Lepage, originally described in J. Comp. Phys. 27, 192(1978). + The current version of the algorithm was described in the Cornell + preprint CLNS-80/447 of March, 1980. + + This code follows most closely the c version by D.R.Yennie, coded + in 1984. + + The input coordinates are x[j], with upper and lower limits xu[j] + and xl[j]. The integration length in the j-th direction is + delx[j]. Each coordinate x[j] is rescaled to a variable y[j] in + the range 0 to 1. The range is divided into bins with boundaries + xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1. The grid + is refined (ie, bins are adjusted) using d[i][j] which is some + variation on the squared sum. A third parameter used in defining + the real coordinate using random numbers is called z. It ranges + from 0 to bins. Its integer part gives the lower index of the bin + into which a call is to be placed, and the remainder gives the + location inside the bin. + + When stratified sampling is used the bins are grouped into boxes, + and the algorithm allocates an equal number of function calls to + each box. + + The variable alpha controls how "stiff" the rebinning algorithm is. + alpha = 0 means never change the grid. Alpha is typically set between + 1 and 2. + + */ + +/* configuration headers */ +#include <config.h> + +/* standard headers */ +#include <math.h> +#include <stdio.h> + +/* gsl headers */ +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_monte_vegas.h> + +/* lib-specific headers */ +#define BINS_MAX 50 /* even integer, will be divided by two */ + +/* A separable grid with coordinates and values */ +#define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)]) +#define NEW_COORD(s,i) ((s)->xin[(i)]) +#define VALUE(s,i,j) ((s)->d[(i)*(s)->dim + (j)]) + +/* predeclare functions */ + +typedef int coord; + +static void init_grid (gsl_monte_vegas_state * s, double xl[], double xu[], + size_t dim); +static void reset_grid_values (gsl_monte_vegas_state * s); +static void init_box_coord (gsl_monte_vegas_state * s, coord box[]); +static int change_box_coord (gsl_monte_vegas_state * s, coord box[]); +static void accumulate_distribution (gsl_monte_vegas_state * s, coord bin[], + double y); +static void random_point (double x[], coord bin[], double *bin_vol, + const coord box[], + const double xl[], const double xu[], + gsl_monte_vegas_state * s, gsl_rng * r); +static void resize_grid (gsl_monte_vegas_state * s, unsigned int bins); +static void refine_grid (gsl_monte_vegas_state * s); + +static void print_lim (gsl_monte_vegas_state * state, + double xl[], double xu[], unsigned long dim); +static void print_head (gsl_monte_vegas_state * state, + unsigned long num_dim, unsigned long calls, + unsigned int it_num, + unsigned int bins, unsigned int boxes); +static void print_res (gsl_monte_vegas_state * state, + unsigned int itr, double res, double err, + double cum_res, double cum_err, + double chi_sq); +static void print_dist (gsl_monte_vegas_state * state, unsigned long dim); +static void print_grid (gsl_monte_vegas_state * state, unsigned long dim); + +int +gsl_monte_vegas_integrate (gsl_monte_function * f, + double xl[], double xu[], + size_t dim, size_t calls, + gsl_rng * r, + gsl_monte_vegas_state * state, + double *result, double *abserr) +{ + double cum_int, cum_sig; + size_t i, k, it; + + if (dim != state->dim) + { + GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL); + } + + for (i = 0; i < dim; i++) + { + if (xu[i] <= xl[i]) + { + GSL_ERROR ("xu must be greater than xl", GSL_EINVAL); + } + + if (xu[i] - xl[i] > GSL_DBL_MAX) + { + GSL_ERROR ("Range of integration is too large, please rescale", + GSL_EINVAL); + } + } + + if (state->stage == 0) + { + init_grid (state, xl, xu, dim); + + if (state->verbose >= 0) + { + print_lim (state, xl, xu, dim); + } + } + + if (state->stage <= 1) + { + state->wtd_int_sum = 0; + state->sum_wgts = 0; + state->chi_sum = 0; + state->it_num = 1; + state->samples = 0; + } + + if (state->stage <= 2) + { + unsigned int bins = state->bins_max; + unsigned int boxes = 1; + + if (state->mode != GSL_VEGAS_MODE_IMPORTANCE_ONLY) + { + /* shooting for 2 calls/box */ + + boxes = floor (pow (calls / 2.0, 1.0 / dim)); + state->mode = GSL_VEGAS_MODE_IMPORTANCE; + + if (2 * boxes >= state->bins_max) + { + /* if bins/box < 2 */ + int box_per_bin = GSL_MAX (boxes / state->bins_max, 1); + + bins = GSL_MIN(boxes / box_per_bin, state->bins_max); + boxes = box_per_bin * bins; + + state->mode = GSL_VEGAS_MODE_STRATIFIED; + } + } + + { + double tot_boxes = pow ((double) boxes, (double) dim); + state->calls_per_box = GSL_MAX (calls / tot_boxes, 2); + calls = state->calls_per_box * tot_boxes; + } + + /* total volume of x-space/(avg num of calls/bin) */ + state->jac = state->vol * pow ((double) bins, (double) dim) / calls; + + state->boxes = boxes; + + /* If the number of bins changes from the previous invocation, bins + are expanded or contracted accordingly, while preserving bin + density */ + + if (bins != state->bins) + { + resize_grid (state, bins); + + if (state->verbose > 1) + { + print_grid (state, dim); + } + } + + if (state->verbose >= 0) + { + print_head (state, + dim, calls, state->it_num, state->bins, state->boxes); + } + } + + state->it_start = state->it_num; + + cum_int = 0.0; + cum_sig = 0.0; + + state->chisq = 0.0; + + for (it = 0; it < state->iterations; it++) + { + double intgrl = 0.0, intgrl_sq = 0.0; + double sig = 0.0; + double wgt; + size_t calls_per_box = state->calls_per_box; + double jacbin = state->jac; + double *x = state->x; + coord *bin = state->bin; + + state->it_num = state->it_start + it; + + reset_grid_values (state); + init_box_coord (state, state->box); + + do + { + double m = 0, q = 0; + double f_sq_sum = 0.0; + + for (k = 0; k < calls_per_box; k++) + { + double fval, bin_vol; + + random_point (x, bin, &bin_vol, state->box, xl, xu, state, r); + + fval = jacbin * bin_vol * GSL_MONTE_FN_EVAL (f, x); + + /* recurrence for mean and variance */ + + { + double d = fval - m; + m += d / (k + 1.0); + q += d * d * (k / (k + 1.0)); + } + + if (state->mode != GSL_VEGAS_MODE_STRATIFIED) + { + double f_sq = fval * fval; + accumulate_distribution (state, bin, f_sq); + } + } + + intgrl += m * calls_per_box; + + f_sq_sum = q * calls_per_box ; + + sig += f_sq_sum ; + + if (state->mode == GSL_VEGAS_MODE_STRATIFIED) + { + accumulate_distribution (state, bin, f_sq_sum); + } + } + while (change_box_coord (state, state->box)); + + /* Compute final results for this iteration */ + + sig = sig / (calls_per_box - 1.0) ; + + if (sig > 0) + { + wgt = 1.0 / sig; + } + else if (state->sum_wgts > 0) + { + wgt = state->sum_wgts / state->samples; + } + else + { + wgt = 0.0; + } + + intgrl_sq = intgrl * intgrl; + + state->result = intgrl; + state->sigma = sqrt(sig); + + if (wgt > 0.0) + { + state->samples++ ; + state->sum_wgts += wgt; + state->wtd_int_sum += intgrl * wgt; + state->chi_sum += intgrl_sq * wgt; + + cum_int = state->wtd_int_sum / state->sum_wgts; + cum_sig = sqrt (1 / state->sum_wgts); + + if (state->samples > 1) + { + state->chisq = (state->chi_sum - state->wtd_int_sum * cum_int) / + (state->samples - 1.0); + } + } + else + { + cum_int += (intgrl - cum_int) / (it + 1.0); + cum_sig = 0.0; + } + + + if (state->verbose >= 0) + { + print_res (state, + state->it_num, intgrl, sqrt (sig), cum_int, cum_sig, + state->chisq); + if (it + 1 == state->iterations && state->verbose > 0) + { + print_grid (state, dim); + } + } + + if (state->verbose > 1) + { + print_dist (state, dim); + } + + refine_grid (state); + + if (state->verbose > 1) + { + print_grid (state, dim); + } + + } + + /* By setting stage to 1 further calls will generate independent + estimates based on the same grid, although it may be rebinned. */ + + state->stage = 1; + + *result = cum_int; + *abserr = cum_sig; + + return GSL_SUCCESS; +} + + + +gsl_monte_vegas_state * +gsl_monte_vegas_alloc (size_t dim) +{ + gsl_monte_vegas_state *s = + (gsl_monte_vegas_state *) malloc (sizeof (gsl_monte_vegas_state)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for vegas state struct", + GSL_ENOMEM, 0); + } + + s->delx = (double *) malloc (dim * sizeof (double)); + + if (s->delx == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for delx", GSL_ENOMEM, 0); + } + + s->d = (double *) malloc (BINS_MAX * dim * sizeof (double)); + + if (s->d == 0) + { + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0); + } + + s->xi = (double *) malloc ((BINS_MAX + 1) * dim * sizeof (double)); + + if (s->xi == 0) + { + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for xi", GSL_ENOMEM, 0); + } + + s->xin = (double *) malloc ((BINS_MAX + 1) * sizeof (double)); + + if (s->xin == 0) + { + free (s->xi); + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0); + } + + s->weight = (double *) malloc (BINS_MAX * sizeof (double)); + + if (s->weight == 0) + { + free (s->xin); + free (s->xi); + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0); + } + + s->box = (coord *) malloc (dim * sizeof (coord)); + + if (s->box == 0) + { + free (s->weight); + free (s->xin); + free (s->xi); + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for box", GSL_ENOMEM, 0); + } + + s->bin = (coord *) malloc (dim * sizeof (coord)); + + if (s->bin == 0) + { + free (s->box); + free (s->weight); + free (s->xin); + free (s->xi); + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for bin", GSL_ENOMEM, 0); + } + + s->x = (double *) malloc (dim * sizeof (double)); + + if (s->x == 0) + { + free (s->bin); + free (s->box); + free (s->weight); + free (s->xin); + free (s->xi); + free (s->d); + free (s->delx); + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->dim = dim; + s->bins_max = BINS_MAX; + + gsl_monte_vegas_init (s); + + return s; +} + +/* Set some default values and whatever */ +int +gsl_monte_vegas_init (gsl_monte_vegas_state * state) +{ + state->stage = 0; + state->alpha = 1.5; + state->verbose = -1; + state->iterations = 5; + state->mode = GSL_VEGAS_MODE_IMPORTANCE; + state->chisq = 0; + state->bins = state->bins_max; + state->ostream = stdout; + + return GSL_SUCCESS; +} + +void +gsl_monte_vegas_free (gsl_monte_vegas_state * s) +{ + free (s->x); + free (s->delx); + free (s->d); + free (s->xi); + free (s->xin); + free (s->weight); + free (s->box); + free (s->bin); + free (s); +} + +static void +init_box_coord (gsl_monte_vegas_state * s, coord box[]) +{ + size_t i; + + size_t dim = s->dim; + + for (i = 0; i < dim; i++) + { + box[i] = 0; + } +} + +/* change_box_coord steps through the box coord like + {0,0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, ... +*/ +static int +change_box_coord (gsl_monte_vegas_state * s, coord box[]) +{ + int j = s->dim - 1; + + int ng = s->boxes; + + while (j >= 0) + { + box[j] = (box[j] + 1) % ng; + + if (box[j] != 0) + { + return 1; + } + + j--; + } + + return 0; +} + +static void +init_grid (gsl_monte_vegas_state * s, double xl[], double xu[], size_t dim) +{ + size_t j; + + double vol = 1.0; + + s->bins = 1; + + for (j = 0; j < dim; j++) + { + double dx = xu[j] - xl[j]; + s->delx[j] = dx; + vol *= dx; + + COORD (s, 0, j) = 0.0; + COORD (s, 1, j) = 1.0; + } + + s->vol = vol; +} + + +static void +reset_grid_values (gsl_monte_vegas_state * s) +{ + size_t i, j; + + size_t dim = s->dim; + size_t bins = s->bins; + + for (i = 0; i < bins; i++) + { + for (j = 0; j < dim; j++) + { + VALUE (s, i, j) = 0.0; + } + } +} + +static void +accumulate_distribution (gsl_monte_vegas_state * s, coord bin[], double y) +{ + size_t j; + size_t dim = s->dim; + + for (j = 0; j < dim; j++) + { + int i = bin[j]; + VALUE (s, i, j) += y; + } +} + +static void +random_point (double x[], coord bin[], double *bin_vol, + const coord box[], const double xl[], const double xu[], + gsl_monte_vegas_state * s, gsl_rng * r) +{ + /* Use the random number generator r to return a random position x + in a given box. The value of bin gives the bin location of the + random position (there may be several bins within a given box) */ + + double vol = 1.0; + + size_t j; + + size_t dim = s->dim; + size_t bins = s->bins; + size_t boxes = s->boxes; + + DISCARD_POINTER(xu); /* prevent warning about unused parameter */ + + for (j = 0; j < dim; ++j) + { + /* box[j] + ran gives the position in the box units, while z + is the position in bin units. */ + + double z = ((box[j] + gsl_rng_uniform_pos (r)) / boxes) * bins; + + int k = z; + + double y, bin_width; + + bin[j] = k; + + if (k == 0) + { + bin_width = COORD (s, 1, j); + y = z * bin_width; + } + else + { + bin_width = COORD (s, k + 1, j) - COORD (s, k, j); + y = COORD (s, k, j) + (z - k) * bin_width; + } + + x[j] = xl[j] + y * s->delx[j]; + + vol *= bin_width; + } + + *bin_vol = vol; +} + + +static void +resize_grid (gsl_monte_vegas_state * s, unsigned int bins) +{ + size_t j, k; + size_t dim = s->dim; + + /* weight is ratio of bin sizes */ + + double pts_per_bin = (double) s->bins / (double) bins; + + for (j = 0; j < dim; j++) + { + double xold; + double xnew = 0; + double dw = 0; + int i = 1; + + for (k = 1; k <= s->bins; k++) + { + dw += 1.0; + xold = xnew; + xnew = COORD (s, k, j); + + for (; dw > pts_per_bin; i++) + { + dw -= pts_per_bin; + NEW_COORD (s, i) = xnew - (xnew - xold) * dw; + } + } + + for (k = 1 ; k < bins; k++) + { + COORD(s, k, j) = NEW_COORD(s, k); + } + + COORD (s, bins, j) = 1; + } + + s->bins = bins; +} + +static void +refine_grid (gsl_monte_vegas_state * s) +{ + size_t i, j, k; + size_t dim = s->dim; + size_t bins = s->bins; + + for (j = 0; j < dim; j++) + { + double grid_tot_j, tot_weight; + double * weight = s->weight; + + double oldg = VALUE (s, 0, j); + double newg = VALUE (s, 1, j); + + VALUE (s, 0, j) = (oldg + newg) / 2; + grid_tot_j = VALUE (s, 0, j); + + /* This implements gs[i][j] = (gs[i-1][j]+gs[i][j]+gs[i+1][j])/3 */ + + for (i = 1; i < bins - 1; i++) + { + double rc = oldg + newg; + oldg = newg; + newg = VALUE (s, i + 1, j); + VALUE (s, i, j) = (rc + newg) / 3; + grid_tot_j += VALUE (s, i, j); + } + VALUE (s, bins - 1, j) = (newg + oldg) / 2; + + grid_tot_j += VALUE (s, bins - 1, j); + + tot_weight = 0; + + for (i = 0; i < bins; i++) + { + weight[i] = 0; + + if (VALUE (s, i, j) > 0) + { + oldg = grid_tot_j / VALUE (s, i, j); + /* damped change */ + weight[i] = pow (((oldg - 1) / oldg / log (oldg)), s->alpha); + } + + tot_weight += weight[i]; + +#ifdef DEBUG + printf("weight[%d] = %g\n", i, weight[i]); +#endif + } + + { + double pts_per_bin = tot_weight / bins; + + double xold; + double xnew = 0; + double dw = 0; + i = 1; + + for (k = 0; k < bins; k++) + { + dw += weight[k]; + xold = xnew; + xnew = COORD (s, k + 1, j); + + for (; dw > pts_per_bin; i++) + { + dw -= pts_per_bin; + NEW_COORD (s, i) = xnew - (xnew - xold) * dw / weight[k]; + } + } + + for (k = 1 ; k < bins ; k++) + { + COORD(s, k, j) = NEW_COORD(s, k); + } + + COORD (s, bins, j) = 1; + } + } +} + + +static void +print_lim (gsl_monte_vegas_state * state, + double xl[], double xu[], unsigned long dim) +{ + unsigned long j; + + fprintf (state->ostream, "The limits of integration are:\n"); + for (j = 0; j < dim; ++j) + fprintf (state->ostream, "\nxl[%lu]=%f xu[%lu]=%f", j, xl[j], j, xu[j]); + fprintf (state->ostream, "\n"); + fflush (state->ostream); +} + +static void +print_head (gsl_monte_vegas_state * state, + unsigned long num_dim, unsigned long calls, + unsigned int it_num, unsigned int bins, unsigned int boxes) +{ + fprintf (state->ostream, + "\nnum_dim=%lu, calls=%lu, it_num=%d, max_it_num=%d ", + num_dim, calls, it_num, state->iterations); + fprintf (state->ostream, + "verb=%d, alph=%.2f,\nmode=%d, bins=%d, boxes=%d\n", + state->verbose, state->alpha, state->mode, bins, boxes); + fprintf (state->ostream, + "\n single.......iteration "); + fprintf (state->ostream, "accumulated......results \n"); + + fprintf (state->ostream, + "iteration integral sigma integral "); + fprintf (state->ostream, " sigma chi-sq/it\n\n"); + fflush (state->ostream); + +} + +static void +print_res (gsl_monte_vegas_state * state, + unsigned int itr, + double res, double err, + double cum_res, double cum_err, + double chi_sq) +{ + fprintf (state->ostream, + "%4d %6.4e %10.2e %6.4e %8.2e %10.2e\n", + itr, res, err, cum_res, cum_err, chi_sq); + fflush (state->ostream); +} + +static void +print_dist (gsl_monte_vegas_state * state, unsigned long dim) +{ + unsigned long i, j; + int p = state->verbose; + if (p < 1) + return; + + for (j = 0; j < dim; ++j) + { + fprintf (state->ostream, "\n axis %lu \n", j); + fprintf (state->ostream, " x g\n"); + for (i = 0; i < state->bins; i++) + { + fprintf (state->ostream, "weight [%11.2e , %11.2e] = ", + COORD (state, i, j), COORD(state,i+1,j)); + fprintf (state->ostream, " %11.2e\n", VALUE (state, i, j)); + + } + fprintf (state->ostream, "\n"); + } + fprintf (state->ostream, "\n"); + fflush (state->ostream); + +} + +static void +print_grid (gsl_monte_vegas_state * state, unsigned long dim) +{ + unsigned long i, j; + int p = state->verbose; + if (p < 1) + return; + + for (j = 0; j < dim; ++j) + { + fprintf (state->ostream, "\n axis %lu \n", j); + fprintf (state->ostream, " x \n"); + for (i = 0; i <= state->bins; i++) + { + fprintf (state->ostream, "%11.2e", COORD (state, i, j)); + if (i % 5 == 4) + fprintf (state->ostream, "\n"); + } + fprintf (state->ostream, "\n"); + } + fprintf (state->ostream, "\n"); + fflush (state->ostream); + +} + |