diff options
Diffstat (limited to 'gsl-1.9/multifit')
30 files changed, 6025 insertions, 0 deletions
diff --git a/gsl-1.9/multifit/ChangeLog b/gsl-1.9/multifit/ChangeLog new file mode 100644 index 0000000..17e4184 --- /dev/null +++ b/gsl-1.9/multifit/ChangeLog @@ -0,0 +1,145 @@ +2007-01-26 Brian Gough <bjg@network-theory.co.uk> + + * fsolver.c (gsl_multifit_fsolver_set): made vector argument x + const + +2006-03-30 Brian Gough <bjg@network-theory.co.uk> + + * fsolver.c (gsl_multifit_fsolver_alloc): minpack algorithms + require n>=p, added check + + * fdfsolver.c (gsl_multifit_fdfsolver_alloc): minpack algorithms + require n>=p, added check + +2006-02-20 Brian Gough <bjg@network-theory.co.uk> + + * multilinear.c (gsl_multifit_linear_est): added linear estimator + + * test_estimator.c (test_estimator): added test for gsl_multifit_linear_est + +2005-07-03 Brian Gough <bjg@network-theory.co.uk> + + * multilinear.c (gsl_multifit_linear_svd): accept a user-specified + tolerance for the SVD cutoff and return effective rank + (gsl_multifit_wlinear_svd): same + +2004-12-23 Brian Gough <bjg@network-theory.co.uk> + + * gsl_multifit_nlin.h: removed unused declaration of + gsl_multifit_fdjacobian + +2004-06-14 Brian Gough <bjg@network-theory.co.uk> + + * lmiterate.c (iterate): handle case where fnorm = 0 to avoid + division by zero + + * covar.c (gsl_multifit_covar): change tolerance test to match + original code, and handle case where tolr = 0 + +2003-03-21 Brian Gough <bjg@network-theory.co.uk> + + * lmset.c (set): removed reference to q, compute QR decomposition + in place + + * lmiterate.c (iterate): removed reference to q, compute QR + decomposition in-place for R + + * lmutil.c: removed compute_qtf + + * lmder.c (lmder_free): removed reference to q + (lmder_alloc): removed reference to q + +Tue Nov 12 22:18:14 2002 Brian Gough <bjg@network-theory.co.uk> + + * lmder.c (lmder_alloc): use GSL_ERROR instead of GSL_ERROR_VAL + for internal alloc functions which return int + +Thu Feb 28 20:15:33 2002 Brian Gough <bjg@network-theory.co.uk> + + * lmiterate.c (iterate): return immediately if evaluation raised + error (Hans E. Plesser) + + * lmpar.c (lmpar): avoid division by zero for w=0 in rank + deficient case + +Mon Oct 8 19:25:55 2001 Brian Gough <bjg@network-theory.co.uk> + + * test.c (main): added extra nist tests + + * lmutil.c (compute_rptdx): fixed bug, permutation in rptdx vector + was incorrectly applied + + * lmpar.c (compute_newton_direction): fixed bug, permutation of + newton direction vector was incorrect (should have been inverse + permutation) + +Mon Jul 30 17:43:21 2001 Brian Gough <bjg@network-theory.co.uk> + + * test_filip.c (test_filip): reduce tolerance on covariance test + slightly for MSVC with /O2 + +Sun Jul 1 22:42:34 2001 Brian Gough <bjg@network-theory.co.uk> + + * multilinear.c: modified to use new-style vector views + + * test_pontius.c: modified to use new-style vector views + + * test_longley.c: modified to use new-style vector views + + * test_fn.c: modified to use new-style vector views + + * test_filip.c: modified to use new-style vector views + +Tue Jun 26 21:41:30 2001 Brian Gough <bjg@network-theory.co.uk> + + * test_filip.c (test_filip): reduce tolerance on covariance test + slightly + +Wed Jun 20 13:11:26 2001 Brian Gough <bjg@network-theory.co.uk> + + * removed unused variable work2 + +Tue Jun 19 23:18:01 2001 Brian Gough <bjg@network-theory.co.uk> + + * multilinear.c: perform column scaling before doing fit to + improve accuracy + (gsl_multifit_linear): use modified Golub-Reinsch SVD for greater + speed + (gsl_multifit_wlinear): use modified Golub-Reinsch SVD for greater + speed + +Wed Jun 6 13:32:22 2001 Brian Gough <bjg@network-theory.co.uk> + + * lmder.c covar.c lmiterate.c lmset.c: updated to use new QR + calling convention (now passes workspace) + +Sat Apr 28 11:46:59 2001 Brian Gough <bjg@network-theory.co.uk> + + * qrsolv.c (qrsolv): removed local declaration of j to avoid + shadowing global j + +Mon Apr 23 13:40:04 2001 Brian Gough <bjg@network-theory.co.uk> + + * qrsolv.c (qrsolv): made function static so it is not exported + +Wed Apr 18 13:39:33 2001 Brian Gough <bjg@network-theory.co.uk> + + * lmpar.c (compute_newton_direction): moved final rescaling inside + loop, as in the original lmpar.f + +Thu Mar 8 15:29:32 2001 Brian Gough <bjg@network-theory.co.uk> + + * lmpar.c (compute_newton_direction): corrected bug that produced + negative index + +Sun Feb 18 16:39:46 2001 Brian Gough <bjg@network-theory.co.uk> + + * fdfsolver.c (gsl_multifit_fdfsolver_alloc): changed so that the + solver _alloc function no longer calls _set, the user must do that + separately. + +Fri Sep 29 19:19:24 2000 Brian Gough <bjg@network-theory.co.uk> + + * Makefile.am multifit/demo.c: removed demo from Makefile since it + is was just a temporary test. + diff --git a/gsl-1.9/multifit/Makefile.am b/gsl-1.9/multifit/Makefile.am new file mode 100644 index 0000000..ab29801 --- /dev/null +++ b/gsl-1.9/multifit/Makefile.am @@ -0,0 +1,20 @@ +noinst_LTLIBRARIES = libgslmultifit.la + +pkginclude_HEADERS = gsl_multifit.h gsl_multifit_nlin.h + +INCLUDES= -I$(top_builddir) + +libgslmultifit_la_SOURCES = multilinear.c work.c lmder.c fsolver.c fdfsolver.c convergence.c gradient.c covar.c + +noinst_HEADERS = lmutil.c lmpar.c lmset.c lmiterate.c qrsolv.c test_brown.c test_enso.c test_filip.c test_fn.c test_hahn1.c test_kirby2.c test_longley.c test_nelson.c test_pontius.c test_estimator.c + +check_PROGRAMS = test #demo + +TESTS = $(check_PROGRAMS) + +test_SOURCES = test.c +test_LDADD = libgslmultifit.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la ../sys/libgslsys.la + +#demo_SOURCES = demo.c +#demo_LDADD = libgslmultifit.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../randist/libgslrandist.la ../rng/libgslrng.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la ../sys/libgslsys.la + diff --git a/gsl-1.9/multifit/Makefile.in b/gsl-1.9/multifit/Makefile.in new file mode 100644 index 0000000..5004e67 --- /dev/null +++ b/gsl-1.9/multifit/Makefile.in @@ -0,0 +1,552 @@ +# 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 = multifit +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) +libgslmultifit_la_LIBADD = +am_libgslmultifit_la_OBJECTS = multilinear.lo work.lo lmder.lo \ + fsolver.lo fdfsolver.lo convergence.lo gradient.lo covar.lo +libgslmultifit_la_OBJECTS = $(am_libgslmultifit_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslmultifit.la ../linalg/libgsllinalg.la \ + ../permutation/libgslpermutation.la ../blas/libgslblas.la \ + ../cblas/libgslcblas.la ../matrix/libgslmatrix.la \ + ../vector/libgslvector.la ../block/libgslblock.la \ + ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la \ + ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.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 = $(libgslmultifit_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslmultifit_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 = libgslmultifit.la +pkginclude_HEADERS = gsl_multifit.h gsl_multifit_nlin.h +INCLUDES = -I$(top_builddir) +libgslmultifit_la_SOURCES = multilinear.c work.c lmder.c fsolver.c fdfsolver.c convergence.c gradient.c covar.c +noinst_HEADERS = lmutil.c lmpar.c lmset.c lmiterate.c qrsolv.c test_brown.c test_enso.c test_filip.c test_fn.c test_hahn1.c test_kirby2.c test_longley.c test_nelson.c test_pontius.c test_estimator.c +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c +test_LDADD = libgslmultifit.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la ../sys/libgslsys.la +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 multifit/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps multifit/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 +libgslmultifit.la: $(libgslmultifit_la_OBJECTS) $(libgslmultifit_la_DEPENDENCIES) + $(LINK) $(libgslmultifit_la_LDFLAGS) $(libgslmultifit_la_OBJECTS) $(libgslmultifit_la_LIBADD) $(LIBS) + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; for p in $$list; do \ + f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f $$p $$f"; \ + rm -f $$p $$f ; \ + done +test$(EXEEXT): $(test_OBJECTS) $(test_DEPENDENCIES) + @rm -f test$(EXEEXT) + $(LINK) $(test_LDFLAGS) $(test_OBJECTS) $(test_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +distclean-libtool: + -rm -f libtool +uninstall-info-am: +install-pkgincludeHEADERS: $(pkginclude_HEADERS) + @$(NORMAL_INSTALL) + test -z "$(pkgincludedir)" || $(mkdir_p) "$(DESTDIR)$(pkgincludedir)" + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \ + f=$(am__strip_dir) \ + echo " $(pkgincludeHEADERS_INSTALL) '$$d$$p' '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + $(pkgincludeHEADERS_INSTALL) "$$d$$p" "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +uninstall-pkgincludeHEADERS: + @$(NORMAL_UNINSTALL) + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + f=$(am__strip_dir) \ + echo " rm -f '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + rm -f "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +check-TESTS: $(TESTS) + @failed=0; all=0; xfail=0; xpass=0; skip=0; \ + srcdir=$(srcdir); export srcdir; \ + list='$(TESTS)'; \ + if test -n "$$list"; then \ + for tst in $$list; do \ + if test -f ./$$tst; then dir=./; \ + elif test -f $$tst; then dir=; \ + else dir="$(srcdir)/"; fi; \ + if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xpass=`expr $$xpass + 1`; \ + failed=`expr $$failed + 1`; \ + echo "XPASS: $$tst"; \ + ;; \ + *) \ + echo "PASS: $$tst"; \ + ;; \ + esac; \ + elif test $$? -ne 77; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xfail=`expr $$xfail + 1`; \ + echo "XFAIL: $$tst"; \ + ;; \ + *) \ + failed=`expr $$failed + 1`; \ + echo "FAIL: $$tst"; \ + ;; \ + esac; \ + else \ + skip=`expr $$skip + 1`; \ + echo "SKIP: $$tst"; \ + fi; \ + done; \ + if test "$$failed" -eq 0; then \ + if test "$$xfail" -eq 0; then \ + banner="All $$all tests passed"; \ + else \ + banner="All $$all tests behaved as expected ($$xfail expected failures)"; \ + fi; \ + else \ + if test "$$xpass" -eq 0; then \ + banner="$$failed of $$all tests failed"; \ + else \ + banner="$$failed of $$all tests did not behave as expected ($$xpass unexpected passes)"; \ + fi; \ + fi; \ + dashes="$$banner"; \ + skipped=""; \ + if test "$$skip" -ne 0; then \ + skipped="($$skip tests were not run)"; \ + test `echo "$$skipped" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$skipped"; \ + fi; \ + report=""; \ + if test "$$failed" -ne 0 && test -n "$(PACKAGE_BUGREPORT)"; then \ + report="Please report to $(PACKAGE_BUGREPORT)"; \ + test `echo "$$report" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$report"; \ + fi; \ + dashes=`echo "$$dashes" | sed s/./=/g`; \ + echo "$$dashes"; \ + echo "$$banner"; \ + test -z "$$skipped" || echo "$$skipped"; \ + test -z "$$report" || echo "$$report"; \ + echo "$$dashes"; \ + test "$$failed" -eq 0; \ + else :; fi + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-TESTS +check: check-am +all-am: Makefile $(LTLIBRARIES) $(HEADERS) +installdirs: + for dir in "$(DESTDIR)$(pkgincludedir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-libtool distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: install-pkgincludeHEADERS + +install-exec-am: + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-info-am uninstall-pkgincludeHEADERS + +.PHONY: CTAGS GTAGS all all-am check check-TESTS check-am clean \ + clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES ctags distclean distclean-compile \ + distclean-generic distclean-libtool distclean-tags distdir dvi \ + dvi-am html html-am info info-am install install-am \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man \ + install-pkgincludeHEADERS install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-info-am \ + uninstall-pkgincludeHEADERS + + +#demo_SOURCES = demo.c +#demo_LDADD = libgslmultifit.la ../linalg/libgsllinalg.la ../permutation/libgslpermutation.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../matrix/libgslmatrix.la ../vector/libgslvector.la ../block/libgslblock.la ../randist/libgslrandist.la ../rng/libgslrng.la ../complex/libgslcomplex.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../utils/libutils.la ../sys/libgslsys.la +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gsl-1.9/multifit/TODO b/gsl-1.9/multifit/TODO new file mode 100644 index 0000000..b67e058 --- /dev/null +++ b/gsl-1.9/multifit/TODO @@ -0,0 +1,7 @@ +The following would also be nice additions to the multifit function suite +(see MS Excel regression output for example): + +1. Produce the correlation coefficient (r) and other statistics. +2. Allow fit variable weighting (not observation weighting). +3. Allow for principal factor computations. + diff --git a/gsl-1.9/multifit/convergence.c b/gsl-1.9/multifit/convergence.c new file mode 100644 index 0000000..73a89a3 --- /dev/null +++ b/gsl-1.9/multifit/convergence.c @@ -0,0 +1,89 @@ +/* multifit/convergence.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_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multifit_nlin.h> + +int +gsl_multifit_test_delta (const gsl_vector * dx, const gsl_vector * x, + double epsabs, double epsrel) +{ + size_t i; + int ok = 1; + const size_t n = x->size ; + + if (epsrel < 0.0) + { + GSL_ERROR ("relative tolerance is negative", GSL_EBADTOL); + } + + for (i = 0 ; i < n ; i++) + { + double xi = gsl_vector_get(x,i); + double dxi = gsl_vector_get(dx,i); + double tolerance = epsabs + epsrel * fabs(xi) ; + + if (fabs(dxi) < tolerance) + { + ok = 1; + } + else + { + ok = 0; + break; + } + } + + if (ok) + return GSL_SUCCESS ; + + return GSL_CONTINUE; +} + +int +gsl_multifit_test_gradient (const gsl_vector * g, double epsabs) +{ + size_t i; + + double residual = 0; + + const size_t n = g->size; + + if (epsabs < 0.0) + { + GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL); + } + + for (i = 0 ; i < n ; i++) + { + double gi = gsl_vector_get(g, i); + + residual += fabs(gi); + } + + + if (residual < epsabs) + { + return GSL_SUCCESS; + } + + return GSL_CONTINUE ; +} diff --git a/gsl-1.9/multifit/covar.c b/gsl-1.9/multifit/covar.c new file mode 100644 index 0000000..05683a7 --- /dev/null +++ b/gsl-1.9/multifit/covar.c @@ -0,0 +1,192 @@ +/* multifit/covar.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 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_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_multifit_nlin.h> + +/* Compute the covariance matrix + + cov = inv (J^T J) + + by QRP^T decomposition of J +*/ + +int +gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar) +{ + double tolr; + + size_t i, j, k; + size_t kmax = 0; + + gsl_matrix * r; + gsl_vector * tau; + gsl_vector * norm; + gsl_permutation * perm; + + size_t m = J->size1, n = J->size2 ; + + if (m < n) + { + GSL_ERROR ("Jacobian be rectangular M x N with M >= N", GSL_EBADLEN); + } + + if (covar->size1 != covar->size2 || covar->size1 != n) + { + GSL_ERROR ("covariance matrix must be square and match second dimension of jacobian", GSL_EBADLEN); + } + + r = gsl_matrix_alloc (m, n); + tau = gsl_vector_alloc (n); + perm = gsl_permutation_alloc (n) ; + norm = gsl_vector_alloc (n) ; + + { + int signum = 0; + gsl_matrix_memcpy (r, J); + gsl_linalg_QRPT_decomp (r, tau, perm, &signum, norm); + } + + + /* Form the inverse of R in the full upper triangle of R */ + + tolr = epsrel * fabs(gsl_matrix_get(r, 0, 0)); + + for (k = 0 ; k < n ; k++) + { + double rkk = gsl_matrix_get(r, k, k); + + if (fabs(rkk) <= tolr) + { + break; + } + + gsl_matrix_set(r, k, k, 1.0/rkk); + + for (j = 0; j < k ; j++) + { + double t = gsl_matrix_get(r, j, k) / rkk; + gsl_matrix_set (r, j, k, 0.0); + + for (i = 0; i <= j; i++) + { + double rik = gsl_matrix_get (r, i, k); + double rij = gsl_matrix_get (r, i, j); + + gsl_matrix_set (r, i, k, rik - t * rij); + } + } + kmax = k; + } + + /* Form the full upper triangle of the inverse of R^T R in the full + upper triangle of R */ + + for (k = 0; k <= kmax ; k++) + { + for (j = 0; j < k; j++) + { + double rjk = gsl_matrix_get (r, j, k); + + for (i = 0; i <= j ; i++) + { + double rij = gsl_matrix_get (r, i, j); + double rik = gsl_matrix_get (r, i, k); + + gsl_matrix_set (r, i, j, rij + rjk * rik); + } + } + + { + double t = gsl_matrix_get (r, k, k); + + for (i = 0; i <= k; i++) + { + double rik = gsl_matrix_get (r, i, k); + + gsl_matrix_set (r, i, k, t * rik); + }; + } + } + + /* Form the full lower triangle of the covariance matrix in the + strict lower triangle of R and in w */ + + for (j = 0 ; j < n ; j++) + { + size_t pj = gsl_permutation_get (perm, j); + + for (i = 0; i <= j; i++) + { + size_t pi = gsl_permutation_get (perm, i); + + double rij; + + if (j > kmax) + { + gsl_matrix_set (r, i, j, 0.0); + rij = 0.0 ; + } + else + { + rij = gsl_matrix_get (r, i, j); + } + + if (pi > pj) + { + gsl_matrix_set (r, pi, pj, rij); + } + else if (pi < pj) + { + gsl_matrix_set (r, pj, pi, rij); + } + + } + + { + double rjj = gsl_matrix_get (r, j, j); + gsl_matrix_set (covar, pj, pj, rjj); + } + } + + + /* symmetrize the covariance matrix */ + + for (j = 0 ; j < n ; j++) + { + for (i = 0; i < j ; i++) + { + double rji = gsl_matrix_get (r, j, i); + + gsl_matrix_set (covar, j, i, rji); + gsl_matrix_set (covar, i, j, rji); + } + } + + gsl_matrix_free (r); + gsl_permutation_free (perm); + gsl_vector_free (tau); + gsl_vector_free (norm); + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multifit/fdfsolver.c b/gsl-1.9/multifit/fdfsolver.c new file mode 100644 index 0000000..4f46a9a --- /dev/null +++ b/gsl-1.9/multifit/fdfsolver.c @@ -0,0 +1,170 @@ +/* multifit/fdfsolver.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 <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multifit_nlin.h> + +gsl_multifit_fdfsolver * +gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T, + size_t n, size_t p) +{ + int status; + + gsl_multifit_fdfsolver * s; + + if (n < p) + { + GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0); + } + + s = (gsl_multifit_fdfsolver *) malloc (sizeof (gsl_multifit_fdfsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for multifit solver struct", + GSL_ENOMEM, 0); + } + + s->x = gsl_vector_calloc (p); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->f = gsl_vector_calloc (n); + + if (s->f == 0) + { + gsl_vector_free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0); + } + + s->J = gsl_matrix_calloc (n,p); + + if (s->J == 0) + { + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0); + } + + s->dx = gsl_vector_calloc (p); + + if (s->dx == 0) + { + gsl_matrix_free (s->J); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0); + } + + s->state = malloc (T->size); + + if (s->state == 0) + { + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for multifit solver state", + GSL_ENOMEM, 0); + } + + s->type = T ; + + status = (s->type->alloc)(s->state, n, p); + + if (status != GSL_SUCCESS) + { + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to set solver", status, 0); + } + + s->fdf = NULL; + + return s; +} + +int +gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s, + gsl_multifit_function_fdf * f, + const gsl_vector * x) +{ + if (s->f->size != f->n) + { + GSL_ERROR ("function size does not match solver", GSL_EBADLEN); + } + + if (s->x->size != x->size) + { + GSL_ERROR ("vector length does not match solver", GSL_EBADLEN); + } + + s->fdf = f; + gsl_vector_memcpy(s->x,x); + + return (s->type->set) (s->state, s->fdf, s->x, s->f, s->J, s->dx); +} + +int +gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s) +{ + return (s->type->iterate) (s->state, s->fdf, s->x, s->f, s->J, s->dx); +} + +void +gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s) +{ + (s->type->free) (s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + gsl_matrix_free (s->J); + free (s); +} + +const char * +gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s) +{ + return s->type->name; +} + +gsl_vector * +gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s) +{ + return s->x; +} + diff --git a/gsl-1.9/multifit/fsolver.c b/gsl-1.9/multifit/fsolver.c new file mode 100644 index 0000000..675f75b --- /dev/null +++ b/gsl-1.9/multifit/fsolver.c @@ -0,0 +1,156 @@ +/* multifit/fsolver.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 <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multifit_nlin.h> + +gsl_multifit_fsolver * +gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T, + size_t n, size_t p) +{ + int status; + + gsl_multifit_fsolver * s; + + if (n < p) + { + GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0); + } + + s = (gsl_multifit_fsolver *) malloc (sizeof (gsl_multifit_fsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for multifit solver struct", + GSL_ENOMEM, 0); + } + + s->x = gsl_vector_calloc (p); + + if (s->x == 0) + { + free (s); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + s->f = gsl_vector_calloc (n); + + if (s->f == 0) + { + gsl_vector_free (s->x); + free (s); + GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0); + } + + s->dx = gsl_vector_calloc (p); + + if (s->dx == 0) + { + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); + GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0); + } + + s->state = malloc (T->size); + + if (s->state == 0) + { + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for multifit solver state", + GSL_ENOMEM, 0); + } + + s->type = T ; + + status = (s->type->alloc)(s->state, n, p); + + if (status != GSL_SUCCESS) + { + (s->type->free)(s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to set solver", status, 0); + } + + s->function = NULL; + + return s; +} + +int +gsl_multifit_fsolver_set (gsl_multifit_fsolver * s, + gsl_multifit_function * f, + const gsl_vector * x) +{ + if (s->f->size != f->n) + { + GSL_ERROR ("function size does not match solver", GSL_EBADLEN); + } + + if (s->x->size != x->size) + { + GSL_ERROR ("vector length does not match solver", GSL_EBADLEN); + } + + s->function = f; + gsl_vector_memcpy(s->x,x); + + return (s->type->set) (s->state, s->function, s->x, s->f, s->dx); +} + +int +gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s) +{ + return (s->type->iterate) (s->state, s->function, s->x, s->f, s->dx); +} + +void +gsl_multifit_fsolver_free (gsl_multifit_fsolver * s) +{ + (s->type->free) (s->state); + free (s->state); + gsl_vector_free (s->dx); + gsl_vector_free (s->x); + gsl_vector_free (s->f); + free (s); +} + +const char * +gsl_multifit_fsolver_name (const gsl_multifit_fsolver * s) +{ + return s->type->name; +} + +gsl_vector * +gsl_multifit_fsolver_position (const gsl_multifit_fsolver * s) +{ + return s->x; +} diff --git a/gsl-1.9/multifit/gradient.c b/gsl-1.9/multifit/gradient.c new file mode 100644 index 0000000..fd31f95 --- /dev/null +++ b/gsl-1.9/multifit/gradient.c @@ -0,0 +1,33 @@ +/* multifit/covar.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_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multifit_nlin.h> +#include <gsl/gsl_blas.h> + +int +gsl_multifit_gradient (const gsl_matrix * J, const gsl_vector * f, + gsl_vector * g) +{ + int status = gsl_blas_dgemv (CblasTrans, 1.0, J, f, 0.0, g); + return status; +} + diff --git a/gsl-1.9/multifit/gsl_multifit.h b/gsl-1.9/multifit/gsl_multifit.h new file mode 100644 index 0000000..e408b0d --- /dev/null +++ b/gsl-1.9/multifit/gsl_multifit.h @@ -0,0 +1,106 @@ +/* multifit/gsl_multifit.h + * + * Copyright (C) 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_MULTIFIT_H__ +#define __GSL_MULTIFIT_H__ + +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.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; /* number of observations */ + size_t p; /* number of parameters */ + gsl_matrix * A; + gsl_matrix * Q; + gsl_matrix * QSI; + gsl_vector * S; + gsl_vector * t; + gsl_vector * xt; + gsl_vector * D; +} +gsl_multifit_linear_workspace; + +gsl_multifit_linear_workspace * +gsl_multifit_linear_alloc (size_t n, size_t p); + +void +gsl_multifit_linear_free (gsl_multifit_linear_workspace * work); + +int +gsl_multifit_linear (const gsl_matrix * X, + const gsl_vector * y, + gsl_vector * c, + gsl_matrix * cov, + double * chisq, + gsl_multifit_linear_workspace * work); + +int +gsl_multifit_linear_svd (const gsl_matrix * X, + const gsl_vector * y, + double tol, + size_t * rank, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, + gsl_multifit_linear_workspace * work); + +int +gsl_multifit_wlinear (const gsl_matrix * X, + const gsl_vector * w, + const gsl_vector * y, + gsl_vector * c, + gsl_matrix * cov, + double * chisq, + gsl_multifit_linear_workspace * work); + +int +gsl_multifit_wlinear_svd (const gsl_matrix * X, + const gsl_vector * w, + const gsl_vector * y, + double tol, + size_t * rank, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, + gsl_multifit_linear_workspace * work); + +int +gsl_multifit_linear_est (const gsl_vector * x, + const gsl_vector * c, + const gsl_matrix * cov, double *y, double *y_err); + + +__END_DECLS + +#endif /* __GSL_MULTIFIT_H__ */ diff --git a/gsl-1.9/multifit/gsl_multifit_nlin.h b/gsl-1.9/multifit/gsl_multifit_nlin.h new file mode 100644 index 0000000..6d37960 --- /dev/null +++ b/gsl-1.9/multifit/gsl_multifit_nlin.h @@ -0,0 +1,172 @@ +/* multifit_nlin/gsl_multifit_nlin.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_MULTIFIT_NLIN_H__ +#define __GSL_MULTIFIT_NLIN_H__ + +#include <stdlib.h> +#include <gsl/gsl_types.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.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 + +int gsl_multifit_gradient (const gsl_matrix * J, const gsl_vector * f, + gsl_vector * g); + +int gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar); + + +/* Definition of vector-valued functions with parameters based on gsl_vector */ + +struct gsl_multifit_function_struct +{ + int (* f) (const gsl_vector * x, void * params, gsl_vector * f); + size_t n; /* number of functions */ + size_t p; /* number of independent variables */ + void * params; +}; + +typedef struct gsl_multifit_function_struct gsl_multifit_function ; + +#define GSL_MULTIFIT_FN_EVAL(F,x,y) (*((F)->f))(x,(F)->params,(y)) + +typedef struct + { + const char *name; + size_t size; + int (*alloc) (void *state, size_t n, size_t p); + int (*set) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); + int (*iterate) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx); + void (*free) (void *state); + } +gsl_multifit_fsolver_type; + +typedef struct + { + const gsl_multifit_fsolver_type * type; + gsl_multifit_function * function ; + gsl_vector * x ; + gsl_vector * f ; + gsl_vector * dx ; + void *state; + } +gsl_multifit_fsolver; + +gsl_multifit_fsolver * +gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T, + size_t n, size_t p); + +void gsl_multifit_fsolver_free (gsl_multifit_fsolver * s); + +int gsl_multifit_fsolver_set (gsl_multifit_fsolver * s, + gsl_multifit_function * f, + const gsl_vector * x); + +int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s); + +const char * gsl_multifit_fsolver_name (const gsl_multifit_fsolver * s); +gsl_vector * gsl_multifit_fsolver_position (const gsl_multifit_fsolver * s); + +/* Definition of vector-valued functions and gradient with parameters + based on gsl_vector */ + +struct gsl_multifit_function_fdf_struct +{ + int (* f) (const gsl_vector * x, void * params, gsl_vector * f); + int (* df) (const gsl_vector * x, void * params, gsl_matrix * df); + int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix *df); + size_t n; /* number of functions */ + size_t p; /* number of independent variables */ + void * params; +}; + +typedef struct gsl_multifit_function_fdf_struct gsl_multifit_function_fdf ; + +#define GSL_MULTIFIT_FN_EVAL_F(F,x,y) ((*((F)->f))(x,(F)->params,(y))) +#define GSL_MULTIFIT_FN_EVAL_DF(F,x,dy) ((*((F)->df))(x,(F)->params,(dy))) +#define GSL_MULTIFIT_FN_EVAL_F_DF(F,x,y,dy) ((*((F)->fdf))(x,(F)->params,(y),(dy))) + +typedef struct + { + const char *name; + size_t size; + int (*alloc) (void *state, size_t n, size_t p); + int (*set) (void *state, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); + int (*iterate) (void *state, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); + void (*free) (void *state); + } +gsl_multifit_fdfsolver_type; + +typedef struct + { + const gsl_multifit_fdfsolver_type * type; + gsl_multifit_function_fdf * fdf ; + gsl_vector * x; + gsl_vector * f; + gsl_matrix * J; + gsl_vector * dx; + void *state; + } +gsl_multifit_fdfsolver; + + +gsl_multifit_fdfsolver * +gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T, + size_t n, size_t p); + +int +gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s, + gsl_multifit_function_fdf * fdf, + const gsl_vector * x); + +int +gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s); + +void +gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s); + +const char * gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s); +gsl_vector * gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s); + +int gsl_multifit_test_delta (const gsl_vector * dx, const gsl_vector * x, + double epsabs, double epsrel); + +int gsl_multifit_test_gradient (const gsl_vector * g, double epsabs); + +/* extern const gsl_multifit_fsolver_type * gsl_multifit_fsolver_gradient; */ + +GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmder; +GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmsder; + + +__END_DECLS + +#endif /* __GSL_MULTIFIT_NLIN_H__ */ diff --git a/gsl-1.9/multifit/lmder.c b/gsl-1.9/multifit/lmder.c new file mode 100644 index 0000000..0dbeca8 --- /dev/null +++ b/gsl-1.9/multifit/lmder.c @@ -0,0 +1,387 @@ +/* multfit/lmder.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 <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_multifit_nlin.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_permutation.h> + + +typedef struct + { + size_t iter; + double xnorm; + double fnorm; + double delta; + double par; + gsl_matrix *r; + gsl_vector *tau; + gsl_vector *diag; + gsl_vector *qtf; + gsl_vector *newton; + gsl_vector *gradient; + gsl_vector *x_trial; + gsl_vector *f_trial; + gsl_vector *df; + gsl_vector *sdiag; + gsl_vector *rptdx; + gsl_vector *w; + gsl_vector *work1; + gsl_permutation * perm; + } +lmder_state_t; + +static int lmder_alloc (void *vstate, size_t n, size_t p); +static int lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static int set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale); +static int lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx); +static void lmder_free (void *vstate); +static int iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale); + +#include "lmutil.c" +#include "lmpar.c" +#include "lmset.c" +#include "lmiterate.c" + + +static int +lmder_alloc (void *vstate, size_t n, size_t p) +{ + lmder_state_t *state = (lmder_state_t *) vstate; + gsl_matrix *r; + gsl_vector *tau, *diag, *qtf, *newton, *gradient, *x_trial, *f_trial, + *df, *sdiag, *rptdx, *w, *work1; + gsl_permutation *perm; + + r = gsl_matrix_calloc (n, p); + + if (r == 0) + { + GSL_ERROR ("failed to allocate space for r", GSL_ENOMEM); + } + + state->r = r; + + tau = gsl_vector_calloc (GSL_MIN(n, p)); + + if (tau == 0) + { + gsl_matrix_free (r); + + GSL_ERROR ("failed to allocate space for tau", GSL_ENOMEM); + } + + state->tau = tau; + + diag = gsl_vector_calloc (p); + + if (diag == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + + GSL_ERROR ("failed to allocate space for diag", GSL_ENOMEM); + } + + state->diag = diag; + + qtf = gsl_vector_calloc (n); + + if (qtf == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + + GSL_ERROR ("failed to allocate space for qtf", GSL_ENOMEM); + } + + state->qtf = qtf; + + newton = gsl_vector_calloc (p); + + if (newton == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + + GSL_ERROR ("failed to allocate space for newton", GSL_ENOMEM); + } + + state->newton = newton; + + gradient = gsl_vector_calloc (p); + + if (gradient == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + + GSL_ERROR ("failed to allocate space for gradient", GSL_ENOMEM); + } + + state->gradient = gradient; + + x_trial = gsl_vector_calloc (p); + + if (x_trial == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + + GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM); + } + + state->x_trial = x_trial; + + f_trial = gsl_vector_calloc (n); + + if (f_trial == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + + GSL_ERROR ("failed to allocate space for f_trial", GSL_ENOMEM); + } + + state->f_trial = f_trial; + + df = gsl_vector_calloc (n); + + if (df == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + + GSL_ERROR ("failed to allocate space for df", GSL_ENOMEM); + } + + state->df = df; + + sdiag = gsl_vector_calloc (p); + + if (sdiag == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + + GSL_ERROR ("failed to allocate space for sdiag", GSL_ENOMEM); + } + + state->sdiag = sdiag; + + + rptdx = gsl_vector_calloc (n); + + if (rptdx == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (sdiag); + + GSL_ERROR ("failed to allocate space for rptdx", GSL_ENOMEM); + } + + state->rptdx = rptdx; + + w = gsl_vector_calloc (n); + + if (w == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (sdiag); + gsl_vector_free (rptdx); + + GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM); + } + + state->w = w; + + work1 = gsl_vector_calloc (p); + + if (work1 == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (sdiag); + gsl_vector_free (rptdx); + gsl_vector_free (w); + + GSL_ERROR ("failed to allocate space for work1", GSL_ENOMEM); + } + + state->work1 = work1; + + perm = gsl_permutation_calloc (p); + + if (perm == 0) + { + gsl_matrix_free (r); + gsl_vector_free (tau); + gsl_vector_free (diag); + gsl_vector_free (qtf); + gsl_vector_free (newton); + gsl_vector_free (gradient); + gsl_vector_free (x_trial); + gsl_vector_free (f_trial); + gsl_vector_free (df); + gsl_vector_free (sdiag); + gsl_vector_free (rptdx); + gsl_vector_free (w); + gsl_vector_free (work1); + + GSL_ERROR ("failed to allocate space for perm", GSL_ENOMEM); + } + + state->perm = perm; + + return GSL_SUCCESS; +} + +static int +lmder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = set (vstate, fdf, x, f, J, dx, 0); + return status ; +} + +static int +lmsder_set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = set (vstate, fdf, x, f, J, dx, 1); + return status ; +} + +static int +lmder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = iterate (vstate, fdf, x, f, J, dx, 0); + return status; +} + +static int +lmsder_iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx) +{ + int status = iterate (vstate, fdf, x, f, J, dx, 1); + return status; +} + +static void +lmder_free (void *vstate) +{ + lmder_state_t *state = (lmder_state_t *) vstate; + + gsl_permutation_free (state->perm); + gsl_vector_free (state->work1); + gsl_vector_free (state->w); + gsl_vector_free (state->rptdx); + gsl_vector_free (state->sdiag); + gsl_vector_free (state->df); + gsl_vector_free (state->f_trial); + gsl_vector_free (state->x_trial); + gsl_vector_free (state->gradient); + gsl_vector_free (state->newton); + gsl_vector_free (state->qtf); + gsl_vector_free (state->diag); + gsl_vector_free (state->tau); + gsl_matrix_free (state->r); +} + +static const gsl_multifit_fdfsolver_type lmder_type = +{ + "lmder", /* name */ + sizeof (lmder_state_t), + &lmder_alloc, + &lmder_set, + &lmder_iterate, + &lmder_free +}; + +static const gsl_multifit_fdfsolver_type lmsder_type = +{ + "lmsder", /* name */ + sizeof (lmder_state_t), + &lmder_alloc, + &lmsder_set, + &lmsder_iterate, + &lmder_free +}; + +const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmder = &lmder_type; +const gsl_multifit_fdfsolver_type *gsl_multifit_fdfsolver_lmsder = &lmsder_type; diff --git a/gsl-1.9/multifit/lmiterate.c b/gsl-1.9/multifit/lmiterate.c new file mode 100644 index 0000000..78d0ea6 --- /dev/null +++ b/gsl-1.9/multifit/lmiterate.c @@ -0,0 +1,222 @@ +static int +iterate (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale) +{ + lmder_state_t *state = (lmder_state_t *) vstate; + + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + gsl_vector *qtf = state->qtf; + gsl_vector *x_trial = state->x_trial; + gsl_vector *f_trial = state->f_trial; + gsl_vector *rptdx = state->rptdx; + gsl_vector *newton = state->newton; + gsl_vector *gradient = state->gradient; + gsl_vector *sdiag = state->sdiag; + gsl_vector *w = state->w; + gsl_vector *work1 = state->work1; + gsl_permutation *perm = state->perm; + + double prered, actred; + double pnorm, fnorm1, fnorm1p, gnorm; + double ratio; + double dirder; + + int iter = 0; + + double p1 = 0.1, p25 = 0.25, p5 = 0.5, p75 = 0.75, p0001 = 0.0001; + + if (state->fnorm == 0.0) + { + return GSL_SUCCESS; + } + + /* Compute qtf = Q^T f */ + + gsl_vector_memcpy (qtf, f); + gsl_linalg_QR_QTvec (r, tau, qtf); + + /* Compute norm of scaled gradient */ + + compute_gradient_direction (r, perm, qtf, diag, gradient); + + { + size_t iamax = gsl_blas_idamax (gradient); + + gnorm = fabs(gsl_vector_get (gradient, iamax) / state->fnorm); + } + + /* Determine the Levenberg-Marquardt parameter */ + +lm_iteration: + + iter++ ; + + { + int status = lmpar (r, perm, qtf, diag, state->delta, &(state->par), newton, gradient, sdiag, dx, w); + if (status) + return status; + } + + /* Take a trial step */ + + gsl_vector_scale (dx, -1.0); /* reverse the step to go downhill */ + + compute_trial_step (x, dx, state->x_trial); + + pnorm = scaled_enorm (diag, dx); + + if (state->iter == 1) + { + if (pnorm < state->delta) + { +#ifdef DEBUG + printf("set delta = pnorm = %g\n" , pnorm); +#endif + state->delta = pnorm; + } + } + + /* Evaluate function at x + p */ + /* return immediately if evaluation raised error */ + { + int status = GSL_MULTIFIT_FN_EVAL_F (fdf, x_trial, f_trial); + if (status) + return status; + } + + fnorm1 = enorm (f_trial); + + /* Compute the scaled actual reduction */ + + actred = compute_actual_reduction (state->fnorm, fnorm1); + +#ifdef DEBUG + printf("lmiterate: fnorm = %g fnorm1 = %g actred = %g\n", state->fnorm, fnorm1, actred); +#endif + + /* Compute rptdx = R P^T dx, noting that |J dx| = |R P^T dx| */ + + compute_rptdx (r, perm, dx, rptdx); + + fnorm1p = enorm (rptdx); + + /* Compute the scaled predicted reduction = |J dx|^2 + 2 par |D dx|^2 */ + + { + double t1 = fnorm1p / state->fnorm; + double t2 = (sqrt(state->par) * pnorm) / state->fnorm; + + prered = t1 * t1 + t2 * t2 / p5; + dirder = -(t1 * t1 + t2 * t2); + } + + /* compute the ratio of the actual to predicted reduction */ + + if (prered > 0) + { + ratio = actred / prered; + } + else + { + ratio = 0; + } + +#ifdef DEBUG + printf("lmiterate: prered = %g dirder = %g ratio = %g\n", prered, dirder,ratio); +#endif + + + /* update the step bound */ + + if (ratio > p25) + { +#ifdef DEBUG + printf("ratio > p25\n"); +#endif + if (state->par == 0 || ratio >= p75) + { + state->delta = pnorm / p5; + state->par *= p5; +#ifdef DEBUG + printf("updated step bounds: delta = %g, par = %g\n", state->delta, state->par); +#endif + } + } + else + { + double temp = (actred >= 0) ? p5 : p5*dirder / (dirder + p5 * actred); + +#ifdef DEBUG + printf("ratio < p25\n"); +#endif + + if (p1 * fnorm1 >= state->fnorm || temp < p1 ) + { + temp = p1; + } + + state->delta = temp * GSL_MIN_DBL (state->delta, pnorm/p1); + + state->par /= temp; +#ifdef DEBUG + printf("updated step bounds: delta = %g, par = %g\n", state->delta, state->par); +#endif + } + + + /* test for successful iteration, termination and stringent tolerances */ + + if (ratio >= p0001) + { + gsl_vector_memcpy (x, x_trial); + gsl_vector_memcpy (f, f_trial); + + /* return immediately if evaluation raised error */ + { + int status = GSL_MULTIFIT_FN_EVAL_DF (fdf, x_trial, J); + if (status) + return status; + } + + /* wa2_j = diag_j * x_j */ + state->xnorm = scaled_enorm(diag, x); + state->fnorm = fnorm1; + state->iter++; + + /* Rescale if necessary */ + + if (scale) + { + update_diag (J, diag); + } + + { + int signum; + gsl_matrix_memcpy (r, J); + gsl_linalg_QRPT_decomp (r, tau, perm, &signum, work1); + } + + return GSL_SUCCESS; + } + else if (fabs(actred) <= GSL_DBL_EPSILON && prered <= GSL_DBL_EPSILON + && p5 * ratio <= 1.0) + { + return GSL_ETOLF ; + } + else if (state->delta <= GSL_DBL_EPSILON * state->xnorm) + { + return GSL_ETOLX; + } + else if (gnorm <= GSL_DBL_EPSILON) + { + return GSL_ETOLG; + } + else if (iter < 10) + { + /* Repeat inner loop if unsuccessful */ + goto lm_iteration; + } + + return GSL_CONTINUE; +} diff --git a/gsl-1.9/multifit/lmpar.c b/gsl-1.9/multifit/lmpar.c new file mode 100644 index 0000000..3f6f105 --- /dev/null +++ b/gsl-1.9/multifit/lmpar.c @@ -0,0 +1,473 @@ +/* multifit/lmpar.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 <gsl/gsl_permute_vector_double.h> + +#include "qrsolv.c" + +static size_t +count_nsing (const gsl_matrix * r) +{ + /* Count the number of nonsingular entries. Returns the index of the + first entry which is singular. */ + + size_t n = r->size2; + size_t i; + + for (i = 0; i < n; i++) + { + double rii = gsl_matrix_get (r, i, i); + + if (rii == 0) + { + break; + } + } + + return i; +} + + +static void +compute_newton_direction (const gsl_matrix * r, const gsl_permutation * perm, + const gsl_vector * qtf, gsl_vector * x) +{ + + /* Compute and store in x the Gauss-Newton direction. If the + Jacobian is rank-deficient then obtain a least squares + solution. */ + + const size_t n = r->size2; + size_t i, j, nsing; + + for (i = 0 ; i < n ; i++) + { + double qtfi = gsl_vector_get (qtf, i); + gsl_vector_set (x, i, qtfi); + } + + nsing = count_nsing (r); + +#ifdef DEBUG + printf("nsing = %d\n", nsing); + printf("r = "); gsl_matrix_fprintf(stdout, r, "%g"); printf("\n"); + printf("qtf = "); gsl_vector_fprintf(stdout, x, "%g"); printf("\n"); +#endif + + for (i = nsing; i < n; i++) + { + gsl_vector_set (x, i, 0.0); + } + + if (nsing > 0) + { + for (j = nsing; j > 0 && j--;) + { + double rjj = gsl_matrix_get (r, j, j); + double temp = gsl_vector_get (x, j) / rjj; + + gsl_vector_set (x, j, temp); + + for (i = 0; i < j; i++) + { + double rij = gsl_matrix_get (r, i, j); + double xi = gsl_vector_get (x, i); + gsl_vector_set (x, i, xi - rij * temp); + } + } + } + + gsl_permute_vector_inverse (perm, x); +} + +static void +compute_newton_correction (const gsl_matrix * r, const gsl_vector * sdiag, + const gsl_permutation * p, gsl_vector * x, + double dxnorm, + const gsl_vector * diag, gsl_vector * w) +{ + size_t n = r->size2; + size_t i, j; + + for (i = 0; i < n; i++) + { + size_t pi = gsl_permutation_get (p, i); + + double dpi = gsl_vector_get (diag, pi); + double xpi = gsl_vector_get (x, pi); + + gsl_vector_set (w, i, dpi * (dpi * xpi) / dxnorm); + } + + for (j = 0; j < n; j++) + { + double sj = gsl_vector_get (sdiag, j); + double wj = gsl_vector_get (w, j); + + double tj = wj / sj; + + gsl_vector_set (w, j, tj); + + for (i = j + 1; i < n; i++) + { + double rij = gsl_matrix_get (r, i, j); + double wi = gsl_vector_get (w, i); + + gsl_vector_set (w, i, wi - rij * tj); + } + } +} + +static void +compute_newton_bound (const gsl_matrix * r, const gsl_vector * x, + double dxnorm, const gsl_permutation * perm, + const gsl_vector * diag, gsl_vector * w) +{ + /* If the jacobian is not rank-deficient then the Newton step + provides a lower bound for the zero of the function. Otherwise + set this bound to zero. */ + + size_t n = r->size2; + + size_t i, j; + + size_t nsing = count_nsing (r); + + if (nsing < n) + { + gsl_vector_set_zero (w); + return; + } + + for (i = 0; i < n; i++) + { + size_t pi = gsl_permutation_get (perm, i); + + double dpi = gsl_vector_get (diag, pi); + double xpi = gsl_vector_get (x, pi); + + gsl_vector_set (w, i, dpi * (dpi * xpi / dxnorm)); + } + + for (j = 0; j < n; j++) + { + double sum = 0; + + for (i = 0; i < j; i++) + { + sum += gsl_matrix_get (r, i, j) * gsl_vector_get (w, i); + } + + { + double rjj = gsl_matrix_get (r, j, j); + double wj = gsl_vector_get (w, j); + + gsl_vector_set (w, j, (wj - sum) / rjj); + } + } +} + +static void +compute_gradient_direction (const gsl_matrix * r, const gsl_permutation * p, + const gsl_vector * qtf, const gsl_vector * diag, + gsl_vector * g) +{ + const size_t n = r->size2; + + size_t i, j; + + for (j = 0; j < n; j++) + { + double sum = 0; + + for (i = 0; i <= j; i++) + { + sum += gsl_matrix_get (r, i, j) * gsl_vector_get (qtf, i); + } + + { + size_t pj = gsl_permutation_get (p, j); + double dpj = gsl_vector_get (diag, pj); + + gsl_vector_set (g, j, sum / dpj); + } + } +} + +static int +lmpar (gsl_matrix * r, const gsl_permutation * perm, const gsl_vector * qtf, + const gsl_vector * diag, double delta, double * par_inout, + gsl_vector * newton, gsl_vector * gradient, gsl_vector * sdiag, + gsl_vector * x, gsl_vector * w) +{ + double dxnorm, gnorm, fp, fp_old, par_lower, par_upper, par_c; + + double par = *par_inout; + + size_t iter = 0; + +#ifdef DEBUG + printf("ENTERING lmpar\n"); +#endif + + + compute_newton_direction (r, perm, qtf, newton); + +#ifdef DEBUG + printf ("newton = "); + gsl_vector_fprintf (stdout, newton, "%g"); + printf ("\n"); + + printf ("diag = "); + gsl_vector_fprintf (stdout, diag, "%g"); + printf ("\n"); +#endif + + /* Evaluate the function at the origin and test for acceptance of + the Gauss-Newton direction. */ + + dxnorm = scaled_enorm (diag, newton); + + fp = dxnorm - delta; + +#ifdef DEBUG + printf ("dxnorm = %g, delta = %g, fp = %g\n", dxnorm, delta, fp); +#endif + + if (fp <= 0.1 * delta) + { + gsl_vector_memcpy (x, newton); +#ifdef DEBUG + printf ("took newton (fp = %g, delta = %g)\n", fp, delta); +#endif + + *par_inout = 0; + + return GSL_SUCCESS; + } + + compute_newton_bound (r, newton, dxnorm, perm, diag, w); + + { + double wnorm = enorm (w); + double phider = wnorm * wnorm; + + /* w == zero if r rank-deficient, + then set lower bound to zero form MINPACK, lmder.f + Hans E. Plesser 2002-02-25 (hans.plesser@itf.nlh.no) */ + if ( wnorm > 0 ) + par_lower = fp / (delta * phider); + else + par_lower = 0.0; + } + +#ifdef DEBUG + printf("par = %g\n", par ); + printf("par_lower = %g\n", par_lower); +#endif + + compute_gradient_direction (r, perm, qtf, diag, gradient); + + gnorm = enorm (gradient); + +#ifdef DEBUG + printf("gradient = "); gsl_vector_fprintf(stdout, gradient, "%g"); printf("\n"); + printf("gnorm = %g\n", gnorm); +#endif + + par_upper = gnorm / delta; + + if (par_upper == 0) + { + par_upper = GSL_DBL_MIN / GSL_MIN_DBL(delta, 0.1); + } + +#ifdef DEBUG + printf("par_upper = %g\n", par_upper); +#endif + + if (par > par_upper) + { +#ifdef DEBUG + printf("set par to par_upper\n"); +#endif + + par = par_upper; + } + else if (par < par_lower) + { +#ifdef DEBUG + printf("set par to par_lower\n"); +#endif + + par = par_lower; + } + + if (par == 0) + { + par = gnorm / dxnorm; +#ifdef DEBUG + printf("set par to gnorm/dxnorm = %g\n", par); +#endif + + } + + /* Beginning of iteration */ + +iteration: + + iter++; + +#ifdef DEBUG + printf("lmpar iteration = %d\n", iter); +#endif + +#ifdef BRIANSFIX + /* Seems like this is described in the paper but not in the MINPACK code */ + + if (par < par_lower || par > par_upper) + { + par = GSL_MAX_DBL (0.001 * par_upper, sqrt(par_lower * par_upper)); + } +#endif + + /* Evaluate the function at the current value of par */ + + if (par == 0) + { + par = GSL_MAX_DBL (0.001 * par_upper, GSL_DBL_MIN); +#ifdef DEBUG + printf("par = 0, set par to = %g\n", par); +#endif + + } + + /* Compute the least squares solution of [ R P x - Q^T f, sqrt(par) D x] + for A = Q R P^T */ + +#ifdef DEBUG + printf ("calling qrsolv with par = %g\n", par); +#endif + + { + double sqrt_par = sqrt(par); + + qrsolv (r, perm, sqrt_par, diag, qtf, x, sdiag, w); + } + + dxnorm = scaled_enorm (diag, x); + + fp_old = fp; + + fp = dxnorm - delta; + +#ifdef DEBUG + printf ("After qrsolv dxnorm = %g, delta = %g, fp = %g\n", dxnorm, delta, fp); + printf ("sdiag = ") ; gsl_vector_fprintf(stdout, sdiag, "%g"); printf("\n"); + printf ("x = ") ; gsl_vector_fprintf(stdout, x, "%g"); printf("\n"); + printf ("r = ") ; gsl_matrix_fprintf(stdout, r, "%g"); printf("\nXXX\n"); +#endif + + /* If the function is small enough, accept the current value of par */ + + if (fabs (fp) <= 0.1 * delta) + goto line220; + + if (par_lower == 0 && fp <= fp_old && fp_old < 0) + goto line220; + + /* Check for maximum number of iterations */ + + if (iter == 10) + goto line220; + + /* Compute the Newton correction */ + + compute_newton_correction (r, sdiag, perm, x, dxnorm, diag, w); + +#ifdef DEBUG + printf ("newton_correction = "); + gsl_vector_fprintf(stdout, w, "%g"); printf("\n"); +#endif + + { + double wnorm = enorm (w); + par_c = fp / (delta * wnorm * wnorm); + } + +#ifdef DEBUG + printf("fp = %g\n", fp); + printf("par_lower = %g\n", par_lower); + printf("par_upper = %g\n", par_upper); + printf("par_c = %g\n", par_c); +#endif + + + /* Depending on the sign of the function, update par_lower or par_upper */ + + if (fp > 0) + { + if (par > par_lower) + { + par_lower = par; +#ifdef DEBUG + printf("fp > 0: set par_lower = par = %g\n", par); +#endif + + } + } + else if (fp < 0) + { + if (par < par_upper) + { +#ifdef DEBUG + printf("fp < 0: set par_upper = par = %g\n", par); +#endif + par_upper = par; + } + } + + /* Compute an improved estimate for par */ + +#ifdef DEBUG + printf("improved estimate par = MAX(%g, %g) \n", par_lower, par+par_c); +#endif + + par = GSL_MAX_DBL (par_lower, par + par_c); + +#ifdef DEBUG + printf("improved estimate par = %g \n", par); +#endif + + + goto iteration; + +line220: + +#ifdef DEBUG + printf("LEAVING lmpar, par = %g\n", par); +#endif + + *par_inout = par; + + return GSL_SUCCESS; +} + + + diff --git a/gsl-1.9/multifit/lmset.c b/gsl-1.9/multifit/lmset.c new file mode 100644 index 0000000..f6a97cc --- /dev/null +++ b/gsl-1.9/multifit/lmset.c @@ -0,0 +1,44 @@ +static int +set (void *vstate, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx, int scale) +{ + lmder_state_t *state = (lmder_state_t *) vstate; + + gsl_matrix *r = state->r; + gsl_vector *tau = state->tau; + gsl_vector *diag = state->diag; + gsl_vector *work1 = state->work1; + gsl_permutation *perm = state->perm; + + int signum; + + GSL_MULTIFIT_FN_EVAL_F_DF (fdf, x, f, J); + + state->par = 0; + state->iter = 1; + state->fnorm = enorm (f); + + gsl_vector_set_all (dx, 0.0); + + /* store column norms in diag */ + + if (scale) + { + compute_diag (J, diag); + } + else + { + gsl_vector_set_all (diag, 1.0); + } + + /* set delta to 100 |D x| or to 100 if |D x| is zero */ + + state->xnorm = scaled_enorm (diag, x); + state->delta = compute_delta (diag, x); + + /* Factorize J into QR decomposition */ + + gsl_matrix_memcpy (r, J); + gsl_linalg_QRPT_decomp (r, tau, perm, &signum, work1); + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multifit/lmutil.c b/gsl-1.9/multifit/lmutil.c new file mode 100644 index 0000000..1f06178 --- /dev/null +++ b/gsl-1.9/multifit/lmutil.c @@ -0,0 +1,134 @@ +static void compute_diag (const gsl_matrix * J, gsl_vector * diag); +static void update_diag (const gsl_matrix * J, gsl_vector * diag); +static double compute_delta (gsl_vector * diag, gsl_vector * x); +static double scaled_enorm (const gsl_vector * d, const gsl_vector * f); +static double enorm (const gsl_vector * f); + +static double +enorm (const gsl_vector * f) +{ + return gsl_blas_dnrm2 (f); +} + +static double +scaled_enorm (const gsl_vector * d, const gsl_vector * f) +{ + double e2 = 0; + size_t i, n = f->size; + for (i = 0; i < n; i++) + { + double fi = gsl_vector_get (f, i); + double di = gsl_vector_get (d, i); + double u = di * fi; + e2 += u * u; + } + return sqrt (e2); +} + +static double +compute_delta (gsl_vector * diag, gsl_vector * x) +{ + double Dx = scaled_enorm (diag, x); + double factor = 100; /* generally recommended value from MINPACK */ + + return (Dx > 0) ? factor * Dx : factor; +} + +static double +compute_actual_reduction (double fnorm, double fnorm1) +{ + double actred; + + if (0.1 * fnorm1 < fnorm) + { + double u = fnorm1 / fnorm; + actred = 1 - u * u; + } + else + { + actred = -1; + } + + return actred; +} + +static void +compute_diag (const gsl_matrix * J, gsl_vector * diag) +{ + size_t i, j, n = J->size1, p = J->size2; + + for (j = 0; j < p; j++) + { + double sum = 0; + + for (i = 0; i < n; i++) + { + double Jij = gsl_matrix_get (J, i, j); + sum += Jij * Jij; + } + if (sum == 0) + sum = 1.0; + + gsl_vector_set (diag, j, sqrt (sum)); + } +} + +static void +update_diag (const gsl_matrix * J, gsl_vector * diag) +{ + size_t i, j, n = diag->size; + + for (j = 0; j < n; j++) + { + double cnorm, diagj, sum = 0; + for (i = 0; i < n; i++) + { + double Jij = gsl_matrix_get (J, i, j); + sum += Jij * Jij; + } + if (sum == 0) + sum = 1.0; + + cnorm = sqrt (sum); + diagj = gsl_vector_get (diag, j); + + if (cnorm > diagj) + gsl_vector_set (diag, j, cnorm); + } +} + +static void +compute_rptdx (const gsl_matrix * r, const gsl_permutation * p, + const gsl_vector * dx, gsl_vector * rptdx) +{ + size_t i, j, N = dx->size; + + for (i = 0; i < N; i++) + { + double sum = 0; + + for (j = i; j < N; j++) + { + size_t pj = gsl_permutation_get (p, j); + + sum += gsl_matrix_get (r, i, j) * gsl_vector_get (dx, pj); + } + + gsl_vector_set (rptdx, i, sum); + } +} + + +static void +compute_trial_step (gsl_vector * x, gsl_vector * dx, gsl_vector * x_trial) +{ + size_t i, N = x->size; + + for (i = 0; i < N; i++) + { + double pi = gsl_vector_get (dx, i); + double xi = gsl_vector_get (x, i); + gsl_vector_set (x_trial, i, xi + pi); + } +} + diff --git a/gsl-1.9/multifit/multilinear.c b/gsl-1.9/multifit/multilinear.c new file mode 100644 index 0000000..6616809 --- /dev/null +++ b/gsl-1.9/multifit/multilinear.c @@ -0,0 +1,432 @@ +/* multifit/multilinear.c + * + * Copyright (C) 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_multifit.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_linalg.h> + +/* Fit + * + * y = X c + * + * where X is an M x N matrix of M observations for N variables. + * + */ + +int +gsl_multifit_linear (const gsl_matrix * X, + const gsl_vector * y, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, gsl_multifit_linear_workspace * work) +{ + size_t rank; + int status = gsl_multifit_linear_svd (X, y, GSL_DBL_EPSILON, &rank, c, + cov, chisq, work); + return status; +} + +/* Handle the general case of the SVD with tolerance and rank */ + +int +gsl_multifit_linear_svd (const gsl_matrix * X, + const gsl_vector * y, + double tol, + size_t * rank, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, gsl_multifit_linear_workspace * work) +{ + if (X->size1 != y->size) + { + GSL_ERROR + ("number of observations in y does not match rows of matrix X", + GSL_EBADLEN); + } + else if (X->size2 != c->size) + { + GSL_ERROR ("number of parameters c does not match columns of matrix X", + GSL_EBADLEN); + } + else if (cov->size1 != cov->size2) + { + GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR); + } + else if (c->size != cov->size1) + { + GSL_ERROR + ("number of parameters does not match size of covariance matrix", + GSL_EBADLEN); + } + else if (X->size1 != work->n || X->size2 != work->p) + { + GSL_ERROR + ("size of workspace does not match size of observation matrix", + GSL_EBADLEN); + } + else if (tol <= 0) + { + GSL_ERROR ("tolerance must be positive", GSL_EINVAL); + } + else + { + const size_t n = X->size1; + const size_t p = X->size2; + + size_t i, j, p_eff; + + gsl_matrix *A = work->A; + gsl_matrix *Q = work->Q; + gsl_matrix *QSI = work->QSI; + gsl_vector *S = work->S; + gsl_vector *xt = work->xt; + gsl_vector *D = work->D; + + /* Copy X to workspace, A <= X */ + + gsl_matrix_memcpy (A, X); + + /* Balance the columns of the matrix A */ + + gsl_linalg_balance_columns (A, D); + + /* Decompose A into U S Q^T */ + + gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt); + + /* Solve y = A c for c */ + + gsl_blas_dgemv (CblasTrans, 1.0, A, y, 0.0, xt); + + /* Scale the matrix Q, Q' = Q S^-1 */ + + gsl_matrix_memcpy (QSI, Q); + + { + double alpha0 = gsl_vector_get (S, 0); + p_eff = 0; + + for (j = 0; j < p; j++) + { + gsl_vector_view column = gsl_matrix_column (QSI, j); + double alpha = gsl_vector_get (S, j); + + if (alpha <= tol * alpha0) { + alpha = 0.0; + } else { + alpha = 1.0 / alpha; + p_eff++; + } + + gsl_vector_scale (&column.vector, alpha); + } + + *rank = p_eff; + } + + gsl_vector_set_zero (c); + + gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c); + + /* Unscale the balancing factors */ + + gsl_vector_div (c, D); + + /* Compute chisq, from residual r = y - X c */ + + { + double s2 = 0, r2 = 0; + + for (i = 0; i < n; i++) + { + double yi = gsl_vector_get (y, i); + gsl_vector_const_view row = gsl_matrix_const_row (X, i); + double y_est, ri; + gsl_blas_ddot (&row.vector, c, &y_est); + ri = yi - y_est; + r2 += ri * ri; + } + + s2 = r2 / (n - p_eff); /* p_eff == rank */ + + *chisq = r2; + + /* Form variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */ + + for (i = 0; i < p; i++) + { + gsl_vector_view row_i = gsl_matrix_row (QSI, i); + double d_i = gsl_vector_get (D, i); + + for (j = i; j < p; j++) + { + gsl_vector_view row_j = gsl_matrix_row (QSI, j); + double d_j = gsl_vector_get (D, j); + double s; + + gsl_blas_ddot (&row_i.vector, &row_j.vector, &s); + + gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j)); + gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j)); + } + } + } + + return GSL_SUCCESS; + } +} + +int +gsl_multifit_wlinear (const gsl_matrix * X, + const gsl_vector * w, + const gsl_vector * y, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, gsl_multifit_linear_workspace * work) +{ + size_t rank; + int status = gsl_multifit_wlinear_svd (X, w, y, GSL_DBL_EPSILON, &rank, c, + cov, chisq, work); + return status; +} + + +int +gsl_multifit_wlinear_svd (const gsl_matrix * X, + const gsl_vector * w, + const gsl_vector * y, + double tol, + size_t * rank, + gsl_vector * c, + gsl_matrix * cov, + double *chisq, gsl_multifit_linear_workspace * work) +{ + if (X->size1 != y->size) + { + GSL_ERROR + ("number of observations in y does not match rows of matrix X", + GSL_EBADLEN); + } + else if (X->size2 != c->size) + { + GSL_ERROR ("number of parameters c does not match columns of matrix X", + GSL_EBADLEN); + } + else if (w->size != y->size) + { + GSL_ERROR ("number of weights does not match number of observations", + GSL_EBADLEN); + } + else if (cov->size1 != cov->size2) + { + GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR); + } + else if (c->size != cov->size1) + { + GSL_ERROR + ("number of parameters does not match size of covariance matrix", + GSL_EBADLEN); + } + else if (X->size1 != work->n || X->size2 != work->p) + { + GSL_ERROR + ("size of workspace does not match size of observation matrix", + GSL_EBADLEN); + } + else + { + const size_t n = X->size1; + const size_t p = X->size2; + + size_t i, j, p_eff; + + gsl_matrix *A = work->A; + gsl_matrix *Q = work->Q; + gsl_matrix *QSI = work->QSI; + gsl_vector *S = work->S; + gsl_vector *t = work->t; + gsl_vector *xt = work->xt; + gsl_vector *D = work->D; + + /* Scale X, A = sqrt(w) X */ + + gsl_matrix_memcpy (A, X); + + for (i = 0; i < n; i++) + { + double wi = gsl_vector_get (w, i); + + if (wi < 0) + wi = 0; + + { + gsl_vector_view row = gsl_matrix_row (A, i); + gsl_vector_scale (&row.vector, sqrt (wi)); + } + } + + /* Balance the columns of the matrix A */ + + gsl_linalg_balance_columns (A, D); + + /* Decompose A into U S Q^T */ + + gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt); + + /* Solve sqrt(w) y = A c for c, by first computing t = sqrt(w) y */ + + for (i = 0; i < n; i++) + { + double wi = gsl_vector_get (w, i); + double yi = gsl_vector_get (y, i); + if (wi < 0) + wi = 0; + gsl_vector_set (t, i, sqrt (wi) * yi); + } + + gsl_blas_dgemv (CblasTrans, 1.0, A, t, 0.0, xt); + + /* Scale the matrix Q, Q' = Q S^-1 */ + + gsl_matrix_memcpy (QSI, Q); + + { + double alpha0 = gsl_vector_get (S, 0); + p_eff = 0; + + for (j = 0; j < p; j++) + { + gsl_vector_view column = gsl_matrix_column (QSI, j); + double alpha = gsl_vector_get (S, j); + + if (alpha <= tol * alpha0) { + alpha = 0.0; + } else { + alpha = 1.0 / alpha; + p_eff++; + } + + gsl_vector_scale (&column.vector, alpha); + } + + *rank = p_eff; + } + + gsl_vector_set_zero (c); + + /* Solution */ + + gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c); + + /* Unscale the balancing factors */ + + gsl_vector_div (c, D); + + /* Form covariance matrix cov = (Q S^-1) (Q S^-1)^T */ + + for (i = 0; i < p; i++) + { + gsl_vector_view row_i = gsl_matrix_row (QSI, i); + double d_i = gsl_vector_get (D, i); + + for (j = i; j < p; j++) + { + gsl_vector_view row_j = gsl_matrix_row (QSI, j); + double d_j = gsl_vector_get (D, j); + double s; + + gsl_blas_ddot (&row_i.vector, &row_j.vector, &s); + + gsl_matrix_set (cov, i, j, s / (d_i * d_j)); + gsl_matrix_set (cov, j, i, s / (d_i * d_j)); + } + } + + /* Compute chisq, from residual r = y - X c */ + + { + double r2 = 0; + + for (i = 0; i < n; i++) + { + double yi = gsl_vector_get (y, i); + double wi = gsl_vector_get (w, i); + gsl_vector_const_view row = gsl_matrix_const_row (X, i); + double y_est, ri; + gsl_blas_ddot (&row.vector, c, &y_est); + ri = yi - y_est; + r2 += wi * ri * ri; + } + + *chisq = r2; + } + + return GSL_SUCCESS; + } +} + + +int +gsl_multifit_linear_est (const gsl_vector * x, + const gsl_vector * c, + const gsl_matrix * cov, double *y, double *y_err) +{ + + if (x->size != c->size) + { + GSL_ERROR ("number of parameters c does not match number of observations x", + GSL_EBADLEN); + } + else if (cov->size1 != cov->size2) + { + GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR); + } + else if (c->size != cov->size1) + { + GSL_ERROR ("number of parameters c does not match size of covariance matrix cov", + GSL_EBADLEN); + } + else + { + size_t i, j; + double var = 0; + + gsl_blas_ddot(x, c, y); /* y = x.c */ + + /* var = x' cov x */ + + for (i = 0; i < x->size; i++) + { + const double xi = gsl_vector_get (x, i); + var += xi * xi * gsl_matrix_get (cov, i, i); + + for (j = 0; j < i; j++) + { + const double xj = gsl_vector_get (x, j); + var += 2 * xi * xj * gsl_matrix_get (cov, i, j); + } + } + + *y_err = sqrt (var); + + return GSL_SUCCESS; + } +} diff --git a/gsl-1.9/multifit/qrsolv.c b/gsl-1.9/multifit/qrsolv.c new file mode 100644 index 0000000..06f0399 --- /dev/null +++ b/gsl-1.9/multifit/qrsolv.c @@ -0,0 +1,218 @@ +/* This function computes the solution to the least squares system + + phi = [ A x = b , lambda D x = 0 ]^2 + + where A is an M by N matrix, D is an N by N diagonal matrix, lambda + is a scalar parameter and b is a vector of length M. + + The function requires the factorization of A into A = Q R P^T, + where Q is an orthogonal matrix, R is an upper triangular matrix + with diagonal elements of non-increasing magnitude and P is a + permuation matrix. The system above is then equivalent to + + [ R z = Q^T b, P^T (lambda D) P z = 0 ] + + where x = P z. If this system does not have full rank then a least + squares solution is obtained. On output the function also provides + an upper triangular matrix S such that + + P^T (A^T A + lambda^2 D^T D) P = S^T S + + Parameters, + + r: On input, contains the full upper triangle of R. On output the + strict lower triangle contains the transpose of the strict upper + triangle of S, and the diagonal of S is stored in sdiag. The full + upper triangle of R is not modified. + + p: the encoded form of the permutation matrix P. column j of P is + column p[j] of the identity matrix. + + lambda, diag: contains the scalar lambda and the diagonal elements + of the matrix D + + qtb: contains the product Q^T b + + x: on output contains the least squares solution of the system + + wa: is a workspace of length N + + */ + +static int +qrsolv (gsl_matrix * r, const gsl_permutation * p, const double lambda, + const gsl_vector * diag, const gsl_vector * qtb, + gsl_vector * x, gsl_vector * sdiag, gsl_vector * wa) +{ + size_t n = r->size2; + + size_t i, j, k, nsing; + + /* Copy r and qtb to preserve input and initialise s. In particular, + save the diagonal elements of r in x */ + + for (j = 0; j < n; j++) + { + double rjj = gsl_matrix_get (r, j, j); + double qtbj = gsl_vector_get (qtb, j); + + for (i = j + 1; i < n; i++) + { + double rji = gsl_matrix_get (r, j, i); + gsl_matrix_set (r, i, j, rji); + } + + gsl_vector_set (x, j, rjj); + gsl_vector_set (wa, j, qtbj); + } + + /* Eliminate the diagonal matrix d using a Givens rotation */ + + for (j = 0; j < n; j++) + { + double qtbpj; + + size_t pj = gsl_permutation_get (p, j); + + double diagpj = lambda * gsl_vector_get (diag, pj); + + if (diagpj == 0) + { + continue; + } + + gsl_vector_set (sdiag, j, diagpj); + + for (k = j + 1; k < n; k++) + { + gsl_vector_set (sdiag, k, 0.0); + } + + /* The transformations to eliminate the row of d modify only a + single element of qtb beyond the first n, which is initially + zero */ + + qtbpj = 0; + + for (k = j; k < n; k++) + { + /* Determine a Givens rotation which eliminates the + appropriate element in the current row of d */ + + double sine, cosine; + + double wak = gsl_vector_get (wa, k); + double rkk = gsl_matrix_get (r, k, k); + double sdiagk = gsl_vector_get (sdiag, k); + + if (sdiagk == 0) + { + continue; + } + + if (fabs (rkk) < fabs (sdiagk)) + { + double cotangent = rkk / sdiagk; + sine = 0.5 / sqrt (0.25 + 0.25 * cotangent * cotangent); + cosine = sine * cotangent; + } + else + { + double tangent = sdiagk / rkk; + cosine = 0.5 / sqrt (0.25 + 0.25 * tangent * tangent); + sine = cosine * tangent; + } + + /* Compute the modified diagonal element of r and the + modified element of [qtb,0] */ + + { + double new_rkk = cosine * rkk + sine * sdiagk; + double new_wak = cosine * wak + sine * qtbpj; + + qtbpj = -sine * wak + cosine * qtbpj; + + gsl_matrix_set(r, k, k, new_rkk); + gsl_vector_set(wa, k, new_wak); + } + + /* Accumulate the transformation in the row of s */ + + for (i = k + 1; i < n; i++) + { + double rik = gsl_matrix_get (r, i, k); + double sdiagi = gsl_vector_get (sdiag, i); + + double new_rik = cosine * rik + sine * sdiagi; + double new_sdiagi = -sine * rik + cosine * sdiagi; + + gsl_matrix_set(r, i, k, new_rik); + gsl_vector_set(sdiag, i, new_sdiagi); + } + } + + /* Store the corresponding diagonal element of s and restore the + corresponding diagonal element of r */ + + { + double rjj = gsl_matrix_get (r, j, j); + double xj = gsl_vector_get(x, j); + + gsl_vector_set (sdiag, j, rjj); + gsl_matrix_set (r, j, j, xj); + } + + } + + /* Solve the triangular system for z. If the system is singular then + obtain a least squares solution */ + + nsing = n; + + for (j = 0; j < n; j++) + { + double sdiagj = gsl_vector_get (sdiag, j); + + if (sdiagj == 0) + { + nsing = j; + break; + } + } + + for (j = nsing; j < n; j++) + { + gsl_vector_set (wa, j, 0.0); + } + + for (k = 0; k < nsing; k++) + { + double sum = 0; + + j = (nsing - 1) - k; + + for (i = j + 1; i < nsing; i++) + { + sum += gsl_matrix_get(r, i, j) * gsl_vector_get(wa, i); + } + + { + double waj = gsl_vector_get (wa, j); + double sdiagj = gsl_vector_get (sdiag, j); + + gsl_vector_set (wa, j, (waj - sum) / sdiagj); + } + } + + /* Permute the components of z back to the components of x */ + + for (j = 0; j < n; j++) + { + size_t pj = gsl_permutation_get (p, j); + double waj = gsl_vector_get (wa, j); + + gsl_vector_set (x, pj, waj); + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multifit/test.c b/gsl-1.9/multifit/test.c new file mode 100644 index 0000000..fede15a --- /dev/null +++ b/gsl-1.9/multifit/test.c @@ -0,0 +1,210 @@ +/* These tests are based on the NIST Statistical Reference Datasets + See http://www.nist.gov/itl/div898/strd/index.html for more + information. */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_multifit.h> +#include <gsl/gsl_multifit_nlin.h> +#include <gsl/gsl_blas.h> + +#include <gsl/gsl_ieee_utils.h> + +#include "test_longley.c" +#include "test_filip.c" +#include "test_pontius.c" +#include "test_brown.c" +#include "test_enso.c" +#include "test_kirby2.c" +#include "test_hahn1.c" +#include "test_nelson.c" +#include "test_fn.c" +#include "test_estimator.c" + +void +test_lmder (gsl_multifit_function_fdf * f, double x0[], + double * X, double F[], double * cov); + +void +test_fdf (const char * name, gsl_multifit_function_fdf * f, + double x0[], double x[], double sumsq, + double sigma[]); + +int +main (void) +{ + gsl_ieee_env_setup(); + + test_longley(); + test_filip(); + test_pontius(); + test_estimator(); + + { + gsl_multifit_function_fdf f = make_fdf (&brown_f, &brown_df, &brown_fdf, + brown_N, brown_P, 0); + + test_lmder(&f, brown_x0, &brown_X[0][0], brown_F, &brown_cov[0][0]); + } + + { + gsl_multifit_function_fdf f = make_fdf (&enso_f, &enso_df, &enso_fdf, + enso_N, enso_P, 0); + + test_fdf("nist-ENSO", &f, enso_x0, enso_x, enso_sumsq, enso_sigma); + } + + { + gsl_multifit_function_fdf f = make_fdf (&kirby2_f, &kirby2_df, &kirby2_fdf, + kirby2_N, kirby2_P, 0); + + test_fdf("nist-kirby2", &f, kirby2_x0, kirby2_x, kirby2_sumsq, kirby2_sigma); + } + + { + gsl_multifit_function_fdf f = make_fdf (&hahn1_f, &hahn1_df, &hahn1_fdf, + hahn1_N, hahn1_P, 0); + + test_fdf("nist-hahn1", &f, hahn1_x0, hahn1_x, hahn1_sumsq, hahn1_sigma); + } + +#ifdef JUNK + { + gsl_multifit_function_fdf f = make_fdf (&nelson_f, &nelson_df, &nelson_fdf, + nelson_N, nelson_P, 0); + + test_fdf("nist-nelson", &f, nelson_x0, nelson_x, nelson_sumsq, nelson_sigma); + } +#endif + + /* now summarize the results */ + + exit (gsl_test_summary ()); +} + + +void +test_lmder (gsl_multifit_function_fdf * f, double x0[], + double * X, double F[], double * cov) +{ + const gsl_multifit_fdfsolver_type *T; + gsl_multifit_fdfsolver *s; + + const size_t n = f->n; + const size_t p = f->p; + + int status; + size_t iter = 0, i; + + gsl_vector_view x = gsl_vector_view_array (x0, p); + + T = gsl_multifit_fdfsolver_lmsder; + s = gsl_multifit_fdfsolver_alloc (T, n, p); + gsl_multifit_fdfsolver_set (s, f, &x.vector); + + do + { + status = gsl_multifit_fdfsolver_iterate (s); + + for (i = 0 ; i < p; i++) + { + gsl_test_rel (gsl_vector_get (s->x, i), X[p*iter+i], 1e-5, + "lmsder, iter=%u, x%u", iter, i); + } + + gsl_test_rel (gsl_blas_dnrm2 (s->f), F[iter], 1e-5, + "lmsder, iter=%u, f", iter); + + iter++; + } + while (iter < 20); + + { + size_t i, j; + gsl_matrix * covar = gsl_matrix_alloc (4, 4); + gsl_multifit_covar (s->J, 0.0, covar); + + for (i = 0; i < 4; i++) + { + for (j = 0; j < 4; j++) + { + gsl_test_rel (gsl_matrix_get(covar,i,j), cov[i*p + j], 1e-7, + "gsl_multifit_covar cov(%d,%d)", i, j) ; + } + } + + gsl_matrix_free (covar); + } + + gsl_multifit_fdfsolver_free (s); + +} + +void +test_fdf (const char * name, gsl_multifit_function_fdf * f, + double x0[], double x_final[], + double f_sumsq, double sigma[]) +{ + const gsl_multifit_fdfsolver_type *T; + gsl_multifit_fdfsolver *s; + + const size_t n = f->n; + const size_t p = f->p; + + int status; + size_t iter = 0; + + gsl_vector_view x = gsl_vector_view_array (x0, p); + + T = gsl_multifit_fdfsolver_lmsder; + s = gsl_multifit_fdfsolver_alloc (T, n, p); + gsl_multifit_fdfsolver_set (s, f, &x.vector); + + do + { + status = gsl_multifit_fdfsolver_iterate (s); + +#ifdef DEBUG + printf("iter = %d status = %d |f| = %.18e x = \n", + iter, status, gsl_blas_dnrm2 (s->f)); + + gsl_vector_fprintf(stdout, s->x, "%.8e"); +#endif + status = gsl_multifit_test_delta (s->dx, s->x, 0.0, 1e-7); + + iter++; + } + while (status == GSL_CONTINUE && iter < 1000); + + { + size_t i; + gsl_matrix * covar = gsl_matrix_alloc (p, p); + gsl_multifit_covar (s->J, 0.0, covar); + + for (i = 0 ; i < p; i++) + { + gsl_test_rel (gsl_vector_get (s->x, i), x_final[i], 1e-5, + "%s, lmsder, x%u", name, i); + } + + + { + double s2 = pow(gsl_blas_dnrm2 (s->f), 2.0); + + gsl_test_rel (s2, f_sumsq, 1e-5, "%s, lmsder, |f|^2", name); + + for (i = 0; i < p; i++) + { + double ei = sqrt(s2/(n-p))*sqrt(gsl_matrix_get(covar,i,i)); + gsl_test_rel (ei, sigma[i], 1e-4, + "%s, sigma(%d)", name, i) ; + } + } + + gsl_matrix_free (covar); + } + + gsl_multifit_fdfsolver_free (s); +} diff --git a/gsl-1.9/multifit/test_brown.c b/gsl-1.9/multifit/test_brown.c new file mode 100644 index 0000000..e1ac937 --- /dev/null +++ b/gsl-1.9/multifit/test_brown.c @@ -0,0 +1,113 @@ +const size_t brown_N = 20; +const size_t brown_P = 4; + +double brown_X[20][4] = { + {24.3485677, 4.71448798, -2.19486633, 2.69405755}, + {22.4116222, 3.93075538, -1.42344852, 2.5233557}, + {17.88886, 2.9290853, .125174936, -3.96823353}, + {17.3237176, 2.99606803, 2.03285653, 2.28992327}, + {17.0906508, 3.02485425, .296995153, .0876226126}, + {16.578006, 3.1036312, -.18617941, .103262914}, + {15.692993, 3.33088442, .0706406887, 1.05923955}, + {14.3232177, 3.85604218, -2.3762839, -3.09486813}, + {14.1279266, 3.97896121, .446109351, 1.40023753}, + {13.6081961, 4.16435075, -1.51250057, -1.52510626}, + {13.4295245, 4.22697223, -.196985195, .532009293}, + {13.0176117, 4.3579261, -.353131208, .301377627}, + {12.2713535, 4.62398535, -.00183585584, .894170703}, + {11.0316144, 5.13967727, -2.38978772, -2.89510064}, + {10.8807981, 5.24558004, .230495952, 1.27315117}, + {10.4029264, 5.41141257, -1.5116632, -1.47615921}, + {10.2574435, 5.46211045, -.299855732, .451893162}, + {9.87863876, 5.57914292, -.368885288, .358086545}, + {9.1894983, 5.82082741, -.230157969, .621476534}, + {8.00589008, 6.27788753, -1.46022815, -1.33468082} +}; + +double brown_F[20] = { + 2474.05541, + 1924.69004, + 1280.63194, + 1244.81867, + 1190.53739, + 1159.34935, + 1108.44426, + 1090.11073, + 1015.92942, + 1002.43533, + 971.221084, + 949.589435, + 911.359899, + 906.522994, + 840.525729, + 833.950164, + 807.557511, + 791.00924, + 761.09598, + 726.787783, +}; + +double brown_cov[4][4] = { + { 1.8893186910e-01, -4.7099989571e-02, 5.2154168404e-01, 1.6608168209e-02}, + {-4.7099989571e-02, 1.1761534388e-02, -1.2987843074e-01, -4.1615942391e-03}, + { 5.2154168404e-01, -1.2987843074e-01, 1.4653936514e+00, 1.5738321686e-02}, + { 1.6608168209e-02, -4.1615942391e-03, 1.5738321686e-02, 4.2348042340e-02}, +}; + +double brown_x0[4] = { 25, 5, -5, -1 }; + +int +brown_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + size_t i; + + for (i = 0; i < 20; i++) + { + double ti = 0.2 * (i + 1); + double ui = x0 + x1 * ti - exp (ti); + double vi = x2 + x3 * sin (ti) - cos (ti); + + gsl_vector_set (f, i, ui * ui + vi * vi); + } + + return GSL_SUCCESS; +} + +int +brown_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double x0 = gsl_vector_get (x, 0); + double x1 = gsl_vector_get (x, 1); + double x2 = gsl_vector_get (x, 2); + double x3 = gsl_vector_get (x, 3); + size_t i; + + for (i = 0; i < 20; i++) + { + double ti = 0.2 * (i + 1); + double ui = x0 + x1 * ti - exp (ti); + double vi = x2 + x3 * sin (ti) - cos (ti); + + gsl_matrix_set (df, i, 0, 2 * ui); + gsl_matrix_set (df, i, 1, 2 * ui * ti); + gsl_matrix_set (df, i, 2, 2 * vi); + gsl_matrix_set (df, i, 3, 2 * vi * sin (ti)); + + } + return GSL_SUCCESS; +} + +int +brown_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + brown_f (x, params, f); + brown_df (x, params, df); + + return GSL_SUCCESS; +} + diff --git a/gsl-1.9/multifit/test_enso.c b/gsl-1.9/multifit/test_enso.c new file mode 100644 index 0000000..6343738 --- /dev/null +++ b/gsl-1.9/multifit/test_enso.c @@ -0,0 +1,279 @@ +const size_t enso_N = 168; +const size_t enso_P = 9; + +double enso_x0[9] = { 10.0, 3.0, 0.5, 44.0, -1.5, 0.5, 26.0, 0.1, 1.5 }; + +double enso_x[9] = { + 1.0510749193E+01, + 3.0762128085E+00, + 5.3280138227E-01, + 4.4311088700E+01, + -1.6231428586E+00, + 5.2554493756E-01, + 2.6887614440E+01, + 2.1232288488E-01, + 1.4966870418E+00 +}; + +double enso_sumsq = 7.8853978668E+02; + +double enso_sigma[9] = { + 1.7488832467E-01, + 2.4310052139E-01, + 2.4354686618E-01, + 9.4408025976E-01, + 2.8078369611E-01, + 4.8073701119E-01, + 4.1612939130E-01, + 5.1460022911E-01, + 2.5434468893E-01 +}; + +double enso_F[168] = { + 12.90000, + 11.30000, + 10.60000, + 11.20000, + 10.90000, + 7.500000, + 7.700000, + 11.70000, + 12.90000, + 14.30000, + 10.90000, + 13.70000, + 17.10000, + 14.00000, + 15.30000, + 8.500000, + 5.700000, + 5.500000, + 7.600000, + 8.600000, + 7.300000, + 7.600000, + 12.70000, + 11.00000, + 12.70000, + 12.90000, + 13.00000, + 10.90000, + 10.400000, + 10.200000, + 8.000000, + 10.90000, + 13.60000, + 10.500000, + 9.200000, + 12.40000, + 12.70000, + 13.30000, + 10.100000, + 7.800000, + 4.800000, + 3.000000, + 2.500000, + 6.300000, + 9.700000, + 11.60000, + 8.600000, + 12.40000, + 10.500000, + 13.30000, + 10.400000, + 8.100000, + 3.700000, + 10.70000, + 5.100000, + 10.400000, + 10.90000, + 11.70000, + 11.40000, + 13.70000, + 14.10000, + 14.00000, + 12.50000, + 6.300000, + 9.600000, + 11.70000, + 5.000000, + 10.80000, + 12.70000, + 10.80000, + 11.80000, + 12.60000, + 15.70000, + 12.60000, + 14.80000, + 7.800000, + 7.100000, + 11.20000, + 8.100000, + 6.400000, + 5.200000, + 12.00000, + 10.200000, + 12.70000, + 10.200000, + 14.70000, + 12.20000, + 7.100000, + 5.700000, + 6.700000, + 3.900000, + 8.500000, + 8.300000, + 10.80000, + 16.70000, + 12.60000, + 12.50000, + 12.50000, + 9.800000, + 7.200000, + 4.100000, + 10.60000, + 10.100000, + 10.100000, + 11.90000, + 13.60000, + 16.30000, + 17.60000, + 15.50000, + 16.00000, + 15.20000, + 11.20000, + 14.30000, + 14.50000, + 8.500000, + 12.00000, + 12.70000, + 11.30000, + 14.50000, + 15.10000, + 10.400000, + 11.50000, + 13.40000, + 7.500000, + 0.6000000, + 0.3000000, + 5.500000, + 5.000000, + 4.600000, + 8.200000, + 9.900000, + 9.200000, + 12.50000, + 10.90000, + 9.900000, + 8.900000, + 7.600000, + 9.500000, + 8.400000, + 10.70000, + 13.60000, + 13.70000, + 13.70000, + 16.50000, + 16.80000, + 17.10000, + 15.40000, + 9.500000, + 6.100000, + 10.100000, + 9.300000, + 5.300000, + 11.20000, + 16.60000, + 15.60000, + 12.00000, + 11.50000, + 8.600000, + 13.80000, + 8.700000, + 8.600000, + 8.600000, + 8.700000, + 12.80000, + 13.20000, + 14.00000, + 13.40000, + 14.80000 +}; + + +int +enso_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double b[9]; + size_t i; + + for (i = 0; i < 9; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 168; i++) + { + double t = (i + 1.0); + double y; + y = b[0]; + y += b[1] * cos(2*M_PI*t/12); + y += b[2] * sin(2*M_PI*t/12); + y += b[4] * cos(2*M_PI*t/b[3]); + y += b[5] * sin(2*M_PI*t/b[3]); + y += b[7] * cos(2*M_PI*t/b[6]); + y += b[8] * sin(2*M_PI*t/b[6]); + + gsl_vector_set (f, i, enso_F[i] - y); + } + + return GSL_SUCCESS; +} + +int +enso_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double b[9]; + size_t i; + + for (i = 0; i < 9; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 168; i++) + { + double t = (i + 1.0); + + gsl_matrix_set (df, i, 0, -1.0); + gsl_matrix_set (df, i, 1, -cos(2*M_PI*t/12)); + gsl_matrix_set (df, i, 2, -sin(2*M_PI*t/12)); + gsl_matrix_set (df, i, 3, + -b[4]*(2*M_PI*t/(b[3]*b[3]))*sin(2*M_PI*t/b[3]) + +b[5]*(2*M_PI*t/(b[3]*b[3]))*cos(2*M_PI*t/b[3])); + gsl_matrix_set (df, i, 4, -cos(2*M_PI*t/b[3])); + gsl_matrix_set (df, i, 5, -sin(2*M_PI*t/b[3])); + gsl_matrix_set (df, i, 6, + -b[7] * (2*M_PI*t/(b[6]*b[6])) * sin(2*M_PI*t/b[6]) + +b[8] * (2*M_PI*t/(b[6]*b[6])) * cos(2*M_PI*t/b[6])); + gsl_matrix_set (df, i, 7, -cos(2*M_PI*t/b[6])); + gsl_matrix_set (df, i, 8, -sin(2*M_PI*t/b[6])); + } + + return GSL_SUCCESS; +} + +int +enso_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + enso_f (x, params, f); + enso_df (x, params, df); + + return GSL_SUCCESS; +} + + + + + diff --git a/gsl-1.9/multifit/test_estimator.c b/gsl-1.9/multifit/test_estimator.c new file mode 100644 index 0000000..51cbaba --- /dev/null +++ b/gsl-1.9/multifit/test_estimator.c @@ -0,0 +1,37 @@ +void +test_estimator () +{ + gsl_vector_view c; + gsl_matrix_view cov; + gsl_vector_view x; + double y, y_err; + + double cov_ij[25] = { + 4.271520, -0.526675, 0.957930, 0.267750, -0.103610, + -0.526675, 5.701680, -0.098080, 0.641845, 0.429780, + 0.957930, -0.098080, 4.584790, 0.375865, 1.510810, + 0.267750, 0.641845, 0.375865, 4.422720, 0.392210, + -0.103610, 0.429780, 1.510810, 0.392210, 5.782750 + + }; + + double c_i[5] = { + -0.627020, 0.848674, 0.216877, -0.057883, 0.596668 + }; + + double x_i[5] = { + 0.99932, 0.23858, 0.19797, 1.44008, -0.15335 + }; + + double y_expected = -5.56037032230000e-01; + double yerr_expected = 3.91891123349318e+00; + + cov = gsl_matrix_view_array(cov_ij, 5, 5); + c = gsl_vector_view_array(c_i, 5); + x = gsl_vector_view_array(x_i, 5); + + gsl_multifit_linear_est(&x.vector , &c.vector, &cov.matrix, &y, &y_err); + + gsl_test_rel (y, y_expected, 256*GSL_DBL_EPSILON, "gsl_multifit_linear_est y"); + gsl_test_rel (y_err, yerr_expected, 256*GSL_DBL_EPSILON, "gsl_multifit_linear_est yerr"); +} diff --git a/gsl-1.9/multifit/test_filip.c b/gsl-1.9/multifit/test_filip.c new file mode 100644 index 0000000..bdb4609 --- /dev/null +++ b/gsl-1.9/multifit/test_filip.c @@ -0,0 +1,202 @@ +size_t filip_n = 82; +size_t filip_p = 11; + +double filip_x[] = { -6.860120914, -4.324130045, -4.358625055, +-4.358426747, -6.955852379, -6.661145254, -6.355462942, -6.118102026, +-7.115148017, -6.815308569, -6.519993057, -6.204119983, -5.853871964, +-6.109523091, -5.79832982, -5.482672118, -5.171791386, -4.851705903, +-4.517126416, -4.143573228, -3.709075441, -3.499489089, -6.300769497, +-5.953504836, -5.642065153, -5.031376979, -4.680685696, -4.329846955, +-3.928486195, -8.56735134, -8.363211311, -8.107682739, -7.823908741, +-7.522878745, -7.218819279, -6.920818754, -6.628932138, -6.323946875, +-5.991399828, -8.781464495, -8.663140179, -8.473531488, -8.247337057, +-7.971428747, -7.676129393, -7.352812702, -7.072065318, -6.774174009, +-6.478861916, -6.159517513, -6.835647144, -6.53165267, -6.224098421, +-5.910094889, -5.598599459, -5.290645224, -4.974284616, -4.64454848, +-4.290560426, -3.885055584, -3.408378962, -3.13200249, -8.726767166, +-8.66695597, -8.511026475, -8.165388579, -7.886056648, -7.588043762, +-7.283412422, -6.995678626, -6.691862621, -6.392544977, -6.067374056, +-6.684029655, -6.378719832, -6.065855188, -5.752272167, -5.132414673, +-4.811352704, -4.098269308, -3.66174277, -3.2644011}; + +double filip_y[] = { 0.8116, 0.9072, 0.9052, 0.9039, 0.8053, 0.8377, +0.8667, 0.8809, 0.7975, 0.8162, 0.8515, 0.8766, 0.8885, 0.8859, +0.8959, 0.8913, 0.8959, 0.8971, 0.9021, 0.909, 0.9139, 0.9199, 0.8692, +0.8872, 0.89, 0.891, 0.8977, 0.9035, 0.9078, 0.7675, 0.7705, 0.7713, +0.7736, 0.7775, 0.7841, 0.7971, 0.8329, 0.8641, 0.8804, 0.7668, +0.7633, 0.7678, 0.7697, 0.77, 0.7749, 0.7796, 0.7897, 0.8131, 0.8498, +0.8741, 0.8061, 0.846, 0.8751, 0.8856, 0.8919, 0.8934, 0.894, 0.8957, +0.9047, 0.9129, 0.9209, 0.9219, 0.7739, 0.7681, 0.7665, 0.7703, +0.7702, 0.7761, 0.7809, 0.7961, 0.8253, 0.8602, 0.8809, 0.8301, +0.8664, 0.8834, 0.8898, 0.8964, 0.8963, 0.9074, 0.9119, 0.9228 } ; + + +void +test_filip () +{ + size_t i, j; + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (filip_n, filip_p); + + gsl_matrix * X = gsl_matrix_alloc (filip_n, filip_p); + gsl_vector_view y = gsl_vector_view_array (filip_y, filip_n); + gsl_vector * c = gsl_vector_alloc (filip_p); + gsl_matrix * cov = gsl_matrix_alloc (filip_p, filip_p); + gsl_vector_view diag; + + double chisq; + + double expected_c[11] = { -1467.48961422980, + -2772.17959193342, + -2316.37108160893, + -1127.97394098372, + -354.478233703349, + -75.1242017393757, + -10.8753180355343, + -1.06221498588947, + -0.670191154593408E-01, + -0.246781078275479E-02, + -0.402962525080404E-04 }; + + double expected_sd[11] = { 298.084530995537, + 559.779865474950, + 466.477572127796, + 227.204274477751, + 71.6478660875927, + 15.2897178747400, + 2.23691159816033, + 0.221624321934227, + 0.142363763154724E-01, + 0.535617408889821E-03, + 0.896632837373868E-05 }; + + double expected_chisq = 0.795851382172941E-03; + + for (i = 0 ; i < filip_n; i++) + { + for (j = 0; j < filip_p; j++) + { + gsl_matrix_set(X, i, j, pow(filip_x[i], j)); + } + } + + gsl_multifit_linear (X, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-7, "filip gsl_fit_multilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-7, "filip gsl_fit_multilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-7, "filip gsl_fit_multilinear c2") ; + gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-7, "filip gsl_fit_multilinear c3") ; + gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-7, "filip gsl_fit_multilinear c4") ; + gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-7, "filip gsl_fit_multilinear c5") ; + gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-7, "filip gsl_fit_multilinear c6") ; + gsl_test_rel (gsl_vector_get(c,7), expected_c[7], 1e-7, "filip gsl_fit_multilinear c7") ; + gsl_test_rel (gsl_vector_get(c,8), expected_c[8], 1e-7, "filip gsl_fit_multilinear c8") ; + gsl_test_rel (gsl_vector_get(c,9), expected_c[9], 1e-7, "filip gsl_fit_multilinear c9") ; + gsl_test_rel (gsl_vector_get(c,10), expected_c[10], 1e-7, "filip gsl_fit_multilinear c10") ; + + diag = gsl_matrix_diagonal (cov); + + gsl_test_rel (gsl_vector_get(&diag.vector,0), pow(expected_sd[0],2.0), 1e-6, "filip gsl_fit_multilinear cov00") ; + gsl_test_rel (gsl_vector_get(&diag.vector,1), pow(expected_sd[1],2.0), 1e-6, "filip gsl_fit_multilinear cov11") ; + gsl_test_rel (gsl_vector_get(&diag.vector,2), pow(expected_sd[2],2.0), 1e-6, "filip gsl_fit_multilinear cov22") ; + gsl_test_rel (gsl_vector_get(&diag.vector,3), pow(expected_sd[3],2.0), 1e-6, "filip gsl_fit_multilinear cov33") ; + gsl_test_rel (gsl_vector_get(&diag.vector,4), pow(expected_sd[4],2.0), 1e-6, "filip gsl_fit_multilinear cov44") ; + gsl_test_rel (gsl_vector_get(&diag.vector,5), pow(expected_sd[5],2.0), 1e-6, "filip gsl_fit_multilinear cov55") ; + gsl_test_rel (gsl_vector_get(&diag.vector,6), pow(expected_sd[6],2.0), 1e-6, "filip gsl_fit_multilinear cov66") ; + gsl_test_rel (gsl_vector_get(&diag.vector,7), pow(expected_sd[7],2.0), 1e-6, "filip gsl_fit_multilinear cov77") ; + gsl_test_rel (gsl_vector_get(&diag.vector,8), pow(expected_sd[8],2.0), 1e-6, "filip gsl_fit_multilinear cov88") ; + gsl_test_rel (gsl_vector_get(&diag.vector,9), pow(expected_sd[9],2.0), 1e-6, "filip gsl_fit_multilinear cov99") ; + gsl_test_rel (gsl_vector_get(&diag.vector,10), pow(expected_sd[10],2.0), 1e-6, "filip gsl_fit_multilinear cov1010") ; + + gsl_test_rel (chisq, expected_chisq, 1e-7, "filip gsl_fit_multilinear chisq") ; + + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_matrix_free(X); + gsl_multifit_linear_free (work); + } + + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (filip_n, filip_p); + + gsl_matrix * X = gsl_matrix_alloc (filip_n, filip_p); + gsl_vector_view y = gsl_vector_view_array (filip_y, filip_n); + gsl_vector * w = gsl_vector_alloc (filip_n); + gsl_vector * c = gsl_vector_alloc (filip_p); + gsl_matrix * cov = gsl_matrix_alloc (filip_p, filip_p); + + double chisq; + + double expected_c[11] = { -1467.48961422980, + -2772.17959193342, + -2316.37108160893, + -1127.97394098372, + -354.478233703349, + -75.1242017393757, + -10.8753180355343, + -1.06221498588947, + -0.670191154593408E-01, + -0.246781078275479E-02, + -0.402962525080404E-04 }; + + /* computed using GNU Calc */ + + double expected_cov[11][11] ={ { 7.9269341767252183262588583867942e9, 1.4880416622254098343441063389706e10, 1.2385811858111487905481427591107e10, 6.0210784406215266653697715794241e9, 1.8936652526181982747116667336389e9, 4.0274900618493109653998118587093e8, 5.8685468011819735806180092394606e7, 5.7873451475721689084330083708901e6, 3.6982719848703747920663262917032e5, 1.3834818802741350637527054170891e4, 2.301758578713219280719633494302e2 }, + { 1.4880416622254098334697515488559e10, 2.7955091668548290835529555438088e10, 2.3286604504243362691678565997033e10, 1.132895006796272983689297219686e10, 3.5657281653312473123348357644683e9, 7.5893300392314445528176646366087e8, 1.1066654886143524811964131660002e8, 1.0921285448484575110763947787775e7, 6.9838139975394769253353547606971e5, 2.6143091775349597218939272614126e4, 4.3523386330348588614289505633539e2 }, + { 1.2385811858111487890788272968677e10, 2.3286604504243362677757802422747e10, 1.9412787917766676553608636489674e10, 9.4516246492862131849077729250098e9, 2.9771226694709917550143152097252e9, 6.3413035086730038062129508949859e8, 9.2536164488309401636559552742339e7, 9.1386304643423333815338760248027e6, 5.8479478338916429826337004060941e5, 2.1905933113294737443808429764554e4, 3.6493161325305557266196635180155e2 }, + { 6.0210784406215266545770691532365e9, 1.1328950067962729823273441573365e10, 9.4516246492862131792040001429636e9, 4.6053152992000107509329772255094e9, 1.4517147860312147098138030287038e9, 3.0944988323328589376402579060072e8, 4.5190223822292688669369522708712e7, 4.4660958693678497534529855690752e6, 2.8599340736122198213681258676423e5, 1.0720394998549386596165641244705e4, 1.7870937745661967319298031044424e2 }, + { 1.8936652526181982701620450132636e9, 3.5657281653312473058825073094524e9, 2.9771226694709917514149924058297e9, 1.451714786031214708936087401632e9, 4.5796563896564815123074920050827e8, 9.7693972414561515534525103622773e7, 1.427717861635658545863942948444e7, 1.4120161287735817621354292900338e6, 9.0484361228623960006818614875557e4, 3.394106783764852373199087455398e3, 5.6617406468519495376287407526295e1 }, + { 4.0274900618493109532650887473599e8, 7.589330039231444534478894935778e8, 6.3413035086730037947153564986653e8, 3.09449883233285893390542947998e8, 9.7693972414561515475770399055121e7, 2.0855726248311948992114244257719e7, 3.0501263034740400533872858749566e6, 3.0187475839310308153394428784224e5, 1.9358204633534233524477930175632e4, 7.2662989867560017077361942813911e2, 1.2129002231061036467607394277965e1 }, + { 5.868546801181973559370854830868e7, 1.1066654886143524778548044386795e8, 9.2536164488309401413296494869777e7, 4.5190223822292688587853853162072e7, 1.4277178616356585441556046753562e7, 3.050126303474040051574715592746e6, 4.4639982579046340884744460329946e5, 4.4212093985989836047285007760238e4, 2.8371395028774486687625333589972e3, 1.0656694507620102300567296504381e2, 1.7799982046359973175080475654123e0 }, + { 5.7873451475721688839974153925406e6, 1.0921285448484575071271480643397e7, 9.1386304643423333540728480344578e6, 4.4660958693678497427674903565664e6, 1.4120161287735817596182229182587e6, 3.0187475839310308117812257613082e5, 4.4212093985989836021482392757677e4, 4.3818874017028389517560906916315e3, 2.813828775753142855163154605027e2, 1.0576188138416671883232607188969e1, 1.7676976288918295012452853715408e-1 }, + { 3.6982719848703747742568351456818e5, 6.9838139975394768959780068745979e5, 5.8479478338916429616547638954781e5, 2.8599340736122198128717796825489e5, 9.0484361228623959793493985226792e4, 1.9358204633534233490579641064343e4, 2.8371395028774486654873647731797e3, 2.8138287757531428535592907878017e2, 1.8081118503579798222896804627964e1, 6.8005074291434681866415478598732e-1, 1.1373581557749643543869665860719e-2 }, + { 1.3834818802741350562839757244708e4, 2.614309177534959709397445440919e4, 2.1905933113294737352721470167247e4, 1.0720394998549386558251721913182e4, 3.3941067837648523632905604575131e3, 7.2662989867560016909534954790835e2, 1.0656694507620102282337905013451e2, 1.0576188138416671871337685672492e1, 6.8005074291434681828743281967838e-1, 2.5593857187900736057022477529078e-2, 4.2831487599116264442963102045936e-4 }, + { 2.3017585787132192669801658674163e2, 4.3523386330348588381716460685124e2, 3.6493161325305557094116270974735e2, 1.7870937745661967246233792737255e2, 5.6617406468519495180024059284629e1, 1.2129002231061036433003571679329e1, 1.7799982046359973135014027410646e0, 1.7676976288918294983059118597214e-1, 1.137358155774964353146460100337e-2, 4.283148759911626442000316269063e-4, 7.172253875245080423800933453952e-6 } }; + + double expected_chisq = 0.795851382172941E-03; + + for (i = 0 ; i < filip_n; i++) + { + for (j = 0; j < filip_p; j++) + { + gsl_matrix_set(X, i, j, pow(filip_x[i], j)); + } + } + + gsl_vector_set_all (w, 1.0); + + gsl_multifit_wlinear (X, w, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-7, "filip gsl_fit_multilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-7, "filip gsl_fit_multilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-7, "filip gsl_fit_multilinear c2") ; + gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-7, "filip gsl_fit_multilinear c3") ; + gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-7, "filip gsl_fit_multilinear c4") ; + gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-7, "filip gsl_fit_multilinear c5") ; + gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-7, "filip gsl_fit_multilinear c6") ; + gsl_test_rel (gsl_vector_get(c,7), expected_c[7], 1e-7, "filip gsl_fit_multilinear c7") ; + gsl_test_rel (gsl_vector_get(c,8), expected_c[8], 1e-7, "filip gsl_fit_multilinear c8") ; + gsl_test_rel (gsl_vector_get(c,9), expected_c[9], 1e-7, "filip gsl_fit_multilinear c9") ; + gsl_test_rel (gsl_vector_get(c,10), expected_c[10], 1e-7, "filip gsl_fit_multilinear c10") ; + + + for (i = 0; i < filip_p; i++) + { + for (j = 0; j < filip_p; j++) + { + gsl_test_rel (gsl_matrix_get(cov,i,j), expected_cov[i][j], 1e-6, + "filip gsl_fit_wmultilinear cov(%d,%d)", i, j) ; + } + } + + gsl_test_rel (chisq, expected_chisq, 1e-7, "filip gsl_fit_multilinear chisq") ; + + gsl_vector_free(w); + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_matrix_free(X); + gsl_multifit_linear_free (work); + } +} diff --git a/gsl-1.9/multifit/test_fn.c b/gsl-1.9/multifit/test_fn.c new file mode 100644 index 0000000..ad190a1 --- /dev/null +++ b/gsl-1.9/multifit/test_fn.c @@ -0,0 +1,25 @@ +gsl_multifit_function_fdf +make_fdf (int (* f) (const gsl_vector *, void *, gsl_vector *), + int (* df) (const gsl_vector *, void *, gsl_matrix *), + int (* fdf) (const gsl_vector *, void *, gsl_vector *, gsl_matrix *), + size_t n, + size_t p, + void * params); + +gsl_multifit_function_fdf +make_fdf (int (* f) (const gsl_vector *, void *, gsl_vector *), + int (* df) (const gsl_vector *, void *, gsl_matrix *), + int (* fdf) (const gsl_vector *, void *, gsl_vector *, gsl_matrix *), + size_t n, + size_t p, + void * params) +{ + gsl_multifit_function_fdf F_new; + F_new.f = f; + F_new.df = df; + F_new.fdf = fdf; + F_new.n = n; + F_new.p = p; + F_new.params = params; + return F_new; +} diff --git a/gsl-1.9/multifit/test_hahn1.c b/gsl-1.9/multifit/test_hahn1.c new file mode 100644 index 0000000..8a2385a --- /dev/null +++ b/gsl-1.9/multifit/test_hahn1.c @@ -0,0 +1,568 @@ +const size_t hahn1_N = 236; +const size_t hahn1_P = 7; + +/* double hahn1_x0[7] = { 10, -1, 0.05, -0.00001, -0.05, 0.001, -0.000001 }; */ + +double hahn1_x0[7] = { 1, -0.1, 0.005, -0.000001, -0.005, 0.0001, -0.0000001}; + +double hahn1_x[7] = { +1.0776351733E+00, +-1.2269296921E-01, +4.0863750610E-03, +-1.4262662514E-06, +-5.7609940901E-03, +2.4053735503E-04, +-1.2314450199E-07 +}; + +double hahn1_sumsq = 1.5324382854E+00; + +double hahn1_sigma[7] = { + 1.7070154742E-01, + 1.2000289189E-02, + 2.2508314937E-04, + 2.7578037666E-07, + 2.4712888219E-04, + 1.0449373768E-05, + 1.3027335327E-08 +}; + +double hahn1_F1[236] = { + .591E0, + 1.547E0, + 2.902E0, + 2.894E0, + 4.703E0, + 6.307E0, + 7.03E0 , + 7.898E0, + 9.470E0, + 9.484E0, + 10.072E0, + 10.163E0, + 11.615E0, + 12.005E0, + 12.478E0, + 12.982E0, + 12.970E0, + 13.926E0, + 14.452E0, + 14.404E0, + 15.190E0, + 15.550E0, + 15.528E0, + 15.499E0, + 16.131E0, + 16.438E0, + 16.387E0, + 16.549E0, + 16.872E0, + 16.830E0, + 16.926E0, + 16.907E0, + 16.966E0, + 17.060E0, + 17.122E0, + 17.311E0, + 17.355E0, + 17.668E0, + 17.767E0, + 17.803E0, + 17.765E0, + 17.768E0, + 17.736E0, + 17.858E0, + 17.877E0, + 17.912E0, + 18.046E0, + 18.085E0, + 18.291E0, + 18.357E0, + 18.426E0, + 18.584E0, + 18.610E0, + 18.870E0, + 18.795E0, + 19.111E0, + .367E0, + .796E0, + 0.892E0, + 1.903E0, + 2.150E0, + 3.697E0, + 5.870E0, + 6.421E0, + 7.422E0, + 9.944E0, + 11.023E0, + 11.87E0 , + 12.786E0, + 14.067E0, + 13.974E0, + 14.462E0, + 14.464E0, + 15.381E0, + 15.483E0, + 15.59E0 , + 16.075E0, + 16.347E0, + 16.181E0, + 16.915E0, + 17.003E0, + 16.978E0, + 17.756E0, + 17.808E0, + 17.868E0, + 18.481E0, + 18.486E0, + 19.090E0, + 16.062E0, + 16.337E0, + 16.345E0, + 16.388E0, + 17.159E0, + 17.116E0, + 17.164E0, + 17.123E0, + 17.979E0, + 17.974E0, + 18.007E0, + 17.993E0, + 18.523E0, + 18.669E0, + 18.617E0, + 19.371E0, + 19.330E0, + 0.080E0, + 0.248E0, + 1.089E0, + 1.418E0, + 2.278E0, + 3.624E0, + 4.574E0, + 5.556E0, + 7.267E0, + 7.695E0, + 9.136E0, + 9.959E0, + 9.957E0, + 11.600E0, + 13.138E0, + 13.564E0, + 13.871E0, + 13.994E0, + 14.947E0, + 15.473E0, + 15.379E0, + 15.455E0, + 15.908E0, + 16.114E0, + 17.071E0, + 17.135E0, + 17.282E0, + 17.368E0, + 17.483E0, + 17.764E0, + 18.185E0, + 18.271E0, + 18.236E0, + 18.237E0, + 18.523E0, + 18.627E0, + 18.665E0, + 19.086E0, + 0.214E0, + 0.943E0, + 1.429E0, + 2.241E0, + 2.951E0, + 3.782E0, + 4.757E0, + 5.602E0, + 7.169E0, + 8.920E0, + 10.055E0, + 12.035E0, + 12.861E0, + 13.436E0, + 14.167E0, + 14.755E0, + 15.168E0, + 15.651E0, + 15.746E0, + 16.216E0, + 16.445E0, + 16.965E0, + 17.121E0, + 17.206E0, + 17.250E0, + 17.339E0, + 17.793E0, + 18.123E0, + 18.49E0 , + 18.566E0, + 18.645E0, + 18.706E0, + 18.924E0, + 19.1E0 , + 0.375E0, + 0.471E0, + 1.504E0, + 2.204E0, + 2.813E0, + 4.765E0, + 9.835E0, + 10.040E0, + 11.946E0, + 12.596E0, + 13.303E0, + 13.922E0, + 14.440E0, + 14.951E0, + 15.627E0, + 15.639E0, + 15.814E0, + 16.315E0, + 16.334E0, + 16.430E0, + 16.423E0, + 17.024E0, + 17.009E0, + 17.165E0, + 17.134E0, + 17.349E0, + 17.576E0, + 17.848E0, + 18.090E0, + 18.276E0, + 18.404E0, + 18.519E0, + 19.133E0, + 19.074E0, + 19.239E0, + 19.280E0, + 19.101E0, + 19.398E0, + 19.252E0, + 19.89E0 , + 20.007E0, + 19.929E0, + 19.268E0, + 19.324E0, + 20.049E0, + 20.107E0, + 20.062E0, + 20.065E0, + 19.286E0, + 19.972E0, + 20.088E0, + 20.743E0, + 20.83E0 , + 20.935E0, + 21.035E0, + 20.93E0 , + 21.074E0, + 21.085E0, + 20.935E0 +}; + + +double hahn1_F0[236] = { + 24.41E0, + 34.82E0, + 44.09E0, + 45.07E0, + 54.98E0, + 65.51E0, + 70.53E0, + 75.70E0, + 89.57E0, + 91.14E0, + 96.40E0, + 97.19E0, + 114.26E0, + 120.25E0, + 127.08E0, + 133.55E0, + 133.61E0, + 158.67E0, + 172.74E0, + 171.31E0, + 202.14E0, + 220.55E0, + 221.05E0, + 221.39E0, + 250.99E0, + 268.99E0, + 271.80E0, + 271.97E0, + 321.31E0, + 321.69E0, + 330.14E0, + 333.03E0, + 333.47E0, + 340.77E0, + 345.65E0, + 373.11E0, + 373.79E0, + 411.82E0, + 419.51E0, + 421.59E0, + 422.02E0, + 422.47E0, + 422.61E0, + 441.75E0, + 447.41E0, + 448.7E0 , + 472.89E0, + 476.69E0, + 522.47E0, + 522.62E0, + 524.43E0, + 546.75E0, + 549.53E0, + 575.29E0, + 576.00E0, + 625.55E0, + 20.15E0, + 28.78E0, + 29.57E0, + 37.41E0, + 39.12E0, + 50.24E0, + 61.38E0, + 66.25E0, + 73.42E0, + 95.52E0, + 107.32E0, + 122.04E0, + 134.03E0, + 163.19E0, + 163.48E0, + 175.70E0, + 179.86E0, + 211.27E0, + 217.78E0, + 219.14E0, + 262.52E0, + 268.01E0, + 268.62E0, + 336.25E0, + 337.23E0, + 339.33E0, + 427.38E0, + 428.58E0, + 432.68E0, + 528.99E0, + 531.08E0, + 628.34E0, + 253.24E0, + 273.13E0, + 273.66E0, + 282.10E0, + 346.62E0, + 347.19E0, + 348.78E0, + 351.18E0, + 450.10E0, + 450.35E0, + 451.92E0, + 455.56E0, + 552.22E0, + 553.56E0, + 555.74E0, + 652.59E0, + 656.20E0, + 14.13E0, + 20.41E0, + 31.30E0, + 33.84E0, + 39.70E0, + 48.83E0, + 54.50E0, + 60.41E0, + 72.77E0, + 75.25E0, + 86.84E0, + 94.88E0, + 96.40E0, + 117.37E0, + 139.08E0, + 147.73E0, + 158.63E0, + 161.84E0, + 192.11E0, + 206.76E0, + 209.07E0, + 213.32E0, + 226.44E0, + 237.12E0, + 330.90E0, + 358.72E0, + 370.77E0, + 372.72E0, + 396.24E0, + 416.59E0, + 484.02E0, + 495.47E0, + 514.78E0, + 515.65E0, + 519.47E0, + 544.47E0, + 560.11E0, + 620.77E0, + 18.97E0, + 28.93E0, + 33.91E0, + 40.03E0, + 44.66E0, + 49.87E0, + 55.16E0, + 60.90E0, + 72.08E0, + 85.15E0, + 97.06E0, + 119.63E0, + 133.27E0, + 143.84E0, + 161.91E0, + 180.67E0, + 198.44E0, + 226.86E0, + 229.65E0, + 258.27E0, + 273.77E0, + 339.15E0, + 350.13E0, + 362.75E0, + 371.03E0, + 393.32E0, + 448.53E0, + 473.78E0, + 511.12E0, + 524.70E0, + 548.75E0, + 551.64E0, + 574.02E0, + 623.86E0, + 21.46E0, + 24.33E0, + 33.43E0, + 39.22E0, + 44.18E0, + 55.02E0, + 94.33E0, + 96.44E0, + 118.82E0, + 128.48E0, + 141.94E0, + 156.92E0, + 171.65E0, + 190.00E0, + 223.26E0, + 223.88E0, + 231.50E0, + 265.05E0, + 269.44E0, + 271.78E0, + 273.46E0, + 334.61E0, + 339.79E0, + 349.52E0, + 358.18E0, + 377.98E0, + 394.77E0, + 429.66E0, + 468.22E0, + 487.27E0, + 519.54E0, + 523.03E0, + 612.99E0, + 638.59E0, + 641.36E0, + 622.05E0, + 631.50E0, + 663.97E0, + 646.9E0 , + 748.29E0, + 749.21E0, + 750.14E0, + 647.04E0, + 646.89E0, + 746.9E0 , + 748.43E0, + 747.35E0, + 749.27E0, + 647.61E0, + 747.78E0, + 750.51E0, + 851.37E0, + 845.97E0, + 847.54E0, + 849.93E0, + 851.61E0, + 849.75E0, + 850.98E0, + 848.23E0 +}; + + +int +hahn1_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double b[7]; + size_t i; + + for (i = 0; i < 7; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 236; i++) + { + double x = hahn1_F0[i]; + double y = ((b[0] + x* (b[1] + x * (b[2] + x * b[3]))) + / (1 + x*(b[4] + x *(b[5] + x*b[6])))); + gsl_vector_set (f, i, hahn1_F1[i] - y); + } + + return GSL_SUCCESS; +} + +int +hahn1_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double b[7]; + size_t i; + + for (i = 0; i < 7; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 236; i++) + { + double x = hahn1_F0[i]; + double u = (b[0] + x*(b[1] + x*(b[2] + x * b[3]))); + double v = (1 + x*(b[4] + x*(b[5] + x*b[6]))); + gsl_matrix_set (df, i, 0, -1/v); + gsl_matrix_set (df, i, 1, -x/v); + gsl_matrix_set (df, i, 2, -x*x/v); + gsl_matrix_set (df, i, 3, -x*x*x/v); + gsl_matrix_set (df, i, 4, x*u/(v*v)); + gsl_matrix_set (df, i, 5, x*x*u/(v*v)); + gsl_matrix_set (df, i, 6, x*x*x*u/(v*v)); + } + + return GSL_SUCCESS; +} + +int +hahn1_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + hahn1_f (x, params, f); + hahn1_df (x, params, df); + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multifit/test_kirby2.c b/gsl-1.9/multifit/test_kirby2.c new file mode 100644 index 0000000..623a488 --- /dev/null +++ b/gsl-1.9/multifit/test_kirby2.c @@ -0,0 +1,392 @@ +const size_t kirby2_N = 151; +const size_t kirby2_P = 5; + +/* double kirby2_x0[5] = { 2, -0.1, 0.003, -0.001, 0.00001 }; */ + +double kirby2_x0[5] = { 1.5, -0.15, 0.0025, -0.0015, 0.00002 }; + +double kirby2_x[5] = { + 1.6745063063E+00, + -1.3927397867E-01, + 2.5961181191E-03, + -1.7241811870E-03, + 2.1664802578E-05 +}; + +double kirby2_sumsq = 3.9050739624E+00; + +double kirby2_sigma[5] = { + 8.7989634338E-02, + 4.1182041386E-03, + 4.1856520458E-05, + 5.8931897355E-05, + 2.0129761919E-07 +}; + +double kirby2_F1[151] = { + 0.0082E0, + 0.0112E0, + 0.0149E0, + 0.0198E0, + 0.0248E0, + 0.0324E0, + 0.0420E0, + 0.0549E0, + 0.0719E0, + 0.0963E0, + 0.1291E0, + 0.1710E0, + 0.2314E0, + 0.3227E0, + 0.4809E0, + 0.7084E0, + 1.0220E0, + 1.4580E0, + 1.9520E0, + 2.5410E0, + 3.2230E0, + 3.9990E0, + 4.8520E0, + 5.7320E0, + 6.7270E0, + 7.8350E0, + 9.0250E0, + 10.2670E0, + 11.5780E0, + 12.9440E0, + 14.3770E0, + 15.8560E0, + 17.3310E0, + 18.8850E0, + 20.5750E0, + 22.3200E0, + 22.3030E0, + 23.4600E0, + 24.0600E0, + 25.2720E0, + 25.8530E0, + 27.1100E0, + 27.6580E0, + 28.9240E0, + 29.5110E0, + 30.7100E0, + 31.3500E0, + 32.5200E0, + 33.2300E0, + 34.3300E0, + 35.0600E0, + 36.1700E0, + 36.8400E0, + 38.0100E0, + 38.6700E0, + 39.8700E0, + 40.0300E0, + 40.5000E0, + 41.3700E0, + 41.6700E0, + 42.3100E0, + 42.7300E0, + 43.4600E0, + 44.1400E0, + 44.5500E0, + 45.2200E0, + 45.9200E0, + 46.3000E0, + 47.0000E0, + 47.6800E0, + 48.0600E0, + 48.7400E0, + 49.4100E0, + 49.7600E0, + 50.4300E0, + 51.1100E0, + 51.5000E0, + 52.1200E0, + 52.7600E0, + 53.1800E0, + 53.7800E0, + 54.4600E0, + 54.8300E0, + 55.4000E0, + 56.4300E0, + 57.0300E0, + 58.0000E0, + 58.6100E0, + 59.5800E0, + 60.1100E0, + 61.1000E0, + 61.6500E0, + 62.5900E0, + 63.1200E0, + 64.0300E0, + 64.6200E0, + 65.4900E0, + 66.0300E0, + 66.8900E0, + 67.4200E0, + 68.2300E0, + 68.7700E0, + 69.5900E0, + 70.1100E0, + 70.8600E0, + 71.4300E0, + 72.1600E0, + 72.7000E0, + 73.4000E0, + 73.9300E0, + 74.6000E0, + 75.1600E0, + 75.8200E0, + 76.3400E0, + 76.9800E0, + 77.4800E0, + 78.0800E0, + 78.6000E0, + 79.1700E0, + 79.6200E0, + 79.8800E0, + 80.1900E0, + 80.6600E0, + 81.2200E0, + 81.6600E0, + 82.1600E0, + 82.5900E0, + 83.1400E0, + 83.5000E0, + 84.0000E0, + 84.4000E0, + 84.8900E0, + 85.2600E0, + 85.7400E0, + 86.0700E0, + 86.5400E0, + 86.8900E0, + 87.3200E0, + 87.6500E0, + 88.1000E0, + 88.4300E0, + 88.8300E0, + 89.1200E0, + 89.5400E0, + 89.8500E0, + 90.2500E0, + 90.5500E0, + 90.9300E0, + 91.2000E0, + 91.5500E0, + 92.2000E0 +}; + + +double kirby2_F0[151] = { + 9.65E0, + 10.74E0, + 11.81E0, + 12.88E0, + 14.06E0, + 15.28E0, + 16.63E0, + 18.19E0, + 19.88E0, + 21.84E0, + 24.00E0, + 26.25E0, + 28.86E0, + 31.85E0, + 35.79E0, + 40.18E0, + 44.74E0, + 49.53E0, + 53.94E0, + 58.29E0, + 62.63E0, + 67.03E0, + 71.25E0, + 75.22E0, + 79.33E0, + 83.56E0, + 87.75E0, + 91.93E0, + 96.10E0, + 100.28E0, + 104.46E0, + 108.66E0, + 112.71E0, + 116.88E0, + 121.33E0, + 125.79E0, + 125.79E0, + 128.74E0, + 130.27E0, + 133.33E0, + 134.79E0, + 137.93E0, + 139.33E0, + 142.46E0, + 143.90E0, + 146.91E0, + 148.51E0, + 151.41E0, + 153.17E0, + 155.97E0, + 157.76E0, + 160.56E0, + 162.30E0, + 165.21E0, + 166.90E0, + 169.92E0, + 170.32E0, + 171.54E0, + 173.79E0, + 174.57E0, + 176.25E0, + 177.34E0, + 179.19E0, + 181.02E0, + 182.08E0, + 183.88E0, + 185.75E0, + 186.80E0, + 188.63E0, + 190.45E0, + 191.48E0, + 193.35E0, + 195.22E0, + 196.23E0, + 198.05E0, + 199.97E0, + 201.06E0, + 202.83E0, + 204.69E0, + 205.86E0, + 207.58E0, + 209.50E0, + 210.65E0, + 212.33E0, + 215.43E0, + 217.16E0, + 220.21E0, + 221.98E0, + 225.06E0, + 226.79E0, + 229.92E0, + 231.69E0, + 234.77E0, + 236.60E0, + 239.63E0, + 241.50E0, + 244.48E0, + 246.40E0, + 249.35E0, + 251.32E0, + 254.22E0, + 256.24E0, + 259.11E0, + 261.18E0, + 264.02E0, + 266.13E0, + 268.94E0, + 271.09E0, + 273.87E0, + 276.08E0, + 278.83E0, + 281.08E0, + 283.81E0, + 286.11E0, + 288.81E0, + 291.08E0, + 293.75E0, + 295.99E0, + 298.64E0, + 300.84E0, + 302.02E0, + 303.48E0, + 305.65E0, + 308.27E0, + 310.41E0, + 313.01E0, + 315.12E0, + 317.71E0, + 319.79E0, + 322.36E0, + 324.42E0, + 326.98E0, + 329.01E0, + 331.56E0, + 333.56E0, + 336.10E0, + 338.08E0, + 340.60E0, + 342.57E0, + 345.08E0, + 347.02E0, + 349.52E0, + 351.44E0, + 353.93E0, + 355.83E0, + 358.32E0, + 360.20E0, + 362.67E0, + 364.53E0, + 367.00E0, + 371.30E0 +}; + + +int +kirby2_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double b[5]; + size_t i; + + for (i = 0; i < 5; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 151; i++) + { + double x = kirby2_F0[i]; + double y = ((b[0] + x* (b[1] + x * b[2])) + / (1 + x*(b[3] + x *b[4]))); + gsl_vector_set (f, i, kirby2_F1[i] - y); + } + + return GSL_SUCCESS; +} + +int +kirby2_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double b[5]; + size_t i; + + for (i = 0; i < 5; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 151; i++) + { + double x = kirby2_F0[i]; + double u = (b[0] + x*(b[1] + x*b[2])); + double v = (1 + x*(b[3] + x*b[4])); + gsl_matrix_set (df, i, 0, -1/v); + gsl_matrix_set (df, i, 1, -x/v); + gsl_matrix_set (df, i, 2, -x*x/v); + gsl_matrix_set (df, i, 3, x*u/(v*v)); + gsl_matrix_set (df, i, 4, x*x*u/(v*v)); + } + + return GSL_SUCCESS; +} + +int +kirby2_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + kirby2_f (x, params, f); + kirby2_df (x, params, df); + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/multifit/test_longley.c b/gsl-1.9/multifit/test_longley.c new file mode 100644 index 0000000..1d2bf03 --- /dev/null +++ b/gsl-1.9/multifit/test_longley.c @@ -0,0 +1,166 @@ +size_t longley_n = 16; +size_t longley_p = 7; + +double longley_x [] = { + 1, 83.0, 234289, 2356, 1590, 107608, 1947, + 1, 88.5, 259426, 2325, 1456, 108632, 1948, + 1, 88.2, 258054, 3682, 1616, 109773, 1949, + 1, 89.5, 284599, 3351, 1650, 110929, 1950, + 1, 96.2, 328975, 2099, 3099, 112075, 1951, + 1, 98.1, 346999, 1932, 3594, 113270, 1952, + 1, 99.0, 365385, 1870, 3547, 115094, 1953, + 1, 100.0, 363112, 3578, 3350, 116219, 1954, + 1, 101.2, 397469, 2904, 3048, 117388, 1955, + 1, 104.6, 419180, 2822, 2857, 118734, 1956, + 1, 108.4, 442769, 2936, 2798, 120445, 1957, + 1, 110.8, 444546, 4681, 2637, 121950, 1958, + 1, 112.6, 482704, 3813, 2552, 123366, 1959, + 1, 114.2, 502601, 3931, 2514, 125368, 1960, + 1, 115.7, 518173, 4806, 2572, 127852, 1961, + 1, 116.9, 554894, 4007, 2827, 130081, 1962 } ; + +double longley_y[] = {60323, 61122, 60171, 61187, 63221, 63639, 64989, 63761, + 66019, 67857, 68169, 66513, 68655, 69564, 69331, 70551}; + + +void +test_longley () +{ + size_t i, j; + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (longley_n, longley_p); + + gsl_matrix_view X = gsl_matrix_view_array (longley_x, longley_n, longley_p); + gsl_vector_view y = gsl_vector_view_array (longley_y, longley_n); + gsl_vector * c = gsl_vector_alloc (longley_p); + gsl_matrix * cov = gsl_matrix_alloc (longley_p, longley_p); + gsl_vector_view diag; + + double chisq; + + double expected_c[7] = { -3482258.63459582, + 15.0618722713733, + -0.358191792925910E-01, + -2.02022980381683, + -1.03322686717359, + -0.511041056535807E-01, + 1829.15146461355 }; + + double expected_sd[7] = { 890420.383607373, + 84.9149257747669, + 0.334910077722432E-01, + 0.488399681651699, + 0.214274163161675, + 0.226073200069370, + 455.478499142212 } ; + + double expected_chisq = 836424.055505915; + + gsl_multifit_linear (&X.matrix, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "longley gsl_fit_multilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "longley gsl_fit_multilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "longley gsl_fit_multilinear c2") ; + gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-10, "longley gsl_fit_multilinear c3") ; + gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-10, "longley gsl_fit_multilinear c4") ; + gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-10, "longley gsl_fit_multilinear c5") ; + gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-10, "longley gsl_fit_multilinear c6") ; + + diag = gsl_matrix_diagonal (cov); + + gsl_test_rel (gsl_vector_get(&diag.vector,0), pow(expected_sd[0],2.0), 1e-10, "longley gsl_fit_multilinear cov00") ; + gsl_test_rel (gsl_vector_get(&diag.vector,1), pow(expected_sd[1],2.0), 1e-10, "longley gsl_fit_multilinear cov11") ; + gsl_test_rel (gsl_vector_get(&diag.vector,2), pow(expected_sd[2],2.0), 1e-10, "longley gsl_fit_multilinear cov22") ; + gsl_test_rel (gsl_vector_get(&diag.vector,3), pow(expected_sd[3],2.0), 1e-10, "longley gsl_fit_multilinear cov33") ; + gsl_test_rel (gsl_vector_get(&diag.vector,4), pow(expected_sd[4],2.0), 1e-10, "longley gsl_fit_multilinear cov44") ; + gsl_test_rel (gsl_vector_get(&diag.vector,5), pow(expected_sd[5],2.0), 1e-10, "longley gsl_fit_multilinear cov55") ; + gsl_test_rel (gsl_vector_get(&diag.vector,6), pow(expected_sd[6],2.0), 1e-10, "longley gsl_fit_multilinear cov66") ; + + gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_multilinear chisq") ; + + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_multifit_linear_free (work); + } + + + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (longley_n, longley_p); + + gsl_matrix_view X = gsl_matrix_view_array (longley_x, longley_n, longley_p); + gsl_vector_view y = gsl_vector_view_array (longley_y, longley_n); + gsl_vector * w = gsl_vector_alloc (longley_n); + gsl_vector * c = gsl_vector_alloc (longley_p); + gsl_matrix * cov = gsl_matrix_alloc (longley_p, longley_p); + + double chisq; + + double expected_c[7] = { -3482258.63459582, + 15.0618722713733, + -0.358191792925910E-01, + -2.02022980381683, + -1.03322686717359, + -0.511041056535807E-01, + 1829.15146461355 }; + + double expected_cov[7][7] = { { 8531122.56783558, +-166.727799925578, 0.261873708176346, 3.91188317230983, +1.1285582054705, -0.889550869422687, -4362.58709870581}, + +{-166.727799925578, 0.0775861253030891, -1.98725210399982e-05, +-0.000247667096727256, -6.82911920718824e-05, 0.000136160797527761, +0.0775255245956248}, + +{0.261873708176346, -1.98725210399982e-05, 1.20690316701888e-08, +1.66429546772984e-07, 3.61843600487847e-08, -6.78805814483582e-08, +-0.00013158719037715}, + +{3.91188317230983, -0.000247667096727256, 1.66429546772984e-07, +2.56665052544717e-06, 6.96541409215597e-07, -9.00858307771567e-07, +-0.00197260370663974}, + +{1.1285582054705, -6.82911920718824e-05, 3.61843600487847e-08, +6.96541409215597e-07, 4.94032602583969e-07, -9.8469143760973e-08, +-0.000576921112208274}, + +{-0.889550869422687, 0.000136160797527761, -6.78805814483582e-08, +-9.00858307771567e-07, -9.8469143760973e-08, 5.49938542664952e-07, +0.000430074434198215}, + +{-4362.58709870581, 0.0775255245956248, -0.00013158719037715, +-0.00197260370663974, -0.000576921112208274, 0.000430074434198215, +2.23229587481535 }} ; + + double expected_chisq = 836424.055505915; + + gsl_vector_set_all (w, 1.0); + + gsl_multifit_wlinear (&X.matrix, w, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "longley gsl_fit_wmultilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "longley gsl_fit_wmultilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "longley gsl_fit_wmultilinear c2") ; + gsl_test_rel (gsl_vector_get(c,3), expected_c[3], 1e-10, "longley gsl_fit_wmultilinear c3") ; + gsl_test_rel (gsl_vector_get(c,4), expected_c[4], 1e-10, "longley gsl_fit_wmultilinear c4") ; + gsl_test_rel (gsl_vector_get(c,5), expected_c[5], 1e-10, "longley gsl_fit_wmultilinear c5") ; + gsl_test_rel (gsl_vector_get(c,6), expected_c[6], 1e-10, "longley gsl_fit_wmultilinear c6") ; + + for (i = 0; i < longley_p; i++) + { + for (j = 0; j < longley_p; j++) + { + gsl_test_rel (gsl_matrix_get(cov,i,j), expected_cov[i][j], 1e-7, + "longley gsl_fit_wmultilinear cov(%d,%d)", i, j) ; + } + } + + gsl_test_rel (chisq, expected_chisq, 1e-10, "longley gsl_fit_wmultilinear chisq") ; + + gsl_vector_free(w); + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_multifit_linear_free (work); + } +} diff --git a/gsl-1.9/multifit/test_nelson.c b/gsl-1.9/multifit/test_nelson.c new file mode 100644 index 0000000..19f594e --- /dev/null +++ b/gsl-1.9/multifit/test_nelson.c @@ -0,0 +1,217 @@ +const size_t nelson_N = 128; +const size_t nelson_P = 3; + +/* double nelson_x0[3] = { 2, 0.0001, -0.01}; */ + +double nelson_x0[3] = { 2.5 , 0.000000005, -0.01}; + +double nelson_x[3] = { + 2.5906836021E+00, + 5.6177717026E-09, + -5.7701013174E-02 +}; + +double nelson_sumsq = 3.7976833176E+00; + +double nelson_sigma[3] = { + 1.9149996413E-02, + 6.1124096540E-09, + 3.9572366543E-03 +}; + +double nelson_F[128][3] = { + { 15.00E0, 1E0, 180E0}, + { 17.00E0, 1E0, 180E0}, + { 15.50E0, 1E0, 180E0}, + { 16.50E0, 1E0, 180E0}, + { 15.50E0, 1E0, 225E0}, + { 15.00E0, 1E0, 225E0}, + { 16.00E0, 1E0, 225E0}, + { 14.50E0, 1E0, 225E0}, + { 15.00E0, 1E0, 250E0}, + { 14.50E0, 1E0, 250E0}, + { 12.50E0, 1E0, 250E0}, + { 11.00E0, 1E0, 250E0}, + { 14.00E0, 1E0, 275E0}, + { 13.00E0, 1E0, 275E0}, + { 14.00E0, 1E0, 275E0}, + { 11.50E0, 1E0, 275E0}, + { 14.00E0, 2E0, 180E0}, + { 16.00E0, 2E0, 180E0}, + { 13.00E0, 2E0, 180E0}, + { 13.50E0, 2E0, 180E0}, + { 13.00E0, 2E0, 225E0}, + { 13.50E0, 2E0, 225E0}, + { 12.50E0, 2E0, 225E0}, + { 12.50E0, 2E0, 225E0}, + { 12.50E0, 2E0, 250E0}, + { 12.00E0, 2E0, 250E0}, + { 11.50E0, 2E0, 250E0}, + { 12.00E0, 2E0, 250E0}, + { 13.00E0, 2E0, 275E0}, + { 11.50E0, 2E0, 275E0}, + { 13.00E0, 2E0, 275E0}, + { 12.50E0, 2E0, 275E0}, + { 13.50E0, 4E0, 180E0}, + { 17.50E0, 4E0, 180E0}, + { 17.50E0, 4E0, 180E0}, + { 13.50E0, 4E0, 180E0}, + { 12.50E0, 4E0, 225E0}, + { 12.50E0, 4E0, 225E0}, + { 15.00E0, 4E0, 225E0}, + { 13.00E0, 4E0, 225E0}, + { 12.00E0, 4E0, 250E0}, + { 13.00E0, 4E0, 250E0}, + { 12.00E0, 4E0, 250E0}, + { 13.50E0, 4E0, 250E0}, + { 10.00E0, 4E0, 275E0}, + { 11.50E0, 4E0, 275E0}, + { 11.00E0, 4E0, 275E0}, + { 9.50E0, 4E0, 275E0}, + { 15.00E0, 8E0, 180E0}, + { 15.00E0, 8E0, 180E0}, + { 15.50E0, 8E0, 180E0}, + { 16.00E0, 8E0, 180E0}, + { 13.00E0, 8E0, 225E0}, + { 10.50E0, 8E0, 225E0}, + { 13.50E0, 8E0, 225E0}, + { 14.00E0, 8E0, 225E0}, + { 12.50E0, 8E0, 250E0}, + { 12.00E0, 8E0, 250E0}, + { 11.50E0, 8E0, 250E0}, + { 11.50E0, 8E0, 250E0}, + { 6.50E0, 8E0, 275E0}, + { 5.50E0, 8E0, 275E0}, + { 6.00E0, 8E0, 275E0}, + { 6.00E0, 8E0, 275E0}, + { 18.50E0, 16E0, 180E0}, + { 17.00E0, 16E0, 180E0}, + { 15.30E0, 16E0, 180E0}, + { 16.00E0, 16E0, 180E0}, + { 13.00E0, 16E0, 225E0}, + { 14.00E0, 16E0, 225E0}, + { 12.50E0, 16E0, 225E0}, + { 11.00E0, 16E0, 225E0}, + { 12.00E0, 16E0, 250E0}, + { 12.00E0, 16E0, 250E0}, + { 11.50E0, 16E0, 250E0}, + { 12.00E0, 16E0, 250E0}, + { 6.00E0, 16E0, 275E0}, + { 6.00E0, 16E0, 275E0}, + { 5.00E0, 16E0, 275E0}, + { 5.50E0, 16E0, 275E0}, + { 12.50E0, 32E0, 180E0}, + { 13.00E0, 32E0, 180E0}, + { 16.00E0, 32E0, 180E0}, + { 12.00E0, 32E0, 180E0}, + { 11.00E0, 32E0, 225E0}, + { 9.50E0, 32E0, 225E0}, + { 11.00E0, 32E0, 225E0}, + { 11.00E0, 32E0, 225E0}, + { 11.00E0, 32E0, 250E0}, + { 10.00E0, 32E0, 250E0}, + { 10.50E0, 32E0, 250E0}, + { 10.50E0, 32E0, 250E0}, + { 2.70E0, 32E0, 275E0}, + { 2.70E0, 32E0, 275E0}, + { 2.50E0, 32E0, 275E0}, + { 2.40E0, 32E0, 275E0}, + { 13.00E0, 48E0, 180E0}, + { 13.50E0, 48E0, 180E0}, + { 16.50E0, 48E0, 180E0}, + { 13.60E0, 48E0, 180E0}, + { 11.50E0, 48E0, 225E0}, + { 10.50E0, 48E0, 225E0}, + { 13.50E0, 48E0, 225E0}, + { 12.00E0, 48E0, 225E0}, + { 7.00E0, 48E0, 250E0}, + { 6.90E0, 48E0, 250E0}, + { 8.80E0, 48E0, 250E0}, + { 7.90E0, 48E0, 250E0}, + { 1.20E0, 48E0, 275E0}, + { 1.50E0, 48E0, 275E0}, + { 1.00E0, 48E0, 275E0}, + { 1.50E0, 48E0, 275E0}, + { 13.00E0, 64E0, 180E0}, + { 12.50E0, 64E0, 180E0}, + { 16.50E0, 64E0, 180E0}, + { 16.00E0, 64E0, 180E0}, + { 11.00E0, 64E0, 225E0}, + { 11.50E0, 64E0, 225E0}, + { 10.50E0, 64E0, 225E0}, + { 10.00E0, 64E0, 225E0}, + { 7.27E0, 64E0, 250E0}, + { 7.50E0, 64E0, 250E0}, + { 6.70E0, 64E0, 250E0}, + { 7.60E0, 64E0, 250E0}, + { 1.50E0, 64E0, 275E0}, + { 1.00E0, 64E0, 275E0}, + { 1.20E0, 64E0, 275E0}, + { 1.20E0, 64E0, 275E0} +}; + + +int +nelson_f (const gsl_vector * x, void *params, gsl_vector * f) +{ + double b[3]; + size_t i; + + for (i = 0; i < 3; i++) + { + b[i] = gsl_vector_get(x, i); + } + + + + for (i = 0; i < 128; i++) + { + double x1 = nelson_F[i][1]; + double x2 = nelson_F[i][2]; + + double y = b[0] - b[1] * x1 * (b[2]*x2 < -100) ? 0.0 : exp(-b[2] * x2); + + gsl_vector_set (f, i, log(nelson_F[i][0]) - y); + } + + return GSL_SUCCESS; +} + +int +nelson_df (const gsl_vector * x, void *params, gsl_matrix * df) +{ + double b[3]; + size_t i; + + for (i = 0; i < 3; i++) + { + b[i] = gsl_vector_get(x, i); + } + + for (i = 0; i < 128; i++) + { + double x1 = nelson_F[i][1]; + double x2 = nelson_F[i][2]; + + gsl_matrix_set (df, i, 0, -1.0); + gsl_matrix_set (df, i, 1, x1 * exp(-b[2] * x2)); + gsl_matrix_set (df, i, 2, -b[1] * x1 * x2 * exp(-b[2] * x2)); + } + + return GSL_SUCCESS; +} + +int +nelson_fdf (const gsl_vector * x, void *params, + gsl_vector * f, gsl_matrix * df) +{ + nelson_f (x, params, f); + nelson_df (x, params, df); + + return GSL_SUCCESS; +} + + + + + diff --git a/gsl-1.9/multifit/test_pontius.c b/gsl-1.9/multifit/test_pontius.c new file mode 100644 index 0000000..69c6905 --- /dev/null +++ b/gsl-1.9/multifit/test_pontius.c @@ -0,0 +1,131 @@ +size_t pontius_n = 40; +size_t pontius_p = 3; + +double pontius_x[] = { 150000, 300000, 450000, 600000, 750000, 900000, +1050000, 1200000, 1350000, 1500000, 1650000, 1800000, 1950000, 2100000, +2250000, 2400000, 2550000, 2700000, 2850000, 3000000, 150000, 300000, +450000, 600000, 750000, 900000, 1050000, 1200000, 1350000, 1500000, +1650000, 1800000, 1950000, 2100000, 2250000, 2400000, 2550000, 2700000, +2850000, 3000000 }; + +double pontius_y[] = { .11019, .21956, .32949, .43899, .54803, .65694, +.76562, .87487, .98292, 1.09146, 1.20001, 1.30822, 1.41599, 1.52399, +1.63194, 1.73947, 1.84646, 1.95392, 2.06128, 2.16844, .11052, .22018, +.32939, .43886, .54798, .65739, .76596, .87474, .98300, 1.09150, +1.20004, 1.30818, 1.41613, 1.52408, 1.63159, 1.73965, 1.84696, +1.95445, 2.06177, 2.16829 }; + +void +test_pontius () +{ + size_t i, j; + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (pontius_n, pontius_p); + + gsl_matrix * X = gsl_matrix_alloc (pontius_n, pontius_p); + gsl_vector_view y = gsl_vector_view_array (pontius_y, pontius_n); + gsl_vector * c = gsl_vector_alloc (pontius_p); + gsl_matrix * cov = gsl_matrix_alloc (pontius_p, pontius_p); + gsl_vector_view diag; + + double chisq; + + double expected_c[3] = { 0.673565789473684E-03, + 0.732059160401003E-06, + -0.316081871345029E-14}; + + double expected_sd[3] = { 0.107938612033077E-03, + 0.157817399981659E-09, + 0.486652849992036E-16 }; + + double expected_chisq = 0.155761768796992E-05; + + for (i = 0 ; i < pontius_n; i++) + { + for (j = 0; j < pontius_p; j++) + { + gsl_matrix_set(X, i, j, pow(pontius_x[i], j)); + } + } + + gsl_multifit_linear (X, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "pontius gsl_fit_multilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "pontius gsl_fit_multilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "pontius gsl_fit_multilinear c2") ; + + diag = gsl_matrix_diagonal (cov); + + gsl_test_rel (gsl_vector_get(&diag.vector,0), pow(expected_sd[0],2.0), 1e-10, "pontius gsl_fit_multilinear cov00") ; + gsl_test_rel (gsl_vector_get(&diag.vector,1), pow(expected_sd[1],2.0), 1e-10, "pontius gsl_fit_multilinear cov11") ; + gsl_test_rel (gsl_vector_get(&diag.vector,2), pow(expected_sd[2],2.0), 1e-10, "pontius gsl_fit_multilinear cov22") ; + + gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_multilinear chisq") ; + + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_matrix_free(X); + gsl_multifit_linear_free (work); + } + + + { + gsl_multifit_linear_workspace * work = + gsl_multifit_linear_alloc (pontius_n, pontius_p); + + gsl_matrix * X = gsl_matrix_alloc (pontius_n, pontius_p); + gsl_vector_view y = gsl_vector_view_array (pontius_y, pontius_n); + gsl_vector * w = gsl_vector_alloc (pontius_n); + gsl_vector * c = gsl_vector_alloc (pontius_p); + gsl_matrix * cov = gsl_matrix_alloc (pontius_p, pontius_p); + + double chisq; + + double expected_c[3] = { 0.673565789473684E-03, + 0.732059160401003E-06, + -0.316081871345029E-14}; + + double expected_chisq = 0.155761768796992E-05; + + double expected_cov[3][3] ={ + {2.76754385964916e-01 , -3.59649122807024e-07, 9.74658869395731e-14}, + {-3.59649122807024e-07, 5.91630591630603e-13, -1.77210703526497e-19}, + {9.74658869395731e-14, -1.77210703526497e-19, 5.62573661988878e-26} }; + + + for (i = 0 ; i < pontius_n; i++) + { + for (j = 0; j < pontius_p; j++) + { + gsl_matrix_set(X, i, j, pow(pontius_x[i], j)); + } + } + + gsl_vector_set_all (w, 1.0); + + gsl_multifit_wlinear (X, w, &y.vector, c, cov, &chisq, work); + + gsl_test_rel (gsl_vector_get(c,0), expected_c[0], 1e-10, "pontius gsl_fit_multilinear c0") ; + gsl_test_rel (gsl_vector_get(c,1), expected_c[1], 1e-10, "pontius gsl_fit_multilinear c1") ; + gsl_test_rel (gsl_vector_get(c,2), expected_c[2], 1e-10, "pontius gsl_fit_multilinear c2") ; + + + for (i = 0; i < pontius_p; i++) + { + for (j = 0; j < pontius_p; j++) + { + gsl_test_rel (gsl_matrix_get(cov,i,j), expected_cov[i][j], 1e-10, + "pontius gsl_fit_wmultilinear cov(%d,%d)", i, j) ; + } + } + + gsl_test_rel (chisq, expected_chisq, 1e-10, "pontius gsl_fit_multilinear chisq") ; + + gsl_vector_free(w); + gsl_vector_free(c); + gsl_matrix_free(cov); + gsl_matrix_free(X); + gsl_multifit_linear_free (work); + } +} diff --git a/gsl-1.9/multifit/work.c b/gsl-1.9/multifit/work.c new file mode 100644 index 0000000..f05e3ff --- /dev/null +++ b/gsl-1.9/multifit/work.c @@ -0,0 +1,133 @@ +/* multifit/work.c + * + * Copyright (C) 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_multifit.h> + +gsl_multifit_linear_workspace * +gsl_multifit_linear_alloc (size_t n, size_t p) +{ + gsl_multifit_linear_workspace *w; + + w = (gsl_multifit_linear_workspace *) + malloc (sizeof (gsl_multifit_linear_workspace)); + + if (w == 0) + { + GSL_ERROR_VAL ("failed to allocate space for multifit_linear struct", + GSL_ENOMEM, 0); + } + + w->n = n; /* number of observations */ + w->p = p; /* number of parameters */ + + w->A = gsl_matrix_alloc (n, p); + + if (w->A == 0) + { + free (w); + GSL_ERROR_VAL ("failed to allocate space for A", GSL_ENOMEM, 0); + } + + w->Q = gsl_matrix_alloc (p, p); + + if (w->Q == 0) + { + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for Q", GSL_ENOMEM, 0); + } + + w->QSI = gsl_matrix_alloc (p, p); + + if (w->QSI == 0) + { + gsl_matrix_free (w->Q); + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for QSI", GSL_ENOMEM, 0); + } + + w->S = gsl_vector_alloc (p); + + if (w->S == 0) + { + gsl_matrix_free (w->QSI); + gsl_matrix_free (w->Q); + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for S", GSL_ENOMEM, 0); + } + + w->t = gsl_vector_alloc (n); + + if (w->t == 0) + { + gsl_vector_free (w->S); + gsl_matrix_free (w->QSI); + gsl_matrix_free (w->Q); + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for t", GSL_ENOMEM, 0); + } + + w->xt = gsl_vector_calloc (p); + + if (w->xt == 0) + { + gsl_vector_free (w->t); + gsl_vector_free (w->S); + gsl_matrix_free (w->QSI); + gsl_matrix_free (w->Q); + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for xt", GSL_ENOMEM, 0); + } + + w->D = gsl_vector_calloc (p); + + if (w->D == 0) + { + gsl_vector_free (w->D); + gsl_vector_free (w->t); + gsl_vector_free (w->S); + gsl_matrix_free (w->QSI); + gsl_matrix_free (w->Q); + gsl_matrix_free (w->A); + free (w); + GSL_ERROR_VAL ("failed to allocate space for xt", GSL_ENOMEM, 0); + } + + return w; +} + +void +gsl_multifit_linear_free (gsl_multifit_linear_workspace * work) +{ + gsl_matrix_free (work->A); + gsl_matrix_free (work->Q); + gsl_matrix_free (work->QSI); + gsl_vector_free (work->S); + gsl_vector_free (work->t); + gsl_vector_free (work->xt); + gsl_vector_free (work->D); + free (work); +} + |