diff options
Diffstat (limited to 'gsl-1.9/histogram')
40 files changed, 6131 insertions, 0 deletions
diff --git a/gsl-1.9/histogram/ChangeLog b/gsl-1.9/histogram/ChangeLog new file mode 100644 index 0000000..18e5ab4 --- /dev/null +++ b/gsl-1.9/histogram/ChangeLog @@ -0,0 +1,156 @@ +2004-11-28 Brian Gough <bjg@network-theory.co.uk> + + * init2d.c (make_uniform): compute uniform range without + cancellation error. + + * init.c (make_uniform): compute uniform range without + cancellation error. + +Tue Aug 27 19:36:43 2002 Brian Gough <bjg@network-theory.co.uk> + + * test2d.c (main): changed test output format from %g to %e for + portability + + * test.c (main): changed test output format from %g to %e for + portability + +Wed Mar 6 22:03:35 2002 Brian Gough <bjg@network-theory.co.uk> + + * test2d.c (main): cleaned up the tests a bit + + * stat2d.c: added checks for wi>0 (Achim Gaedke) + +Sat Jan 26 17:09:10 2002 Brian Gough <bjg@network-theory.co.uk> + + * stat2d.c: added <math.h> include for sqrt + +Fri Jan 18 21:45:35 2002 Brian Gough <bjg@network-theory.co.uk> + + * stat2d.c: functions to compute statistics of 2d histograms + (Achim Gaedke) + +Mon Jan 14 19:34:31 2002 Brian Gough <bjg@network-theory.co.uk> + + * stat.c (gsl_histogram_sum): new function to sum bins (Achim + Gaedke) + + * maxval2d.c (gsl_histogram2d_sum): new function to sum bins + (Achim Gaedke) + +Thu Oct 18 14:48:07 2001 Brian Gough <bjg@network-theory.co.uk> + + * pdf2d.c (gsl_histogram2d_pdf_alloc): changed the definition of + the pdf alloc function to be consistent with the rest of the + library + + * pdf.c (gsl_histogram_pdf_alloc): changed the definition of the + pdf alloc function to be consistent with the rest of the library + + * init2d.c (gsl_histogram2d_alloc): added an alloc function for + consistency + + * init.c (gsl_histogram_alloc): added an alloc function for + consistency + +Wed Sep 12 13:38:40 2001 Brian Gough <bjg@network-theory.co.uk> + + * stat.c (gsl_histogram_mean): fixed calculation of mean/sigma and + made it part of the library + +Sun Aug 19 13:31:35 2001 Brian Gough <bjg@network-theory.co.uk> + + * test_gsl_histogram.sh: moved to top-level directory + +Sat Aug 18 22:21:26 2001 Brian Gough <bjg@network-theory.co.uk> + + * gsl-histogram.c: moved to top-level directory + +Mon Jun 25 17:41:42 2001 Brian Gough <bjg@network-theory.co.uk> + + * init2d.c (gsl_histogram2d_set_ranges_uniform): added range + initialization functions suggested by Achim Gaedke + + * init.c (gsl_histogram_set_ranges_uniform): added range + initialization functions suggested by Achim Gaedke + +Tue Apr 17 22:13:05 2001 Brian Gough <bjg@network-theory.co.uk> + + * get2d.c: include "find.c" + +Mon Apr 16 20:13:45 2001 Brian Gough <bjg@network-theory.co.uk> + + * get2d.c (gsl_histogram2d_get): removed unnecessary reference to + find2d + +Mon Jan 22 13:55:13 2001 Brian Gough <bjg@network-theory.co.uk> + + * find.c (find): optimize for the linear case, include own binary + search for speed + + * add.c (gsl_histogram_accumulate): fix check of array bound for + index + +Sun May 28 12:23:46 2000 Brian Gough <bjg@network-theory.co.uk> + + * test2d.c (main): use binary mode "b" when reading and writing + binary files + + * test.c (main): use binary mode "b" when reading and writing + binary files + +Wed Apr 26 15:09:22 2000 Brian Gough <bjg@network-theory.co.uk> + + * oper2d.c (gsl_histogram2d_shift): added function for shifting + histogram by a constant offset + + * oper.c (gsl_histogram_shift): added function for shifting + histogram by a constant offset + +Wed Apr 19 17:27:44 2000 Brian Gough <bjg@network-theory.co.uk> + + * added numerous extensions from Simone Piccardi + +2000-04-01 Mark Galassi <rosalia@lanl.gov> + + * *.c: changed 0 -> GSL_SUCCESS where appropriate; THANKS to Dave + Morrison. + +Fri Nov 19 15:31:51 1999 Brian Gough <bjg@network-theory.co.uk> + + * gsl-histogram.c (main): free memory before exit, eliminates + warning from checkergcc + + * test_gsl_histogram.sh: added a test for the gsl-histogram + program + +Fri Oct 1 15:47:01 1999 Brian Gough <bjg@network-theory.co.uk> + + * file.c file2d.c: converted to use new block i/o functions + +Wed Aug 18 11:41:40 1999 Brian Gough <bjg@network-theory.co.uk> + + * eliminated obvious memory leaks from the tests, so that we can + check that the _free functions work correctly + + * gsl-histogram.c (main): removed unused variable + +Fri Aug 6 11:19:37 1999 Brian Gough <bjg@network-theory.co.uk> + + * removed dependence on rand() and RAND_MAX + +1999-08-05 Mark Galassi <rosalia@lanl.gov> + + * gsl-histogram.c (main): fixed a simple logic bug. Thanks to + Barak Pearlmutter (bap@cs.unm.edu) for the patch. + +1998-11-06 <bjg@ancho.lanl.gov> + + * used a cast of (int) when attempting to print size_t variables + with %d + +Wed Sep 16 15:08:59 1998 Brian Gough <bjg@vvv.lanl.gov> + + * gsl-histogram.c (main): made the number of bins + optional. If you don't specify it then it uses bins of width 1. + + diff --git a/gsl-1.9/histogram/Makefile.am b/gsl-1.9/histogram/Makefile.am new file mode 100644 index 0000000..0adb316 --- /dev/null +++ b/gsl-1.9/histogram/Makefile.am @@ -0,0 +1,20 @@ +noinst_LTLIBRARIES = libgslhistogram.la + +pkginclude_HEADERS = gsl_histogram.h gsl_histogram2d.h + +INCLUDES= -I$(top_builddir) + +libgslhistogram_la_SOURCES = add.c get.c init.c params.c reset.c file.c pdf.c gsl_histogram.h add2d.c get2d.c init2d.c params2d.c reset2d.c file2d.c pdf2d.c gsl_histogram2d.h calloc_range.c calloc_range2d.c copy.c copy2d.c maxval.c maxval2d.c oper.c oper2d.c stat.c stat2d.c + +noinst_HEADERS = urand.c find.c find2d.c + +check_PROGRAMS = test +TESTS = $(check_PROGRAMS) + +EXTRA_DIST = urand.c + +test_SOURCES = test.c test1d.c test2d.c test1d_resample.c test2d_resample.c test1d_trap.c test2d_trap.c +test_LDADD = libgslhistogram.la ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la + +CLEANFILES = test.txt test.dat + diff --git a/gsl-1.9/histogram/Makefile.in b/gsl-1.9/histogram/Makefile.in new file mode 100644 index 0000000..3ebd2b6 --- /dev/null +++ b/gsl-1.9/histogram/Makefile.in @@ -0,0 +1,553 @@ +# 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 = histogram +DIST_COMMON = $(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) +libgslhistogram_la_LIBADD = +am_libgslhistogram_la_OBJECTS = add.lo get.lo init.lo params.lo \ + reset.lo file.lo pdf.lo add2d.lo get2d.lo init2d.lo \ + params2d.lo reset2d.lo file2d.lo pdf2d.lo calloc_range.lo \ + calloc_range2d.lo copy.lo copy2d.lo maxval.lo maxval2d.lo \ + oper.lo oper2d.lo stat.lo stat2d.lo +libgslhistogram_la_OBJECTS = $(am_libgslhistogram_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) test1d.$(OBJEXT) test2d.$(OBJEXT) \ + test1d_resample.$(OBJEXT) test2d_resample.$(OBJEXT) \ + test1d_trap.$(OBJEXT) test2d_trap.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslhistogram.la ../block/libgslblock.la \ + ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la \ + ../test/libgsltest.la ../sys/libgslsys.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 = $(libgslhistogram_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslhistogram_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 = libgslhistogram.la +pkginclude_HEADERS = gsl_histogram.h gsl_histogram2d.h +INCLUDES = -I$(top_builddir) +libgslhistogram_la_SOURCES = add.c get.c init.c params.c reset.c file.c pdf.c gsl_histogram.h add2d.c get2d.c init2d.c params2d.c reset2d.c file2d.c pdf2d.c gsl_histogram2d.h calloc_range.c calloc_range2d.c copy.c copy2d.c maxval.c maxval2d.c oper.c oper2d.c stat.c stat2d.c +noinst_HEADERS = urand.c find.c find2d.c +TESTS = $(check_PROGRAMS) +EXTRA_DIST = urand.c +test_SOURCES = test.c test1d.c test2d.c test1d_resample.c test2d_resample.c test1d_trap.c test2d_trap.c +test_LDADD = libgslhistogram.la ../block/libgslblock.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la +CLEANFILES = test.txt test.dat +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 histogram/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps histogram/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 +libgslhistogram.la: $(libgslhistogram_la_OBJECTS) $(libgslhistogram_la_DEPENDENCIES) + $(LINK) $(libgslhistogram_la_LDFLAGS) $(libgslhistogram_la_OBJECTS) $(libgslhistogram_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: + -test -z "$(CLEANFILES)" || rm -f $(CLEANFILES) + +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 + +# 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/histogram/TODO b/gsl-1.9/histogram/TODO new file mode 100644 index 0000000..1e5b07a --- /dev/null +++ b/gsl-1.9/histogram/TODO @@ -0,0 +1,3 @@ +* Implement N-d histograms (Simone Piccardi <piccardi@fi.infn.it> is +working on something here). + diff --git a/gsl-1.9/histogram/add.c b/gsl-1.9/histogram/add.c new file mode 100644 index 0000000..1697d41 --- /dev/null +++ b/gsl-1.9/histogram/add.c @@ -0,0 +1,55 @@ +/* histogram/add.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +#include "find.c" + +int +gsl_histogram_increment (gsl_histogram * h, double x) +{ + int status = gsl_histogram_accumulate (h, x, 1.0); + return status; +} + +int +gsl_histogram_accumulate (gsl_histogram * h, double x, double weight) +{ + const size_t n = h->n; + size_t index = 0; + + int status = find (h->n, h->range, x, &index); + + if (status) + { + return GSL_EDOM; + } + + if (index >= n) + { + GSL_ERROR ("index lies outside valid range of 0 .. n - 1", + GSL_ESANITY); + } + + h->bin[index] += weight; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/add2d.c b/gsl-1.9/histogram/add2d.c new file mode 100644 index 0000000..066119d --- /dev/null +++ b/gsl-1.9/histogram/add2d.c @@ -0,0 +1,66 @@ +/* histogram/add2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> +#include <gsl/gsl_histogram2d.h> + +#include "find2d.c" + +int +gsl_histogram2d_increment (gsl_histogram2d * h, double x, double y) +{ + int status = gsl_histogram2d_accumulate (h, x, y, 1.0); + return status; +} + +int +gsl_histogram2d_accumulate (gsl_histogram2d * h, + double x, double y, double weight) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + + size_t i = 0, j = 0; + + int status = find2d (h->nx, h->xrange, h->ny, h->yrange, x, y, &i, &j); + + if (status) + { + return GSL_EDOM; + } + + if (i >= nx) + { + GSL_ERROR ("index lies outside valid range of 0 .. nx - 1", + GSL_ESANITY); + } + + if (j >= ny) + { + GSL_ERROR ("index lies outside valid range of 0 .. ny - 1", + GSL_ESANITY); + } + + h->bin[i * ny + j] += weight; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/calloc_range.c b/gsl-1.9/histogram/calloc_range.c new file mode 100644 index 0000000..cb6e346 --- /dev/null +++ b/gsl-1.9/histogram/calloc_range.c @@ -0,0 +1,108 @@ +/* gsl_histogram_calloc_range.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram_calloc_range.c: + * Routine to create a variable binning histogram providing + * an input range vector. Need GSL library and header. + * Do range check and allocate the histogram data. + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +gsl_histogram * +gsl_histogram_calloc_range (size_t n, double *range) +{ + size_t i; + gsl_histogram *h; + + /* check arguments */ + + if (n == 0) + { + GSL_ERROR_VAL ("histogram length n must be positive integer", + GSL_EDOM, 0); + } + + /* check ranges */ + + for (i = 0; i < n; i++) + { + if (range[i] >= range[i + 1]) + { + GSL_ERROR_VAL ("histogram bin extremes must be " + "in increasing order", GSL_EDOM, 0); + } + } + + /* Allocate histogram */ + + h = (gsl_histogram *) malloc (sizeof (gsl_histogram)); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram struct", + GSL_ENOMEM, 0); + } + + h->range = (double *) malloc ((n + 1) * sizeof (double)); + + if (h->range == 0) + { + /* exception in constructor, avoid memory leak */ + free (h); + GSL_ERROR_VAL ("failed to allocate space for histogram ranges", + GSL_ENOMEM, 0); + } + + h->bin = (double *) malloc (n * sizeof (double)); + + if (h->bin == 0) + { + /* exception in constructor, avoid memory leak */ + free (h->range); + free (h); + GSL_ERROR_VAL ("failed to allocate space for histogram bins", + GSL_ENOMEM, 0); + } + + /* initialize ranges */ + + for (i = 0; i <= n; i++) + { + h->range[i] = range[i]; + } + + /* clear contents */ + + for (i = 0; i < n; i++) + { + h->bin[i] = 0; + } + + h->n = n; + + return h; +} diff --git a/gsl-1.9/histogram/calloc_range2d.c b/gsl-1.9/histogram/calloc_range2d.c new file mode 100644 index 0000000..98d0caf --- /dev/null +++ b/gsl-1.9/histogram/calloc_range2d.c @@ -0,0 +1,153 @@ +/* gsl_histogram2d_calloc_range.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram2d_calloc_range.c: + * Routine to create a variable binning 2D histogram providing + * the input range vectors. Need GSL library and header. + * Do range check and allocate the histogram data. + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> +/* + * Routine that create a 2D histogram using the given + * values for X and Y ranges + */ +gsl_histogram2d * +gsl_histogram2d_calloc_range (size_t nx, size_t ny, + double *xrange, + double *yrange) +{ + size_t i, j; + gsl_histogram2d *h; + + /* check arguments */ + + if (nx == 0) + { + GSL_ERROR_VAL ("histogram length nx must be positive integer", + GSL_EDOM, 0); + } + + if (ny == 0) + { + GSL_ERROR_VAL ("histogram length ny must be positive integer", + GSL_EDOM, 0); + } + + /* init ranges */ + + for (i = 0; i < nx; i++) + { + if (xrange[i] >= xrange[i + 1]) + { + GSL_ERROR_VAL ("histogram xrange not in increasing order", + GSL_EDOM, 0); + } + } + + for (j = 0; j < ny; j++) + { + if (yrange[j] >= yrange[j + 1]) + { + GSL_ERROR_VAL ("histogram yrange not in increasing order" + ,GSL_EDOM, 0); + } + } + + /* Allocate histogram */ + + h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d)); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram struct", + GSL_ENOMEM, 0); + } + + h->xrange = (double *) malloc ((nx + 1) * sizeof (double)); + + if (h->xrange == 0) + { + /* exception in constructor, avoid memory leak */ + free (h); + + GSL_ERROR_VAL ("failed to allocate space for histogram xrange", + GSL_ENOMEM, 0); + } + + h->yrange = (double *) malloc ((ny + 1) * sizeof (double)); + + if (h->yrange == 0) + { + /* exception in constructor, avoid memory leak */ + free (h); + + GSL_ERROR_VAL ("failed to allocate space for histogram yrange", + GSL_ENOMEM, 0); + } + + h->bin = (double *) malloc (nx * ny * sizeof (double)); + + if (h->bin == 0) + { + /* exception in constructor, avoid memory leak */ + free (h->xrange); + free (h->yrange); + free (h); + + GSL_ERROR_VAL ("failed to allocate space for histogram bins", + GSL_ENOMEM, 0); + } + + /* init histogram */ + + /* init ranges */ + + for (i = 0; i <= nx; i++) + { + h->xrange[i] = xrange[i]; + } + + for (j = 0; j <= ny; j++) + { + h->yrange[j] = yrange[j]; + } + + /* clear contents */ + + for (i = 0; i < nx; i++) + { + for (j = 0; j < ny; j++) + { + h->bin[i * ny + j] = 0; + } + } + + h->nx = nx; + h->ny = ny; + + return h; +} diff --git a/gsl-1.9/histogram/copy.c b/gsl-1.9/histogram/copy.c new file mode 100644 index 0000000..848dcf1 --- /dev/null +++ b/gsl-1.9/histogram/copy.c @@ -0,0 +1,91 @@ +/* gsl_histogram_copy.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram_copy.c: + * Routine to copy an histogram. + * Need GSL library and headers. + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +/* + * gsl_histogram_copy: + * copy the contents of an histogram into another + */ + +int +gsl_histogram_memcpy (gsl_histogram * dest, const gsl_histogram * src) +{ + size_t n = src->n; + size_t i; + + if (dest->n != src->n) + { + GSL_ERROR ("histograms have different sizes, cannot copy", + GSL_EINVAL); + } + + for (i = 0; i <= n; i++) + { + dest->range[i] = src->range[i]; + } + + for (i = 0; i < n; i++) + { + dest->bin[i] = src->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram_duplicate: + * duplicate an histogram creating + * an identical new one + */ + +gsl_histogram * +gsl_histogram_clone (const gsl_histogram * src) +{ + size_t n = src->n; + size_t i; + gsl_histogram *h; + + h = gsl_histogram_calloc_range (n, src->range); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram struct", + GSL_ENOMEM, 0); + } + + for (i = 0; i < n; i++) + { + h->bin[i] = src->bin[i]; + } + + return h; +} diff --git a/gsl-1.9/histogram/copy2d.c b/gsl-1.9/histogram/copy2d.c new file mode 100644 index 0000000..94fa2d6 --- /dev/null +++ b/gsl-1.9/histogram/copy2d.c @@ -0,0 +1,96 @@ +/* gsl_histogram2d_copy.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram2d_copy.c: + * Routine to copy a 2D histogram. + * Need GSL library and header. + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +/* + * gsl_histogram2d_copy: + * copy the contents of an histogram into another + */ +int +gsl_histogram2d_memcpy (gsl_histogram2d * dest, const gsl_histogram2d * src) +{ + size_t nx = src->nx; + size_t ny = src->ny; + size_t i; + if (dest->nx != src->nx || dest->ny != src->ny) + { + GSL_ERROR ("histograms have different sizes, cannot copy", + GSL_EINVAL); + } + + for (i = 0; i <= nx; i++) + { + dest->xrange[i] = src->xrange[i]; + } + + for (i = 0; i <= ny; i++) + { + dest->yrange[i] = src->yrange[i]; + } + + for (i = 0; i < nx * ny; i++) + { + dest->bin[i] = src->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_duplicate: + * duplicate an histogram creating + * an identical new one + */ + +gsl_histogram2d * +gsl_histogram2d_clone (const gsl_histogram2d * src) +{ + size_t nx = src->nx; + size_t ny = src->ny; + size_t i; + gsl_histogram2d *h; + + h = gsl_histogram2d_calloc_range (nx, ny, src->xrange, src->yrange); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram struct", + GSL_ENOMEM, 0); + } + + for (i = 0; i < nx * ny; i++) + { + h->bin[i] = src->bin[i]; + } + + return h; +} diff --git a/gsl-1.9/histogram/file.c b/gsl-1.9/histogram/file.c new file mode 100644 index 0000000..7013dc0 --- /dev/null +++ b/gsl-1.9/histogram/file.c @@ -0,0 +1,127 @@ +/* histogram/file.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_block.h> +#include <gsl/gsl_histogram.h> + +int +gsl_histogram_fread (FILE * stream, gsl_histogram * h) +{ + int status = gsl_block_raw_fread (stream, h->range, h->n + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fread (stream, h->bin, h->n, 1); + return status; +} + +int +gsl_histogram_fwrite (FILE * stream, const gsl_histogram * h) +{ + int status = gsl_block_raw_fwrite (stream, h->range, h->n + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fwrite (stream, h->bin, h->n, 1); + return status; +} + +int +gsl_histogram_fprintf (FILE * stream, const gsl_histogram * h, + const char *range_format, const char *bin_format) +{ + size_t i; + const size_t n = h->n; + + for (i = 0; i < n; i++) + { + int status = fprintf (stream, range_format, h->range[i]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, range_format, h->range[i + 1]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, bin_format, h->bin[i]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc ('\n', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + } + + return GSL_SUCCESS; +} + +int +gsl_histogram_fscanf (FILE * stream, gsl_histogram * h) +{ + size_t i; + const size_t n = h->n; + double upper; + + for (i = 0; i < n; i++) + { + int status = fscanf (stream, + "%lg %lg %lg", h->range + i, &upper, + h->bin + i); + + if (status != 3) + { + GSL_ERROR ("fscanf failed", GSL_EFAILED); + } + } + + h->range[n] = upper; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/file2d.c b/gsl-1.9/histogram/file2d.c new file mode 100644 index 0000000..2bf1b6e --- /dev/null +++ b/gsl-1.9/histogram/file2d.c @@ -0,0 +1,184 @@ +/* histogram/file2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_block.h> +#include <gsl/gsl_histogram2d.h> + +int +gsl_histogram2d_fread (FILE * stream, gsl_histogram2d * h) +{ + int status = gsl_block_raw_fread (stream, h->xrange, h->nx + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fread (stream, h->yrange, h->ny + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fread (stream, h->bin, h->nx * h->ny, 1); + + return status; +} + +int +gsl_histogram2d_fwrite (FILE * stream, const gsl_histogram2d * h) +{ + int status = gsl_block_raw_fwrite (stream, h->xrange, h->nx + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fwrite (stream, h->yrange, h->ny + 1, 1); + + if (status) + return status; + + status = gsl_block_raw_fwrite (stream, h->bin, h->nx * h->ny, 1); + return status; +} + +int +gsl_histogram2d_fprintf (FILE * stream, const gsl_histogram2d * h, + const char *range_format, const char *bin_format) +{ + size_t i, j; + const size_t nx = h->nx; + const size_t ny = h->ny; + int status; + + for (i = 0; i < nx; i++) + { + for (j = 0; j < ny; j++) + { + status = fprintf (stream, range_format, h->xrange[i]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, range_format, h->xrange[i + 1]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, range_format, h->yrange[j]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, range_format, h->yrange[j + 1]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc (' ', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + + status = fprintf (stream, bin_format, h->bin[i * ny + j]); + + if (status < 0) + { + GSL_ERROR ("fprintf failed", GSL_EFAILED); + } + + status = putc ('\n', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + } + status = putc ('\n', stream); + + if (status == EOF) + { + GSL_ERROR ("putc failed", GSL_EFAILED); + } + } + + return GSL_SUCCESS; +} + +int +gsl_histogram2d_fscanf (FILE * stream, gsl_histogram2d * h) +{ + size_t i, j; + const size_t nx = h->nx; + const size_t ny = h->ny; + double xupper, yupper; + + for (i = 0; i < nx; i++) + { + for (j = 0; j < ny; j++) + { + int status = fscanf (stream, + "%lg %lg %lg %lg %lg", + h->xrange + i, &xupper, + h->yrange + j, &yupper, + h->bin + i * ny + j); + + if (status != 5) + { + GSL_ERROR ("fscanf failed", GSL_EFAILED); + } + } + h->yrange[ny] = yupper; + } + + h->xrange[nx] = xupper; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/find.c b/gsl-1.9/histogram/find.c new file mode 100644 index 0000000..b252a36 --- /dev/null +++ b/gsl-1.9/histogram/find.c @@ -0,0 +1,86 @@ +/* histogram/find.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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. + */ + +/* determines whether to optimize for linear ranges */ +#define LINEAR_OPT 1 + +static int find (const size_t n, const double range[], + const double x, size_t * i); + +static int +find (const size_t n, const double range[], const double x, size_t * i) +{ + size_t i_linear, lower, upper, mid; + + if (x < range[0]) + { + return -1; + } + + if (x >= range[n]) + { + return +1; + } + + /* optimize for linear case */ + +#ifdef LINEAR_OPT + { + double u = (x - range[0]) / (range[n] - range[0]); + i_linear = (size_t) (u * n); + } + + if (x >= range[i_linear] && x < range[i_linear + 1]) + { + *i = i_linear; + return 0; + } +#endif + + /* perform binary search */ + + upper = n ; + lower = 0 ; + + while (upper - lower > 1) + { + mid = (upper + lower) / 2 ; + + if (x >= range[mid]) + { + lower = mid ; + } + else + { + upper = mid ; + } + } + + *i = lower ; + + /* sanity check the result */ + + if (x < range[lower] || x >= range[lower + 1]) + { + GSL_ERROR ("x not found in range", GSL_ESANITY); + } + + return 0; +} + diff --git a/gsl-1.9/histogram/find2d.c b/gsl-1.9/histogram/find2d.c new file mode 100644 index 0000000..1b8e5b5 --- /dev/null +++ b/gsl-1.9/histogram/find2d.c @@ -0,0 +1,50 @@ +/* histogram/find2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 "find.c" + +static int +find2d (const size_t nx, const double xrange[], + const size_t ny, const double yrange[], + const double x, const double y, + size_t * i, size_t * j); + + +static int +find2d (const size_t nx, const double xrange[], + const size_t ny, const double yrange[], + const double x, const double y, + size_t * i, size_t * j) +{ + int status = find (nx, xrange, x, i); + + if (status) + { + return status; + } + + status = find (ny, yrange, y, j); + + if (status) + { + return status; + } + + return 0; +} diff --git a/gsl-1.9/histogram/get.c b/gsl-1.9/histogram/get.c new file mode 100644 index 0000000..6e832a5 --- /dev/null +++ b/gsl-1.9/histogram/get.c @@ -0,0 +1,69 @@ +/* histogram/get.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +#include "find.c" + +double +gsl_histogram_get (const gsl_histogram * h, size_t i) +{ + const size_t n = h->n; + + if (i >= n) + { + GSL_ERROR_VAL ("index lies outside valid range of 0 .. n - 1", + GSL_EDOM, 0); + } + + return h->bin[i]; +} + +int +gsl_histogram_get_range (const gsl_histogram * h, size_t i, + double *lower, double *upper) +{ + const size_t n = h->n; + + if (i >= n) + { + GSL_ERROR ("index lies outside valid range of 0 .. n - 1", GSL_EDOM); + } + + *lower = h->range[i]; + *upper = h->range[i + 1]; + + return GSL_SUCCESS; +} + +int +gsl_histogram_find (const gsl_histogram * h, + const double x, size_t * i) +{ + int status = find (h->n, h->range, x, i); + + if (status) + { + GSL_ERROR ("x not found in range of h", GSL_EDOM); + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/get2d.c b/gsl-1.9/histogram/get2d.c new file mode 100644 index 0000000..593e2de --- /dev/null +++ b/gsl-1.9/histogram/get2d.c @@ -0,0 +1,101 @@ +/* histogram/get2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +#include "find.c" + +double +gsl_histogram2d_get (const gsl_histogram2d * h, const size_t i, const size_t j) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + + if (i >= nx) + { + GSL_ERROR_VAL ("index i lies outside valid range of 0 .. nx - 1", + GSL_EDOM, 0); + } + + if (j >= ny) + { + GSL_ERROR_VAL ("index j lies outside valid range of 0 .. ny - 1", + GSL_EDOM, 0); + } + + return h->bin[i * ny + j]; +} + +int +gsl_histogram2d_get_xrange (const gsl_histogram2d * h, const size_t i, + double *xlower, double *xupper) +{ + const size_t nx = h->nx; + + if (i >= nx) + { + GSL_ERROR ("index i lies outside valid range of 0 .. nx - 1", GSL_EDOM); + } + + *xlower = h->xrange[i]; + *xupper = h->xrange[i + 1]; + + return GSL_SUCCESS; +} + +int +gsl_histogram2d_get_yrange (const gsl_histogram2d * h, const size_t j, + double *ylower, double *yupper) +{ + const size_t ny = h->ny; + + if (j >= ny) + { + GSL_ERROR ("index j lies outside valid range of 0 .. ny - 1", GSL_EDOM); + } + + *ylower = h->yrange[j]; + *yupper = h->yrange[j + 1]; + + return GSL_SUCCESS; +} + +int +gsl_histogram2d_find (const gsl_histogram2d * h, + const double x, const double y, + size_t * i, size_t * j) +{ + int status = find (h->nx, h->xrange, x, i); + + if (status) + { + GSL_ERROR ("x not found in range of h", GSL_EDOM); + } + + status = find (h->ny, h->yrange, y, j); + + if (status) + { + GSL_ERROR ("y not found in range of h", GSL_EDOM); + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/gsl_histogram.h b/gsl-1.9/histogram/gsl_histogram.h new file mode 100644 index 0000000..e246e3a --- /dev/null +++ b/gsl-1.9/histogram/gsl_histogram.h @@ -0,0 +1,134 @@ +/* histogram/gsl_histogram.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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. + */ + +#ifndef __GSL_HISTOGRAM_H__ +#define __GSL_HISTOGRAM_H__ + +#include <stdlib.h> +#include <stdio.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 n ; + double * range ; + double * bin ; +} gsl_histogram ; + +typedef struct { + size_t n ; + double * range ; + double * sum ; +} gsl_histogram_pdf ; + +gsl_histogram * gsl_histogram_alloc (size_t n); + +gsl_histogram * gsl_histogram_calloc (size_t n); +gsl_histogram * gsl_histogram_calloc_uniform (const size_t n, const double xmin, const double xmax); +void gsl_histogram_free (gsl_histogram * h); +int gsl_histogram_increment (gsl_histogram * h, double x); +int gsl_histogram_accumulate (gsl_histogram * h, double x, double weight); +int gsl_histogram_find (const gsl_histogram * h, + const double x, size_t * i); + +double gsl_histogram_get (const gsl_histogram * h, size_t i); +int gsl_histogram_get_range (const gsl_histogram * h, size_t i, + double * lower, double * upper); + +double gsl_histogram_max (const gsl_histogram * h); +double gsl_histogram_min (const gsl_histogram * h); +size_t gsl_histogram_bins (const gsl_histogram * h); + +void gsl_histogram_reset (gsl_histogram * h); + +gsl_histogram * gsl_histogram_calloc_range(size_t n, double * range); + +int +gsl_histogram_set_ranges (gsl_histogram * h, const double range[], size_t size); +int +gsl_histogram_set_ranges_uniform (gsl_histogram * h, double xmin, double xmax); + + + +int +gsl_histogram_memcpy(gsl_histogram * dest, const gsl_histogram * source); + +gsl_histogram * +gsl_histogram_clone(const gsl_histogram * source); + +double gsl_histogram_max_val (const gsl_histogram * h); + +size_t gsl_histogram_max_bin (const gsl_histogram * h); + +double gsl_histogram_min_val (const gsl_histogram * h); + +size_t gsl_histogram_min_bin (const gsl_histogram * h); + +int +gsl_histogram_equal_bins_p(const gsl_histogram *h1, const gsl_histogram *h2); + +int +gsl_histogram_add(gsl_histogram *h1, const gsl_histogram *h2); + +int +gsl_histogram_sub(gsl_histogram *h1, const gsl_histogram *h2); + +int +gsl_histogram_mul(gsl_histogram *h1, const gsl_histogram *h2); + +int +gsl_histogram_div(gsl_histogram *h1, const gsl_histogram *h2); + +int +gsl_histogram_scale(gsl_histogram *h, double scale); + +int +gsl_histogram_shift (gsl_histogram * h, double shift); + + +double gsl_histogram_sigma (const gsl_histogram * h); + +double gsl_histogram_mean (const gsl_histogram * h); + +double gsl_histogram_sum (const gsl_histogram * h); + +int gsl_histogram_fwrite (FILE * stream, const gsl_histogram * h) ; +int gsl_histogram_fread (FILE * stream, gsl_histogram * h); +int gsl_histogram_fprintf (FILE * stream, const gsl_histogram * h, + const char * range_format, const char * bin_format); +int gsl_histogram_fscanf (FILE * stream, gsl_histogram * h); + +gsl_histogram_pdf * gsl_histogram_pdf_alloc (const size_t n); +int gsl_histogram_pdf_init (gsl_histogram_pdf * p, const gsl_histogram * h); +void gsl_histogram_pdf_free (gsl_histogram_pdf * p); +double gsl_histogram_pdf_sample (const gsl_histogram_pdf * p, double r); + +__END_DECLS + +#endif /* __GSL_HISTOGRAM_H__ */ diff --git a/gsl-1.9/histogram/gsl_histogram2d.h b/gsl-1.9/histogram/gsl_histogram2d.h new file mode 100644 index 0000000..b4628aa --- /dev/null +++ b/gsl-1.9/histogram/gsl_histogram2d.h @@ -0,0 +1,172 @@ +/* histogram/gsl_histogram2d.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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. + */ + +#ifndef __GSL_HISTOGRAM2D_H__ +#define __GSL_HISTOGRAM2D_H__ + +#include <stdlib.h> +#include <stdio.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 nx, ny ; + double * xrange ; + double * yrange ; + double * bin ; +} gsl_histogram2d ; + +typedef struct { + size_t nx, ny ; + double * xrange ; + double * yrange ; + double * sum ; +} gsl_histogram2d_pdf ; + +gsl_histogram2d * gsl_histogram2d_alloc (const size_t nx, const size_t ny); +gsl_histogram2d * gsl_histogram2d_calloc (const size_t nx, const size_t ny); +gsl_histogram2d * gsl_histogram2d_calloc_uniform (const size_t nx, const size_t ny, + const double xmin, const double xmax, + const double ymin, const double ymax); + +void gsl_histogram2d_free (gsl_histogram2d * h); + +int gsl_histogram2d_increment (gsl_histogram2d * h, double x, double y); +int gsl_histogram2d_accumulate (gsl_histogram2d * h, + double x, double y, double weight); +int gsl_histogram2d_find (const gsl_histogram2d * h, + const double x, const double y, size_t * i, size_t * j); + +double gsl_histogram2d_get (const gsl_histogram2d * h, const size_t i, const size_t j); +int gsl_histogram2d_get_xrange (const gsl_histogram2d * h, const size_t i, + double * xlower, double * xupper); +int gsl_histogram2d_get_yrange (const gsl_histogram2d * h, const size_t j, + double * ylower, double * yupper); + + +double gsl_histogram2d_xmax (const gsl_histogram2d * h); +double gsl_histogram2d_xmin (const gsl_histogram2d * h); +size_t gsl_histogram2d_nx (const gsl_histogram2d * h); + +double gsl_histogram2d_ymax (const gsl_histogram2d * h); +double gsl_histogram2d_ymin (const gsl_histogram2d * h); +size_t gsl_histogram2d_ny (const gsl_histogram2d * h); + +void gsl_histogram2d_reset (gsl_histogram2d * h); + +gsl_histogram2d * +gsl_histogram2d_calloc_range(size_t nx, size_t ny, + double *xrange, double *yrange); + +int +gsl_histogram2d_set_ranges_uniform (gsl_histogram2d * h, + double xmin, double xmax, + double ymin, double ymax); + +int +gsl_histogram2d_set_ranges (gsl_histogram2d * h, + const double xrange[], size_t xsize, + const double yrange[], size_t ysize); + +int +gsl_histogram2d_memcpy(gsl_histogram2d *dest, const gsl_histogram2d *source); + +gsl_histogram2d * +gsl_histogram2d_clone(const gsl_histogram2d * source); + +double +gsl_histogram2d_max_val(const gsl_histogram2d *h); + +void +gsl_histogram2d_max_bin (const gsl_histogram2d *h, size_t *i, size_t *j); + +double +gsl_histogram2d_min_val(const gsl_histogram2d *h); + +void +gsl_histogram2d_min_bin (const gsl_histogram2d *h, size_t *i, size_t *j); + +double +gsl_histogram2d_xmean (const gsl_histogram2d * h); + +double +gsl_histogram2d_ymean (const gsl_histogram2d * h); + +double +gsl_histogram2d_xsigma (const gsl_histogram2d * h); + +double +gsl_histogram2d_ysigma (const gsl_histogram2d * h); + +double +gsl_histogram2d_cov (const gsl_histogram2d * h); + +double +gsl_histogram2d_sum (const gsl_histogram2d *h); + +int +gsl_histogram2d_equal_bins_p(const gsl_histogram2d *h1, + const gsl_histogram2d *h2) ; + +int +gsl_histogram2d_add(gsl_histogram2d *h1, const gsl_histogram2d *h2); + +int +gsl_histogram2d_sub(gsl_histogram2d *h1, const gsl_histogram2d *h2); + +int +gsl_histogram2d_mul(gsl_histogram2d *h1, const gsl_histogram2d *h2); + +int +gsl_histogram2d_div(gsl_histogram2d *h1, const gsl_histogram2d *h2); + +int +gsl_histogram2d_scale(gsl_histogram2d *h, double scale); + +int +gsl_histogram2d_shift(gsl_histogram2d *h, double shift); + +int gsl_histogram2d_fwrite (FILE * stream, const gsl_histogram2d * h) ; +int gsl_histogram2d_fread (FILE * stream, gsl_histogram2d * h); +int gsl_histogram2d_fprintf (FILE * stream, const gsl_histogram2d * h, + const char * range_format, + const char * bin_format); +int gsl_histogram2d_fscanf (FILE * stream, gsl_histogram2d * h); + +gsl_histogram2d_pdf * gsl_histogram2d_pdf_alloc (const size_t nx, const size_t ny); +int gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * p, const gsl_histogram2d * h); +void gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * p); +int gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * p, + double r1, double r2, + double * x, double * y); + +__END_DECLS + +#endif /* __GSL_HISTOGRAM2D_H__ */ + diff --git a/gsl-1.9/histogram/init.c b/gsl-1.9/histogram/init.c new file mode 100644 index 0000000..6ff47f5 --- /dev/null +++ b/gsl-1.9/histogram/init.c @@ -0,0 +1,198 @@ +/* histogram/init.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +gsl_histogram * +gsl_histogram_alloc (size_t n) +{ + gsl_histogram *h; + + if (n == 0) + { + GSL_ERROR_VAL ("histogram length n must be positive integer", + GSL_EDOM, 0); + } + + h = (gsl_histogram *) malloc (sizeof (gsl_histogram)); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram struct", + GSL_ENOMEM, 0); + } + + h->range = (double *) malloc ((n + 1) * sizeof (double)); + + if (h->range == 0) + { + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram ranges", + GSL_ENOMEM, 0); + } + + h->bin = (double *) malloc (n * sizeof (double)); + + if (h->bin == 0) + { + free (h->range); + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram bins", + GSL_ENOMEM, 0); + } + + h->n = n; + + return h; +} + +static void +make_uniform (double range[], size_t n, double xmin, double xmax) +{ + size_t i; + + for (i = 0; i <= n; i++) + { + double f1 = ((double) (n-i) / (double) n); + double f2 = ((double) i / (double) n); + range[i] = f1 * xmin + f2 * xmax; + } +} + +gsl_histogram * +gsl_histogram_calloc_uniform (const size_t n, const double xmin, + const double xmax) +{ + gsl_histogram *h; + + if (xmin >= xmax) + { + GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0); + } + + h = gsl_histogram_calloc (n); + + if (h == 0) + { + return h; + } + + make_uniform (h->range, n, xmin, xmax); + + return h; +} + +gsl_histogram * +gsl_histogram_calloc (size_t n) +{ + gsl_histogram * h = gsl_histogram_alloc (n); + + if (h == 0) + { + return h; + } + + { + size_t i; + + for (i = 0; i < n + 1; i++) + { + h->range[i] = i; + } + + for (i = 0; i < n; i++) + { + h->bin[i] = 0; + } + } + + h->n = n; + + return h; +} + + +void +gsl_histogram_free (gsl_histogram * h) +{ + free (h->range); + free (h->bin); + free (h); +} + +/* These initialization functions suggested by Achim Gaedke */ + +int +gsl_histogram_set_ranges_uniform (gsl_histogram * h, double xmin, double xmax) +{ + size_t i; + const size_t n = h->n; + + if (xmin >= xmax) + { + GSL_ERROR ("xmin must be less than xmax", GSL_EINVAL); + } + + /* initialize ranges */ + + make_uniform (h->range, n, xmin, xmax); + + /* clear contents */ + + for (i = 0; i < n; i++) + { + h->bin[i] = 0; + } + + return GSL_SUCCESS; +} + +int +gsl_histogram_set_ranges (gsl_histogram * h, const double range[], size_t size) +{ + size_t i; + const size_t n = h->n; + + if (size != (n+1)) + { + GSL_ERROR ("size of range must match size of histogram", GSL_EINVAL); + } + + /* initialize ranges */ + + for (i = 0; i <= n; i++) + { + h->range[i] = range[i]; + } + + /* clear contents */ + + for (i = 0; i < n; i++) + { + h->bin[i] = 0; + } + + return GSL_SUCCESS; +} + diff --git a/gsl-1.9/histogram/init2d.c b/gsl-1.9/histogram/init2d.c new file mode 100644 index 0000000..bed068f --- /dev/null +++ b/gsl-1.9/histogram/init2d.c @@ -0,0 +1,300 @@ +/* histogram/init2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +gsl_histogram2d * +gsl_histogram2d_alloc (const size_t nx, const size_t ny) +{ + gsl_histogram2d *h; + + if (nx == 0) + { + GSL_ERROR_VAL ("histogram2d length nx must be positive integer", + GSL_EDOM, 0); + } + + if (ny == 0) + { + GSL_ERROR_VAL ("histogram2d length ny must be positive integer", + GSL_EDOM, 0); + } + + h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d)); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram2d struct", + GSL_ENOMEM, 0); + } + + h->xrange = (double *) malloc ((nx + 1) * sizeof (double)); + + if (h->xrange == 0) + { + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d x ranges", + GSL_ENOMEM, 0); + } + + h->yrange = (double *) malloc ((ny + 1) * sizeof (double)); + + if (h->yrange == 0) + { + free (h->xrange); + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d y ranges", + GSL_ENOMEM, 0); + } + + h->bin = (double *) malloc (nx * ny * sizeof (double)); + + if (h->bin == 0) + { + free (h->xrange); + free (h->yrange); + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram bins", + GSL_ENOMEM, 0); + } + + h->nx = nx; + h->ny = ny; + + return h; +} + +static void +make_uniform (double range[], size_t n, double xmin, double xmax) +{ + size_t i; + + for (i = 0; i <= n; i++) + { + double f1 = ((double) (n-i) / (double) n); + double f2 = ((double) i / (double) n); + range[i] = f1 * xmin + f2 * xmax; + } +} + +gsl_histogram2d * +gsl_histogram2d_calloc_uniform (const size_t nx, const size_t ny, + const double xmin, const double xmax, + const double ymin, const double ymax) +{ + gsl_histogram2d *h; + + if (xmin >= xmax) + { + GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0); + } + + if (ymin >= ymax) + { + GSL_ERROR_VAL ("ymin must be less than ymax", GSL_EINVAL, 0); + } + + h = gsl_histogram2d_calloc (nx, ny); + + if (h == 0) + { + return h; + } + + make_uniform (h->xrange, nx, xmin, xmax); + make_uniform (h->yrange, ny, ymin, ymax); + + return h; +} + +gsl_histogram2d * +gsl_histogram2d_calloc (const size_t nx, const size_t ny) +{ + gsl_histogram2d *h; + + if (nx == 0) + { + GSL_ERROR_VAL ("histogram2d length nx must be positive integer", + GSL_EDOM, 0); + } + + if (ny == 0) + { + GSL_ERROR_VAL ("histogram2d length ny must be positive integer", + GSL_EDOM, 0); + } + + h = (gsl_histogram2d *) malloc (sizeof (gsl_histogram2d)); + + if (h == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram2d struct", + GSL_ENOMEM, 0); + } + + h->xrange = (double *) malloc ((nx + 1) * sizeof (double)); + + if (h->xrange == 0) + { + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d x ranges", + GSL_ENOMEM, 0); + } + + h->yrange = (double *) malloc ((ny + 1) * sizeof (double)); + + if (h->yrange == 0) + { + free (h->xrange); + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d y ranges", + GSL_ENOMEM, 0); + } + + h->bin = (double *) malloc (nx * ny * sizeof (double)); + + if (h->bin == 0) + { + free (h->xrange); + free (h->yrange); + free (h); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram bins", + GSL_ENOMEM, 0); + } + + { + size_t i; + + for (i = 0; i < nx + 1; i++) + { + h->xrange[i] = i; + } + + for (i = 0; i < ny + 1; i++) + { + h->yrange[i] = i; + } + + for (i = 0; i < nx * ny; i++) + { + h->bin[i] = 0; + } + } + + h->nx = nx; + h->ny = ny; + + return h; +} + + +void +gsl_histogram2d_free (gsl_histogram2d * h) +{ + free (h->xrange); + free (h->yrange); + free (h->bin); + free (h); +} + + +int +gsl_histogram2d_set_ranges_uniform (gsl_histogram2d * h, + double xmin, double xmax, + double ymin, double ymax) +{ + size_t i; + const size_t nx = h->nx, ny = h->ny; + + if (xmin >= xmax) + { + GSL_ERROR_VAL ("xmin must be less than xmax", GSL_EINVAL, 0); + } + + if (ymin >= ymax) + { + GSL_ERROR_VAL ("ymin must be less than ymax", GSL_EINVAL, 0); + } + + /* initialize ranges */ + + make_uniform (h->xrange, nx, xmin, xmax); + make_uniform (h->yrange, ny, ymin, ymax); + + /* clear contents */ + + for (i = 0; i < nx * ny; i++) + { + h->bin[i] = 0; + } + + return GSL_SUCCESS; +} + +int +gsl_histogram2d_set_ranges (gsl_histogram2d * h, + const double xrange[], size_t xsize, + const double yrange[], size_t ysize) +{ + size_t i; + const size_t nx = h->nx, ny = h->ny; + + if (xsize != (nx + 1)) + { + GSL_ERROR_VAL ("size of xrange must match size of histogram", + GSL_EINVAL, 0); + } + + if (ysize != (ny + 1)) + { + GSL_ERROR_VAL ("size of yrange must match size of histogram", + GSL_EINVAL, 0); + } + + /* initialize ranges */ + + for (i = 0; i <= nx; i++) + { + h->xrange[i] = xrange[i]; + } + + for (i = 0; i <= ny; i++) + { + h->yrange[i] = yrange[i]; + } + + /* clear contents */ + + for (i = 0; i < nx * ny; i++) + { + h->bin[i] = 0; + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/maxval.c b/gsl-1.9/histogram/maxval.c new file mode 100644 index 0000000..2031537 --- /dev/null +++ b/gsl-1.9/histogram/maxval.c @@ -0,0 +1,101 @@ +/* gsl_histogram_maxval.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram_maxval.c: + * Routine to find maximum and minumum content of a hisogram. + * Need GSL library and header. + * Contains the routines: + * gsl_histogram_max_val find max content values + * gsl_histogram_min_val find min content values + * gsl_histogram_bin_max find coordinates of max contents bin + * gsl_histogram_bin_min find coordinates of min contents bin + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +double +gsl_histogram_max_val (const gsl_histogram * h) +{ + const size_t n = h->n; + size_t i; + double max = h->bin[0]; + for (i = 0; i < n; i++) + { + if (h->bin[i] > max) + { + max = h->bin[i]; + } + } + return max; +} + +size_t +gsl_histogram_max_bin (const gsl_histogram * h) +{ + size_t i; + size_t imax = 0; + double max = h->bin[0]; + for (i = 0; i < h->n; i++) + { + if (h->bin[i] > max) + { + max = h->bin[i]; + imax = i; + } + } + return imax; +} + +double +gsl_histogram_min_val (const gsl_histogram * h) +{ + size_t i; + double min = h->bin[0]; + for (i = 0; i < h->n; i++) + { + if (h->bin[i] < min) + { + min = h->bin[i]; + } + } + return min; +} + +size_t +gsl_histogram_min_bin (const gsl_histogram * h) +{ + size_t i; + size_t imin = 0; + double min = h->bin[0]; + for (i = 0; i < h->n; i++) + { + if (h->bin[i] < min) + { + min = h->bin[i]; + imin = i; + } + } + return imin; +} diff --git a/gsl-1.9/histogram/maxval2d.c b/gsl-1.9/histogram/maxval2d.c new file mode 100644 index 0000000..916fc85 --- /dev/null +++ b/gsl-1.9/histogram/maxval2d.c @@ -0,0 +1,140 @@ +/* gsl_histogram2d_maxval.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram2d_maxval.c: + * Routine to find maximum and minumum content of a 2D hisogram. + * Need GSL library and header. + * Contains the routines: + * gsl_histogram2d_max_val find max content values + * gsl_histogram2d_min_val find min content values + * gsl_histogram2d_bin_max find coordinates of max contents bin + * gsl_histogram2d_bin_min find coordinates of min contents bin + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> +/* + * Return the maximum contents value of a 2D histogram + */ +double +gsl_histogram2d_max_val (const gsl_histogram2d * h) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + double max = h->bin[0 * ny + 0]; + + for (i = 0; i < nx * ny; i++) + { + if (h->bin[i] > max) + { + max = h->bin[i]; + } + } + return max; +} + +/* + * Find the bin index for maximum value of a 2D histogram + */ + +void +gsl_histogram2d_max_bin (const gsl_histogram2d * h, size_t * imax_out, size_t * jmax_out) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t imax = 0, jmax = 0, i, j; + double max = h->bin[0 * ny + 0]; + + for (i = 0; i < nx; i++) + { + for (j = 0; j < ny; j++) + { + double x = h->bin[i * ny + j]; + + if (x > max) + { + max = x; + imax = i; + jmax = j; + } + } + } + + *imax_out = imax; + *jmax_out = jmax; +} + +/* + * Return the minimum contents value of a 2D histogram + */ + +double +gsl_histogram2d_min_val (const gsl_histogram2d * h) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + double min = h->bin[0 * ny + 0]; + + for (i = 0; i < nx * ny; i++) + { + if (h->bin[i] < min) + { + min = h->bin[i]; + } + } + + return min; +} +/* + * Find the bin index for minimum value of a 2D histogram + */ +void +gsl_histogram2d_min_bin (const gsl_histogram2d * h, size_t * imin_out, size_t * jmin_out) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t imin = 0, jmin = 0, i, j; + double min = h->bin[0 * ny + 0]; + + for (i = 0; i < nx; i++) + { + for (j = 0; j < ny; j++) + { + double x = h->bin[i * ny + j]; + + if (x < min) + { + min = x; + imin = i; + jmin = j; + } + } + } + + *imin_out = imin; + *jmin_out = jmin; +} + diff --git a/gsl-1.9/histogram/oper.c b/gsl-1.9/histogram/oper.c new file mode 100644 index 0000000..1d96a6d --- /dev/null +++ b/gsl-1.9/histogram/oper.c @@ -0,0 +1,195 @@ +/* gsl_histogram_oper.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram_oper.c: + * Routine to make operation on histograms. + * Need GSL library and header. + * Contains the routines: + * gsl_histogram_same_binning check if two histograms have the same binning + * gsl_histogram_add add two histograms + * gsl_histogram_sub subctract two histograms + * gsl_histogram_mult multiply two histograms + * gsl_histogram_div divide two histograms + * gsl_histogram_scale scale histogram contents + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +/* + * gsl_histogram_same_binning: + * control if two histograms have the + * same binning + */ + +int +gsl_histogram_equal_bins_p (const gsl_histogram * h1, const gsl_histogram * h2) +{ + if (h1->n != h2->n) + { + return 0; + } + + { + size_t i; + /* init ranges */ + + for (i = 0; i <= h1->n; i++) + { + if (h1->range[i] != h2->range[i]) + { + return 0; + } + } + } + + return 1; +} + +/* + * gsl_histogram_add: + * add two histograms + */ +int +gsl_histogram_add (gsl_histogram * h1, const gsl_histogram * h2) +{ + size_t i; + + if (!gsl_histogram_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < h1->n; i++) + { + h1->bin[i] += h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram_sub: + * subtract two histograms + */ + +int +gsl_histogram_sub (gsl_histogram * h1, const gsl_histogram * h2) +{ + size_t i; + + if (!gsl_histogram_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < h1->n; i++) + { + h1->bin[i] -= h2->bin[i]; + } + + return GSL_SUCCESS; + +} + +/* + * gsl_histogram_mult: + * multiply two histograms + */ + +int +gsl_histogram_mul (gsl_histogram * h1, const gsl_histogram * h2) +{ + size_t i; + + if (!gsl_histogram_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < h1->n; i++) + { + h1->bin[i] *= h2->bin[i]; + } + + return GSL_SUCCESS; +} +/* + * gsl_histogram_div: + * divide two histograms + */ +int +gsl_histogram_div (gsl_histogram * h1, const gsl_histogram * h2) +{ + size_t i; + + if (!gsl_histogram_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < h1->n; i++) + { + h1->bin[i] /= h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram_scale: + * scale a histogram by a numeric factor + */ + +int +gsl_histogram_scale (gsl_histogram * h, double scale) +{ + size_t i; + + for (i = 0; i < h->n; i++) + { + h->bin[i] *= scale; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram_shift: + * shift a histogram by a numeric offset + */ + +int +gsl_histogram_shift (gsl_histogram * h, double shift) +{ + size_t i; + + for (i = 0; i < h->n; i++) + { + h->bin[i] += shift; + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/histogram/oper2d.c b/gsl-1.9/histogram/oper2d.c new file mode 100644 index 0000000..054f78f --- /dev/null +++ b/gsl-1.9/histogram/oper2d.c @@ -0,0 +1,203 @@ +/* gsl_histogram2d_oper.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram2d_oper.c: + * Routine to make operation on 2D histograms. + * Need GSL library and header. + * Contains the routines: + * gsl_histogram2d_same_binning check if two histograms have the same binning + * gsl_histogram2d_add add two histogram + * gsl_histogram2d_sub subctract two histogram + * gsl_histogram2d_mult multiply two histogram + * gsl_histogram2d_div divide two histogram + * gsl_histogram2d_scale scale histogram contents + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +/* + * gsl_histogram2d_same_binning: + * control if two histogram have the + * same binning + */ +int +gsl_histogram2d_equal_bins_p (const gsl_histogram2d * h1, + const gsl_histogram2d * h2) +{ + + if ((h1->nx != h2->nx) || (h1->ny != h2->ny)) + { + return 0; + } + { + size_t i; + /* init ranges */ + for (i = 0; i <= (h1->nx); i++) + { + if (h1->xrange[i] != h2->xrange[i]) + { + return 0; + } + } + for (i = 0; i <= (h1->ny); i++) + { + if (h1->yrange[i] != h2->yrange[i]) + { + return 0; + } + } + } + return 1; +} + +/* + * gsl_histogram2d_add: + * add two histogram + */ + +int +gsl_histogram2d_add (gsl_histogram2d * h1, const gsl_histogram2d * h2) +{ + size_t i; + + if (!gsl_histogram2d_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < (h1->nx) * (h1->ny); i++) + { + h1->bin[i] += h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_sub: + * subtract two histogram + */ + +int +gsl_histogram2d_sub (gsl_histogram2d * h1, const gsl_histogram2d * h2) +{ + size_t i; + + if (!gsl_histogram2d_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < (h1->nx) * (h1->ny); i++) + { + h1->bin[i] -= h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_mult: + * multiply two histogram + */ + +int +gsl_histogram2d_mul (gsl_histogram2d * h1, const gsl_histogram2d * h2) +{ + size_t i; + + if (!gsl_histogram2d_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < (h1->nx) * (h1->ny); i++) + { + h1->bin[i] *= h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_div: + * divide two histogram + */ + +int +gsl_histogram2d_div (gsl_histogram2d * h1, const gsl_histogram2d * h2) +{ + size_t i; + + if (!gsl_histogram2d_equal_bins_p (h1, h2)) + { + GSL_ERROR ("histograms have different binning", GSL_EINVAL); + } + + for (i = 0; i < (h1->nx) * (h1->ny); i++) + { + h1->bin[i] /= h2->bin[i]; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_scale: + * scale a histogram by a numeric factor + */ + +int +gsl_histogram2d_scale (gsl_histogram2d * h, double scale) +{ + size_t i; + + for (i = 0; i < (h->nx) * (h->ny); i++) + { + h->bin[i] *= scale; + } + + return GSL_SUCCESS; +} + +/* + * gsl_histogram2d_shift: + * shift a histogram by a numeric offset + */ + +int +gsl_histogram2d_shift (gsl_histogram2d * h, double shift) +{ + size_t i; + + for (i = 0; i < (h->nx) * (h->ny); i++) + { + h->bin[i] += shift; + } + + return GSL_SUCCESS; +} + diff --git a/gsl-1.9/histogram/params.c b/gsl-1.9/histogram/params.c new file mode 100644 index 0000000..9e2ca12 --- /dev/null +++ b/gsl-1.9/histogram/params.c @@ -0,0 +1,42 @@ +/* histogram/params.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +double +gsl_histogram_max (const gsl_histogram * h) +{ + const int n = h->n; + + return h->range[n]; +} + +double +gsl_histogram_min (const gsl_histogram * h) +{ + return h->range[0]; +} + +size_t +gsl_histogram_bins (const gsl_histogram * h) +{ + return h->n; +} diff --git a/gsl-1.9/histogram/params2d.c b/gsl-1.9/histogram/params2d.c new file mode 100644 index 0000000..3ca946e --- /dev/null +++ b/gsl-1.9/histogram/params2d.c @@ -0,0 +1,60 @@ +/* histogram/params2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +double +gsl_histogram2d_xmax (const gsl_histogram2d * h) +{ + const int nx = h->nx; + return h->xrange[nx]; +} + +double +gsl_histogram2d_xmin (const gsl_histogram2d * h) +{ + return h->xrange[0]; +} + +double +gsl_histogram2d_ymax (const gsl_histogram2d * h) +{ + const int ny = h->ny; + return h->yrange[ny]; +} + +double +gsl_histogram2d_ymin (const gsl_histogram2d * h) +{ + return h->yrange[0]; +} + +size_t +gsl_histogram2d_nx (const gsl_histogram2d * h) +{ + return h->nx; +} + +size_t +gsl_histogram2d_ny (const gsl_histogram2d * h) +{ + return h->ny; +} diff --git a/gsl-1.9/histogram/pdf.c b/gsl-1.9/histogram/pdf.c new file mode 100644 index 0000000..78f9a22 --- /dev/null +++ b/gsl-1.9/histogram/pdf.c @@ -0,0 +1,153 @@ +/* histogram/pdf.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +#include "find.c" + +double +gsl_histogram_pdf_sample (const gsl_histogram_pdf * p, double r) +{ + size_t i; + int status; + +/* Wrap the exclusive top of the bin down to the inclusive bottom of + the bin. Since this is a single point it should not affect the + distribution. */ + + if (r == 1.0) + { + r = 0.0; + } + + status = find (p->n, p->sum, r, &i); + + if (status) + { + GSL_ERROR_VAL ("cannot find r in cumulative pdf", GSL_EDOM, 0); + } + else + { + double delta = (r - p->sum[i]) / (p->sum[i + 1] - p->sum[i]); + double x = p->range[i] + delta * (p->range[i + 1] - p->range[i]); + return x; + } +} + +gsl_histogram_pdf * +gsl_histogram_pdf_alloc (const size_t n) +{ + gsl_histogram_pdf *p; + + if (n == 0) + { + GSL_ERROR_VAL ("histogram pdf length n must be positive integer", + GSL_EDOM, 0); + } + + p = (gsl_histogram_pdf *) malloc (sizeof (gsl_histogram_pdf)); + + if (p == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram pdf struct", + GSL_ENOMEM, 0); + } + + p->range = (double *) malloc ((n + 1) * sizeof (double)); + + if (p->range == 0) + { + free (p); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram pdf ranges", + GSL_ENOMEM, 0); + } + + p->sum = (double *) malloc ((n + 1) * sizeof (double)); + + if (p->sum == 0) + { + free (p->range); + free (p); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram pdf sums", + GSL_ENOMEM, 0); + } + + p->n = n; + + return p; +} + +int +gsl_histogram_pdf_init (gsl_histogram_pdf * p, const gsl_histogram * h) +{ + size_t i; + size_t n = p->n; + + if (n != h->n) + { + GSL_ERROR ("histogram length must match pdf length", GSL_EINVAL); + } + + for (i = 0; i < n; i++) + { + if (h->bin[i] < 0) + { + GSL_ERROR ("histogram bins must be non-negative to compute" + "a probability distribution", GSL_EDOM); + } + } + + for (i = 0; i < n + 1; i++) + { + p->range[i] = h->range[i]; + } + + { + double mean = 0, sum = 0; + + for (i = 0; i < n; i++) + { + mean += (h->bin[i] - mean) / ((double) (i + 1)); + } + + p->sum[0] = 0; + + for (i = 0; i < n; i++) + { + sum += (h->bin[i] / mean) / n; + p->sum[i + 1] = sum; + } + } + + return GSL_SUCCESS; +} + + +void +gsl_histogram_pdf_free (gsl_histogram_pdf * p) +{ + free (p->range); + free (p->sum); + free (p); +} diff --git a/gsl-1.9/histogram/pdf2d.c b/gsl-1.9/histogram/pdf2d.c new file mode 100644 index 0000000..f62dea5 --- /dev/null +++ b/gsl-1.9/histogram/pdf2d.c @@ -0,0 +1,186 @@ +/* histogram/pdf2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> +#include <gsl/gsl_histogram2d.h> + +#include "find.c" + +int +gsl_histogram2d_pdf_sample (const gsl_histogram2d_pdf * p, + double r1, double r2, + double *x, double *y) +{ + size_t k; + int status; + +/* Wrap the exclusive top of the bin down to the inclusive bottom of + the bin. Since this is a single point it should not affect the + distribution. */ + + if (r2 == 1.0) + { + r2 = 0.0; + } + if (r1 == 1.0) + { + r1 = 0.0; + } + + status = find (p->nx * p->ny, p->sum, r1, &k); + + if (status) + { + GSL_ERROR ("cannot find r1 in cumulative pdf", GSL_EDOM); + } + else + { + size_t i = k / p->ny; + size_t j = k - (i * p->ny); + double delta = (r1 - p->sum[k]) / (p->sum[k + 1] - p->sum[k]); + *x = p->xrange[i] + delta * (p->xrange[i + 1] - p->xrange[i]); + *y = p->yrange[j] + r2 * (p->yrange[j + 1] - p->yrange[j]); + return GSL_SUCCESS; + } +} + +gsl_histogram2d_pdf * +gsl_histogram2d_pdf_alloc (const size_t nx, const size_t ny) +{ + const size_t n = nx * ny; + + gsl_histogram2d_pdf *p; + + if (n == 0) + { + GSL_ERROR_VAL ("histogram2d pdf length n must be positive integer", + GSL_EDOM, 0); + } + + p = (gsl_histogram2d_pdf *) malloc (sizeof (gsl_histogram2d_pdf)); + + if (p == 0) + { + GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf struct", + GSL_ENOMEM, 0); + } + + p->xrange = (double *) malloc ((nx + 1) * sizeof (double)); + + if (p->xrange == 0) + { + free (p); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf xranges", + GSL_ENOMEM, 0); + } + + p->yrange = (double *) malloc ((ny + 1) * sizeof (double)); + + if (p->yrange == 0) + { + free (p->xrange); + free (p); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf yranges", + GSL_ENOMEM, 0); + } + + p->sum = (double *) malloc ((n + 1) * sizeof (double)); + + if (p->sum == 0) + { + free (p->yrange); + free (p->xrange); + free (p); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for histogram2d pdf sums", + GSL_ENOMEM, 0); + } + + p->nx = nx; + p->ny = ny; + + return p; +} + +int +gsl_histogram2d_pdf_init (gsl_histogram2d_pdf * p, const gsl_histogram2d * h) +{ + size_t i; + const size_t nx = p->nx; + const size_t ny = p->ny; + const size_t n = nx * ny; + + if (nx != h->nx || ny != h->ny) + { + GSL_ERROR ("histogram2d size must match pdf size", GSL_EDOM); + } + + for (i = 0; i < n; i++) + { + if (h->bin[i] < 0) + { + GSL_ERROR ("histogram bins must be non-negative to compute" + "a probability distribution", GSL_EDOM); + } + } + + for (i = 0; i < nx + 1; i++) + { + p->xrange[i] = h->xrange[i]; + } + + for (i = 0; i < ny + 1; i++) + { + p->yrange[i] = h->yrange[i]; + } + + { + double mean = 0, sum = 0; + + for (i = 0; i < n; i++) + { + mean += (h->bin[i] - mean) / ((double) (i + 1)); + } + + p->sum[0] = 0; + + for (i = 0; i < n; i++) + { + sum += (h->bin[i] / mean) / n; + p->sum[i + 1] = sum; + } + } + + return GSL_SUCCESS; +} + + +void +gsl_histogram2d_pdf_free (gsl_histogram2d_pdf * p) +{ + free (p->xrange); + free (p->yrange); + free (p->sum); + free (p); +} diff --git a/gsl-1.9/histogram/reset.c b/gsl-1.9/histogram/reset.c new file mode 100644 index 0000000..1206122 --- /dev/null +++ b/gsl-1.9/histogram/reset.c @@ -0,0 +1,34 @@ +/* histogram/reset.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +void +gsl_histogram_reset (gsl_histogram * h) +{ + size_t i; + const size_t n = h->n; + + for (i = 0; i < n; i++) + { + h->bin[i] = 0; + } +} diff --git a/gsl-1.9/histogram/reset2d.c b/gsl-1.9/histogram/reset2d.c new file mode 100644 index 0000000..160b1da --- /dev/null +++ b/gsl-1.9/histogram/reset2d.c @@ -0,0 +1,35 @@ +/* histogram/reset2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +void +gsl_histogram2d_reset (gsl_histogram2d * h) +{ + size_t i; + const size_t nx = h->nx; + const size_t ny = h->ny; + + for (i = 0; i < nx * ny; i++) + { + h->bin[i] = 0; + } +} diff --git a/gsl-1.9/histogram/stat.c b/gsl-1.9/histogram/stat.c new file mode 100644 index 0000000..5ee9e9c --- /dev/null +++ b/gsl-1.9/histogram/stat.c @@ -0,0 +1,140 @@ +/* gsl_histogram_stat.c + * Copyright (C) 2000 Simone Piccardi + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +/*************************************************************** + * + * File gsl_histogram_stat.c: + * Routines for statisticalcomputations on histograms. + * Need GSL library and header. + * Contains the routines: + * gsl_histogram_mean compute histogram mean + * gsl_histogram_sigma compute histogram sigma + * + * Author: S. Piccardi + * Jan. 2000 + * + ***************************************************************/ +#include <config.h> +#include <stdlib.h> +#include <math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram.h> + +/* FIXME: We skip negative values in the histogram h->bin[i] < 0, + since those correspond to negative weights (BJG) */ + +double +gsl_histogram_mean (const gsl_histogram * h) +{ + const size_t n = h->n; + size_t i; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wmean = 0; + long double W = 0; + + for (i = 0; i < n; i++) + { + double xi = (h->range[i + 1] + h->range[i]) / 2; + double wi = h->bin[i]; + + if (wi > 0) + { + W += wi; + wmean += (xi - wmean) * (wi / W); + } + } + + return wmean; +} + +double +gsl_histogram_sigma (const gsl_histogram * h) +{ + const size_t n = h->n; + size_t i; + + long double wvariance = 0 ; + long double wmean = 0; + long double W = 0; + + /* FIXME: should use a single pass formula here, as given in + N.J.Higham 'Accuracy and Stability of Numerical Methods', p.12 */ + + /* Compute the mean */ + + for (i = 0; i < n; i++) + { + double xi = (h->range[i + 1] + h->range[i]) / 2; + double wi = h->bin[i]; + + if (wi > 0) + { + W += wi; + wmean += (xi - wmean) * (wi / W); + } + } + + /* Compute the variance */ + + W = 0.0; + + for (i = 0; i < n; i++) + { + double xi = ((h->range[i + 1]) + (h->range[i])) / 2; + double wi = h->bin[i]; + + if (wi > 0) { + const long double delta = (xi - wmean); + W += wi ; + wvariance += (delta * delta - wvariance) * (wi / W); + } + } + + { + double sigma = sqrt (wvariance) ; + return sigma; + } +} + + +/* + sum up all bins of histogram + */ + +double +gsl_histogram_sum(const gsl_histogram * h) +{ + double sum=0; + size_t i=0; + size_t n; + n=h->n; + + while(i < n) + sum += h->bin[i++]; + + return sum; +} + diff --git a/gsl-1.9/histogram/stat2d.c b/gsl-1.9/histogram/stat2d.c new file mode 100644 index 0000000..5ad7adf --- /dev/null +++ b/gsl-1.9/histogram/stat2d.c @@ -0,0 +1,266 @@ +/* histogram/stat2d.c + * Copyright (C) 2002 Achim Gaedke + * + * This library 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 library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +/*************************************************************** + * + * File histogram/stat2d.c: + * Routine to return statistical values of the content of a 2D hisogram. + * + * Contains the routines: + * gsl_histogram2d_sum sum up all bin values + * gsl_histogram2d_xmean determine mean of x values + * gsl_histogram2d_ymean determine mean of y values + * + * Author: Achim Gaedke Achim.Gaedke@zpr.uni-koeln.de + * Jan. 2002 + * + ***************************************************************/ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_histogram2d.h> + +/* + sum up all bins of histogram2d + */ + +double +gsl_histogram2d_sum (const gsl_histogram2d * h) +{ + const size_t n = h->nx * h->ny; + double sum = 0; + size_t i = 0; + + while (i < n) + sum += h->bin[i++]; + + return sum; +} + +double +gsl_histogram2d_xmean (const gsl_histogram2d * h) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + size_t j; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wmean = 0; + long double W = 0; + + for (i = 0; i < nx; i++) + { + double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0; + double wi = 0; + + for (j = 0; j < ny; j++) + { + double wij = h->bin[i * ny + j]; + if (wij > 0) + wi += wij; + } + if (wi > 0) + { + W += wi; + wmean += (xi - wmean) * (wi / W); + } + } + + return wmean; +} + +double +gsl_histogram2d_ymean (const gsl_histogram2d * h) +{ + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + size_t j; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wmean = 0; + long double W = 0; + + for (j = 0; j < ny; j++) + { + double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0; + double wj = 0; + + for (i = 0; i < nx; i++) + { + double wij = h->bin[i * ny + j]; + if (wij > 0) + wj += wij; + } + + if (wj > 0) + { + W += wj; + wmean += (yj - wmean) * (wj / W); + } + } + + return wmean; +} + +double +gsl_histogram2d_xsigma (const gsl_histogram2d * h) +{ + const double xmean = gsl_histogram2d_xmean (h); + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + size_t j; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wvariance = 0; + long double W = 0; + + for (i = 0; i < nx; i++) + { + double xi = (h->xrange[i + 1] + h->xrange[i]) / 2 - xmean; + double wi = 0; + + for (j = 0; j < ny; j++) + { + double wij = h->bin[i * ny + j]; + if (wij > 0) + wi += wij; + } + + if (wi > 0) + { + W += wi; + wvariance += ((xi * xi) - wvariance) * (wi / W); + } + } + + { + double xsigma = sqrt (wvariance); + return xsigma; + } +} + +double +gsl_histogram2d_ysigma (const gsl_histogram2d * h) +{ + const double ymean = gsl_histogram2d_ymean (h); + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + size_t j; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wvariance = 0; + long double W = 0; + + for (j = 0; j < ny; j++) + { + double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean; + double wj = 0; + + for (i = 0; i < nx; i++) + { + double wij = h->bin[i * ny + j]; + if (wij > 0) + wj += wij; + } + if (wj > 0) + { + W += wj; + wvariance += ((yj * yj) - wvariance) * (wj / W); + } + } + + { + double ysigma = sqrt (wvariance); + return ysigma; + } +} + +double +gsl_histogram2d_cov (const gsl_histogram2d * h) +{ + const double xmean = gsl_histogram2d_xmean (h); + const double ymean = gsl_histogram2d_ymean (h); + const size_t nx = h->nx; + const size_t ny = h->ny; + size_t i; + size_t j; + + /* Compute the bin-weighted arithmetic mean M of a histogram using the + recurrence relation + + M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) + W(n) = W(n-1) + w(n) + + */ + + long double wcovariance = 0; + long double W = 0; + + for (j = 0; j < ny; j++) + { + for (i = 0; i < nx; i++) + { + double xi = (h->xrange[i + 1] + h->xrange[i]) / 2.0 - xmean; + double yj = (h->yrange[j + 1] + h->yrange[j]) / 2.0 - ymean; + double wij = h->bin[i * ny + j]; + + if (wij > 0) + { + W += wij; + wcovariance += ((xi * yj) - wcovariance) * (wij / W); + } + } + } + + return wcovariance; + +} diff --git a/gsl-1.9/histogram/test.c b/gsl-1.9/histogram/test.c new file mode 100644 index 0000000..c170de8 --- /dev/null +++ b/gsl-1.9/histogram/test.c @@ -0,0 +1,45 @@ +/* histogram/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_histogram.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +void test1d (void); +void test2d (void); +void test1d_resample (void); +void test2d_resample (void); +void test1d_trap (void); +void test2d_trap (void); + +int +main (void) +{ + test1d(); + test2d(); + test1d_resample(); + test2d_resample(); + test1d_trap(); + test2d_trap(); + + exit (gsl_test_summary ()); +} diff --git a/gsl-1.9/histogram/test1d.c b/gsl-1.9/histogram/test1d.c new file mode 100644 index 0000000..ff2d89a --- /dev/null +++ b/gsl-1.9/histogram/test1d.c @@ -0,0 +1,483 @@ +/* histogram/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_histogram.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +#define N 397 +#define NR 10 + +void +test1d (void) +{ + double xr[NR + 1] = + {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}; + + gsl_histogram *h, *h1, *hr, *g; + size_t i, j; + + gsl_ieee_env_setup (); + + h = gsl_histogram_calloc (N); + h1 = gsl_histogram_calloc (N); + g = gsl_histogram_calloc (N); + + gsl_test (h->range == 0, "gsl_histogram_alloc returns valid range pointer"); + gsl_test (h->bin == 0, "gsl_histogram_alloc returns valid bin pointer"); + gsl_test (h->n != N, "gsl_histogram_alloc returns valid size"); + + + hr = gsl_histogram_calloc_range (NR, xr); + + gsl_test (hr->range == 0, "gsl_histogram_calloc_range returns valid range pointer"); + gsl_test (hr->bin == 0, "gsl_histogram_calloc_range returns valid bin pointer"); + gsl_test (hr->n != NR, "gsl_histogram_calloc_range returns valid size"); + + { + int status = 0; + for (i = 0; i <= NR; i++) + { + if (hr->range[i] != xr[i]) + { + status = 1; + } + }; + + gsl_test (status, "gsl_histogram_calloc_range creates range"); + } + + for (i = 0; i <= NR; i++) + { + hr->range[i] = 0.0; + } + + { + int status = gsl_histogram_set_ranges (hr, xr, NR+1); + + for (i = 0; i <= NR; i++) + { + if (hr->range[i] != xr[i]) + { + status = 1; + } + }; + + gsl_test (status, "gsl_histogram_set_range sets range"); + } + + + for (i = 0; i < N; i++) + { + gsl_histogram_accumulate (h, (double) i, (double) i); + }; + + { + int status = 0; + + for (i = 0; i < N; i++) + { + if (h->bin[i] != (double) i) + { + status = 1; + } + }; + + gsl_test (status, "gsl_histogram_accumulate writes into array"); + } + + + { + int status = 0; + + for (i = 0; i < N; i++) + { + if (gsl_histogram_get (h, i) != i) + status = 1; + }; + gsl_test (status, "gsl_histogram_get reads from array"); + } + + for (i = 0; i <= N; i++) + { + h1->range[i] = 100.0 + i; + } + + gsl_histogram_memcpy (h1, h); + + { + int status = 0; + for (i = 0; i <= N; i++) + { + if (h1->range[i] != h->range[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram_memcpy copies bin ranges"); + } + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (gsl_histogram_get (h1, i) != gsl_histogram_get (h, i)) + status = 1; + }; + gsl_test (status, "gsl_histogram_memcpy copies bin values"); + } + + gsl_histogram_free (h1); + + h1 = gsl_histogram_clone (h); + + { + int status = 0; + for (i = 0; i <= N; i++) + { + if (h1->range[i] != h->range[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram_clone copies bin ranges"); + } + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (gsl_histogram_get (h1, i) != gsl_histogram_get (h, i)) + status = 1; + }; + gsl_test (status, "gsl_histogram_clone copies bin values"); + } + + gsl_histogram_reset (h); + + { + int status = 0; + + for (i = 0; i < N; i++) + { + if (h->bin[i] != 0) + status = 1; + } + gsl_test (status, "gsl_histogram_reset zeros array"); + } + + + { + int status = 0; + + for (i = 0; i < N; i++) + { + gsl_histogram_increment (h, (double) i); + + for (j = 0; j <= i; j++) + { + if (h->bin[j] != 1) + { + status = 1; + } + } + + for (j = i + 1; j < N; j++) + { + if (h->bin[j] != 0) + { + status = 1; + } + } + } + + gsl_test (status, "gsl_histogram_increment increases bin value"); + } + + { + int status = 0; + for (i = 0; i < N; i++) + { + double x0 = 0, x1 = 0; + + gsl_histogram_get_range (h, i, &x0, &x1); + + if (x0 != i || x1 != i + 1) + { + status = 1; + } + } + gsl_test (status, "gsl_histogram_getbinrange returns bin range"); + } + + { + int status = 0; + if (gsl_histogram_max (h) != N) + status = 1; + gsl_test (status, "gsl_histogram_max returns maximum"); + } + + { + int status = 0; + if (gsl_histogram_min (h) != 0) + status = 1; + gsl_test (status, "gsl_histogram_min returns minimum"); + } + + { + int status = 0; + if (gsl_histogram_bins (h) != N) + status = 1; + gsl_test (status, "gsl_histogram_bins returns number of bins"); + } + + h->bin[2] = 123456.0; + h->bin[4] = -654321; + + { + double max = gsl_histogram_max_val (h); + gsl_test (max != 123456.0, "gsl_histogram_max_val finds maximum value"); + } + + { + double min = gsl_histogram_min_val (h); + gsl_test (min != -654321.0, "gsl_histogram_min_val finds minimum value"); + } + + { + size_t imax = gsl_histogram_max_bin (h); + gsl_test (imax != 2, "gsl_histogram_max_bin finds maximum value bin"); + } + + { + size_t imin = gsl_histogram_min_bin (h); + gsl_test (imin != 4, "gsl_histogram_min_bin find minimum value bin"); + } + + for (i = 0; i < N; i++) + { + h->bin[i] = i + 27; + g->bin[i] = (i + 27) * (i + 1); + } + + { + double sum=gsl_histogram_sum (h); + gsl_test(sum != N*27+((N-1)*N)/2, "gsl_histogram_sum sums all bin values"); + } + + gsl_histogram_memcpy (h1, g); + gsl_histogram_add (h1, h); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != g->bin[i] + h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_add histogram addition"); + } + + gsl_histogram_memcpy (h1, g); + gsl_histogram_sub (h1, h); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != g->bin[i] - h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_sub histogram subtraction"); + } + + + gsl_histogram_memcpy (h1, g); + gsl_histogram_mul (h1, h); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != g->bin[i] * h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_mul histogram multiplication"); + } + + + gsl_histogram_memcpy (h1, g); + gsl_histogram_div (h1, h); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != g->bin[i] / h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_div histogram division"); + } + + gsl_histogram_memcpy (h1, g); + gsl_histogram_scale (h1, 0.5); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != 0.5 * g->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_scale histogram scaling"); + } + + gsl_histogram_memcpy (h1, g); + gsl_histogram_shift (h1, 0.25); + + { + int status = 0; + for (i = 0; i < N; i++) + { + if (h1->bin[i] != 0.25 + g->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram_shift histogram shift"); + } + + + gsl_histogram_free (h); /* free whatever is in h */ + + h = gsl_histogram_calloc_uniform (N, 0.0, 1.0); + + gsl_test (h->range == 0, + "gsl_histogram_calloc_uniform returns valid range pointer"); + gsl_test (h->bin == 0, + "gsl_histogram_calloc_uniform returns valid bin pointer"); + gsl_test (h->n != N, + "gsl_histogram_calloc_uniform returns valid size"); + + gsl_histogram_accumulate (h, 0.0, 1.0); + gsl_histogram_accumulate (h, 0.1, 2.0); + gsl_histogram_accumulate (h, 0.2, 3.0); + gsl_histogram_accumulate (h, 0.3, 4.0); + + { + size_t i1, i2, i3, i4; + double expected; + int status = gsl_histogram_find (h, 0.0, &i1); + status = gsl_histogram_find (h, 0.1, &i2); + status = gsl_histogram_find (h, 0.2, &i3); + status = gsl_histogram_find (h, 0.3, &i4); + + for (i = 0; i < N; i++) + { + if (i == i1) + { + expected = 1.0; + } + else if (i == i2) + { + expected = 2.0; + } + else if (i == i3) + { + expected = 3.0; + } + else if (i == i4) + { + expected = 4.0; + } + else + { + expected = 0.0; + } + + if (h->bin[i] != expected) + { + status = 1; + } + + } + gsl_test (status, "gsl_histogram_find returns index"); + } + + + { + FILE *f = fopen ("test.txt", "w"); + gsl_histogram_fprintf (f, h, "%.19e", "%.19e"); + fclose (f); + } + + { + FILE *f = fopen ("test.txt", "r"); + gsl_histogram *hh = gsl_histogram_calloc (N); + int status = 0; + + gsl_histogram_fscanf (f, hh); + + for (i = 0; i < N; i++) + { + if (h->range[i] != hh->range[i]) + status = 1; + if (h->bin[i] != hh->bin[i]) + status = 1; + } + if (h->range[N] != hh->range[N]) + status = 1; + + gsl_test (status, "gsl_histogram_fprintf and fscanf"); + + gsl_histogram_free (hh); + fclose (f); + } + + { + FILE *f = fopen ("test.dat", "wb"); + gsl_histogram_fwrite (f, h); + fclose (f); + } + + { + FILE *f = fopen ("test.dat", "rb"); + gsl_histogram *hh = gsl_histogram_calloc (N); + int status = 0; + + gsl_histogram_fread (f, hh); + + for (i = 0; i < N; i++) + { + if (h->range[i] != hh->range[i]) + status = 1; + if (h->bin[i] != hh->bin[i]) + status = 1; + } + if (h->range[N] != hh->range[N]) + status = 1; + + gsl_test (status, "gsl_histogram_fwrite and fread"); + + gsl_histogram_free (hh); + fclose (f); + } + + gsl_histogram_free (h); + gsl_histogram_free (g); + gsl_histogram_free (h1); + gsl_histogram_free (hr); +} diff --git a/gsl-1.9/histogram/test1d_resample.c b/gsl-1.9/histogram/test1d_resample.c new file mode 100644 index 0000000..525c126 --- /dev/null +++ b/gsl-1.9/histogram/test1d_resample.c @@ -0,0 +1,99 @@ +/* histogram/test1d_resample.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_histogram.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +#include "urand.c" + +void +test1d_resample (void) +{ + size_t i; + int status = 0; + + gsl_histogram *h; + + gsl_ieee_env_setup (); + + h = gsl_histogram_calloc_uniform (10, 0.0, 1.0); + + gsl_histogram_increment (h, 0.1); + gsl_histogram_increment (h, 0.2); + gsl_histogram_increment (h, 0.2); + gsl_histogram_increment (h, 0.3); + + { + gsl_histogram_pdf *p = gsl_histogram_pdf_alloc (10); + + gsl_histogram *hh = gsl_histogram_calloc_uniform (100, 0.0, 1.0); + + gsl_histogram_pdf_init (p, h); + + for (i = 0; i < 100000; i++) + { + double u = urand(); + double x = gsl_histogram_pdf_sample (p, u); + gsl_histogram_increment (hh, x); + } + + for (i = 0; i < 100; i++) + { + double y = gsl_histogram_get (hh, i) / 2500; + double x, xmax; + size_t k; + double ya; + + gsl_histogram_get_range (hh, i, &x, &xmax); + + gsl_histogram_find (h, x, &k); + ya = gsl_histogram_get (h, k); + + if (ya == 0) + { + if (y != 0) + { + printf ("%d: %g vs %g\n", (int) i, y, ya); + status = 1; + } + } + else + { + double err = 1 / sqrt (gsl_histogram_get (hh, i)); + double sigma = fabs ((y - ya) / (ya * err)); + if (sigma > 3) + { + status = 1; + printf ("%g vs %g err=%g sigma=%g\n", y, ya, err, sigma); + } + } + } + + gsl_histogram_pdf_free (p) ; + gsl_histogram_free (hh); + + gsl_test (status, "gsl_histogram_pdf_sample within statistical errors"); + } + + gsl_histogram_free (h); +} diff --git a/gsl-1.9/histogram/test1d_trap.c b/gsl-1.9/histogram/test1d_trap.c new file mode 100644 index 0000000..511cb1f --- /dev/null +++ b/gsl-1.9/histogram/test1d_trap.c @@ -0,0 +1,130 @@ +/* histogram/test_trap.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_histogram.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_ieee_utils.h> + +#define N 397 + +static void my_error_handler (const char *reason, const char *file, + int line, int err); +static int status = 0; + +void +test1d_trap (void) +{ + gsl_histogram *h; + double result, lower, upper; + size_t i; + + gsl_set_error_handler (&my_error_handler); + + gsl_ieee_env_setup (); + + status = 0; + h = gsl_histogram_calloc (0); + gsl_test (!status, "gsl_histogram_calloc traps zero-length histogram"); + gsl_test (h != 0, + "gsl_histogram_calloc returns NULL for zero-length histogram"); + + status = 0; + h = gsl_histogram_calloc_uniform (0, 0.0, 1.0); + gsl_test (!status, + "gsl_histogram_calloc_uniform traps zero-length histogram"); + gsl_test (h != 0, + "gsl_histogram_calloc_uniform returns NULL for zero-length histogram"); + + status = 0; + h = gsl_histogram_calloc_uniform (10, 1.0, 1.0); + gsl_test (!status, + "gsl_histogram_calloc_uniform traps equal endpoints"); + gsl_test (h != 0, + "gsl_histogram_calloc_uniform returns NULL for equal endpoints"); + + status = 0; + h = gsl_histogram_calloc_uniform (10, 2.0, 1.0); + gsl_test (!status, + "gsl_histogram_calloc_uniform traps invalid range"); + gsl_test (h != 0, + "gsl_histogram_calloc_uniform returns NULL for invalid range"); + + h = gsl_histogram_calloc_uniform (N, 0.0, 1.0); + + status = gsl_histogram_accumulate (h, 1.0, 10.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x at xmax"); + + status = gsl_histogram_accumulate (h, 2.0, 100.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x above xmax"); + + status = gsl_histogram_accumulate (h, -1.0, 1000.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_accumulate traps x below xmin"); + + status = gsl_histogram_increment (h, 1.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x at xmax"); + + status = gsl_histogram_increment (h, 2.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x above xmax"); + + status = gsl_histogram_increment (h, -1.0); + gsl_test (status != GSL_EDOM, "gsl_histogram_increment traps x below xmin"); + + + result = gsl_histogram_get (h, N); + gsl_test (result != 0, "gsl_histogram_get traps index at n"); + + result = gsl_histogram_get (h, N + 1); + gsl_test (result != 0, "gsl_histogram_get traps index above n"); + + status = gsl_histogram_get_range (h, N, &lower, &upper); + gsl_test (status != GSL_EDOM, + "gsl_histogram_get_range traps index at n"); + + status = gsl_histogram_get_range (h, N + 1, &lower, &upper); + gsl_test (status != GSL_EDOM, + "gsl_histogram_get_range traps index above n"); + + + status = 0; + gsl_histogram_find (h, -0.01, &i); + gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x below xmin"); + + status = 0; + gsl_histogram_find (h, 1.0, &i); + gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x at xmax"); + + status = 0; + gsl_histogram_find (h, 1.1, &i); + gsl_test (status != GSL_EDOM, "gsl_histogram_find traps x above xmax"); + + gsl_histogram_free (h); +} + + +static 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); + status = 1; +} diff --git a/gsl-1.9/histogram/test2d.c b/gsl-1.9/histogram/test2d.c new file mode 100644 index 0000000..df09155 --- /dev/null +++ b/gsl-1.9/histogram/test2d.c @@ -0,0 +1,753 @@ +/* histogram/test2d.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_machine.h> +#include <gsl/gsl_histogram2d.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +#define M 107 +#define N 239 +#define M1 17 +#define N1 23 +#define MR 10 +#define NR 5 + +void +test2d (void) +{ + double xr[MR + 1] = + { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 }; + + double yr[NR + 1] = { 90.0, 91.0, 92.0, 93.0, 94.0, 95.0 }; + + gsl_histogram2d *h, *h1, *g, *hr; + size_t i, j, k; + + gsl_ieee_env_setup (); + + h = gsl_histogram2d_calloc (M, N); + h1 = gsl_histogram2d_calloc (M, N); + g = gsl_histogram2d_calloc (M, N); + + gsl_test (h->xrange == 0, + "gsl_histogram2d_calloc returns valid xrange pointer"); + gsl_test (h->yrange == 0, + "gsl_histogram2d_calloc returns valid yrange pointer"); + gsl_test (h->bin == 0, "gsl_histogram2d_calloc returns valid bin pointer"); + + gsl_test (h->nx != M, "gsl_histogram2d_calloc returns valid nx"); + gsl_test (h->ny != N, "gsl_histogram2d_calloc returns valid ny"); + + hr = gsl_histogram2d_calloc_range (MR, NR, xr, yr); + + gsl_test (hr->xrange == 0, + "gsl_histogram2d_calloc_range returns valid xrange pointer"); + gsl_test (hr->yrange == 0, + "gsl_histogram2d_calloc_range returns valid yrange pointer"); + gsl_test (hr->bin == 0, + "gsl_histogram2d_calloc_range returns valid bin pointer"); + + gsl_test (hr->nx != MR, "gsl_histogram2d_calloc_range returns valid nx"); + gsl_test (hr->ny != NR, "gsl_histogram2d_calloc_range returns valid ny"); + + { + int status = 0; + for (i = 0; i <= MR; i++) + { + if (hr->xrange[i] != xr[i]) + { + status = 1; + } + }; + + gsl_test (status, + "gsl_histogram2d_calloc_range creates xrange"); + } + + { + int status = 0; + for (i = 0; i <= NR; i++) + { + if (hr->yrange[i] != yr[i]) + { + status = 1; + } + }; + + gsl_test (status, + "gsl_histogram2d_calloc_range creates yrange"); + } + + for (i = 0; i <= MR; i++) + { + hr->xrange[i] = 0.0; + } + + for (i = 0; i <= NR; i++) + { + hr->yrange[i] = 0.0; + } + + { + int status = gsl_histogram2d_set_ranges (hr, xr, MR + 1, yr, NR + 1); + + for (i = 0; i <= MR; i++) + { + if (hr->xrange[i] != xr[i]) + { + status = 1; + } + }; + + gsl_test (status, "gsl_histogram2d_set_ranges sets xrange"); + } + + { + int status = 0; + for (i = 0; i <= NR; i++) + { + if (hr->yrange[i] != yr[i]) + { + status = 1; + } + }; + + gsl_test (status, "gsl_histogram2d_set_ranges sets yrange"); + } + + + k = 0; + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + k++; + gsl_histogram2d_accumulate (h, (double) i, (double) j, (double) k); + }; + } + + { + int status = 0; + k = 0; + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + k++; + if (h->bin[i * N + j] != (double) k) + { + status = 1; + } + } + } + gsl_test (status, + "gsl_histogram2d_accumulate writes into array"); + } + + { + int status = 0; + k = 0; + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + k++; + if (gsl_histogram2d_get (h, i, j) != (double) k) + status = 1; + }; + } + gsl_test (status, "gsl_histogram2d_get reads from array"); + } + + for (i = 0; i <= M; i++) + { + h1->xrange[i] = 100.0 + i; + } + + for (i = 0; i <= N; i++) + { + h1->yrange[i] = 900.0 + i * i; + } + + gsl_histogram2d_memcpy (h1, h); + + { + int status = 0; + for (i = 0; i <= M; i++) + { + if (h1->xrange[i] != h->xrange[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram2d_memcpy copies bin xranges"); + } + + { + int status = 0; + for (i = 0; i <= N; i++) + { + if (h1->yrange[i] != h->yrange[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram2d_memcpy copies bin yranges"); + } + + { + int status = 0; + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + if (gsl_histogram2d_get (h1, i, j) != + gsl_histogram2d_get (h, i, j)) + status = 1; + } + } + gsl_test (status, "gsl_histogram2d_memcpy copies bin values"); + } + + gsl_histogram2d_free (h1); + + h1 = gsl_histogram2d_clone (h); + + { + int status = 0; + for (i = 0; i <= M; i++) + { + if (h1->xrange[i] != h->xrange[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram2d_clone copies bin xranges"); + } + + { + int status = 0; + for (i = 0; i <= N; i++) + { + if (h1->yrange[i] != h->yrange[i]) + status = 1; + }; + gsl_test (status, "gsl_histogram2d_clone copies bin yranges"); + } + + { + int status = 0; + for (i = 0; i < M; i++) + { + for (j = 0; j < N; j++) + { + if (gsl_histogram2d_get (h1, i, j) != + gsl_histogram2d_get (h, i, j)) + status = 1; + } + } + gsl_test (status, "gsl_histogram2d_clone copies bin values"); + } + + + gsl_histogram2d_reset (h); + + { + int status = 0; + + for (i = 0; i < M * N; i++) + { + if (h->bin[i] != 0) + status = 1; + } + gsl_test (status, "gsl_histogram2d_reset zeros array"); + } + + gsl_histogram2d_free (h); + h = gsl_histogram2d_calloc (M1, N1); + + { + + int status = 0; + + for (i = 0; i < M1; i++) + { + for (j = 0; j < N1; j++) + { + gsl_histogram2d_increment (h, (double) i, (double) j); + + for (k = 0; k <= i * N1 + j; k++) + { + if (h->bin[k] != 1) + { + status = 1; + } + } + + for (k = i * N1 + j + 1; k < M1 * N1; k++) + { + if (h->bin[k] != 0) + { + status = 1; + } + } + } + } + gsl_test (status, "gsl_histogram2d_increment increases bin value"); + } + + gsl_histogram2d_free (h); + h = gsl_histogram2d_calloc (M, N); + + { + int status = 0; + for (i = 0; i < M; i++) + { + double x0 = 0, x1 = 0; + gsl_histogram2d_get_xrange (h, i, &x0, &x1); + + if (x0 != i || x1 != i + 1) + { + status = 1; + } + } + gsl_test (status, + "gsl_histogram2d_get_xlowerlimit and xupperlimit"); + } + + + { + int status = 0; + for (i = 0; i < N; i++) + { + double y0 = 0, y1 = 0; + gsl_histogram2d_get_yrange (h, i, &y0, &y1); + + if (y0 != i || y1 != i + 1) + { + status = 1; + } + } + gsl_test (status, + "gsl_histogram2d_get_ylowerlimit and yupperlimit"); + } + + + { + int status = 0; + if (gsl_histogram2d_xmax (h) != M) + status = 1; + gsl_test (status, "gsl_histogram2d_xmax"); + } + + { + int status = 0; + if (gsl_histogram2d_xmin (h) != 0) + status = 1; + gsl_test (status, "gsl_histogram2d_xmin"); + } + + { + int status = 0; + if (gsl_histogram2d_nx (h) != M) + status = 1; + gsl_test (status, "gsl_histogram2d_nx"); + } + + { + int status = 0; + if (gsl_histogram2d_ymax (h) != N) + status = 1; + gsl_test (status, "gsl_histogram2d_ymax"); + } + + { + int status = 0; + if (gsl_histogram2d_ymin (h) != 0) + status = 1; + gsl_test (status, "gsl_histogram2d_ymin"); + } + + { + int status = 0; + if (gsl_histogram2d_ny (h) != N) + status = 1; + gsl_test (status, "gsl_histogram2d_ny"); + } + + h->bin[3 * N + 2] = 123456.0; + h->bin[4 * N + 3] = -654321; + + { + double max = gsl_histogram2d_max_val (h); + gsl_test (max != 123456.0, "gsl_histogram2d_max_val finds maximum value"); + } + + { + double min = gsl_histogram2d_min_val (h); + gsl_test (min != -654321.0, + "gsl_histogram2d_min_val finds minimum value"); + } + + { + size_t imax, jmax; + gsl_histogram2d_max_bin (h, &imax, &jmax); + gsl_test (imax != 3 + || jmax != 2, + "gsl_histogram2d_max_bin finds maximum value bin"); + } + + { + size_t imin, jmin; + gsl_histogram2d_min_bin (h, &imin, &jmin); + gsl_test (imin != 4 + || jmin != 3, "gsl_histogram2d_min_bin find minimum value bin"); + } + + for (i = 0; i < M * N; i++) + { + h->bin[i] = i + 27; + g->bin[i] = (i + 27) * (i + 1); + } + + { + double sum = gsl_histogram2d_sum (h); + gsl_test (sum != N * M * 27 + ((N * M - 1) * N * M) / 2, + "gsl_histogram2d_sum sums all bin values"); + } + + { + /* first test... */ + const double xpos = 0.6; + const double ypos = 0.85; + double xmean; + double ymean; + size_t xbin; + size_t ybin; + gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N); + gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1); + gsl_histogram2d_increment (h3, xpos, ypos); + gsl_histogram2d_find (h3, xpos, ypos, &xbin, &ybin); + xmean = gsl_histogram2d_xmean (h3); + ymean = gsl_histogram2d_ymean (h3); + + { + double expected_xmean = (h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0; + double expected_ymean = (h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0; + gsl_test_abs (xmean, expected_xmean, 100.0 * GSL_DBL_EPSILON, + "gsl_histogram2d_xmean"); + gsl_test_abs (ymean, expected_ymean, 100.0 * GSL_DBL_EPSILON, + "gsl_histogram2d_ymean"); + }; + gsl_histogram2d_free (h3); + } + + { + /* test it with bivariate normal distribution */ + const double xmean = 0.7; + const double ymean = 0.7; + const double xsigma = 0.1; + const double ysigma = 0.1; + const double correl = 0.5; + const double norm = + 10.0 / M_PI / xsigma / ysigma / sqrt (1.0 - correl * correl); + size_t xbin; + size_t ybin; + gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N); + gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1); + /* initialize with 2d gauss pdf in two directions */ + for (xbin = 0; xbin < M; xbin++) + { + double xi = + ((h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0 - xmean) / xsigma; + for (ybin = 0; ybin < N; ybin++) + { + double yi = + ((h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0 - + ymean) / ysigma; + double prob = + norm * exp (-(xi * xi - 2.0 * correl * xi * yi + yi * yi) / + 2.0 / (1 - correl * correl)); + h3->bin[xbin * N + ybin] = prob; + } + } + { + double xs = gsl_histogram2d_xsigma (h3); + double ys = gsl_histogram2d_ysigma (h3); + /* evaluate results and compare with parameters */ + + gsl_test_abs (gsl_histogram2d_xmean (h3), xmean, 2.0/M, + "gsl_histogram2d_xmean histogram mean(x)"); + gsl_test_abs (gsl_histogram2d_ymean (h3), ymean, 2.0/N, + "gsl_histogram2d_ymean histogram mean(y)"); + gsl_test_abs (xs, xsigma, 2.0/M, + "gsl_histogram2d_xsigma histogram stdev(x)"); + gsl_test_abs (ys, ysigma, 2.0/N, + "gsl_histogram2d_ysigma histogram stdev(y)"); + gsl_test_abs (gsl_histogram2d_cov (h3) / xs / ys, correl, + 2.0/((M < N) ? M : N), + "gsl_histogram2d_cov histogram covariance"); + } + gsl_histogram2d_free (h3); + } + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_add (h1, h); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != g->bin[i] + h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_add histogram addition"); + } + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_sub (h1, h); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != g->bin[i] - h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_sub histogram subtraction"); + } + + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_mul (h1, h); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != g->bin[i] * h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_mul histogram multiplication"); + } + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_div (h1, h); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != g->bin[i] / h->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_div histogram division"); + } + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_scale (h1, 0.5); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != 0.5 * g->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_scale histogram scaling"); + } + + gsl_histogram2d_memcpy (h1, g); + gsl_histogram2d_shift (h1, 0.25); + + { + int status = 0; + for (i = 0; i < M * N; i++) + { + if (h1->bin[i] != 0.25 + g->bin[i]) + status = 1; + } + gsl_test (status, "gsl_histogram2d_shift histogram shift"); + } + + gsl_histogram2d_free (h); /* free whatever is in h */ + + h = gsl_histogram2d_calloc_uniform (M1, N1, 0.0, 5.0, 0.0, 5.0); + + gsl_test (h->xrange == 0, + "gsl_histogram2d_calloc_uniform returns valid range pointer"); + gsl_test (h->yrange == 0, + "gsl_histogram2d_calloc_uniform returns valid range pointer"); + gsl_test (h->bin == 0, + "gsl_histogram2d_calloc_uniform returns valid bin pointer"); + gsl_test (h->nx != M1, "gsl_histogram2d_calloc_uniform returns valid nx"); + gsl_test (h->ny != N1, "gsl_histogram2d_calloc_uniform returns valid ny"); + + gsl_histogram2d_accumulate (h, 0.0, 3.01, 1.0); + gsl_histogram2d_accumulate (h, 0.1, 2.01, 2.0); + gsl_histogram2d_accumulate (h, 0.2, 1.01, 3.0); + gsl_histogram2d_accumulate (h, 0.3, 0.01, 4.0); + + { + size_t i1, i2, i3, i4; + size_t j1, j2, j3, j4; + double expected; + int status; + status = gsl_histogram2d_find (h, 0.0, 3.01, &i1, &j1); + status = gsl_histogram2d_find (h, 0.1, 2.01, &i2, &j2); + status = gsl_histogram2d_find (h, 0.2, 1.01, &i3, &j3); + status = gsl_histogram2d_find (h, 0.3, 0.01, &i4, &j4); + + for (i = 0; i < M1; i++) + { + for (j = 0; j < N1; j++) + { + if (i == i1 && j == j1) + { + expected = 1.0; + } + else if (i == i2 && j == j2) + { + expected = 2.0; + } + else if (i == i3 && j == j3) + { + expected = 3.0; + } + else if (i == i4 && j == j4) + { + expected = 4.0; + } + else + { + expected = 0.0; + } + + if (h->bin[i * N1 + j] != expected) + { + status = 1; + } + } + } + gsl_test (status, "gsl_histogram2d_find returns index"); + } + + { + FILE *f = fopen ("test.txt", "w"); + gsl_histogram2d_fprintf (f, h, "%.19e", "%.19e"); + fclose (f); + } + + { + FILE *f = fopen ("test.txt", "r"); + gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1); + int status = 0; + + gsl_histogram2d_fscanf (f, hh); + + for (i = 0; i <= M1; i++) + { + if (h->xrange[i] != hh->xrange[i]) + { + printf ("xrange[%d] : %g orig vs %g\n", + (int) i, h->xrange[i], hh->xrange[i]); + status = 1; + } + } + + for (j = 0; j <= N1; j++) + { + if (h->yrange[j] != hh->yrange[j]) + { + printf ("yrange[%d] : %g orig vs %g\n", + (int) j, h->yrange[j], hh->yrange[j]); + status = 1; + } + } + + for (i = 0; i < M1 * N1; i++) + { + if (h->bin[i] != hh->bin[i]) + { + printf ("bin[%d] : %g orig vs %g\n", + (int) i, h->bin[i], hh->bin[i]); + status = 1; + } + } + + gsl_test (status, "gsl_histogram2d_fprintf and fscanf"); + + gsl_histogram2d_free (hh); + fclose (f); + } + + { + FILE *f = fopen ("test.dat", "wb"); + gsl_histogram2d_fwrite (f, h); + fclose (f); + } + + { + FILE *f = fopen ("test.dat", "rb"); + gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1); + int status = 0; + + gsl_histogram2d_fread (f, hh); + + for (i = 0; i <= M1; i++) + { + if (h->xrange[i] != hh->xrange[i]) + { + printf ("xrange[%d] : %g orig vs %g\n", + (int) i, h->xrange[i], hh->xrange[i]); + status = 1; + } + } + + for (j = 0; j <= N1; j++) + { + if (h->yrange[j] != hh->yrange[j]) + { + printf ("yrange[%d] : %g orig vs %g\n", + (int) j, h->yrange[j], hh->yrange[j]); + status = 1; + } + } + + for (i = 0; i < M1 * N1; i++) + { + if (h->bin[i] != hh->bin[i]) + { + printf ("bin[%d] : %g orig vs %g\n", + (int) i, h->bin[i], hh->bin[i]); + status = 1; + } + } + + gsl_test (status, "gsl_histogram2d_fwrite and fread"); + + gsl_histogram2d_free (hh); + fclose (f); + } + + gsl_histogram2d_free (h); + gsl_histogram2d_free (h1); + gsl_histogram2d_free (g); + gsl_histogram2d_free (hr); +} diff --git a/gsl-1.9/histogram/test2d_resample.c b/gsl-1.9/histogram/test2d_resample.c new file mode 100644 index 0000000..4df207f --- /dev/null +++ b/gsl-1.9/histogram/test2d_resample.c @@ -0,0 +1,115 @@ +/* histogram/test2d_resample.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <gsl/gsl_histogram2d.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_ieee_utils.h> + +#include "urand.c" + +void +test2d_resample (void) +{ + size_t i, j; + int status = 0; + double total = 0; + size_t N = 200000; + + gsl_histogram2d *h; + + gsl_ieee_env_setup (); + + h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 0.0, 1.0); + + for (i = 0; i < 10; i++) + { + for (j = 0; j < 10; j++) + { + double w = 10.0 * i + j; + total += w; + gsl_histogram2d_accumulate (h, 0.1 * i, 0.1 * i, w); + } + } + + { + gsl_histogram2d_pdf *p = gsl_histogram2d_pdf_alloc (10,10); + + gsl_histogram2d *hh = gsl_histogram2d_calloc_uniform (20, 20, + 0.0, 1.0, + 0.0, 1.0); + + gsl_histogram2d_pdf_init (p, h); + + for (i = 0; i < N; i++) + { + double u = urand(); + double v = urand(); + double x, y; + status = gsl_histogram2d_pdf_sample (p, u, v, &x, &y); + status = gsl_histogram2d_increment (hh, x, y); + } + + status = 0; + for (i = 0; i < 20; i++) + { + for (j = 0; j < 20; j++) + { + double z = 4 * total * gsl_histogram2d_get (hh, i, j) / (double) N; + size_t k1, k2; + double ya; + double x, xmax, y, ymax; + + gsl_histogram2d_get_xrange (hh, i, &x, &xmax); + gsl_histogram2d_get_yrange (hh, j, &y, &ymax); + + gsl_histogram2d_find (h, x, y, &k1, &k2); + ya = gsl_histogram2d_get (h, k1, k2); + + if (ya == 0) + { + if (z != 0) + { + status = 1; + printf ("(%d,%d): %g vs %g\n", (int)i, (int)j, z, ya); + } + } + else + { + double err = 1 / sqrt (gsl_histogram2d_get (hh, i, j)); + double sigma = fabs ((z - ya) / (ya * err)); + if (sigma > 3) + { + status = 1; + printf ("%g vs %g err=%g sigma=%g\n", z, ya, err, sigma); + } + } + } + } + + gsl_histogram2d_pdf_free (p) ; + gsl_histogram2d_free (hh) ; + + gsl_test (status, "gsl_histogram2d_pdf_sample within statistical errors"); + } + + gsl_histogram2d_free (h) ; +} diff --git a/gsl-1.9/histogram/test2d_trap.c b/gsl-1.9/histogram/test2d_trap.c new file mode 100644 index 0000000..b597ae5 --- /dev/null +++ b/gsl-1.9/histogram/test2d_trap.c @@ -0,0 +1,203 @@ +/* histogram/test2d_trap.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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 <stdio.h> +#include <gsl/gsl_histogram2d.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_ieee_utils.h> + +#define N 107 +#define M 239 + +static void my_error_handler (const char *reason, const char *file, + int line, int err); +static int status = 0; + +void +test2d_trap (void) +{ + gsl_histogram2d *h; + double result, lower, upper; + size_t i, j; + + gsl_set_error_handler (&my_error_handler); + + gsl_ieee_env_setup (); + + status = 0; + h = gsl_histogram2d_calloc (0, 10); + gsl_test (!status, "gsl_histogram_calloc traps zero-width histogram"); + gsl_test (h != 0, + "gsl_histogram2d_calloc returns NULL for zero-width histogram"); + + + status = 0; + h = gsl_histogram2d_calloc (10, 0); + gsl_test (!status, "gsl_histogram_calloc traps zero-length histogram"); + gsl_test (h != 0, + "gsl_histogram2d_calloc returns NULL for zero-length histogram"); + + + status = 0; + h = gsl_histogram2d_calloc_uniform (0, 10, 0.0, 1.0, 0.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps zero-width histogram"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for zero-width histogram"); + + status = 0; + h = gsl_histogram2d_calloc_uniform (10, 0, 0.0, 1.0, 0.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps zero-length histogram"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for zero-length histogram"); + + status = 0; + h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 1.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps equal endpoints"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for equal endpoints"); + + status = 0; + h = gsl_histogram2d_calloc_uniform (10, 10, 1.0, 1.0, 0.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps equal endpoints"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for equal endpoints"); + + status = 0; + h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 2.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps invalid range"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for invalid range"); + + status = 0; + h = gsl_histogram2d_calloc_uniform (10, 10, 2.0, 1.0, 0.0, 1.0); + gsl_test (!status, + "gsl_histogram2d_calloc_uniform traps invalid range"); + gsl_test (h != 0, + "gsl_histogram2d_calloc_uniform returns NULL for invalid range"); + + h = gsl_histogram2d_calloc_uniform (N, M, 0.0, 1.0, 0.0, 1.0); + + status = gsl_histogram2d_accumulate (h, 1.0, 0.0, 10.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x at xmax"); + + status = gsl_histogram2d_accumulate (h, 2.0, 0.0, 100.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x above xmax"); + + status = gsl_histogram2d_accumulate (h, -1.0, 0.0, 1000.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps x below xmin"); + + status = gsl_histogram2d_accumulate (h, 0.0, 1.0, 10.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y at ymax"); + + status = gsl_histogram2d_accumulate (h, 0.0, 2.0, 100.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y above ymax"); + + status = gsl_histogram2d_accumulate (h, 0.0, -1.0, 1000.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_accumulate traps y below ymin"); + + + status = gsl_histogram2d_increment (h, 1.0, 0.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x at xmax"); + + status = gsl_histogram2d_increment (h, 2.0, 0.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x above xmax"); + + status = gsl_histogram2d_increment (h, -1.0, 0.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps x below xmin"); + + + status = gsl_histogram2d_increment (h, 0.0, 1.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y at ymax"); + + status = gsl_histogram2d_increment (h, 0.0, 2.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y above ymax"); + + status = gsl_histogram2d_increment (h, 0.0, -1.0); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_increment traps y below ymin"); + + + result = gsl_histogram2d_get (h, N, 0); + gsl_test (result != 0, "gsl_histogram2d_get traps x index at nx"); + + result = gsl_histogram2d_get (h, N + 1, 0); + gsl_test (result != 0, "gsl_histogram2d_get traps x index above nx"); + + result = gsl_histogram2d_get (h, 0, M); + gsl_test (result != 0, "gsl_histogram2d_get traps y index at ny"); + + result = gsl_histogram2d_get (h, 0, M + 1); + gsl_test (result != 0, "gsl_histogram2d_get traps y index above ny"); + + + status = gsl_histogram2d_get_xrange (h, N, &lower, &upper); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_xrange traps index at nx"); + + status = gsl_histogram2d_get_xrange (h, N + 1, &lower, &upper); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_xrange traps index above nx"); + + status = gsl_histogram2d_get_yrange (h, M, &lower, &upper); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_yrange traps index at ny"); + + status = gsl_histogram2d_get_yrange (h, M + 1, &lower, &upper); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_get_yrange traps index above ny"); + + status = 0; + gsl_histogram2d_find (h, -0.01, 0.0, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x below xmin"); + + status = 0; + gsl_histogram2d_find (h, 1.0, 0.0, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x at xmax"); + + status = 0; + gsl_histogram2d_find (h, 1.1, 0.0, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps x above xmax"); + + + status = 0; + gsl_histogram2d_find (h, 0.0, -0.01, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y below ymin"); + + status = 0; + gsl_histogram2d_find (h, 0.0, 1.0, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y at ymax"); + + status = 0; + gsl_histogram2d_find (h, 0.0, 1.1, &i, &j); + gsl_test (status != GSL_EDOM, "gsl_histogram2d_find traps y above ymax"); + + gsl_histogram2d_free (h); +} + + +static 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); + status = 1; +} diff --git a/gsl-1.9/histogram/urand.c b/gsl-1.9/histogram/urand.c new file mode 100644 index 0000000..9c6450f --- /dev/null +++ b/gsl-1.9/histogram/urand.c @@ -0,0 +1,26 @@ +/* histogram/urand.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * 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. + */ + +static double urand (void); + +static double urand (void) { + static unsigned long int x = 1; + x = (1103515245 * x + 12345) & 0x7fffffffUL ; + return x / 2147483648.0 ; +} |