Index: trunk/tests/test_SpecialFunctions.cc =================================================================== --- trunk/tests/test_SpecialFunctions.cc (revision 816) +++ trunk/tests/test_SpecialFunctions.cc (revision 817) @@ -1,71 +1,72 @@ #include #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/nm/SpecialFunctions.hh" #define SQRT2L 1.414213562373095048801689L #define SQRPI 1.77245385090551602729816748334L using namespace npstat; using namespace std; static long double inverseErf(const long double fval) { long double x = npstat::inverseGaussCdf((fval + 1.0L)/2.0L)/SQRT2L; for (unsigned i=0; i<2; ++i) { const long double guessed = erfl(x); const long double deri = 2.0/SQRPI*expl(-x*x); x += (fval - guessed)/deri; } return x; } namespace { TEST(localInverseErf) { for (unsigned i=0; i<100; ++i) { long double x = -1 + 2.0L/100*(i + 0.5L); long double erfval = erfl(x); long double inverse = inverseErf(erfval); CHECK_CLOSE(x, inverse, 1.0e-19L); } } TEST(incompleteGamma) { const long double eps = 1.0e-9; CHECK_CLOSE(1.0-exp(-5.0), incompleteGamma(1.0, 5.0), eps); CHECK_CLOSE(1.0-0.0404276819945128, incompleteGamma(2.0, 5.0), eps); + CHECK_CLOSE(0.87534798051691886, incompleteGamma(3.0, 5.0), eps); CHECK_CLOSE(5.0, inverseIncompleteGamma(1.0, 1.0-exp(-5.0)), eps); CHECK_CLOSE(5.0, inverseIncompleteGamma(2.0, 1.0-0.0404276819945128), eps); CHECK_CLOSE(5.0, inverseIncompleteGammaC(1.0, exp(-5.0)), eps); CHECK_CLOSE(5.0, inverseIncompleteGammaC(2.0, 0.0404276819945128), eps); } TEST(normalDensityDerivative) { const long double eps = 1.0e-20; // f[x_] := Exp[-x^2/2]/Sqrt[2 Pi] // N[D[f[x], {x, 0}]/.x->1/10, 25] CHECK_CLOSE(0.3969525474770117655105297L, normalDensityDerivative(0, 0.1L), eps); // N[D[f[x], {x, 1}]/.x->2/10, 25] CHECK_CLOSE(-0.07820853879509117560175475L, normalDensityDerivative(1, 0.2L), eps); // N[D[f[x], {x, 2}]/.x->3/10, 25] CHECK_CLOSE(-0.3470629120690769179033078L, normalDensityDerivative(2, 0.3L), eps); // N[D[f[x], {x, 3}]/.x->4/10, 25] CHECK_CLOSE(0.4183548793845752775973313L, normalDensityDerivative(3, 0.4L), eps); // N[D[f[x], {x, 4}]/.x->5/10, 25] CHECK_CLOSE(0.5501020730692179340229382L, normalDensityDerivative(4, 0.5L), eps); } } Index: trunk/tests/test_DistributionsND.cc =================================================================== --- trunk/tests/test_DistributionsND.cc (revision 816) +++ trunk/tests/test_DistributionsND.cc (revision 817) @@ -1,62 +1,148 @@ #include "UnitTest++.h" #include "test_utils.hh" +#include "geners/binaryIO.hh" + #include "npstat/nm/lapack_double.h" #include "npstat/nm/SimpleFunctors.hh" +#include "npstat/stat/Distributions1D.hh" #include "npstat/stat/DistributionsND.hh" #include "npstat/stat/ScalableGaussND.hh" +#include "npstat/stat/EllipticalDistribution.hh" +#include "npstat/stat/MultivariateSumsqAccumulator.hh" +#include "npstat/stat/InMemoryNtuple.hh" +#include "npstat/stat/DistributionNDReader.hh" using namespace npstat; using namespace std; namespace { + class GaussND2 : public EllipticalDistribution + { + public: + typedef EllipticalDistribution Base; + + inline GaussND2(const double* location, const unsigned dim, + const Matrix& covmat) + : EllipticalDistribution(location, dim, + covmat.symFcn(SquareRoot()), + Exponential1D(0.0, 2.0), + Gamma1D(0.0, 2.0, dim/2.0)) {} + + inline virtual ~GaussND2() {} + + inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} + inline virtual bool write(std::ostream& os) const + {return Base::classId().write(os) && Base::write(os);} + + static inline const char* classname() {return "test::GaussND2";} + static inline unsigned version() {return 1;} + static inline GaussND2* read(const gs::ClassId& id, std::istream& in) + { + static const gs::ClassId current(gs::ClassId::makeId()); + current.ensureSameId(id); + CPP11_auto_ptr pb = gs::read_obj(in); + return new GaussND2(*pb); + } + + private: + inline GaussND2(const Base& b) : EllipticalDistribution(b) {} + }; + void io_test(const AbsDistributionND& d) { std::ostringstream os; CHECK(d.classId().write(os)); CHECK(d.write(os)); std::istringstream is(os.str()); gs::ClassId id(is, 1); AbsDistributionND* phoenix = AbsDistributionND::read(id, is); CHECK(*phoenix == d); delete phoenix; } TEST(LinTransformedDistroND) { const double eps = 1.0e-10; const lapack_double covdata[2][2] = {{1.0, 0.3}, {0.3, 0.25}}; const Matrix covmat(2, 2, &covdata[0][0]); const double shifts[2] = {0.5, 0.2}; const double zeros[2] = {0.0, 0.0}; const double ones[2] = {1.0, 1.0}; const ScalableGaussND g0(zeros, ones, 2); const GaussND g1(shifts, 2, covmat); const Matrix covSqr = covmat.symFcn(SquareRoot()); const LinTransformedDistroND g2(g0, shifts, 2, covSqr); double x[2], r[2], y[2]; for (unsigned i=0; i<100; ++i) { x[0] = (test_rng()-0.5)*5.0; x[1] = test_rng()-0.5; const double d0 = g1.density(x, 2); const double d1 = g2.density(x, 2); CHECK_CLOSE(d0, d1, eps); r[0] = test_rng(); r[1] = test_rng(); g1.unitMap(r, 2, x); g2.unitMap(r, 2, y); CHECK_CLOSE(x[0], y[0], eps); CHECK_CLOSE(x[1], y[1], eps); } io_test(g0); io_test(g1); io_test(g2); } + + TEST(EllipticalDistribution) + { + const double eps = 1.0e-10; + const unsigned maxdim = 10; + const unsigned ntries = 20; + + StaticDistributionNDReader::registerClass(); + + for (unsigned dim=1; dim mu(dim); + for (unsigned i=0; i nt(simpleColumnNames(dim)); + const unsigned nfills = 5U*dim; + std::vector rnd(dim); + for (unsigned ifill=0; ifill sums; + nt.cycleOverRows(sums); + MultivariateSumsqAccumulator<> sumsqs(sums); + nt.cycleOverRows(sumsqs); + const Matrix covmat(sumsqs.covMat()); + + // Now we can construct a reference normal + // and an elliptical normal + const GaussND g1(&mu[0], dim, covmat); + const GaussND2 g2(&mu[0], dim, covmat); + + io_test(g1); + io_test(g2); + + for (unsigned itry=0; itry 1, - 'sinclude' => 1, - 'm4_define' => 1, - 'AM_AUTOMAKE_VERSION' => 1, - 'AC_FC_PP_DEFINE' => 1, - '_AM_COND_IF' => 1, - 'AM_NLS' => 1, - 'AC_CANONICAL_TARGET' => 1, - 'AM_ENABLE_MULTILIB' => 1, - 'AM_PROG_CXX_C_O' => 1, - 'include' => 1, - '_AM_COND_ENDIF' => 1, - 'm4_include' => 1, - 'define' => 1, - 'AC_CONFIG_MACRO_DIR_TRACE' => 1, - 'AC_FC_FREEFORM' => 1, 'AC_CANONICAL_BUILD' => 1, - 'AM_PROG_FC_C_O' => 1, - 'm4_pattern_allow' => 1, + 'AC_CANONICAL_HOST' => 1, + 'AM_XGETTEXT_OPTION' => 1, + 'AM_PATH_GUILE' => 1, + 'AM_MAKEFILE_INCLUDE' => 1, + 'm4_include' => 1, + 'AC_FC_SRCEXT' => 1, + 'AC_CANONICAL_TARGET' => 1, 'm4_pattern_forbid' => 1, - 'AM_PROG_MKDIR_P' => 1, + '_AM_COND_ELSE' => 1, + 'AC_PROG_LIBTOOL' => 1, + 'AC_CONFIG_FILES' => 1, 'AC_SUBST' => 1, - 'AM_GNU_GETTEXT_INTL_SUBDIR' => 1, - 'AM_PROG_MOC' => 1, + 'AM_EXTRA_RECURSIVE_TARGETS' => 1, + 'include' => 1, + 'AC_FC_PP_DEFINE' => 1, + 'AM_PROG_MKDIR_P' => 1, + 'AM_PROG_LIBTOOL' => 1, 'AC_CANONICAL_SYSTEM' => 1, - 'AM_MAINTAINER_MODE' => 1, + 'AM_AUTOMAKE_VERSION' => 1, + 'AM_SILENT_RULES' => 1, + 'AC_CONFIG_LINKS' => 1, + 'AM_GNU_GETTEXT_INTL_SUBDIR' => 1, + 'AM_PROG_CC_C_O' => 1, + 'AC_DEFINE_TRACE_LITERAL' => 1, + 'AC_INIT' => 1, '_AM_SUBST_NOTMAKE' => 1, + 'AM_GNU_GETTEXT' => 1, + 'AM_PROG_MOC' => 1, + 'AC_CONFIG_LIBOBJ_DIR' => 1, + 'LT_INIT' => 1, + 'AC_CONFIG_MACRO_DIR_TRACE' => 1, + 'm4_define' => 1, + 'AC_FC_FREEFORM' => 1, + 'AM_CONDITIONAL' => 1, + '_LT_AC_TAGCONFIG' => 1, 'm4_sinclude' => 1, 'AC_FC_PP_SRCEXT' => 1, - 'AC_CONFIG_HEADERS' => 1, - 'AC_LIBSOURCE' => 1, - 'AC_CONFIG_FILES' => 1, 'AC_CONFIG_SUBDIRS' => 1, - 'GTK_DOC_CHECK' => 1, - 'AC_CANONICAL_HOST' => 1, - 'AC_CONFIG_LINKS' => 1, - 'AC_REQUIRE_AUX_FILE' => 1, - '_AM_COND_ELSE' => 1, - 'LT_INIT' => 1, - '_LT_AC_TAGCONFIG' => 1, - 'AU_DEFINE' => 1, + 'AC_CONFIG_AUX_DIR' => 1, + '_AM_MAKEFILE_INCLUDE' => 1, 'AH_OUTPUT' => 1, - 'AC_PROG_LIBTOOL' => 1, + 'AC_LIBSOURCE' => 1, + 'AC_CONFIG_HEADERS' => 1, 'AM_POT_TOOLS' => 1, - 'AM_SILENT_RULES' => 1, - 'AM_GNU_GETTEXT' => 1, - 'AM_PATH_GUILE' => 1, - 'AM_PROG_F77_C_O' => 1, - 'AM_PROG_AR' => 1, - 'AC_SUBST_TRACE' => 1, - 'LT_SUPPORTED_TAG' => 1, - 'AC_DEFINE_TRACE_LITERAL' => 1, - 'AC_CONFIG_AUX_DIR' => 1, - 'AM_CONDITIONAL' => 1, - '_m4_warn' => 1, - 'AC_FC_SRCEXT' => 1, - 'AM_PROG_LIBTOOL' => 1, 'IT_PROG_INTLTOOL' => 1, - 'AC_CONFIG_LIBOBJ_DIR' => 1, - 'AM_INIT_AUTOMAKE' => 1, + '_m4_warn' => 1, + 'AM_PROG_AR' => 1, + 'm4_pattern_allow' => 1, + 'AM_MAINTAINER_MODE' => 1, + 'sinclude' => 1, + '_AM_COND_ENDIF' => 1, + 'AM_NLS' => 1, 'LT_CONFIG_LTDL_DIR' => 1, - 'AM_MAKEFILE_INCLUDE' => 1, - 'AM_PROG_CC_C_O' => 1, - 'AM_EXTRA_RECURSIVE_TARGETS' => 1, - '_AM_MAKEFILE_INCLUDE' => 1, - 'AM_XGETTEXT_OPTION' => 1 + 'define' => 1, + 'AM_PROG_FC_C_O' => 1, + 'GTK_DOC_CHECK' => 1, + 'AM_ENABLE_MULTILIB' => 1, + 'AC_REQUIRE_AUX_FILE' => 1, + 'AU_DEFINE' => 1, + 'AM_PROG_CXX_C_O' => 1, + 'AM_INIT_AUTOMAKE' => 1, + 'AC_SUBST_TRACE' => 1, + 'AM_PROG_F77_C_O' => 1, + '_AM_COND_IF' => 1, + 'LT_SUPPORTED_TAG' => 1 } ], 'Autom4te::Request' ), bless( [ '1', 1, [ '/usr/share/autoconf' ], [ '/usr/share/autoconf/autoconf/autoconf.m4f', '/usr/share/aclocal-1.16/internal/ac-config-macro-dirs.m4', '/usr/share/aclocal/ltargz.m4', '/usr/share/aclocal/ltdl.m4', '/usr/share/aclocal/pkg.m4', '/usr/share/aclocal-1.16/amversion.m4', '/usr/share/aclocal-1.16/auxdir.m4', '/usr/share/aclocal-1.16/cond.m4', '/usr/share/aclocal-1.16/depend.m4', '/usr/share/aclocal-1.16/depout.m4', '/usr/share/aclocal-1.16/extra-recurs.m4', '/usr/share/aclocal-1.16/init.m4', '/usr/share/aclocal-1.16/install-sh.m4', '/usr/share/aclocal-1.16/lead-dot.m4', '/usr/share/aclocal-1.16/make.m4', '/usr/share/aclocal-1.16/missing.m4', '/usr/share/aclocal-1.16/options.m4', '/usr/share/aclocal-1.16/prog-cc-c-o.m4', '/usr/share/aclocal-1.16/runlog.m4', '/usr/share/aclocal-1.16/sanity.m4', '/usr/share/aclocal-1.16/silent.m4', '/usr/share/aclocal-1.16/strip.m4', '/usr/share/aclocal-1.16/substnot.m4', '/usr/share/aclocal-1.16/tar.m4', 'm4/libtool.m4', 'm4/ltoptions.m4', 'm4/ltsugar.m4', 'm4/ltversion.m4', 'm4/lt~obsolete.m4', 'configure.ac' ], { - '_AM_CONFIG_MACRO_DIRS' => 1, + 'AC_DEFUN' => 1, + 'AC_ENABLE_STATIC' => 1, + 'AC_LTDL_SHLIBPATH' => 1, + 'AC_LIBTOOL_PICMODE' => 1, + 'AC_PROG_LD' => 1, + 'PKG_CHECK_EXISTS' => 1, + 'AM_PROG_CC_C_O' => 1, 'AM_PROG_INSTALL_STRIP' => 1, - '_AC_AM_CONFIG_HEADER_HOOK' => 1, + 'LT_LANG' => 1, 'AC_DISABLE_STATIC' => 1, - 'AC_DISABLE_FAST_INSTALL' => 1, - 'AM_PROG_NM' => 1, - 'AC_LTDL_DLSYM_USCORE' => 1, - 'LT_SYS_DLOPEN_SELF' => 1, - 'AM_ENABLE_STATIC' => 1, - 'm4_pattern_allow' => 1, - 'LTDL_INSTALLABLE' => 1, - '_AM_AUTOCONF_VERSION' => 1, + 'PKG_CHECK_MODULES_STATIC' => 1, + 'LTOPTIONS_VERSION' => 1, + '_AM_DEPENDENCIES' => 1, + 'AC_DEPLIBS_CHECK_METHOD' => 1, + 'AC_LIBTOOL_SYS_GLOBAL_SYMBOL_PIPE' => 1, + 'AM_PROG_LIBTOOL' => 1, + 'AC_PATH_MAGIC' => 1, + 'LT_CMD_MAX_LEN' => 1, + 'm4_include' => 1, + 'AC_ENABLE_SHARED' => 1, + 'AC_LTDL_SYMBOL_USCORE' => 1, + 'AC_LIBTOOL_SETUP' => 1, 'AM_OUTPUT_DEPENDENCY_COMMANDS' => 1, - '_LT_DLL_DEF_P' => 1, + '_AM_AUTOCONF_VERSION' => 1, + 'AC_LIBTOOL_CONFIG' => 1, + 'AC_LIBTOOL_OBJDIR' => 1, + 'PKG_NOARCH_INSTALLDIR' => 1, + 'AC_LIBTOOL_SYS_LIB_STRIP' => 1, + 'AM_PROG_NM' => 1, + 'AC_LIB_LTDL' => 1, + '_LT_PREPARE_SED_QUOTE_VARS' => 1, + '_LT_AC_PROG_CXXCPP' => 1, 'AM_DISABLE_SHARED' => 1, 'LT_AC_PROG_SED' => 1, - '_AM_DEPENDENCIES' => 1, - 'AC_LIBTOOL_DLOPEN' => 1, - '_AM_SUBST_NOTMAKE' => 1, - '_LT_AC_LANG_GCJ_CONFIG' => 1, - '_LT_AC_SHELL_INIT' => 1, - '_LT_CC_BASENAME' => 1, - 'AC_LIBTOOL_PROG_COMPILER_NO_RTTI' => 1, - 'LT_FUNC_ARGZ' => 1, - 'LT_SYS_DLSEARCH_PATH' => 1, - 'AC_LIBTOOL_LANG_GCJ_CONFIG' => 1, - '_LT_PREPARE_SED_QUOTE_VARS' => 1, - '_AM_MANGLE_OPTION' => 1, - 'LT_AC_PROG_EGREP' => 1, - 'LT_SYS_DLOPEN_DEPLIBS' => 1, - 'AC_LTDL_ENABLE_INSTALL' => 1, - 'AC_PROG_LIBTOOL' => 1, - '_LT_PROG_ECHO_BACKSLASH' => 1, - '_LT_LINKER_OPTION' => 1, - '_LT_REQUIRED_DARWIN_CHECKS' => 1, + 'LT_SYS_DLOPEN_SELF' => 1, + 'PKG_INSTALLDIR' => 1, + 'LT_SYS_MODULE_PATH' => 1, + 'PKG_CHECK_VAR' => 1, + 'LT_CONFIG_LTDL_DIR' => 1, + '_LT_COMPILER_OPTION' => 1, + 'PKG_PROG_PKG_CONFIG' => 1, + 'LT_AC_PROG_GCJ' => 1, + 'LT_FUNC_DLSYM_USCORE' => 1, 'LTOBSOLETE_VERSION' => 1, - 'LTOPTIONS_VERSION' => 1, - 'PKG_CHECK_MODULES_STATIC' => 1, - 'AC_DEFUN_ONCE' => 1, - '_LT_AC_PROG_CXXCPP' => 1, - 'AM_INIT_AUTOMAKE' => 1, - 'AC_ENABLE_FAST_INSTALL' => 1, - 'LT_LANG' => 1, + 'LT_SYS_MODULE_EXT' => 1, + 'AC_LIBTOOL_PROG_COMPILER_NO_RTTI' => 1, + '_LT_AC_PROG_ECHO_BACKSLASH' => 1, + 'LT_LIB_DLLOAD' => 1, + 'AC_LIBTOOL_POSTDEP_PREDEP' => 1, 'AC_LIBTOOL_LANG_CXX_CONFIG' => 1, - 'LTDL_CONVENIENCE' => 1, - 'LT_SYS_SYMBOL_USCORE' => 1, - 'LT_LIB_M' => 1, - '_LT_AC_LANG_CXX' => 1, - 'AC_LIBTOOL_SYS_DYNAMIC_LINKER' => 1, - '_m4_warn' => 1, - 'LT_AC_PROG_GCJ' => 1, - 'AC_LTDL_SYSSEARCHPATH' => 1, + 'AC_LIBTOOL_DLOPEN_SELF' => 1, + 'AC_LIBTOOL_LINKER_OPTION' => 1, + '_AM_MANGLE_OPTION' => 1, + 'LT_SYS_DLSEARCH_PATH' => 1, + '_LT_AC_TAGCONFIG' => 1, + 'AC_CONFIG_MACRO_DIR' => 1, + 'AC_LIBTOOL_PROG_LD_SHLIBS' => 1, + '_AM_OUTPUT_DEPENDENCY_COMMANDS' => 1, + 'AC_LIBLTDL_CONVENIENCE' => 1, + 'AC_LIBTOOL_WIN32_DLL' => 1, + 'AM_SUBST_NOTMAKE' => 1, + 'AM_SILENT_RULES' => 1, + '_AC_PROG_LIBTOOL' => 1, + 'AC_CHECK_LIBM' => 1, + 'AC_LIBTOOL_DLOPEN' => 1, + 'AC_LIBTOOL_PROG_CC_C_O' => 1, + 'AM_EXTRA_RECURSIVE_TARGETS' => 1, + '_LT_AC_TRY_DLOPEN_SELF' => 1, + '_LT_AC_LANG_F77_CONFIG' => 1, + '_AM_PROG_TAR' => 1, 'LT_PATH_NM' => 1, - 'm4_include' => 1, - 'AC_LIBTOOL_OBJDIR' => 1, - 'PKG_CHECK_VAR' => 1, - '_LT_LINKER_BOILERPLATE' => 1, - '_AM_SET_OPTIONS' => 1, - '_LT_AC_CHECK_DLFCN' => 1, - 'PKG_INSTALLDIR' => 1, - 'AC_LIBTOOL_CXX' => 1, - 'AC_LIBTOOL_LANG_C_CONFIG' => 1, + '_LT_AC_LANG_GCJ_CONFIG' => 1, 'LT_PROG_GO' => 1, - 'AM_SUBST_NOTMAKE' => 1, - 'AU_DEFUN' => 1, - '_LT_AC_FILE_LTDLL_C' => 1, - 'AM_DISABLE_STATIC' => 1, + 'LT_AC_PROG_RC' => 1, + '_LT_AC_LANG_GCJ' => 1, + '_LT_AC_LANG_RC_CONFIG' => 1, + 'AC_LTDL_SYS_DLOPEN_DEPLIBS' => 1, + 'AC_DISABLE_FAST_INSTALL' => 1, + '_LT_AC_LANG_F77' => 1, + 'AC_LIBTOOL_SYS_DYNAMIC_LINKER' => 1, + 'LT_SUPPORTED_TAG' => 1, '_LTDL_SETUP' => 1, - 'AC_LIB_LTDL' => 1, - 'AC_LIBTOOL_PICMODE' => 1, - '_LT_PROG_FC' => 1, - 'AC_LIBTOOL_PROG_LD_SHLIBS' => 1, - 'AC_PROG_NM' => 1, + '_LT_AC_LANG_C_CONFIG' => 1, + 'AC_LTDL_SHLIBEXT' => 1, + 'LTDL_INIT' => 1, '_LT_PATH_TOOL_PREFIX' => 1, - 'AC_ENABLE_STATIC' => 1, - 'AC_LIBTOOL_LANG_RC_CONFIG' => 1, + 'LTVERSION_VERSION' => 1, + '_LT_AC_SYS_COMPILER' => 1, + 'AC_LTDL_ENABLE_INSTALL' => 1, + 'AC_LIBTOOL_COMPILER_OPTION' => 1, + 'LT_SYS_SYMBOL_USCORE' => 1, 'AC_LIBTOOL_SYS_OLD_ARCHIVE' => 1, - 'AC_LIBTOOL_GCJ' => 1, - 'AC_LTDL_DLLIB' => 1, - 'AC_LTDL_PREOPEN' => 1, - 'LT_SYS_MODULE_EXT' => 1, - 'AC_DISABLE_SHARED' => 1, - 'AM_SET_DEPDIR' => 1, - 'AM_SILENT_RULES' => 1, - '_LT_AC_LANG_F77' => 1, - '_LT_AC_LANG_C_CONFIG' => 1, - 'AC_PROG_LD_RELOAD_FLAG' => 1, 'AM_SANITY_CHECK' => 1, - 'LT_OUTPUT' => 1, - 'PKG_NOARCH_INSTALLDIR' => 1, - 'AC_LIBTOOL_SETUP' => 1, - 'LTVERSION_VERSION' => 1, - 'AM_SET_LEADING_DOT' => 1, - 'LT_LIB_DLLOAD' => 1, - 'AC_CHECK_LIBM' => 1, - 'AM_PROG_LIBTOOL' => 1, - 'PKG_CHECK_MODULES' => 1, - 'AM_PROG_CC_C_O' => 1, - 'AM_MAKE_INCLUDE' => 1, - 'AC_PATH_MAGIC' => 1, - 'AM_PROG_INSTALL_SH' => 1, - 'AC_CONFIG_MACRO_DIR' => 1, - 'LT_PROG_GCJ' => 1, - 'AC_LIBTOOL_POSTDEP_PREDEP' => 1, - '_AM_PROG_CC_C_O' => 1, - 'AM_MISSING_PROG' => 1, - 'AC_LIBTOOL_LINKER_OPTION' => 1, - 'AC_PROG_LD_GNU' => 1, - '_LT_AC_LANG_GCJ' => 1, - '_LT_COMPILER_OPTION' => 1, - '_AM_SET_OPTION' => 1, - 'AC_LIBTOOL_PROG_CC_C_O' => 1, - 'AC_LTDL_SHLIBPATH' => 1, - 'LT_PATH_LD' => 1, - 'AC_PATH_TOOL_PREFIX' => 1, - 'include' => 1, - 'AC_LIBTOOL_FC' => 1, - 'AC_DEFUN' => 1, - 'AC_LTDL_SYS_DLOPEN_DEPLIBS' => 1, - 'PKG_PROG_PKG_CONFIG' => 1, - '_LT_LIBOBJ' => 1, + 'LTSUGAR_VERSION' => 1, + 'AM_AUX_DIR_EXPAND' => 1, + 'AC_LTDL_OBJDIR' => 1, + 'LT_LIB_M' => 1, + '_LT_LINKER_OPTION' => 1, 'AC_LIBTOOL_RC' => 1, - 'm4_pattern_forbid' => 1, - 'AC_LIBLTDL_CONVENIENCE' => 1, - 'LT_PROG_RC' => 1, - 'AM_MISSING_HAS_RUN' => 1, - 'AC_PROG_LD' => 1, + '_LT_AC_TAGVAR' => 1, + 'AC_LIBTOOL_LANG_F77_CONFIG' => 1, + 'LT_SYS_DLOPEN_DEPLIBS' => 1, + '_AM_SUBST_NOTMAKE' => 1, + '_AC_AM_CONFIG_HEADER_HOOK' => 1, + 'AC_CONFIG_MACRO_DIR_TRACE' => 1, + '_LT_AC_SHELL_INIT' => 1, + 'LT_INIT' => 1, + '_AM_SET_OPTIONS' => 1, + 'include' => 1, + '_LT_PROG_CXX' => 1, + 'AM_SET_LEADING_DOT' => 1, 'AM_ENABLE_SHARED' => 1, + 'AM_AUTOMAKE_VERSION' => 1, 'LT_WITH_LTDL' => 1, - 'AC_LIBTOOL_PROG_COMPILER_PIC' => 1, '_AM_IF_OPTION' => 1, - 'AC_LIBTOOL_DLOPEN_SELF' => 1, - 'AC_LTDL_SHLIBEXT' => 1, - '_LT_AC_LANG_CXX_CONFIG' => 1, + 'AM_DISABLE_STATIC' => 1, + '_AM_SET_OPTION' => 1, + '_AM_CONFIG_MACRO_DIRS' => 1, + 'AC_PROG_LIBTOOL' => 1, + 'm4_pattern_forbid' => 1, + 'AM_ENABLE_STATIC' => 1, + 'AC_DEFUN_ONCE' => 1, + 'AC_PROG_LD_RELOAD_FLAG' => 1, '_LT_PROG_F77' => 1, - 'AC_LIBTOOL_COMPILER_OPTION' => 1, - '_PKG_SHORT_ERRORS_SUPPORTED' => 1, - 'LT_CONFIG_LTDL_DIR' => 1, - '_AM_PROG_TAR' => 1, - 'AM_EXTRA_RECURSIVE_TARGETS' => 1, - '_LT_PROG_LTMAIN' => 1, - 'AC_LTDL_SYMBOL_USCORE' => 1, - '_LT_AC_PROG_ECHO_BACKSLASH' => 1, - '_AM_OUTPUT_DEPENDENCY_COMMANDS' => 1, - '_LT_AC_TAGVAR' => 1, - 'AC_WITH_LTDL' => 1, - 'AC_LIBLTDL_INSTALLABLE' => 1, - 'AC_LIBTOOL_SYS_LIB_STRIP' => 1, - 'AC_CONFIG_MACRO_DIR_TRACE' => 1, - '_LT_AC_TRY_DLOPEN_SELF' => 1, - 'AC_LIBTOOL_F77' => 1, - 'AM_AUTOMAKE_VERSION' => 1, - 'AC_DEPLIBS_CHECK_METHOD' => 1, - 'AC_LTDL_OBJDIR' => 1, 'AM_PROG_LD' => 1, - 'AC_ENABLE_SHARED' => 1, - '_LT_PROG_CXX' => 1, - 'AC_LIBTOOL_LANG_F77_CONFIG' => 1, '_LT_AC_LOCK' => 1, - 'AC_LIBTOOL_WIN32_DLL' => 1, + '_LT_LIBOBJ' => 1, + 'AC_LIBTOOL_F77' => 1, + '_PKG_SHORT_ERRORS_SUPPORTED' => 1, + 'AM_INIT_AUTOMAKE' => 1, + 'AM_SET_CURRENT_AUTOMAKE_VERSION' => 1, + 'AC_PROG_LD_GNU' => 1, + 'AC_DISABLE_SHARED' => 1, + '_LT_PROG_ECHO_BACKSLASH' => 1, + 'AC_LTDL_PREOPEN' => 1, 'AC_LIBTOOL_SYS_MAX_CMD_LEN' => 1, - 'PKG_CHECK_EXISTS' => 1, - 'AM_AUX_DIR_EXPAND' => 1, - 'LT_SYS_MODULE_PATH' => 1, - '_LT_AC_LANG_RC_CONFIG' => 1, - 'AC_LIBTOOL_CONFIG' => 1, + 'LT_PROG_RC' => 1, + '_LT_LINKER_BOILERPLATE' => 1, + '_m4_warn' => 1, + 'AC_PROG_EGREP' => 1, + 'AM_SET_DEPDIR' => 1, + 'AC_PROG_NM' => 1, + '_LT_DLL_DEF_P' => 1, + 'AC_LIBTOOL_LANG_C_CONFIG' => 1, + 'AM_PROG_INSTALL_SH' => 1, + 'AM_CONDITIONAL' => 1, + 'AC_ENABLE_FAST_INSTALL' => 1, + 'AC_LTDL_DLSYM_USCORE' => 1, + '_LT_COMPILER_BOILERPLATE' => 1, 'AC_LIBTOOL_SYS_HARD_LINK_LOCKS' => 1, + '_LT_AC_FILE_LTDLL_C' => 1, + '_LT_PROG_FC' => 1, + 'AC_LIBTOOL_PROG_COMPILER_PIC' => 1, + 'AM_RUN_LOG' => 1, + 'AM_MAKE_INCLUDE' => 1, + 'AC_LIBTOOL_LANG_GCJ_CONFIG' => 1, + 'LTDL_INSTALLABLE' => 1, + 'LT_OUTPUT' => 1, + 'AC_LIBTOOL_FC' => 1, + 'AC_WITH_LTDL' => 1, + '_AM_PROG_CC_C_O' => 1, + 'AU_DEFUN' => 1, '_LT_WITH_SYSROOT' => 1, + 'AC_LIBLTDL_INSTALLABLE' => 1, + 'AM_MISSING_PROG' => 1, + 'AC_LIBTOOL_LANG_RC_CONFIG' => 1, + 'AM_MISSING_HAS_RUN' => 1, + '_LT_AC_CHECK_DLFCN' => 1, + 'LT_FUNC_ARGZ' => 1, + 'AC_LTDL_DLLIB' => 1, + 'LTDL_CONVENIENCE' => 1, + 'AC_LIBTOOL_CXX' => 1, + 'AM_DEP_TRACK' => 1, '_LT_AC_SYS_LIBPATH_AIX' => 1, - 'LT_FUNC_DLSYM_USCORE' => 1, - 'AC_PROG_EGREP' => 1, - 'AM_RUN_LOG' => 1, - 'LT_CMD_MAX_LEN' => 1, - 'LTSUGAR_VERSION' => 1, - '_LT_AC_SYS_COMPILER' => 1, - 'LT_INIT' => 1, - '_LT_AC_TAGCONFIG' => 1, - 'LTDL_INIT' => 1, + '_LT_AC_LANG_CXX' => 1, + '_LT_PROG_LTMAIN' => 1, 'AC_LIBTOOL_PROG_LD_HARDCODE_LIBPATH' => 1, - 'LT_AC_PROG_RC' => 1, - 'AC_LIBTOOL_SYS_GLOBAL_SYMBOL_PIPE' => 1, - '_LT_AC_LANG_F77_CONFIG' => 1, - '_LT_COMPILER_BOILERPLATE' => 1, - '_AC_PROG_LIBTOOL' => 1, - 'AM_DEP_TRACK' => 1, - 'LT_SUPPORTED_TAG' => 1, - 'AM_CONDITIONAL' => 1, - 'AM_SET_CURRENT_AUTOMAKE_VERSION' => 1 + 'AC_LIBTOOL_GCJ' => 1, + 'LT_AC_PROG_EGREP' => 1, + 'LT_PATH_LD' => 1, + '_LT_AC_LANG_CXX_CONFIG' => 1, + 'PKG_CHECK_MODULES' => 1, + 'm4_pattern_allow' => 1, + '_LT_REQUIRED_DARWIN_CHECKS' => 1, + 'LT_PROG_GCJ' => 1, + '_LT_CC_BASENAME' => 1, + 'AC_LTDL_SYSSEARCHPATH' => 1, + 'AC_PATH_TOOL_PREFIX' => 1 } ], 'Autom4te::Request' ), bless( [ '2', 1, [ '/usr/share/autoconf' ], [ '/usr/share/autoconf/autoconf/autoconf.m4f', 'aclocal.m4', 'configure.ac' ], { - 'AC_CONFIG_LINKS' => 1, - 'AC_REQUIRE_AUX_FILE' => 1, - 'LT_INIT' => 1, + 'AC_FC_SRCEXT' => 1, + 'm4_include' => 1, + 'AC_PROG_LIBTOOL' => 1, + 'm4_pattern_forbid' => 1, '_AM_COND_ELSE' => 1, + 'AC_CANONICAL_TARGET' => 1, + 'AC_CONFIG_FILES' => 1, + 'AM_EXTRA_RECURSIVE_TARGETS' => 1, + 'AC_SUBST' => 1, + 'AC_CANONICAL_BUILD' => 1, + 'AM_XGETTEXT_OPTION' => 1, + 'AC_CANONICAL_HOST' => 1, + 'AM_MAKEFILE_INCLUDE' => 1, + 'AM_PATH_GUILE' => 1, + 'AM_GNU_GETTEXT' => 1, + 'AM_PROG_MOC' => 1, + '_AM_SUBST_NOTMAKE' => 1, + 'AC_INIT' => 1, + 'AC_CONFIG_MACRO_DIR_TRACE' => 1, + 'LT_INIT' => 1, + 'AC_CONFIG_LIBOBJ_DIR' => 1, + 'AC_FC_PP_SRCEXT' => 1, + 'm4_sinclude' => 1, '_LT_AC_TAGCONFIG' => 1, + 'AM_CONDITIONAL' => 1, + 'AC_FC_FREEFORM' => 1, + 'AC_CONFIG_AUX_DIR' => 1, + 'AC_CONFIG_SUBDIRS' => 1, + 'AC_FC_PP_DEFINE' => 1, + 'AM_PROG_MKDIR_P' => 1, + 'include' => 1, + 'AC_CONFIG_LINKS' => 1, + 'AM_AUTOMAKE_VERSION' => 1, + 'AM_SILENT_RULES' => 1, + 'AM_PROG_LIBTOOL' => 1, + 'AC_CANONICAL_SYSTEM' => 1, + 'AM_PROG_CC_C_O' => 1, + 'AM_GNU_GETTEXT_INTL_SUBDIR' => 1, + 'AC_DEFINE_TRACE_LITERAL' => 1, + 'AM_PROG_AR' => 1, + 'm4_pattern_allow' => 1, + 'AC_LIBSOURCE' => 1, + '_AM_MAKEFILE_INCLUDE' => 1, 'AH_OUTPUT' => 1, 'AM_POT_TOOLS' => 1, - 'AC_PROG_LIBTOOL' => 1, - 'AM_SILENT_RULES' => 1, - 'AM_GNU_GETTEXT' => 1, - 'AM_PATH_GUILE' => 1, + 'AC_CONFIG_HEADERS' => 1, + 'IT_PROG_INTLTOOL' => 1, + '_m4_warn' => 1, 'AM_PROG_F77_C_O' => 1, - 'AM_PROG_AR' => 1, 'AC_SUBST_TRACE' => 1, - 'AC_DEFINE_TRACE_LITERAL' => 1, - 'LT_SUPPORTED_TAG' => 1, - 'AC_CONFIG_AUX_DIR' => 1, - 'AM_CONDITIONAL' => 1, - '_m4_warn' => 1, - 'AC_FC_SRCEXT' => 1, - 'AM_PROG_LIBTOOL' => 1, - 'IT_PROG_INTLTOOL' => 1, - 'AC_CONFIG_LIBOBJ_DIR' => 1, 'AM_INIT_AUTOMAKE' => 1, - 'LT_CONFIG_LTDL_DIR' => 1, - 'AM_PROG_CC_C_O' => 1, - 'AM_MAKEFILE_INCLUDE' => 1, - 'AM_EXTRA_RECURSIVE_TARGETS' => 1, - '_AM_MAKEFILE_INCLUDE' => 1, - 'AM_XGETTEXT_OPTION' => 1, - 'AC_INIT' => 1, - 'sinclude' => 1, - 'AC_FC_PP_DEFINE' => 1, - 'AM_AUTOMAKE_VERSION' => 1, - '_AM_COND_IF' => 1, - 'AM_NLS' => 1, - 'AC_CANONICAL_TARGET' => 1, 'AM_PROG_CXX_C_O' => 1, - 'AM_ENABLE_MULTILIB' => 1, - 'include' => 1, + '_AM_COND_IF' => 1, + 'LT_SUPPORTED_TAG' => 1, '_AM_COND_ENDIF' => 1, - 'm4_include' => 1, - 'AC_CONFIG_MACRO_DIR_TRACE' => 1, - 'AC_FC_FREEFORM' => 1, - 'AC_CANONICAL_BUILD' => 1, - 'AM_PROG_FC_C_O' => 1, - 'm4_pattern_allow' => 1, - 'm4_pattern_forbid' => 1, - 'AM_PROG_MKDIR_P' => 1, - 'AC_SUBST' => 1, - 'AM_GNU_GETTEXT_INTL_SUBDIR' => 1, - 'AM_PROG_MOC' => 1, - 'AC_CANONICAL_SYSTEM' => 1, - 'm4_sinclude' => 1, + 'sinclude' => 1, 'AM_MAINTAINER_MODE' => 1, - '_AM_SUBST_NOTMAKE' => 1, - 'AC_FC_PP_SRCEXT' => 1, - 'AC_CONFIG_HEADERS' => 1, - 'AC_LIBSOURCE' => 1, - 'AC_CONFIG_FILES' => 1, - 'AC_CONFIG_SUBDIRS' => 1, - 'AC_CANONICAL_HOST' => 1, - 'GTK_DOC_CHECK' => 1 + 'LT_CONFIG_LTDL_DIR' => 1, + 'AM_NLS' => 1, + 'AC_REQUIRE_AUX_FILE' => 1, + 'AM_ENABLE_MULTILIB' => 1, + 'GTK_DOC_CHECK' => 1, + 'AM_PROG_FC_C_O' => 1 } ], 'Autom4te::Request' ), bless( [ '3', 1, [ '/usr/share/autoconf' ], [ '/usr/share/autoconf/autoconf/autoconf.m4f', 'aclocal.m4', '/usr/share/autoconf/autoconf/trailer.m4', 'configure.ac' ], { - 'AC_CONFIG_HEADERS' => 1, + 'AM_CONDITIONAL' => 1, + 'AC_FC_FREEFORM' => 1, 'AC_FC_PP_SRCEXT' => 1, - '_AM_SUBST_NOTMAKE' => 1, - 'AM_MAINTAINER_MODE' => 1, + '_LT_AC_TAGCONFIG' => 1, 'm4_sinclude' => 1, - 'AC_CANONICAL_HOST' => 1, - 'GTK_DOC_CHECK' => 1, 'AC_CONFIG_SUBDIRS' => 1, - 'AC_CONFIG_FILES' => 1, - 'AC_LIBSOURCE' => 1, + 'AC_CONFIG_AUX_DIR' => 1, + '_AM_SUBST_NOTMAKE' => 1, + 'AC_INIT' => 1, 'AM_PROG_MOC' => 1, + 'AM_GNU_GETTEXT' => 1, + 'AC_CONFIG_LIBOBJ_DIR' => 1, + 'LT_INIT' => 1, + 'AC_CONFIG_MACRO_DIR_TRACE' => 1, 'AM_GNU_GETTEXT_INTL_SUBDIR' => 1, - 'AC_SUBST' => 1, + 'AM_PROG_CC_C_O' => 1, + 'AC_DEFINE_TRACE_LITERAL' => 1, + 'include' => 1, 'AM_PROG_MKDIR_P' => 1, - 'm4_pattern_forbid' => 1, - 'AC_CANONICAL_SYSTEM' => 1, - 'm4_include' => 1, - 'm4_pattern_allow' => 1, - 'AM_PROG_FC_C_O' => 1, - 'AC_CANONICAL_BUILD' => 1, - 'AC_CONFIG_MACRO_DIR_TRACE' => 1, - 'AC_FC_FREEFORM' => 1, - 'AM_NLS' => 1, - 'AC_CANONICAL_TARGET' => 1, - '_AM_COND_IF' => 1, 'AC_FC_PP_DEFINE' => 1, 'AM_AUTOMAKE_VERSION' => 1, - 'sinclude' => 1, - 'AC_INIT' => 1, - '_AM_COND_ENDIF' => 1, - 'include' => 1, - 'AM_PROG_CXX_C_O' => 1, - 'AM_ENABLE_MULTILIB' => 1, - 'AC_CONFIG_LIBOBJ_DIR' => 1, - 'IT_PROG_INTLTOOL' => 1, + 'AM_SILENT_RULES' => 1, + 'AC_CANONICAL_SYSTEM' => 1, 'AM_PROG_LIBTOOL' => 1, + 'AC_CONFIG_LINKS' => 1, + 'AC_CONFIG_FILES' => 1, + 'AC_SUBST' => 1, + 'AM_EXTRA_RECURSIVE_TARGETS' => 1, 'AC_FC_SRCEXT' => 1, - 'AM_XGETTEXT_OPTION' => 1, - '_AM_MAKEFILE_INCLUDE' => 1, + 'm4_include' => 1, + '_AM_COND_ELSE' => 1, + 'm4_pattern_forbid' => 1, + 'AC_CANONICAL_TARGET' => 1, + 'AC_PROG_LIBTOOL' => 1, 'AM_MAKEFILE_INCLUDE' => 1, - 'AM_PROG_CC_C_O' => 1, - 'AM_EXTRA_RECURSIVE_TARGETS' => 1, - 'LT_CONFIG_LTDL_DIR' => 1, + 'AM_PATH_GUILE' => 1, + 'AC_CANONICAL_BUILD' => 1, + 'AC_CANONICAL_HOST' => 1, + 'AM_XGETTEXT_OPTION' => 1, + '_AM_COND_IF' => 1, + 'LT_SUPPORTED_TAG' => 1, 'AM_INIT_AUTOMAKE' => 1, + 'AM_PROG_CXX_C_O' => 1, + 'AM_PROG_F77_C_O' => 1, 'AC_SUBST_TRACE' => 1, + 'AM_PROG_FC_C_O' => 1, + 'GTK_DOC_CHECK' => 1, + 'AM_ENABLE_MULTILIB' => 1, + 'AC_REQUIRE_AUX_FILE' => 1, + 'AM_MAINTAINER_MODE' => 1, + '_AM_COND_ENDIF' => 1, + 'sinclude' => 1, + 'AM_NLS' => 1, + 'LT_CONFIG_LTDL_DIR' => 1, + 'm4_pattern_allow' => 1, 'AM_PROG_AR' => 1, + 'IT_PROG_INTLTOOL' => 1, '_m4_warn' => 1, - 'AM_CONDITIONAL' => 1, - 'AC_CONFIG_AUX_DIR' => 1, - 'AC_DEFINE_TRACE_LITERAL' => 1, - 'LT_SUPPORTED_TAG' => 1, - 'AM_GNU_GETTEXT' => 1, - 'AM_SILENT_RULES' => 1, - 'AM_POT_TOOLS' => 1, - 'AC_PROG_LIBTOOL' => 1, 'AH_OUTPUT' => 1, - 'AM_PROG_F77_C_O' => 1, - 'AM_PATH_GUILE' => 1, - 'AC_REQUIRE_AUX_FILE' => 1, - 'AC_CONFIG_LINKS' => 1, - '_LT_AC_TAGCONFIG' => 1, - 'LT_INIT' => 1, - '_AM_COND_ELSE' => 1 + '_AM_MAKEFILE_INCLUDE' => 1, + 'AC_LIBSOURCE' => 1, + 'AC_CONFIG_HEADERS' => 1, + 'AM_POT_TOOLS' => 1 } ], 'Autom4te::Request' ) ); Index: trunk/npstat/nm/SpecialFunctions.cc =================================================================== --- trunk/npstat/nm/SpecialFunctions.cc (revision 816) +++ trunk/npstat/nm/SpecialFunctions.cc (revision 817) @@ -1,2135 +1,2137 @@ // Most of the code below is cherry-picked from the Cephes library // written by Stephen L. Moshier #include #include #include #include #include #include #include "npstat/nm/SpecialFunctions.hh" #include "npstat/nm/MathUtils.hh" #include "npstat/nm/Matrix.hh" #include "npstat/rng/permutation.hh" #define SQRTPIL 1.772453850905516027298167L #define SQRTTWOPIL 2.506628274631000502415765L #define SQTPI 2.50662827463100050242 #define LOGPI 1.14472988584940017414 #define MAXGAM 171.624376956302725 #define MAXSTIR 143.01608 #define MAXLOG 709.782712893383996732224 #define MINLOG (-708.396418532264106224) #define MACHEP 1.2e-16 #define MAXLGM 2.556348e305 #define MAXNUM 1.79769313486231570815e308 #define LS2PI 0.91893853320467274178 #define GAUSS_MAX_SIGMA 37.615 static double igam(double a, double x); static double invgauss(double x); static double polevl(const double x, const double *p, int i) { double ans = *p++; do { ans = ans * x + *p++; } while (--i); return ans; } static double p1evl(double x, const double *p, const int N) { double ans = x + *p++; int i = N - 1; do { ans = ans * x + *p++; } while (--i); return ans; } /* Gamma function computed by Stirling's formula. * The polynomial STIR is valid for 33 <= x <= 172. */ static double stirf(const double x) { static const double STIR[5] = { 7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3, 3.47222221605458667310E-3, 8.33333333333482257126E-2, }; double y, w; w = 1.0/x; w = 1.0 + w * polevl( w, STIR, 4 ); y = exp(x); if( x > MAXSTIR ) { // Avoid overflow in pow() const double v = pow( x, 0.5 * x - 0.25 ); y = v * (v / y); } else { y = pow( x, x - 0.5 ) / y; } return SQTPI * y * w; } // Logarithm of the gamma function static double lgam(double x, int* sgngam_out=0) { int sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; double p, q, u, w, z; int i; static const double A[] = { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4, -2.77777777730099687205E-3, 8.33333333333331927722E-2 }; static const double B[] = { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5, -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }; static const double C[] = { -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5, -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }; if( x < -34.0 ) { q = -x; w = lgam(q, &sgngam); if (sgngam_out) *sgngam_out = sgngam; p = floor(q); if( p == q ) goto loverf; i = static_cast(p); if( (i & 1) == 0 ) sgngam = -1; else sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; z = q - p; if( z > 0.5 ) { p += 1.0; z = p - q; } z = q * sin( M_PI * z ); if( z == 0.0 ) goto loverf; // z = log(PI) - log( z ) - w; z = LOGPI - log( z ) - w; return( z ); } if( x < 13.0 ) { z = 1.0; p = 0.0; u = x; while( u >= 3.0 ) { p -= 1.0; u = x + p; z *= u; } while( u < 2.0 ) { if( u == 0.0 ) goto loverf; z /= u; p += 1.0; u = x + p; } if( z < 0.0 ) { sgngam = -1; z = -z; } else sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; if( u == 2.0 ) return( log(z) ); p -= 2.0; x = x + p; p = x * polevl( x, B, 5 ) / p1evl( x, C, 6); return( log(z) + p ); } if( x > MAXLGM ) { loverf: assert(!"Overflow in lgam"); return( sgngam * MAXNUM ); } q = ( x - 0.5 ) * log(x) - x + LS2PI; if( x > 1.0e8 ) return( q ); p = 1.0/(x*x); if( x >= 1000.0 ) q += (( 7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) *p + 0.0833333333333333333333) / x; else q += polevl( p, A, 4 ) / x; return( q ); } // Complementary incomplete gamma static double igamc(double a, double x) { const double big = 4.503599627370496e15; const double biginv = 2.22044604925031308085e-16; double ans, ax, c, yc, r, t, y, z; double pk, pkm1, pkm2, qk, qkm1, qkm2; if( (x <= 0) || ( a <= 0) ) return( 1.0 ); if( (x < 1.0) || (x < a) ) return( 1.0 - igam(a,x) ); // Compute x**a * exp(-x) / gamma(a) ax = a * log(x) - x - lgam(a); if( ax < -MAXLOG ) { // mtherr( "igamc", UNDERFLOW ); // assert(!"Underflow in igamc"); return 0.0; } ax = exp(ax); // continued fraction y = 1.0 - a; z = x + y + 1.0; c = 0.0; pkm2 = 1.0; qkm2 = x; pkm1 = x + 1.0; qkm1 = z * x; ans = pkm1/qkm1; do { c += 1.0; y += 1.0; z += 2.0; yc = y * c; pk = pkm1 * z - pkm2 * yc; qk = qkm1 * z - qkm2 * yc; if( qk != 0 ) { r = pk/qk; t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( fabs(pk) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } } while( t > MACHEP ); return( ans * ax ); } // Incomplete gamma static double igam(double a, double x) { double ans, ax, c, r; if( (x <= 0) || ( a <= 0) ) return( 0.0 ); if( (x > 1.0) && (x > a ) ) return( 1.0 - igamc(a,x) ); // Compute x**a * exp(-x) / gamma(a) ax = a * log(x) - x - lgam(a); if( ax < -MAXLOG ) { // mtherr( "igam", UNDERFLOW ); assert(!"Underflow in igam"); return( 0.0 ); } ax = exp(ax); // power series r = a; c = 1.0; ans = 1.0; do { r += 1.0; c *= x/r; ans += c; } while( c/ans > MACHEP ); return( ans * ax/a ); } // Inverse of complemented incomplete gamma integral. // Work only for y0 <= 0.5. static double igami(double a, double y0) { double x0, x1, x, yl, yh, y, d, lgm, dithresh; int i, dir; assert(y0 <= 0.5); /* bound the solution */ x0 = MAXNUM; yl = 0; x1 = 0; yh = 1.0; dithresh = 5.0 * MACHEP; /* approximation to inverse function */ d = 1.0/(9.0*a); y = ( 1.0 - d - invgauss(y0) * sqrt(d) ); x = a * y * y * y; lgm = lgam(a); for( i=0; i<10; i++ ) { if( x > x0 || x < x1 ) goto ihalve; y = igamc(a,x); if( y < yl || y > yh ) goto ihalve; if( y < y0 ) { x0 = x; yl = y; } else { x1 = x; yh = y; } /* compute the derivative of the function at this point */ d = (a - 1.0) * log(x) - x - lgm; if( d < -MAXLOG ) goto ihalve; d = -exp(d); /* compute the step to the next approximation of x */ d = (y - y0)/d; if( fabs(d/x) < MACHEP ) goto done; x = x - d; } /* Resort to interval halving if Newton iteration did not converge. */ ihalve: d = 0.0625; if( x0 == MAXNUM ) { if( x <= 0.0 ) x = 1.0; while( x0 == MAXNUM ) { x = (1.0 + d) * x; y = igamc( a, x ); if( y < y0 ) { x0 = x; yl = y; break; } d = d + d; } } d = 0.5; dir = 0; for( i=0; i<400; i++ ) { x = x1 + d * (x0 - x1); y = igamc( a, x ); lgm = (x0 - x1)/(x1 + x0); if( fabs(lgm) < dithresh ) break; lgm = (y - y0)/y0; if( fabs(lgm) < dithresh ) break; if( x <= 0.0 ) break; if( y >= y0 ) { x1 = x; yh = y; if( dir < 0 ) { dir = 0; d = 0.5; } else if( dir > 1 ) d = 0.5 * d + 0.5; else d = (y0 - yl)/(yh - yl); dir += 1; } else { x0 = x; yl = y; if( dir > 0 ) { dir = 0; d = 0.5; } else if( dir < -1 ) d = 0.5 * d; else d = (y0 - yl)/(yh - yl); dir -= 1; } } if( x == 0.0 ) // mtherr( "igami", UNDERFLOW ); assert(!"Underflow in igami"); done: return( x ); } static double invgauss(const double P) { assert(P > 0.0); assert(P < 1.0); // Translated from PPND16 algorithm of Wichura (originally in Fortran). // This was not taken from Cephes. static const double ZERO = 0., ONE = 1., HALF = 0.5, SPLIT1 = 0.425, SPLIT2 = 5., CONST1 = 0.180625, CONST2 = 1.6; static const double A0 = 3.3871328727963666080, A1 = 1.3314166789178437745E+2, A2 = 1.9715909503065514427E+3, A3 = 1.3731693765509461125E+4, A4 = 4.5921953931549871457E+4, A5 = 6.7265770927008700853E+4, A6 = 3.3430575583588128105E+4, A7 = 2.5090809287301226727E+3, B1 = 4.2313330701600911252E+1, B2 = 6.8718700749205790830E+2, B3 = 5.3941960214247511077E+3, B4 = 2.1213794301586595867E+4, B5 = 3.9307895800092710610E+4, B6 = 2.8729085735721942674E+4, B7 = 5.2264952788528545610E+3; static const double C0 = 1.42343711074968357734, C1 = 4.63033784615654529590, C2 = 5.76949722146069140550, C3 = 3.64784832476320460504, C4 = 1.27045825245236838258, C5 = 2.41780725177450611770E-1, C6 = 2.27238449892691845833E-2, C7 = 7.74545014278341407640E-4, D1 = 2.05319162663775882187, D2 = 1.67638483018380384940, D3 = 6.89767334985100004550E-1, D4 = 1.48103976427480074590E-1, D5 = 1.51986665636164571966E-2, D6 = 5.47593808499534494600E-4, D7 = 1.05075007164441684324E-9; static const double E0 = 6.65790464350110377720, E1 = 5.46378491116411436990, E2 = 1.78482653991729133580, E3 = 2.96560571828504891230E-1, E4 = 2.65321895265761230930E-2, E5 = 1.24266094738807843860E-3, E6 = 2.71155556874348757815E-5, E7 = 2.01033439929228813265E-7, F1 = 5.99832206555887937690E-1, F2 = 1.36929880922735805310E-1, F3 = 1.48753612908506148525E-2, F4 = 7.86869131145613259100E-4, F5 = 1.84631831751005468180E-5, F6 = 1.42151175831644588870E-7, F7 = 2.04426310338993978564E-15; const double Q = P - HALF; double R, PPND16; if (fabs(Q) <= SPLIT1) { R = CONST1 - Q * Q; PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) * R + A2) * R + A1) * R + A0) / (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) * R + B2) * R + B1) * R + ONE); } else { if (Q < ZERO) R = P; else R = ONE - P; R = sqrt(-log(R)); if (R <= SPLIT2) { R = R - CONST2; PPND16 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3) * R + C2) * R + C1) * R + C0) / (((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) * R + D2) * R + D1) * R + ONE); } else { R = R - SPLIT2 ; PPND16 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3) * R + E2) * R + E1) * R + E0) / (((((((F7 * R + F6) * R + F5) * R + F4) * R + F3) * R + F2) * R + F1) * R + ONE); } if (Q < ZERO) PPND16 = -PPND16; } return PPND16; } /* Continued fraction expansion #1 * for incomplete beta integral */ static double incbcf( double a, double b, double x ) { static const double big = 4.503599627370496e15; static const double biginv = 2.22044604925031308085e-16; double xk, pk, pkm1, pkm2, qk, qkm1, qkm2; double k1, k2, k3, k4, k5, k6, k7, k8; double r, t, ans, thresh; int n; k1 = a; k2 = a + b; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = b - 1.0; k7 = k4; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; ans = 1.0; r = 1.0; n = 0; thresh = 3.0 * MACHEP; do { xk = -( x * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( x * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; if( t < thresh ) goto cdone; k1 += 1.0; k2 += 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 -= 1.0; k7 += 2.0; k8 += 2.0; if( (fabs(qk) + fabs(pk)) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } if( (fabs(qk) < biginv) || (fabs(pk) < biginv) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 300 ); cdone: return(ans); } /* Continued fraction expansion #2 * for incomplete beta integral */ static double incbd( double a, double b, double x ) { static const double big = 4.503599627370496e15; static const double biginv = 2.22044604925031308085e-16; double xk, pk, pkm1, pkm2, qk, qkm1, qkm2; double k1, k2, k3, k4, k5, k6, k7, k8; double r, t, ans, z, thresh; int n; k1 = a; k2 = b - 1.0; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = a + b; k7 = a + 1.0;; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; z = x / (1.0-x); ans = 1.0; r = 1.0; n = 0; thresh = 3.0 * MACHEP; do { xk = -( z * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( z * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; if( t < thresh ) goto cdone; k1 += 1.0; k2 -= 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 += 1.0; k7 += 2.0; k8 += 2.0; if( (fabs(qk) + fabs(pk)) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } if( (fabs(qk) < biginv) || (fabs(pk) < biginv) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 300 ); cdone: return(ans); } /* Power series for incomplete beta integral. Use when b*x is small and x not too close to 1. */ static double pseries( double a, double b, double x ) { double s, t, u, v, n, t1, z, ai; ai = 1.0 / a; u = (1.0 - b) * x; v = u / (a + 1.0); t1 = v; t = u; n = 2.0; s = 0.0; z = MACHEP * ai; while( fabs(v) > z ) { u = (n - b) * x / n; t *= u; v = t / (a + n); s += v; n += 1.0; } s += t1; s += ai; u = a * log(x); if( (a+b) < MAXGAM && fabs(u) < MAXLOG ) { t = npstat::Gamma(a+b)/(npstat::Gamma(a)*npstat::Gamma(b)); s = s * t * pow(x,a); } else { t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s); if( t < MINLOG ) s = 0.0; else s = exp(t); } return(s); } /* Incomplete beta function */ static double incbet( double aa, double bb, double xx ) { double a, b, t, x, xc, w, y; int flag; if( aa <= 0.0 || bb <= 0.0 ) goto domerr; if( (xx <= 0.0) || ( xx >= 1.0) ) { if( xx == 0.0 ) return(0.0); if( xx == 1.0 ) return( 1.0 ); domerr: // mtherr( "incbet", DOMAIN ); assert(!"Domain error in incbet"); return( 0.0 ); } flag = 0; if( (bb * xx) <= 1.0 && xx <= 0.95) { t = pseries(aa, bb, xx); goto done; } w = 1.0 - xx; /* Reverse a and b if x is greater than the mean. */ if( xx > (aa/(aa+bb)) ) { flag = 1; a = bb; b = aa; xc = xx; x = w; } else { a = aa; b = bb; xc = w; x = xx; } if( flag == 1 && (b * x) <= 1.0 && x <= 0.95) { t = pseries(a, b, x); goto done; } /* Choose expansion for better convergence. */ y = x * (a+b-2.0) - (a-1.0); if( y < 0.0 ) w = incbcf( a, b, x ); else w = incbd( a, b, x ) / xc; /* Multiply w by the factor a b _ _ _ x (1-x) | (a+b) / ( a | (a) | (b) ) . */ y = a * log(x); t = b * log(xc); if( (a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG ) { t = pow(xc,b); t *= pow(x,a); t /= a; t *= w; t *= npstat::Gamma(a+b) / (npstat::Gamma(a) * npstat::Gamma(b)); goto done; } /* Resort to logarithms. */ y += t + lgam(a+b) - lgam(a) - lgam(b); y += log(w/a); if( y < MINLOG ) t = 0.0; else t = exp(y); done: if( flag == 1 ) { if( t <= MACHEP ) t = 1.0 - MACHEP; else t = 1.0 - t; } return( t ); } /* Inverse incomplete beta function */ static double incbi( double aa, double bb, double yy0 ) { double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; int i, rflg, dir, nflg; i = 0; if( yy0 <= 0.0 ) return(0.0); if( yy0 >= 1.0 ) return(1.0); x0 = 0.0; yl = 0.0; x1 = 1.0; yh = 1.0; nflg = 0; if( aa <= 1.0 || bb <= 1.0 ) { dithresh = 1.0e-6; rflg = 0; a = aa; b = bb; y0 = yy0; x = a/(a+b); y = incbet( a, b, x ); goto ihalve; } else { dithresh = 1.0e-4; } /* approximation to inverse function */ yp = -npstat::inverseGaussCdf(yy0); if( yy0 > 0.5 ) { rflg = 1; a = bb; b = aa; y0 = 1.0 - yy0; yp = -yp; } else { rflg = 0; a = aa; b = bb; y0 = yy0; } lgm = (yp * yp - 3.0)/6.0; x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) ); d = yp * sqrt( x + lgm ) / x - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) ) * (lgm + 5.0/6.0 - 2.0/(3.0*x)); d = 2.0 * d; if( d < MINLOG ) { x = 1.0; goto under; } x = a/( a + b * exp(d) ); y = incbet( a, b, x ); yp = (y - y0)/y0; if( fabs(yp) < 0.2 ) goto newt; /* Resort to interval halving if not close enough. */ ihalve: dir = 0; di = 0.5; for( i=0; i<100; i++ ) { if( i != 0 ) { x = x0 + di * (x1 - x0); if( x == 1.0 ) x = 1.0 - MACHEP; if( x == 0.0 ) { di = 0.5; x = x0 + di * (x1 - x0); if( x == 0.0 ) goto under; } y = incbet( a, b, x ); yp = (x1 - x0)/(x1 + x0); if( fabs(yp) < dithresh ) goto newt; yp = (y-y0)/y0; if( fabs(yp) < dithresh ) goto newt; } if( y < y0 ) { x0 = x; yl = y; if( dir < 0 ) { dir = 0; di = 0.5; } else if( dir > 3 ) di = 1.0 - (1.0 - di) * (1.0 - di); else if( dir > 1 ) di = 0.5 * di + 0.5; else di = (y0 - y)/(yh - yl); dir += 1; if( x0 > 0.75 ) { if( rflg == 1 ) { rflg = 0; a = aa; b = bb; y0 = yy0; } else { rflg = 1; a = bb; b = aa; y0 = 1.0 - yy0; } x = 1.0 - x; y = incbet( a, b, x ); x0 = 0.0; yl = 0.0; x1 = 1.0; yh = 1.0; goto ihalve; } } else { x1 = x; if( rflg == 1 && x1 < MACHEP ) { x = 0.0; goto done; } yh = y; if( dir > 0 ) { dir = 0; di = 0.5; } else if( dir < -3 ) di = di * di; else if( dir < -1 ) di = 0.5 * di; else di = (y - y0)/(yh - yl); dir -= 1; } } // Partial loss of precision. Hm... Will ignore for now - igv // mtherr( "incbi", PLOSS ); if( x0 >= 1.0 ) { x = 1.0 - MACHEP; goto done; } if( x <= 0.0 ) { under: // mtherr( "incbi", UNDERFLOW ); assert(!"Underflow in incbi"); x = 0.0; goto done; } newt: if( nflg ) goto done; nflg = 1; lgm = lgam(a+b) - lgam(a) - lgam(b); for( i=0; i<8; i++ ) { /* Compute the function at this point. */ if( i != 0 ) y = incbet(a,b,x); if( y < yl ) { x = x0; y = yl; } else if( y > yh ) { x = x1; y = yh; } else if( y < y0 ) { x0 = x; yl = y; } else { x1 = x; yh = y; } if( x == 1.0 || x == 0.0 ) break; /* Compute the derivative of the function at this point. */ d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm; if( d < MINLOG ) goto done; if( d > MAXLOG ) break; d = exp(d); /* Compute the step to the next approximation of x. */ d = (y - y0)/d; xt = x - d; if( xt <= x0 ) { y = (x - x0) / (x1 - x0); xt = x0 + 0.5 * y * (x - x0); if( xt <= 0.0 ) break; } if( xt >= x1 ) { y = (x1 - x) / (x1 - x0); xt = x1 - 0.5 * y * (x1 - x); if( xt >= 1.0 ) break; } x = xt; if( fabs(d/x) < 128.0 * MACHEP ) goto done; } /* Did not converge. */ dithresh = 256.0 * MACHEP; goto ihalve; done: if( rflg ) { if( x <= MACHEP ) x = 1.0 - MACHEP; else x = 1.0 - x; } return( x ); } static void buildGaussDeriCoeffs(const unsigned order, std::vector* coeffs) { assert(coeffs); const unsigned dim = order + 1U; coeffs->resize(dim); if (order == 0U) { (*coeffs)[0] = 1LL; return; } npstat::Matrix deriOp(dim, dim, 0); for (unsigned row=0; row(row+1U); for (unsigned row=1; row& result = deriOp.pow(order); for (unsigned row=0; row= 0 // px = Phi( x ) (cumulative std. normal) // pxs = Phi( x * sqrt( ( 1 - rho ) / ( 1 + rho ) ) ) // Phi2diag() should not be used directly but Phi2() instead. static double Phi2diag( const double x, const double a, const double px, const double pxs ) { if( a <= 0.0 ) return px; // rho == 1 if( a >= 1.0 ) return px * px; // rho == 0 double b = 2.0 - a, sqrt_ab = sqrt( a * b ); // b == 1 + rho const double c1 = 6.36619772367581343e-1; // == 2 / PI const double c2 = 1.25331413731550025; // == sqrt( PI / 2 ) const double c3 = 1.57079632679489662; // == PI / 2 const double c4 = 1.591549430918953358e-1; // == 1 / ( 2 * PI ) // avoid asin() for values close to 1 (vertical tangent) double asr = ( a > 0.1 ? asin( 1.0 - a ) : acos( sqrt_ab ) ); // return upper bound if absolute error within double precision double comp = px * pxs; if( comp * ( 1.0 - a - c1 * asr ) < 5e-17 ) return b * comp; // initialize odd and even terms for recursion double tmp = c2 * x; double alpha = a * x * x / b; double a_even = -tmp * a; double a_odd = -sqrt_ab * alpha; double beta = x * x; double b_even = tmp * sqrt_ab; double b_odd = sqrt_ab * beta; double delta = 2.0 * x * x / b; double d_even = ( 1.0 - a ) * c3 - asr; double d_odd = tmp * ( sqrt_ab - a ); // recursion; update odd and even terms in each step double res = 0.0, res_new = d_even + d_odd; int k = 2; while( res != res_new ) { d_even = ( a_odd + b_odd + delta * d_even ) / k; a_even *= alpha / k; b_even *= beta / k; k++; a_odd *= alpha / k; b_odd *= beta / k; d_odd = ( a_even + b_even + delta * d_odd ) / k; k++; res = res_new; res_new += d_even + d_odd; } // transform; check against upper and lower bounds res *= exp( -x * x / b ) * c4; return std::max( ( 1.0 + c1 * asr ) * comp, b * comp - std::max( 0.0, res ) ); } static double squarethis(const double v) { return v * v; } // Auxiliary function Phi2help() will be called (twice, with x and y interchanged) by Phi2() // In particular, certain pathological cases will have been dealt with in Phi2() before // Axis is transformed into diagonal, and Phi2diag() is called static double Phi2help( const double x, const double y, const double rho ) { // note: case x == y == 0.0 is dealt with in Phi2() if( fabs(x) <= DBL_MIN ) return ( y >= 0.0 ? 0.0 : 0.5 ); double s = sqrt( ( 1.0 - rho ) * ( 1.0 + rho ) ); double a = 0.0, b1 = -fabs( x ), b2 = 0.0; // avoid cancellation by treating the cases rho -> +-1 // separately (cutoff at |0.99| might be optimized); if( rho > 0.99 ) { double tmp = sqrt( ( 1.0 - rho ) / ( 1.0 + rho ) ); b2 = -fabs( ( x - y ) / s - x * tmp ); a = squarethis( ( x - y ) / x / s - tmp ); } else if( rho < -0.99 ) { double tmp = sqrt( ( 1.0 + rho ) / ( 1.0 - rho ) ); b2 = -fabs( ( x + y ) / s - x * tmp ); a = squarethis( ( x + y ) / x / s - tmp ); } else { b2 = -fabs( rho * x - y ) / s; a = squarethis( b2 / x ); } // Phi(): cumulative std. normal (has to be defined elsewhere) double p1 = Phi( b1 ), p2 = Phi( b2 ); // reduction to the diagonal double q = 0.0; if( a <= 1.0 ) q = 0.5 * Phi2diag( b1, 2.0 * a / ( 1.0 + a ), p1, p2 ); else q = p1 * p2 - 0.5 * Phi2diag( b2, 2.0 / ( 1.0 + a ), p2, p1 ); // finally, transformation according to conditions on x, y, rho int c1 = ( y / x >= rho ), c2 = ( x < 0.0 ), c3 = c2 && ( y >= 0.0 ); return ( c1 && c3 ? q - 0.5 : c1 && c2 ? q : c1 ? 0.5 - p1 + q : c3 ? p1 - q - 0.5 : c2 ? p1 - q : 0.5 - q ); } // Evaluation of the cumulative bivariate normal distribution // for arbitrary (x,y) and correlation rho. This function is // described in the article "Recursive Numerical Evaluation of // the Cumulative Bivariate Normal Distribution" by Christian Meyer, // Journal of Statistical Software, Vol. 52, Issue 10, Mar 2013. // // It appears not to work very well for large rho (for example, // Phi2(-0, 0.5, -0.98865243049242746) is wrong), so for large // values of rho implementation by Genz is used (which has its // own problems). // static double Phi2( const double x, const double y, const double rho ) { // special case |rho| == 1 => lower or upper Frechet bound // Phi(): cumulative std. normal (has to be defined elsewhere) if( ( 1.0 - rho ) * ( 1.0 + rho ) <= 0.0 ) { if( rho > 0.0 ) return Phi( std::min( x, y ) ); else return std::max( 0.0, std::min( 1.0, Phi( x ) + Phi( y ) - 1.0 ) ); } // special case x == y == 0, avoids complications in Phi2help() if( x == 0.0 && y == 0.0 ) { if( rho > 0.0 ) return Phi2diag( 0.0, 1.0 - rho, 0.5, 0.5 ); else return 0.5 - Phi2diag( 0.0, 1.0 + rho, 0.5, 0.5 ); } // standard case by reduction formula const double result = std::max( 0.0, std::min( 1.0, Phi2help( x, y, rho ) + Phi2help( y, x, rho ) ) ); // std::cout.precision(17); // std::cout << "Phi2(" << x << ", " << y << ", " << rho << ") = " << result << std::endl; return result; } /* * A function for computing bivariate normal probabilities. * * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenz@wsu.edu * * This function is based on the method described by * Drezner, Z and G.O. Wesolowsky, (1989), * On the computation of the bivariate normal integral, * Journal of Statist. Comput. Simul. 35, pp. 101-107, * with major modifications for double precision, and for |R| close to 1. * * BVND calculates the probability that X > DH and Y > DK. * Note: Prob( X < DH, Y < DK ) = BVND( -DH, -DK, R ). * * Parameters * * DH DOUBLE PRECISION, integration limit * DK DOUBLE PRECISION, integration limit * R DOUBLE PRECISION, correlation coefficient * * Translated from Fortran by igv (July 2014). */ static double bvnd(const double DH, const double DK, const double R) { const double TWOPI = 2.0*M_PI; // Gauss-Legendre points static const double X[11][4] = { {0., 0. , 0. , 0. }, {0., -0.9324695142031522, -0.9815606342467191, -0.9931285991850949}, {0., -0.6612093864662647, -0.9041172563704750, -0.9639719272779138}, {0., -0.2386191860831970, -0.7699026741943050, -0.9122344282513259}, {0., 0. , -0.5873179542866171, -0.8391169718222188}, {0., 0. , -0.3678314989981802, -0.7463319064601508}, {0., 0. , -0.1252334085114692, -0.6360536807265150}, {0., 0. , 0. , -0.5108670019508271}, {0., 0. , 0. , -0.3737060887154196}, {0., 0. , 0. , -0.2277858511416451}, {0., 0. , 0. , -0.7652652113349733} }; // Gauss-Legendre weights static const double W[11][4] = { {0., 0. , 0. , 0. }, {0., 0.1713244923791705, 0.04717533638651177, 0.01761400713915212}, {0., 0.3607615730481384, 0.1069393259953183 , 0.04060142980038694}, {0., 0.4679139345726904, 0.1600783285433464 , 0.06267204833410906}, {0., 0. , 0.2031674267230659 , 0.08327674157670475}, {0., 0. , 0.2334925365383547 , 0.1019301198172404 }, {0., 0. , 0.2491470458134029 , 0.1181945319615184 }, {0., 0. , 0. , 0.1316886384491766 }, {0., 0. , 0. , 0.1420961093183821 }, {0., 0. , 0. , 0.1491729864726037 }, {0., 0. , 0. , 0.1527533871307259 } }; int I, IS, LG, NG; double AS, A, B, C, D, RS, XS, BVN; double SN, ASR, H, K, BS, HS, HK; if ( fabs(R) < 0.3 ) { NG = 1; LG = 3; } else if ( fabs(R) < 0.75 ) { NG = 2; LG = 6; } else { NG = 3; LG = 10; } H = DH; K = DK; HK = H*K; BVN = 0.0; if ( fabs(R) < 0.925 ) { if ( fabs(R) > 0.0 ) { HS = ( H*H + K*K )/2.0; ASR = asin(R); for (I = 1; I <= LG; ++I) { for (IS = -1; IS <= 1; IS += 2) { SN = sin( ASR*( IS*X[I][NG] + 1.0 )/2.0 ); BVN += W[I][NG]*exp( ( SN*HK-HS )/( 1.0-SN*SN ) ); } } BVN = BVN*ASR/( 2.0*TWOPI ); } BVN += Phi(-H)*Phi(-K); } else { if ( R < 0.0 ) { K = -K; HK = -HK; } if ( fabs(R) < 1.0 ) { AS = ( 1 - R )*( 1 + R ); A = sqrt(AS); BS = pow(H - K, 2); C = ( 4 - HK )/8.0; D = ( 12 - HK )/16.0; ASR = -( BS/AS + HK )/2.0; if ( ASR > -100.0 ) BVN = A*exp(ASR)*( 1 - C*( BS - AS )*( 1 - D*BS/5 )/3 + C*D*AS*AS/5 ); if ( -HK < 100.0 ) { B = sqrt(BS); BVN -= exp( -HK/2 )*SQTPI*Phi(-B/A)*B *( 1 - C*BS*( 1 - D*BS/5 )/3 ); } A /= 2.0; for (I = 1; I <= LG; ++I) { for (IS = -1; IS <= 1; IS += 2) { XS = pow( A*( IS*X[I][NG] + 1 ), 2); RS = sqrt( 1 - XS ); ASR = -( BS/XS + HK )/2; if ( ASR > -100 ) BVN += A*W[I][NG]*exp( ASR ) *( exp( -HK*XS/( 2*pow(1 + RS, 2) ) )/RS - ( 1 + C*XS*( 1 + D*XS ) ) ); } } BVN = -BVN/TWOPI; } if ( R > 0.0 ) BVN += Phi( -std::max( H, K ) ); else { BVN = -BVN; if ( K > H ) { if ( H < 0.0 ) BVN += Phi(K) - Phi(H); else BVN += Phi(-H) - Phi(-K); } } } return BVN; } namespace npstat { double inverseGaussCdf(const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseGaussCdf: argument outside of [0, 1] interval"); if (x == 1.0) return GAUSS_MAX_SIGMA; else if (x == 0.0) return -GAUSS_MAX_SIGMA; else return invgauss(x); } double incompleteBeta(const double a, const double b, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::incompleteBeta: argument outside of [0, 1] interval"); if (!(a > 0.0 && b > 0.0)) throw std::invalid_argument( "In npstat::incompleteBeta: parameters must be positive"); return incbet(a, b, x); } double inverseIncompleteBeta(const double a, const double b, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteBeta: " "argument outside of [0, 1] interval"); if (!(a > 0.0 && b > 0.0)) throw std::invalid_argument( "In npstat::inverseIncompleteBeta: parameters must be positive"); return incbi(a, b, x); } double Gamma(double x) { if (x <= 0.0) throw std::domain_error( "In npstat::Gamma: argument must be positive"); // The Stirling formula overflows for x >= MAXGAM if (x >= MAXGAM) throw std::overflow_error( "In npstat::Gamma: argument is too large"); const unsigned ix = x; if (ix && x == static_cast(ix)) return ldfactorial(ix - 1U); if (x == static_cast(ix) + 0.5) return ldfactorial(2U*ix)/ldfactorial(ix)/powl(4.0L,ix)*SQRTPIL; static const double P[] = { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }; static const double Q[] = { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2, 3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 }; double p, z = 1.0, q = x; if (q > 33.0) return stirf(q); while( x >= 3.0 ) { x -= 1.0; z *= x; } while( x < 2.0 ) { if( x < 1.e-9 ) return( z/((1.0 + 0.5772156649015329 * x) * x) ); z /= x; x += 1.0; } if( x == 2.0 ) return(z); x -= 2.0; p = polevl( x, P, 6 ); q = polevl( x, Q, 7 ); return( z * p / q ); } double incompleteGamma(const double a, const double x) { return igam(a, x); } double incompleteGammaC(const double a, const double x) { return igamc(a, x); } double inverseIncompleteGammaC(const double a, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteGammaC: " "argument outside of [0, 1] interval"); if (x <= 0.5) return igami(a, x); else return inverseIncompleteGamma(a, 1.0 - x); } double inverseIncompleteGamma(const double a, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteGamma: " "argument outside of [0, 1] interval"); - if (x >= 0.5) + if (x == 0.0) + return 0.0; + else if (x >= 0.5) return igami(a, 1.0 - x); else { const double targetEps = 8.0*MACHEP; double xmin = 0.0; double xmax = igami(a, 0.5); while ((xmax - xmin)/(xmax + xmin + 1.0) > targetEps) { const double xtry = (xmax + xmin)/2.0; if (igam(a, xtry) > x) xmax = xtry; else xmin = xtry; } return (xmax + xmin)/2.0; } } long double dawsonIntegral(const long double x_in) { const unsigned deg = 20U; const long double x = fabsl(x_in); long double result = 0.0L; if (x == 0.0L) ; else if (x <= 1.0L) { const long double epsilon = LDBL_EPSILON/2.0L; const long double twoxsq = -2.0L*x*x; long double term = 1.0L; long double dfact = 1.0L; long double xn = 1.0L; unsigned y = 0; do { result += term; y += 1U; dfact *= (2U*y + 1U); xn *= twoxsq; term = xn/dfact; } while (fabsl(term) > epsilon * fabsl(result)); result *= x; } else if (x <= 1.75L) { static const long double coeffs[deg + 1U] = { +4.563960711239483142081e-1L, -9.268566100670767619861e-2L, -7.334392170021420220239e-3L, +3.379523740404396755124e-3L, -3.085898448678595090813e-4L, -1.519846724619319512311e-5L, +4.903955822454009397182e-6L, -2.106910538629224721838e-7L, -2.930676220603996192089e-8L, +3.326790071774057337673e-9L, +3.335593457695769191326e-11L, -2.279104036721012221982e-11L, +7.877561633156348806091e-13L, +9.173158167107974472228e-14L, -7.341175636102869400671e-15L, -1.763370444125849029511e-16L, +3.792946298506435014290e-17L, -4.251969162435936250171e-19L, -1.358295820818448686821e-19L, +5.268740962820224108235e-21L, +3.414939674304748094484e-22L }; const long double midpoint = (1.75L + 1.0L)/2.0L; const long double halflen = (1.75L - 1.0L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 2.5L) { static const long double coeffs[deg + 1U] = { +2.843711194548592808550e-1L, -6.791774139166808940530e-2L, +6.955211059379384327814e-3L, -2.726366582146839486784e-4L, -6.516682485087925163874e-5L, +1.404387911504935155228e-5L, -1.103288540946056915318e-6L, -1.422154597293404846081e-8L, +1.102714664312839585330e-8L, -8.659211557383544255053e-10L, -8.048589443963965285748e-12L, +6.092061709996351761426e-12L, -3.580977611213519234324e-13L, -1.085173558590137965737e-14L, +2.411707924175380740802e-15L, -7.760751294610276598631e-17L, -6.701490147030045891595e-18L, +6.350145841254563572100e-19L, -2.034625734538917052251e-21L, -2.260543651146274653910e-21L, +9.782419961387425633151e-23L }; const long double midpoint = (2.5L + 1.75L)/2.0L; const long double halflen = (2.5L - 1.75L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 3.25L) { static const long double coeffs[deg + 1U] = { +1.901351274204578126827e-1L, -3.000575522193632460118e-2L, +2.672138524890489432579e-3L, -2.498237548675235150519e-4L, +2.013483163459701593271e-5L, -8.454663603108548182962e-7L, -8.036589636334016432368e-8L, +2.055498509671357933537e-8L, -2.052151324060186596995e-9L, +8.584315967075483822464e-11L, +5.062689357469596748991e-12L, -1.038671167196342609090e-12L, +6.367962851860231236238e-14L, +3.084688422647419767229e-16L, -3.417946142546575188490e-16L, +2.311567730100119302160e-17L, -6.170132546983726244716e-20L, -9.133176920944950460847e-20L, +5.712092431423316128728e-21L, +1.269641078369737220790e-23L, -2.072659711527711312699e-23L }; const long double midpoint = (3.25L + 2.5L)/2.0L; const long double halflen = (3.25L - 2.5L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 4.25L) { static const long double coeffs[deg + 1U] = { +1.402884974484995678749e-1L, -2.053975371995937033959e-2L, +1.595388628922920119352e-3L, -1.336894584910985998203e-4L, +1.224903774178156286300e-5L, -1.206856028658387948773e-6L, +1.187997233269528945503e-7L, -1.012936061496824448259e-8L, +5.244408240062370605664e-10L, +2.901444759022254846562e-11L, -1.168987502493903926906e-11L, +1.640096995420504465839e-12L, -1.339190668554209618318e-13L, +3.643815972666851044790e-15L, +6.922486581126169160232e-16L, -1.158761251467106749752e-16L, +8.164320395639210093180e-18L, -5.397918405779863087588e-20L, -5.052069908100339242896e-20L, +5.322512674746973445361e-21L, -1.869294542789169825747e-22L }; const long double midpoint = (4.25L + 3.25L)/2.0L; const long double halflen = (4.25L - 3.25L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 5.5L) { static const long double coeffs[deg + 1U] = { +1.058610209741581514157e-1L, -1.429297757627935191694e-2L, +9.911301703835545472874e-4L, -7.079903107876049846509e-5L, +5.229587914675267516134e-6L, -4.016071345964089296212e-7L, +3.231734714422926453741e-8L, -2.752870944370338482109e-9L, +2.503059741885009530630e-10L, -2.418699000594890423278e-11L, +2.410158905786160001792e-12L, -2.327254341132174000949e-13L, +1.958284411563056492727e-14L, -1.099893145048991004460e-15L, -2.959085292526991317697e-17L, +1.966366179276295203082e-17L, -3.314408783993662492621e-18L, +3.635520318133814622089e-19L, -2.550826919215104648800e-20L, +3.830090587178262542288e-22L, +1.836693763159216122739e-22L }; const long double midpoint = (5.5L + 4.25L)/2.0L; const long double halflen = (5.5L - 4.25L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 7.25L) { static const long double coeffs[deg + 1U] = { +8.024637207807814739314e-2L, -1.136614891549306029413e-2L, +8.164249750628661856014e-4L, -5.951964778701328943018e-5L, +4.407349502747483429390e-6L, -3.317746826184531133862e-7L, +2.541483569880571680365e-8L, -1.983391157250772649001e-9L, +1.579050614491277335581e-10L, -1.284592098551537518322e-11L, +1.070070857004674207604e-12L, -9.151832297362522251950e-14L, +8.065447314948125338081e-15L, -7.360105847607056315915e-16L, +6.995966000187407197283e-17L, -6.964349343411584120055e-18L, +7.268789359189778223225e-19L, -7.885125241947769024019e-20L, +8.689022564130615225208e-21L, -9.353211304381231554634e-22L +9.218280404899298404756e-23L }; const long double midpoint = (7.25L + 5.5L)/2.0L; const long double halflen = (7.25L - 5.5L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else { const int nterms = 50; long double term[nterms + 1]; const long double epsilon = LDBL_EPSILON/2.0L; const long double xsq = x*x; const long double twox = 2.0L*x; const long double twoxsq = 2.0L*xsq; long double xn = twoxsq; long double factorial = 1.0L; int n = 2; term[0] = 1.0L; term[1] = 1.0L/xn; for (; n <= nterms; ++n) { xn *= twoxsq; factorial *= (2*n - 1); term[n] = factorial/xn; if (term[n] < epsilon) break; } if (n > nterms) n = nterms; for (; n >= 0; --n) result += term[n]; result /= twox; } if (x_in < 0.0L) result *= -1.0L; return result; } long double inverseExpsqIntegral(const long double x_in) { const long double smallx = sqrtl(LDBL_EPSILON); const long double x = fabsl(x_in); long double result = 0.0L; // At the very beginning, the forward function behaves as x + x^3/3 if (x < smallx) result = x; // Are we below the integral upper limit of 1/2? else if (x < 0.5449871041836222236624201L) { static const long double coeffs[] = { 0.2579094020165510878448047L, 0.2509616007789004064162848L, -0.008037441097819308338583395L, -0.0009645198158737501684728918L, 0.0001298963476188684667851777L, 2.83303393956705648285913e-6L, -1.879031091504004970753414e-6L, 8.89033414496187588540693e-8L, 2.190598580039940838615595e-8L, -2.958785324082878640818902e-9L, -1.393055660359914508583854e-10L, 5.934164278634234871603216e-11L, -2.038318546964313720058418e-12L, -8.718325681035840557959458e-13L, 1.013108738230624177815487e-13L, 7.80188419892238096141015e-15L, -2.410959710954351931245264e-15L, 4.28768887009677018630627e-17L, 4.075320326600366577580808e-17L, -3.925657400094226595458262e-18L, -4.52032825617882087470387e-19L, 1.08172986581534678569357e-19L, 7.619258918621063777867205e-23L, -2.043372106982501109979208e-21L, 1.577535083736277024466253e-22L, 2.659861973971217145144534e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 0.2724935520918111118312101L; const long double halflen = midpoint; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } // Are we below the integral upper limit of 1? else if (x < 1.462651745907181608804049L) { static const long double coeffs[] = { 0.7736025700148903044157852L, 0.2483094370727645101600738L, -0.0236047661681907541554642L, 0.001718250905865626423590199L, -3.383723505479159579142638e-6L, -0.00002706856684823251983606182L, 5.564585822987557167324481e-6L, -6.297400054716145999859568e-7L, 1.793804726508794465230363e-8L, 9.974525492565834910508807e-9L, -2.630092849889704961303573e-9L, 3.584362156719671871930982e-10L, -1.848479157680947543606944e-11L, -4.505786185834443101923114e-12L, 1.497497786049082262992397e-12L, -2.345976552921671025035597e-13L, 1.675613965321671764579391e-14L, 2.080526674977816958132577e-15L, -9.222206647061674196609463e-16L, 1.635927767802667303380293e-16L, -1.458481834781130100146786e-17L, -8.477734754986502993701378e-19L, 5.883853143077372799032853e-19L, -1.178516243865358739330946e-19L, 1.244597734746846080964198e-20L, 1.834734630829857873717171e-22L, -3.79697806729613468087694e-22L, 8.640525246432963515541941e-23L, -1.072975125060893395507772e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 1.003819425045401916233234L; const long double halflen = 0.4588323208617796925708142L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } // Are we below the integral upper limit of 3/2? else if (x < 4.063114058624186262070998L) { static const long double coeffs[] = { 1.237474264255216149826922L, 0.2515376732927722040966062L, 0.01255809348578646846677752L, -0.001555984193633101089113648L, -0.00003266175488449684464535135L, 0.00001861132693718087179337154L, 3.038788024730725430182981e-7L, -3.055587383693358717339231e-7L, 2.195404333387651282823557e-10L, 5.223891523975639259623118e-9L, -8.74172529522376286448331e-11L, -9.289565568740836123320397e-11L, 3.037929091913979095346417e-12L, 1.697081673433425594636218e-12L, -8.376226831184981721677363e-14L, -3.144033465184008496066638e-14L, 2.107731553034119316083898e-15L, 5.865893869124205701237769e-16L, -5.061884931584969465772662e-17L, -1.096017981587409945712935e-17L, 1.182058217010476282761648e-18L, 2.040660311333339268713157e-19L, -2.708723954709022960299832e-20L, -3.767175977793052434129249e-21L, 6.118697880934383791686931e-22L, 6.976999259292224769625502e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 0.9003422806239746400989134L; const long double halflen = 0.2836972833794905219117775L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } // Are we below the integral upper limit of 2? else if (x < 16.4526277655072302247364L) { static const long double coeffs[] = { 1.750318675095801084701066L, 0.2503081949080952607456335L, -0.0003683988977021577213047602L, -0.0003053217430231520892093692L, 0.00004984978059867808464521542L, -2.910067173452788393366537e-6L, -1.234138038866370318487004e-7L, 3.698515524218906567339791e-8L, -2.598716127513287147666407e-9L, -8.028249755820916232073154e-11L, 3.38719730193993051252875e-11L, -2.805688393204576054583932e-12L, -4.636644559366711367183408e-14L, 3.430651608591744447051587e-14L, -3.234144260654980846019301e-15L, -1.479571272220063709849795e-17L, 3.663824718035758023592578e-17L, -3.866514289799821609911763e-18L, 2.280975526810673800594954e-20L, 4.039168081650367776987e-20L, -4.73516370177886271900388e-21L, 7.25430770939280219423293e-23L, 4.526173229162084391018431e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 1.428752297062238501405865L; const long double halflen = 0.2447127330587733393951737L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } // Are we below the integral upper limit of 3? else if (x < 1444.545122892714154713760L) { static const long double coeffs[] = { 2.50294099696623892772079L, 0.4994357077985400736034567L, -0.00291622179390211582491955L, 0.0005756988194441929049790467L, -0.00002810720142885454120222505L, -0.00001091749276776510487675747L, 3.29093808250468757948933e-6L, -4.888188874032568897843031e-7L, 4.172102909667302951211426e-8L, -4.599575541801073701575834e-10L, -6.053857271985015243693571e-10L, 1.509946372242848844836431e-10L, -2.455367299247951048634247e-11L, 2.672752200212154657233337e-12L, -9.103858678297231810743428e-14L, -3.718269817594061105825935e-14L, 1.069438018054794309965648e-14L, -1.744407969349348513403047e-15L, 1.883808777693390076608131e-16L, -7.095219918870126016632479e-18L, -2.561647683564104671111204e-18L, 7.953272251856515822155603e-19L, -1.363875890657281736925731e-19L, 1.519215969460881098898382e-20L, -6.014662801208487581697383e-22L, -2.017882982607174277058989e-22L, 6.264772567011851344575421e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 2.185393865913887927395372L; const long double halflen = 0.5119288357928760865943334L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } else { // At this point, sqrtl(logl(x)) looks more or less like // a straight line. We can do Newton-Raphson efficiently // in this transformed variable. const long double target = sqrtl(logl(x)); long double xtry = target, xold = 0.0L; while (fabsl((xtry - xold)/xtry) > 2.0*LDBL_EPSILON) { xold = xtry; const long double daws = dawsonIntegral(xtry); const long double sqrval = sqrtl(xtry*xtry + logl(daws)); const long double deriv = 0.5L/sqrval/daws; xtry += (target - sqrval)/deriv; } result = xtry; } if (x_in < 0.0L) result *= -1.0L; return result; } long double normalDensityDerivative(const unsigned order, const long double x) { static std::vector > coeffs; const unsigned ncoeffs = coeffs.size(); if (order >= ncoeffs) coeffs.resize(order + 1); unsigned nc = coeffs[order].size(); if (nc == 0U) { buildGaussDeriCoeffs(order, &coeffs[order]); nc = coeffs[order].size(); } assert(nc == order + 1U); return polySeriesSum(&coeffs[order][0], order, x)* expl(-x*x/2.0L)/SQRTTWOPIL; } double bivariateNormalIntegral(const double rho, const double x, const double y) { const double absrho = fabs(rho); if (absrho > 1.0) throw std::domain_error( "In npstat::bivariateNormalIntegral: " "correlation coefficient is outside of the [-1, 1] interval"); if (absrho < 0.925) return Phi2(-x, -y, rho); else return bvnd(x, y, rho); } double besselK(const double nu, const double x) { if (nu < 0.0 || x <= 0.0) throw std::domain_error( "In npstat::besselK: arguments outside of the allowed range"); const unsigned maxIter = 10000U; const double fpMin = DBL_MIN/DBL_EPSILON; const double xmin = 2.0; const int nl = std::floor(nu + 0.5); const double xmu = nu - nl; assert(std::abs(xmu) <= 0.5); const double xmu2 = xmu*xmu; const double xi = 1.0/x; const double xi2 = 2.0*xi; double h = nu*xi; if (h < fpMin) h = fpMin; double b = xi2*nu; double d = 0.0; double c = h; bool converged = false; for (unsigned i=0; i= 0; --l) { const double ritemp = fact*ril + ripl; fact -= xi; ripl = fact*ritemp + ril; ril = ritemp; } double rkmu, rk1; if (x < xmin) { const double xover2 = x/2.0; const long double ln2overx = -logl(xover2); const long double sigma = xmu*ln2overx; const double sqrEps = sqrt(DBL_EPSILON); const double pimu = M_PI*xmu; const double prefact = std::abs(pimu) < 0.1*sqrEps ? 1.0 : pimu/sin(pimu); const long double sqrlEps = sqrtl(LDBL_EPSILON); const long double fact2 = std::abs(sigma) < 0.1*sqrlEps ? 1.0L : sinhl(sigma)/sigma; const long double xover2gam = powl(xover2, xmu); const double gamplus = Gamma(1.0 + xmu); const double gamminus = Gamma(1.0 - xmu); const double gam2 = 0.5/gamplus + 0.5/gamminus; double gam1; if (std::abs(xmu) > 0.01) gam1 = (0.5/gamminus - 0.5/gamplus)/xmu; else { static const long double coeffs[] = { -0.577215664901532860606512L, 0.0420026350340952355290039L, 0.0421977345555443367482083L, 0.00721894324666309954239501L, 0.000215241674114950972815730L, 0.0000201348547807882386556894L }; gam1 = polySeriesSum(coeffs, 5U, xmu2); } long double p = 0.5/xover2gam*gamplus; long double q = 0.5*xover2gam*gamminus; long double ff = prefact*(gam1*coshl(sigma) + gam2*fact2*ln2overx); long double sum = ff; long double sum1 = p; long double ck = 1.0L; const double xover2sq = xover2*xover2; converged = false; for (unsigned i=1; i<=maxIter && !converged; ++i) { ff = (i*ff + p + q)/(i*i - xmu2); ck *= (xover2sq/i); p /= (i - xmu); q /= (i + xmu); const long double del = ck*ff; sum += del; const long double del1 = ck*(p - i*ff); sum1 += del1; if (std::abs(del) < std::abs(sum)*DBL_EPSILON) converged = true; } if (!converged) throw std::runtime_error( "In npstat::besselK: series failed to converge"); rkmu = sum; rk1 = sum1*xi2; } else { b = 2.0*(1.0 + x); d = 1.0/b; h = d; double delh = d; double q1 = 0.0; double q2 = 1.0; const double a1 = 0.25 - xmu2; double q = a1; c = a1; double a = -a1; double s = 1.0 + q*delh; converged = false; for (unsigned i=1; i<=maxIter && !converged; ++i) { a -= 2*i; c = -a*c/(i+1.0); const double qnew = (q1 - b*q2)/a; q1 = q2; q2 = qnew; q += c*qnew; b += 2.0; d = 1.0/(b+a*d); delh = (b*d-1.0)*delh; h += delh; const double dels = q*delh; s += dels; if (std::abs(dels/s) <= DBL_EPSILON) converged = true; } if (!converged) throw std::runtime_error( "In npstat::besselK: failed to converge in cf2"); h = a1*h; rkmu = std::sqrt(M_PI/(2.0*x))*std::exp(-x)/s; rk1 = rkmu*(xmu + x + 0.5 - h)*xi; } for (int i=1; i<=nl; ++i) { const double rktemp = (xmu + i)*xi2*rk1 + rkmu; rkmu = rk1; rk1 = rktemp; } return rkmu; } } Index: trunk/npstat/stat/Distributions1D.cc =================================================================== --- trunk/npstat/stat/Distributions1D.cc (revision 816) +++ trunk/npstat/stat/Distributions1D.cc (revision 817) @@ -1,2849 +1,2859 @@ #include #include #include #include #include "geners/binaryIO.hh" #include "npstat/nm/MathUtils.hh" #include "npstat/nm/SpecialFunctions.hh" #include "npstat/nm/interpolate.hh" #include "npstat/stat/Distributions1D.hh" #include "npstat/stat/StatUtils.hh" #include "npstat/stat/distributionReadError.hh" #define SQR2PI 2.5066282746310005 #define SQRT2 1.41421356237309505 #define SQRPI 1.77245385090551603 #define SQRT2L 1.414213562373095048801689L #define SQRPIL 1.77245385090551602729816748L #define TWOPIL 6.28318530717958647692528676656L static long double inverseErf(const long double fval) { long double x = npstat::inverseGaussCdf((fval + 1.0L)/2.0L)/SQRT2L; for (unsigned i=0; i<2; ++i) { const long double guessed = erfl(x); const long double deri = 2.0L/SQRPIL*expl(-x*x); x += (fval - guessed)/deri; } return x; } static unsigned improved_random(npstat::AbsRandomGenerator& g, long double* generatedRandom) { const long double extra = sqrt(DBL_EPSILON); long double u = 0.0L; unsigned callcount = 0; while (u <= 0.0L || u >= 1.0L) { u = g()*(1.0L + extra) - extra/2.0L; u += (g() - 0.5L)*extra; callcount += 2U; } *generatedRandom = u; return callcount; } static unsigned gauss_random(const double mean, const double sigma, npstat::AbsRandomGenerator& g, double* generatedRandom) { assert(generatedRandom); long double r1 = 0.0L, r2 = 0.0L; const unsigned calls = improved_random(g, &r1) + improved_random(g, &r2); *generatedRandom = mean + sigma*sqrtl(-2.0L*logl(r1))*sinl(TWOPIL*(r2-0.5L)); return calls; } // static unsigned gauss_random(const double mean, const double sigma, // npstat::AbsRandomGenerator& g, // double* generatedRandom) // { // assert(generatedRandom); // long double r1 = 0.0L; // const unsigned count = improved_random(g, &r1); // *generatedRandom = mean + sigma*SQRT2*inverseErf(2.0L*r1 - 1.0L); // return count; // } namespace npstat { bool SymmetricBeta1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, n_); return !os.fail(); } SymmetricBeta1D* SymmetricBeta1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { double n; gs::read_pod(in, &n); if (!in.fail()) return new SymmetricBeta1D(location, scale, n); } distributionReadError(in, classname()); return 0; } bool SymmetricBeta1D::isEqual(const AbsDistribution1D& otherBase) const { const SymmetricBeta1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && n_ == r.n_; } SymmetricBeta1D::SymmetricBeta1D(const double location, const double scale, const double power) : AbsScalableDistribution1D(location, scale), n_(power), intpow_(-1) { norm_ = calculateNorm(); } SymmetricBeta1D::SymmetricBeta1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), n_(params[0]), intpow_(-1) { norm_ = calculateNorm(); } double SymmetricBeta1D::calculateNorm() { static const double normcoeffs[11] = { 0.5, 0.75, 0.9375, 1.09375, 1.23046875, 1.353515625, 1.46630859375, 1.571044921875, 1.6692352294921875, 1.76197052001953125, 1.85006904602050781}; if (n_ <= -1.0) throw std::invalid_argument( "In npstat::SymmetricBeta1D::calculateNorm: " "invalid power parameter"); const int intpow = static_cast(floor(n_)); if (static_cast(intpow) == n_ && intpow >= 0 && intpow <= 10) { intpow_ = intpow; return normcoeffs[intpow]; } else return Gamma(1.5 + n_)/sqrt(M_PI)/Gamma(1.0 + n_); } double SymmetricBeta1D::unscaledDensity(const double x) const { const double oneminusrsq = 1.0 - x*x; if (oneminusrsq <= 0.0) return 0.0; else { double fcn; switch (intpow_) { case 0: fcn = 1.0; break; case 1: fcn = oneminusrsq; break; case 2: fcn = oneminusrsq*oneminusrsq; break; case 3: fcn = oneminusrsq*oneminusrsq*oneminusrsq; break; case 4: { const double tmp2 = oneminusrsq*oneminusrsq; fcn = tmp2*tmp2; } break; case 5: { const double tmp2 = oneminusrsq*oneminusrsq; fcn = tmp2*tmp2*oneminusrsq; } break; case 6: { const double tmp2 = oneminusrsq*oneminusrsq; fcn = tmp2*tmp2*tmp2; } break; default: fcn = pow(oneminusrsq, n_); break; } return norm_*fcn; } } double SymmetricBeta1D::unscaledCdf(const double x) const { if (x >= 1.0) return 1.0; else if (x <= -1.0) return 0.0; else { double cdf; if (n_ == 0.0) cdf = (x + 1.0)/2.0; else if (n_ == 1.0) cdf = (x*(3.0 - x*x) + 2.0)/4.0; else if (n_ == 2.0) { const double tmp1 = 1.0 + x; cdf = tmp1*tmp1*tmp1*(8.0 + 3.0*x*(x - 3.0))/16.0; } else if (n_ == 3.0) { const double xsq = x*x; cdf = (16. + x*(35. + xsq*(xsq*(21. - 5.*xsq) - 35.0)))/32.0; } else if (n_ == 4.0) { const double tmp1 = 1.0 + x; const double tmp1sq = tmp1*tmp1; const double tmp2 = 128. + x*(-325. + x*(345. + x*(35.*x - 175.))); cdf = tmp1sq*tmp1sq*tmp1*tmp2/256.0; } else cdf = incompleteBeta(n_+1.0, n_+1.0, (x + 1.0)/2.0); if (cdf < 0.0) cdf = 0.0; if (cdf > 1.0) cdf = 1.0; return cdf; } } double SymmetricBeta1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::SymmetricBeta1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return -1.0; else if (r1 == 1.0) return 1.0; else { double r; if (n_ == 0.0) r = r1*2.0 - 1.0; else r = 2.0*inverseIncompleteBeta(n_+1.0, n_+1.0, r1) - 1.0; if (r < -1.0) r = -1.0; else if (r > 1.0) r = 1.0; return r; } } bool Beta1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, alpha_); gs::write_pod(os, beta_); return !os.fail(); } Beta1D* Beta1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, a, b; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &a); gs::read_pod(in, &b); if (!in.fail()) return new Beta1D(location, scale, a, b); } distributionReadError(in, classname()); return 0; } bool Beta1D::isEqual(const AbsDistribution1D& otherBase) const { const Beta1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && alpha_ == r.alpha_ && beta_ == r.beta_; } Beta1D::Beta1D(const double location, const double scale, const double pa, const double pb) : AbsScalableDistribution1D(location, scale), alpha_(pa), beta_(pb) { norm_ = calculateNorm(); } Beta1D::Beta1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), alpha_(params[0]), beta_(params[1]) { norm_ = calculateNorm(); } double Beta1D::calculateNorm() const { if (!(alpha_ > 0.0 && beta_ > 0.0)) throw std::invalid_argument( "In npstat::Beta1D::calculateNorm: invalid power parameters"); return Gamma(alpha_ + beta_)/Gamma(alpha_)/Gamma(beta_); } double Beta1D::unscaledDensity(const double x) const { if (x <= 0.0 || x >= 1.0) return 0.0; else if (alpha_ == 1.0 && beta_ == 1.0) return 1.0; else return norm_*pow(x, alpha_-1.0)*pow(1.0-x, beta_-1.0); } double Beta1D::unscaledCdf(const double x) const { if (x >= 1.0) return 1.0; else if (x <= 0.0) return 0.0; else if (alpha_ == 1.0 && beta_ == 1.0) return x; else return incompleteBeta(alpha_, beta_, x); } double Beta1D::unscaledExceedance(const double x) const { return 1.0 - unscaledCdf(x); } double Beta1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Beta1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 0.0; else if (r1 == 1.0) return 1.0; else if (alpha_ == 1.0 && beta_ == 1.0) return r1; else return inverseIncompleteBeta(alpha_, beta_, r1); } Gamma1D::Gamma1D(const double location, const double scale, const double a) : AbsScalableDistribution1D(location, scale), alpha_(a) { initialize(); } Gamma1D::Gamma1D(double location, double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), alpha_(params[0]) { initialize(); } void Gamma1D::initialize() { if (!(alpha_ > 0.0)) throw std::invalid_argument( "In npstat::Gamma1D::initialize: invalid power parameter"); norm_ = 1.0/Gamma(alpha_); + uplim_ = -log(DBL_MIN); } bool Gamma1D::isEqual(const AbsDistribution1D& otherBase) const { const Gamma1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && alpha_ == r.alpha_; } double Gamma1D::unscaledDensity(const double x) const { - if (x > 0.0) + if (x > 0.0 && x <= uplim_) return norm_*pow(x, alpha_-1.0)*exp(-x); else return 0.0; } double Gamma1D::unscaledCdf(const double x) const { - if (x > 0.0) - return incompleteGamma(alpha_, x); - else + if (x <= 0.0) return 0.0; + else if (x >= uplim_) + return 1.0; + else + return incompleteGamma(alpha_, x); } double Gamma1D::unscaledExceedance(const double x) const { - if (x > 0.0) - return incompleteGammaC(alpha_, x); - else + if (x <= 0.0) return 1.0; + else if (x >= uplim_) + return 0.0; + else + return incompleteGammaC(alpha_, x); } double Gamma1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Gamma1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); - return inverseIncompleteGamma(alpha_, r1); + if (r1 == 0.0) + return 0.0; + else if (r1 == 1.0) + return uplim_; + else + return inverseIncompleteGamma(alpha_, r1); } bool Gamma1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, alpha_); return !os.fail(); } Gamma1D* Gamma1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, a; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &a); if (!in.fail()) return new Gamma1D(location, scale, a); } distributionReadError(in, classname()); return 0; } Gauss1D::Gauss1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale), xmin_(inverseGaussCdf(0.0)), xmax_(inverseGaussCdf(1.0)) { } Gauss1D::Gauss1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale), xmin_(inverseGaussCdf(0.0)), xmax_(inverseGaussCdf(1.0)) { } double Gauss1D::unscaledDensity(const double x) const { if (x < xmin_ || x > xmax_) return 0.0; else return exp(-x*x/2.0)/SQR2PI; } unsigned Gauss1D::random(AbsRandomGenerator& g, double* generatedRandom) const { return gauss_random(location(), scale(), g, generatedRandom); } bool Gauss1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Gauss1D* Gauss1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Gauss1D(location, scale); } double Gauss1D::unscaledCdf(const double x) const { if (x <= xmin_) return 0.0; if (x >= xmax_) return 1.0; if (x < 0.0) return erfc(-x/SQRT2)/2.0; else return (1.0 + erf(x/SQRT2))/2.0; } double Gauss1D::unscaledExceedance(const double x) const { if (x <= xmin_) return 1.0; if (x >= xmax_) return 0.0; if (x > 0.0) return erfc(x/SQRT2)/2.0; else return (1.0 - erf(x/SQRT2))/2.0; } double Gauss1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Gauss1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); return inverseGaussCdf(r1); } double Uniform1D::unscaledDensity(const double x) const { if (x >= 0.0 && x <= 1.0) return 1.0; else return 0.0; } bool Uniform1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Uniform1D* Uniform1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Uniform1D(location, scale); } double Uniform1D::unscaledCdf(const double x) const { if (x <= 0.0) return 0.0; else if (x >= 1.0) return 1.0; else return x; } double Uniform1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Uniform1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); return r1; } double IsoscelesTriangle1D::unscaledDensity(const double x) const { if (x > -1.0 && x < 1.0) return 1.0 - fabs(x); else return 0.0; } bool IsoscelesTriangle1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } IsoscelesTriangle1D* IsoscelesTriangle1D::read( const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new IsoscelesTriangle1D(location, scale); } double IsoscelesTriangle1D::unscaledCdf(const double x) const { if (x <= -1.0) return 0.0; else if (x >= 1.0) return 1.0; else if (x <= 0.0) { const double tmp = 1.0 + x; return 0.5*tmp*tmp; } else { const double tmp = 1.0 - x; return 1.0 - 0.5*tmp*tmp; } } double IsoscelesTriangle1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::IsoscelesTriangle1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return -1.0; else if (r1 == 1.0) return 1.0; else if (r1 <= 0.5) return sqrt(2.0*r1) - 1.0; else return 1.0 - sqrt((1.0 - r1)*2.0); } bool Exponential1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Exponential1D* Exponential1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Exponential1D(location, scale); } double Exponential1D::unscaledDensity(const double x) const { if (x < 0.0) return 0.0; const double eval = exp(-x); return eval < DBL_MIN ? 0.0 : eval; } double Exponential1D::unscaledCdf(const double x) const { return x > 0.0 ? 1.0 - exp(-x) : 0.0; } double Exponential1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Exponential1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 1.0) return -log(DBL_MIN); else return -log(1.0 - r1); } double Exponential1D::unscaledExceedance(const double x) const { if (x < 0.0) return 1.0; const double eval = exp(-x); return eval < DBL_MIN ? 0.0 : eval; } bool Logistic1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Logistic1D* Logistic1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Logistic1D(location, scale); } double Logistic1D::unscaledDensity(const double x) const { const double eval = exp(-x); if (eval < DBL_MIN) return 0.0; else { const double tmp = 1.0 + eval; return eval/tmp/tmp; } } double Logistic1D::unscaledCdf(const double x) const { const double lmax = -log(DBL_MIN); if (x <= -lmax) return 0.0; else if (x >= lmax) return 1.0; else return 1.0/(1.0 + exp(-x)); } double Logistic1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Logistic1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); const double lmax = -log(DBL_MIN); if (r1 == 0.0) return -lmax; else if (r1 == 1.0) return lmax; else return log(r1/(1.0 - r1)); } double Logistic1D::unscaledExceedance(const double x) const { const double lmax = -log(DBL_MIN); if (x >= lmax) return 0.0; else if (x <= -lmax) return 1.0; else { const double eval = exp(-x); return eval/(1.0 + eval); } } bool Quadratic1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, a_); gs::write_pod(os, b_); return !os.fail(); } Quadratic1D* Quadratic1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { double a, b; gs::read_pod(in, &a); gs::read_pod(in, &b); if (!in.fail()) return new Quadratic1D(location, scale, a, b); } distributionReadError(in, classname()); return 0; } bool Quadratic1D::isEqual(const AbsDistribution1D& otherBase) const { const Quadratic1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && a_ == r.a_ && b_ == r.b_; } Quadratic1D::Quadratic1D(const double location, const double scale, const double a, const double b) : AbsScalableDistribution1D(location, scale), a_(a), b_(b) { verifyNonNegative(); } Quadratic1D::Quadratic1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), a_(params[0]), b_(params[1]) { verifyNonNegative(); } void Quadratic1D::verifyNonNegative() { const double a = a_; const double b = b_; if (b == 0.0) { if (fabs(a) > 1.0) throw std::invalid_argument( "In npstat::Quadratic1D::verifyNonNegative:" " invalid distribution parameters"); } else { double x1 = 0.0, x2 = 0.0; const double sixb = 6*b; if (solveQuadratic((2*a-sixb)/sixb, (1-a+b)/sixb, &x1, &x2)) { if (!(fabs(x1 - 0.5) >= 0.5 && fabs(x2 - 0.5) >= 0.5)) throw std::invalid_argument( "In npstat::Quadratic1D::verifyNonNegative:" " invalid distribution parameters"); } } } double Quadratic1D::unscaledDensity(const double x) const { if (x < 0.0 || x > 1.0) return 0.0; else return 1.0 + b_ - a_ + x*(2.0*a_ + 6.0*b_*(x - 1.0)); } double Quadratic1D::unscaledCdf(const double x) const { if (x <= 0.0) return 0.0; else if (x >= 1.0) return 1.0; else return x*(1.0 + b_ - a_ + x*(a_ - 3.0*b_ + 2.0*b_*x)); } double Quadratic1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Quadratic1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 0.0; else if (r1 == 1.0) return 1.0; else { const double a = a_; const double b = b_; if (b == 0.0) { if (a == 0.0) return r1; else { double x0 = 0.0, x1 = 0.0; const unsigned n = solveQuadratic( (1.0 - a)/a, -r1/a, &x0, &x1); if (!n) throw std::runtime_error( "In npstat::Quadratic1D::unscaledQuantile: " "no solutions"); if (fabs(x0 - 0.5) < fabs(x1 - 0.5)) return x0; else return x1; } } else { const double twob = 2*b; double x[3] = {0.0}; const unsigned n = solveCubic( (a - 3*b)/twob, (1 - a + b)/twob, -r1/twob, x); if (n == 1U) return x[0]; else { unsigned ibest = 0; double dbest = fabs(x[0] - 0.5); for (unsigned i=1; i()); current.ensureSameId(id); double location, scale, a, b; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &a); gs::read_pod(in, &b); if (!in.fail()) return new LogQuadratic1D(location, scale, a, b); } distributionReadError(in, classname()); return 0; } bool LogQuadratic1D::isEqual(const AbsDistribution1D& otherBase) const { const LogQuadratic1D& r = static_cast( otherBase); return AbsScalableDistribution1D::isEqual(r) && a_ == r.a_ && b_ == r.b_; } LogQuadratic1D::LogQuadratic1D(const double location, const double scale, const double a, const double b) : AbsScalableDistribution1D(location, scale), a_(a), b_(b) { normalize(); } LogQuadratic1D::LogQuadratic1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), a_(params[0]), b_(params[1]) { normalize(); } inline long double LogQuadratic1D::quadInteg(const long double x) const { return dawsonIntegral(x)*expl(x*x); } void LogQuadratic1D::normalize() { ref_ = 0.0L; range_ = 1.0L; k_ = 0.0; s_ = 0.0; norm_ = 1.0; const double b = b_; const double a = a_; if (b > DBL_EPSILON) { k_ = sqrt(6.0*b); s_ = 0.5 - a/(6.0*b); ref_ = quadInteg(k_*s_); range_ = ref_ - quadInteg(k_*(s_ - 1.0)); norm_ = k_/range_; } else if (b < -DBL_EPSILON) { k_ = sqrt(-6.0*b); s_ = 0.5 - a/(6.0*b); ref_ = erfl(k_*s_); range_ = ref_ - erfl(k_*(s_ - 1.0)); norm_ = 2.0*k_/SQRPI/range_; } else if (fabs(a) > DBL_EPSILON) { range_ = expl(2.0L*a) - 1.0L; if (fabs(a) > 1.e-10) norm_ = a/sinh(a); } } double LogQuadratic1D::unscaledDensity(const double x) const { if (x < 0.0 || x > 1.0) return 0.0; const double b = b_; if (fabs(b) > DBL_EPSILON) { const double delta = x - s_; return norm_*exp(6.0*b*delta*delta); } const double a = a_; if (fabs(a) > DBL_EPSILON) return norm_*exp((2.0*x - 1.0)*a); else return 1.0; } double LogQuadratic1D::unscaledCdf(const double x) const { if (x <= 0.0) return 0.0; else if (x >= 1.0) return 1.0; else { const double b = b_; const double a = a_; if (b > DBL_EPSILON) return (ref_ - quadInteg(k_*(s_ - x)))/range_; else if (b < -DBL_EPSILON) return (ref_ - erfl(k_*(s_ - x)))/range_; else if (fabs(a) > DBL_EPSILON) return (expl(2.0L*a*x) - 1.0L)/range_; else return x; } } double LogQuadratic1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::LogQuadratic1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 0.0; else if (r1 == 1.0) return 1.0; else { const double b = b_; const double a = a_; double q = 0.0; if (b > DBL_EPSILON) q = s_ - inverseExpsqIntegral(ref_ - r1*range_)/k_; else if (b < -DBL_EPSILON) q = s_ - inverseErf(ref_ - r1*range_)/k_; else if (fabs(a) > DBL_EPSILON) q = logl(r1*range_ + 1.0L)/2.0/a; else q = r1; if (q < 0.0) q = 0.0; else if (q > 1.0) q = 1.0; return q; } } bool TruncatedGauss1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, nsigma_); return !os.fail(); } TruncatedGauss1D* TruncatedGauss1D::read( const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, nsig; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &nsig); if (!in.fail()) return new TruncatedGauss1D(location, scale, nsig); } distributionReadError(in, classname()); return 0; } bool TruncatedGauss1D::isEqual(const AbsDistribution1D& otherBase) const { const TruncatedGauss1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && nsigma_ == r.nsigma_; } void TruncatedGauss1D::initialize() { if (nsigma_ <= 0.0) throw std::invalid_argument( "In npstat::TruncatedGauss1D::initialize: " "invalid truncation parameter"); const double maxSig = inverseGaussCdf(1.0); if (nsigma_ >= maxSig) { nsigma_ = maxSig; cdf0_ = 0.0; norm_ = 1.0; } else { cdf0_ = erfc(nsigma_/SQRT2)/2.0; const double u = (1.0 + erf(nsigma_/SQRT2))/2.0; norm_ = 1.0/(u - cdf0_); } } TruncatedGauss1D::TruncatedGauss1D(const double location, const double scale, const double i_nsigma) : AbsScalableDistribution1D(location, scale), nsigma_(i_nsigma) { initialize(); } TruncatedGauss1D::TruncatedGauss1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), nsigma_(params[0]) { initialize(); } double TruncatedGauss1D::unscaledDensity(const double x) const { if (fabs(x) > nsigma_) return 0.0; else return norm_*exp(-x*x/2.0)/SQR2PI; } unsigned TruncatedGauss1D::random(AbsRandomGenerator& g, double* generatedRandom) const { const double m = location(); const double s = scale(); unsigned cnt = gauss_random(m, s, g, generatedRandom); while (fabs(*generatedRandom - m) > nsigma_*s) cnt += gauss_random(m, s, g, generatedRandom); return cnt; } double TruncatedGauss1D::unscaledCdf(const double x) const { if (x <= -nsigma_) return 0.0; else if (x >= nsigma_) return 1.0; else if (x < 0.0) return (erfc(-x/SQRT2)/2.0 - cdf0_)*norm_; else return ((1.0 + erf(x/SQRT2))/2.0 - cdf0_)*norm_; } double TruncatedGauss1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::TruncatedGauss1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return -nsigma_; else if (r1 == 1.0) return nsigma_; else return inverseGaussCdf(r1/norm_ + cdf0_); } bool MirroredGauss1D::isEqual(const AbsDistribution1D& otherBase) const { const MirroredGauss1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && mu0_ == r.mu0_ && sigma0_ == r.sigma0_; } MirroredGauss1D::MirroredGauss1D(const double location, const double scale, const double mean, double const sigma) : AbsScalableDistribution1D(location, scale), mu0_(mean), sigma0_(sigma) { validate(); } MirroredGauss1D::MirroredGauss1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), mu0_(params[0]), sigma0_(params[1]) { validate(); } bool MirroredGauss1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, mu0_); gs::write_pod(os, sigma0_); return !os.fail(); } MirroredGauss1D* MirroredGauss1D::read( const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, mu, sig; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &mu); gs::read_pod(in, &sig); if (!in.fail()) return new MirroredGauss1D(location, scale, mu, sig); } distributionReadError(in, classname()); return 0; } double MirroredGauss1D::unscaledDensity(const double x) const { if (x < 0.0 || x > 1.0) return 0.0; Gauss1D g(x, sigma0_); long double acc = g.density(mu0_)*1.0L + g.density(-mu0_); for (unsigned k=1; ; ++k) { const long double old = acc; acc += g.density(2.0*k + mu0_); acc += g.density(2.0*k - mu0_); acc += g.density(-2.0*k + mu0_); acc += g.density(-2.0*k - mu0_); if (old == acc) break; } return acc; } long double MirroredGauss1D::ldCdf(const double x) const { if (x <= 0.0) return 0.0L; else if (x >= 1.0) return 1.0L; else { long double acc = 0.0L; { Gauss1D g(mu0_, sigma0_); acc += (g.cdf(x) - g.cdf(0.0)); } { Gauss1D g(-mu0_, sigma0_); acc += (g.cdf(x) - g.cdf(0.0)); } for (unsigned k=1; ; ++k) { const long double old = acc; { Gauss1D g(2.0*k + mu0_, sigma0_); acc += (g.cdf(x) - g.cdf(0.0)); } { Gauss1D g(2.0*k - mu0_, sigma0_); acc += (g.cdf(x) - g.cdf(0.0)); } { Gauss1D g(-2.0*k + mu0_, sigma0_); acc -= (g.exceedance(x) - g.exceedance(0.0)); } { Gauss1D g(-2.0*k - mu0_, sigma0_); acc -= (g.exceedance(x) - g.exceedance(0.0)); } if (old == acc) break; } return acc; } } double MirroredGauss1D::unscaledCdf(const double x) const { return ldCdf(x); } double MirroredGauss1D::unscaledExceedance(const double x) const { return 1.0L - ldCdf(x); } double MirroredGauss1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::MirroredGauss1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 0.0; else if (r1 == 1.0) return 1.0; else { const long double ldr1 = r1; double xmin = 0.0; double xmax = 1.0; for (unsigned i=0; i<1000; ++i) { if ((xmax - xmin)/xmax <= 2.0*DBL_EPSILON) break; const double xtry = (xmin + xmax)/2.0; const long double ld = ldCdf(xtry); if (ld == ldr1) return xtry; else if (ld > ldr1) xmax = xtry; else xmin = xtry; } return (xmin + xmax)/2.0; } } void MirroredGauss1D::validate() { if (mu0_ < 0.0 || mu0_ > 1.0) throw std::invalid_argument( "In MirroredGauss1D::validate: interval mean must be within [0, 1]"); if (sigma0_ <= 0.0) throw std::invalid_argument( "In MirroredGauss1D::validate: interval sigma must be positive"); } BifurcatedGauss1D::BifurcatedGauss1D( const double location, const double scale, const double i_leftSigmaFraction, const double i_nSigmasLeft, const double i_nSigmasRight) : AbsScalableDistribution1D(location, scale), leftSigma_(i_leftSigmaFraction*2.0), rightSigma_(2.0 - leftSigma_), nSigmasLeft_(i_nSigmasLeft), nSigmasRight_(i_nSigmasRight) { initialize(); } BifurcatedGauss1D::BifurcatedGauss1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), leftSigma_(params[0]*2.0), rightSigma_(2.0 - leftSigma_), nSigmasLeft_(params[1]), nSigmasRight_(params[2]) { initialize(); } void BifurcatedGauss1D::initialize() { if (leftSigma_ < 0.0 || rightSigma_ < 0.0) throw std::invalid_argument( "In npstat::BifurcatedGauss1D::initialize: " "invalid left sigma fraction"); if (nSigmasLeft_ < 0.0) throw std::invalid_argument( "In npstat::BifurcatedGauss1D::initialize: " "invalid left truncation parameter"); if (nSigmasRight_ < 0.0) throw std::invalid_argument( "In npstat::BifurcatedGauss1D::initialize: " "invalid right truncation parameter"); if (nSigmasLeft_ + nSigmasRight_ == 0.0) throw std::invalid_argument( "In npstat::BifurcatedGauss1D::initialize: " "both truncation parameters can not be 0"); const double maxNSigma = inverseGaussCdf(1.0); if (nSigmasRight_ > maxNSigma) nSigmasRight_ = maxNSigma; if (nSigmasLeft_ > maxNSigma) nSigmasLeft_ = maxNSigma; cdf0Left_ = erfc(nSigmasLeft_/SQRT2)/2.0; cdf0Right_ = (1.0 + erf(nSigmasRight_/SQRT2))/2.0; assert(cdf0Right_ > cdf0Left_); const double leftArea = (0.5 - cdf0Left_)*leftSigma_; const double rightArea = (cdf0Right_ - 0.5)*rightSigma_; norm_ = 1.0/(leftArea + rightArea); leftCdfFrac_ = leftArea/(leftArea + rightArea); } double BifurcatedGauss1D::unscaledDensity(const double x) const { if (x == 0.0) return norm_/SQR2PI; else if (x > 0.0) { if (x > rightSigma_*nSigmasRight_) return 0.0; else { const double dx = x/rightSigma_; return norm_*exp(-dx*dx/2.0)/SQR2PI; } } else { if (x < -leftSigma_*nSigmasLeft_) return 0.0; else { const double dx = x/leftSigma_; return norm_*exp(-dx*dx/2.0)/SQR2PI; } } } double BifurcatedGauss1D::unscaledCdf(const double x) const { if (x == 0.0) return leftCdfFrac_; else if (x > 0.0) { if (x > rightSigma_*nSigmasRight_) return 1.0; else { const double dx = x/rightSigma_; const double cdfDelta = (1.0 + erf(dx/SQRT2))/2.0 - cdf0Right_; return 1.0 + cdfDelta*(1.0 - leftCdfFrac_)/(cdf0Right_ - 0.5); } } else { if (x < -leftSigma_*nSigmasLeft_) return 0.0; else { const double dx = x/leftSigma_; const double cdfDelta = erfc(-dx/SQRT2)/2.0 - cdf0Left_; return cdfDelta*leftCdfFrac_/(0.5 - cdf0Left_); } } } double BifurcatedGauss1D::unscaledExceedance(const double x) const { return 1.0 - unscaledCdf(x); } double BifurcatedGauss1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::BifurcatedGauss1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return -nSigmasLeft_*leftSigma_; else if (r1 == 1.0) return nSigmasRight_*rightSigma_; else { if (r1 == leftCdfFrac_) return 0.0; else if (r1 < leftCdfFrac_) { // Map 0 into cdf0Left_ and leftCdfFrac_ into 0.5 const double arg = r1/leftCdfFrac_*(0.5-cdf0Left_) + cdf0Left_; return leftSigma_*inverseGaussCdf(arg); } else { // Map leftCdfFrac_ into 0.5 and 1.0 into cdf0Right_ const double d = (r1 - leftCdfFrac_)/(1.0 - leftCdfFrac_); const double arg = 0.5 + d*(cdf0Right_ - 0.5); return rightSigma_*inverseGaussCdf(arg); } } } bool BifurcatedGauss1D::isEqual(const AbsDistribution1D& otherBase) const { const BifurcatedGauss1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && leftSigma_ == r.leftSigma_ && rightSigma_ == r.rightSigma_ && nSigmasLeft_ == r.nSigmasLeft_ && nSigmasRight_ == r.nSigmasRight_ && norm_ == r.norm_ && leftCdfFrac_ == r.leftCdfFrac_ && cdf0Left_ == r.cdf0Left_ && cdf0Right_ == r.cdf0Right_; } bool BifurcatedGauss1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, leftSigma_); gs::write_pod(os, rightSigma_); gs::write_pod(os, nSigmasLeft_); gs::write_pod(os, nSigmasRight_); gs::write_pod(os, norm_); gs::write_pod(os, leftCdfFrac_); gs::write_pod(os, cdf0Left_); gs::write_pod(os, cdf0Right_); return !os.fail(); } BifurcatedGauss1D::BifurcatedGauss1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale) { } BifurcatedGauss1D* BifurcatedGauss1D::read( const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { BifurcatedGauss1D* ptr = new BifurcatedGauss1D(location, scale); gs::read_pod(in, &ptr->leftSigma_); gs::read_pod(in, &ptr->rightSigma_); gs::read_pod(in, &ptr->nSigmasLeft_); gs::read_pod(in, &ptr->nSigmasRight_); gs::read_pod(in, &ptr->norm_); gs::read_pod(in, &ptr->leftCdfFrac_); gs::read_pod(in, &ptr->cdf0Left_); gs::read_pod(in, &ptr->cdf0Right_); if (!in.fail() && ptr->leftSigma_ >= 0.0 && ptr->rightSigma_ >= 0.0 && ptr->leftSigma_ + ptr->rightSigma_ > 0.0 && ptr->nSigmasLeft_ >= 0.0 && ptr->nSigmasRight_ >= 0.0 && ptr->nSigmasLeft_ + ptr->nSigmasRight_ > 0.0 && ptr->norm_ > 0.0 && ptr->leftCdfFrac_ >= 0.0 && ptr->cdf0Left_ >= 0.0 && ptr->cdf0Right_ > ptr->cdf0Left_) return ptr; else delete ptr; } distributionReadError(in, classname()); return 0; } Cauchy1D::Cauchy1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale), support_(sqrt(DBL_MAX/M_PI)) { } Cauchy1D::Cauchy1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale), support_(sqrt(DBL_MAX/M_PI)) { } bool Cauchy1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Cauchy1D* Cauchy1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Cauchy1D(location, scale); } double Cauchy1D::unscaledDensity(const double x) const { if (fabs(x) < support_) return 1.0/M_PI/(1.0 + x*x); else return 0.0; } double Cauchy1D::unscaledCdf(const double x) const { if (x < -support_) return 0.0; else if (x > support_) return 1.0; else return atan(x)/M_PI + 0.5; } double Cauchy1D::unscaledQuantile(const double x) const { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::Cauchy1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (x == 0.0) return -support_; else if (x == 1.0) return support_; else return tan(M_PI*(x - 0.5)); } bool LogNormal::isEqual(const AbsDistribution1D& otherBase) const { const LogNormal& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && skew_ == r.skew_; } bool LogNormal::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, skew_); return !os.fail(); } LogNormal* LogNormal::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, skew; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &skew); if (!in.fail()) return new LogNormal(location, scale, skew); } distributionReadError(in, classname()); return 0; } void LogNormal::initialize() { logw_ = 0.0; s_ = 0.0; xi_ = 0.0; emgamovd_ = 0.0; if (skew_) { const double b1 = skew_*skew_; const double tmp = pow((2.0+b1+sqrt(b1*(4.0+b1)))/2.0, 1.0/3.0); const double w = tmp + 1.0/tmp - 1.0; logw_ = log(w); if (logw_ > 0.0) { s_ = sqrt(logw_); emgamovd_ = 1.0/sqrt(w*(w-1.0)); xi_ = -emgamovd_*sqrt(w); } else { // This is not different from a Gaussian within // the numerical precision of our calculations logw_ = 0.0; skew_ = 0.0; } } } LogNormal::LogNormal(const double mean, const double stdev, const double skewness) : AbsScalableDistribution1D(mean, stdev), skew_(skewness) { initialize(); } LogNormal::LogNormal(const double mean, const double stdev, const std::vector& params) : AbsScalableDistribution1D(mean, stdev), skew_(params[0]) { initialize(); } double LogNormal::unscaledDensity(const double x) const { if (skew_) { const double diff = skew_ > 0.0 ? x - xi_ : -x - xi_; if (diff <= 0.0) return 0.0; else { const double lg = log(diff/emgamovd_); return exp(-lg*lg/2.0/logw_)/s_/SQR2PI/diff; } } else { // This is a Gaussian return exp(-x*x/2.0)/SQR2PI; } } double LogNormal::unscaledCdf(const double x) const { if (skew_) { const double diff = skew_ > 0.0 ? x - xi_ : -x - xi_; double posCdf = 0.0; if (diff > 0.0) posCdf = (1.0 + erf(log(diff/emgamovd_)/s_/SQRT2))/2.0; return skew_ > 0.0 ? posCdf : 1.0 - posCdf; } else return (1.0 + erf(x/SQRT2))/2.0; } double LogNormal::unscaledExceedance(const double x) const { return 1.0 - unscaledCdf(x); } double LogNormal::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::LogNormal::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); const double g = inverseGaussCdf(skew_ >= 0.0 ? r1 : 1.0 - r1); if (skew_) { const double v = emgamovd_*exp(s_*g) + xi_; return skew_ > 0.0 ? v : -v; } else return g; } Moyal1D::Moyal1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale), xmax_(-2.0*log(DBL_MIN*SQR2PI)), xmin_(-log(xmax_)) { } Moyal1D::Moyal1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale), xmax_(-2.0*log(DBL_MIN*SQR2PI)), xmin_(-log(xmax_)) { } double Moyal1D::unscaledDensity(const double x) const { if (x <= xmin_ || x >= xmax_) return 0.0; else return exp(-0.5*(x + exp(-x)))/SQR2PI; } double Moyal1D::unscaledCdf(const double x) const { if (x <= xmin_) return 0.0; else if (x >= xmax_) return 1.0; else return incompleteGammaC(0.5, 0.5*exp(-x)); } double Moyal1D::unscaledExceedance(const double x) const { if (x <= xmin_) return 1.0; else if (x >= xmax_) return 0.0; else return incompleteGamma(0.5, 0.5*exp(-x)); } double Moyal1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Moyal1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return xmin_; else if (r1 == 1.0) return xmax_; else { const double d = inverseIncompleteGammaC(0.5, r1); return -log(2.0*d); } } bool Moyal1D::write(std::ostream& os) const { return AbsScalableDistribution1D::write(os); } Moyal1D* Moyal1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (!AbsScalableDistribution1D::read(in, &location, &scale)) { distributionReadError(in, classname()); return 0; } return new Moyal1D(location, scale); } bool Pareto1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, c_); return !os.fail(); } Pareto1D* Pareto1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, c; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &c); if (!in.fail()) return new Pareto1D(location, scale, c); } distributionReadError(in, classname()); return 0; } bool Pareto1D::isEqual(const AbsDistribution1D& otherBase) const { const Pareto1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && c_ == r.c_; } void Pareto1D::initialize() { if (c_ <= 0.0) throw std::invalid_argument( "In npstat::Pareto1D::initialize: power parameter must be positive"); if (c_ > 1.0) support_ = pow(1.0/DBL_MIN, 1.0/c_); else support_ = 1.0/DBL_MIN; } Pareto1D::Pareto1D(const double location, const double scale, const double powerParam) : AbsScalableDistribution1D(location, scale), c_(powerParam) { initialize(); } Pareto1D::Pareto1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), c_(params[0]) { initialize(); } double Pareto1D::unscaledDensity(const double x) const { if (x < 1.0 || x > support_) return 0.0; else return c_/pow(x, c_ + 1.0); } double Pareto1D::unscaledCdf(const double x) const { if (x <= 1.0) return 0.0; else if (x >= support_) return 1.0; else return 1.0 - 1.0/pow(x, c_); } double Pareto1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Pareto1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 1.0; else if (r1 == 1.0) return support_; else return pow(1.0 - r1, -1.0/c_); } double Pareto1D::unscaledExceedance(const double x) const { if (x <= 1.0) return 1.0; else if (x >= support_) return 0.0; else return 1.0/pow(x, c_); } bool UniPareto1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, c_); return !os.fail(); } UniPareto1D* UniPareto1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, c; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &c); if (!in.fail()) return new UniPareto1D(location, scale, c); } distributionReadError(in, classname()); return 0; } bool UniPareto1D::isEqual(const AbsDistribution1D& otherBase) const { const UniPareto1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && c_ == r.c_; } void UniPareto1D::initialize() { if (c_ <= 0.0) throw std::invalid_argument( "In npstat::UniPareto1D::initialize: power parameter must be positive"); if (c_ > 1.0) support_ = pow(1.0/DBL_MIN, 1.0/c_); else support_ = 1.0/DBL_MIN; amplitude_ = c_/(c_ + 1.0); } UniPareto1D::UniPareto1D(const double location, const double scale, const double powerParam) : AbsScalableDistribution1D(location, scale), c_(powerParam) { initialize(); } UniPareto1D::UniPareto1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), c_(params[0]) { initialize(); } double UniPareto1D::unscaledDensity(const double x) const { if (x < 0.0 || x > support_) return 0.0; else if (x <= 1.0) return amplitude_; else return amplitude_/pow(x, c_ + 1.0); } double UniPareto1D::unscaledCdf(const double x) const { if (x <= 0.0) return 0.0; else if (x >= support_) return 1.0; else if (x <= 1.0) return x*amplitude_; else return amplitude_ + (1.0 - amplitude_)*(1.0 - 1.0/pow(x, c_)); } double UniPareto1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::UniPareto1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 <= amplitude_) return r1/amplitude_; else if (r1 == 1.0) return support_; else return pow(1.0 - (r1 - amplitude_)/(1.0 - amplitude_), -1.0/c_); } double UniPareto1D::unscaledExceedance(const double x) const { if (x > 1.0) return (1.0 - amplitude_)/pow(x, c_); else return 1.0 - unscaledCdf(x); } bool Huber1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, tailWeight_); return !os.fail(); } Huber1D* Huber1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale, tailWeight; if (AbsScalableDistribution1D::read(in, &location, &scale)) { gs::read_pod(in, &tailWeight); if (!in.fail()) return new Huber1D(location, scale, tailWeight); } distributionReadError(in, classname()); return 0; } bool Huber1D::isEqual(const AbsDistribution1D& otherBase) const { const Huber1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && tailWeight_ == r.tailWeight_; } void Huber1D::initialize() { if (!(tailWeight_ >= 0.0 && tailWeight_ < 1.0)) throw std::invalid_argument( "In npstat::Huber1D::initialize: " "tail weight not inside [0, 1) interval"); if (tailWeight_ == 0.0) { // Pure Gaussian a_ = DBL_MAX; normfactor_ = 1.0/sqrt(2.0*M_PI); support_ = inverseGaussCdf(1.0); cdf0_ = 0.0; } else { // Solve the equation for "a" by bisection const double eps = 2.0*DBL_EPSILON; double c = -2.0*log(tailWeight_); assert(c > 0.0); assert(weight(c) <= tailWeight_); double b = 0.0; while ((c - b)/(c + b) > eps) { const double half = (c + b)/2.0; if (weight(half) >= tailWeight_) b = half; else c = half; } a_ = (c + b)/2.0; normfactor_ = 0.5/(exp(-a_*a_/2.0)/a_ + sqrt(M_PI/2.0)*erf(a_/SQRT2)); support_ = a_/2.0 - log(DBL_MIN)/a_; cdf0_ = (1.0 + erf(-a_/SQRT2))/2.0; } } Huber1D::Huber1D(const double location, const double scale, const double tailWeight) : AbsScalableDistribution1D(location, scale), tailWeight_(tailWeight) { initialize(); } Huber1D::Huber1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), tailWeight_(params[0]) { initialize(); } double Huber1D::unscaledDensity(const double x) const { const double absx = fabs(x); if (absx <= a_) return normfactor_*exp(-x*x/2.0); else return normfactor_*exp(a_*(a_/2.0 - absx)); } double Huber1D::unscaledCdf(const double x) const { if (tailWeight_ == 0.0) return (1.0 + erf(x/SQRT2))/2.0; if (x < -a_) return normfactor_*exp((a_*(a_ + 2*x))/2.0)/a_; else if (x <= a_) { static const double sq1 = sqrt(M_PI/2.0); static const double sq2 = sqrt(2.0); return normfactor_*sq1*(erf(a_/sq2) + erf(x/sq2)) + tailWeight_/2; } else return 1.0 - normfactor_*exp((a_*(a_ - 2*x))/2.0)/a_; } double Huber1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Huber1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (tailWeight_ == 0.0) return inverseGaussCdf(r1); if (r1 == 0.0) return -support_; else if (r1 == 1.0) return support_; else if (r1 <= tailWeight_/2.0) return log(r1*a_/normfactor_)/a_ - 0.5*a_; else if (r1 < 1.0 - tailWeight_/2.0) { const double t = (r1 - tailWeight_/2.0)/normfactor_/SQR2PI + cdf0_; return inverseGaussCdf(t); } else return 0.5*a_ - log((1.0 - r1)*a_/normfactor_)/a_; } double Huber1D::weight(const double a) const { static const double sq1 = sqrt(M_PI/2.0); return 1.0/(1.0 + a*exp(a*a/2.0)*sq1*erf(a/SQRT2)); } bool Tabulated1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, deg_); gs::write_pod_vector(os, table_); return !os.fail(); } Tabulated1D* Tabulated1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { unsigned deg = 0; gs::read_pod(in, °); std::vector table; gs::read_pod_vector(in, &table); if (!in.fail() && table.size()) return new Tabulated1D(location, scale, &table[0], table.size(), deg); } distributionReadError(in, classname()); return 0; } bool Tabulated1D::isEqual(const AbsDistribution1D& otherBase) const { const Tabulated1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && table_ == r.table_ && deg_ == r.deg_; } Tabulated1D::Tabulated1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale) { const unsigned npara = params.size(); if (npara) initialize(¶ms[0], npara, std::min(3U, npara-1U)); else initialize(static_cast(0), 0U, 0U); } void Tabulated1D::normalize() { cdf_.clear(); cdf_.reserve(len_); cdf_.push_back(0.0); long double sum = 0.0L; for (unsigned i=0; i(sum)); } double norm = cdf_[len_ - 1]; assert(norm > 0.0); for (unsigned i=0; i0; --i) { sum += intervalInteg(i-1); exceed_[i-1] = static_cast(sum); } norm = exceed_[0]; for (unsigned i=0; i 1.0) return 0.0; if (x == 0.0) return table_[0]; if (x == 1.0) return table_[len_ - 1]; unsigned idx = static_cast(x/step_); if (idx >= len_ - 1) idx = len_ - 2; const double dx = x/step_ - idx; switch (deg_) { case 0: if (dx < 0.5) return table_[idx]; else return table_[idx + 1]; case 1: return interpolate_linear(dx, table_[idx], table_[idx + 1]); case 2: if (idx == 0) return interpolate_quadratic(dx, table_[idx], table_[idx + 1], table_[idx + 2]); else if (idx == len_ - 2) return interpolate_quadratic(dx+1.0, table_[idx - 1], table_[idx], table_[idx + 1]); else { const double v0 = interpolate_quadratic( dx, table_[idx], table_[idx + 1], table_[idx + 2]); const double v1 = interpolate_quadratic( dx+1.0, table_[idx - 1], table_[idx], table_[idx + 1]); return (v0 + v1)/2.0; } case 3: if (idx == 0) return interpolate_cubic(dx, table_[idx], table_[idx + 1], table_[idx + 2], table_[idx + 3]); else if (idx == len_ - 2) return interpolate_cubic(dx+2.0, table_[idx-2], table_[idx-1], table_[idx], table_[idx + 1]); else return interpolate_cubic(dx+1.0, table_[idx - 1], table_[idx], table_[idx + 1], table_[idx + 2]); default: assert(0); return 0.0; } } double Tabulated1D::unscaledDensity(const double x) const { const double v = this->interpolate(x); if (v >= 0.0) return v; else return 0.0; } double Tabulated1D::unscaledCdf(double x) const { if (x <= 0.0) return 0.0; if (x >= 1.0) return 1.0; unsigned idx = static_cast(x/step_); if (idx >= len_ - 1) idx = len_ - 2; double v; switch (deg_) { case 0: { const double dx = x/step_ - idx; if (dx < 0.5) v = table_[idx]*dx*step_; else v = (table_[idx]*0.5 + table_[idx + 1]*(dx - 0.5))*step_; } break; default: v = interpIntegral(step_*idx, x); } return cdf_[idx] + v; } double Tabulated1D::unscaledExceedance(double x) const { if (x <= 0.0) return 1.0; if (x >= 1.0) return 0.0; unsigned idx = static_cast(x/step_); if (idx >= len_ - 1) idx = len_ - 2; double v; switch (deg_) { case 0: { const double dx = x/step_ - idx; if (dx < 0.5) v = table_[idx]*dx*step_; else v = (table_[idx]*0.5 + table_[idx + 1]*(dx - 0.5))*step_; } break; default: v = interpIntegral(step_*idx, x); } return exceed_[idx] - v; } double Tabulated1D::interpIntegral(const double a, const double b) const { static const double legendreRootOver2 = 0.5/sqrt(3.0); const double x0 = (b + a)/2.0; const double step = b - a; const double v0 = unscaledDensity(x0 - legendreRootOver2*step); const double v1 = unscaledDensity(x0 + legendreRootOver2*step); return step*(v0 + v1)/2.0; } double Tabulated1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::Tabulated1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0) return 0.0; if (r1 == 1.0) return 1.0; unsigned idx = std::lower_bound(cdf_.begin(), cdf_.end(), r1) - cdf_.begin() - 1U; double xlo = step_*idx; const double dcdf = r1 - cdf_[idx]; assert(dcdf > 0.0); switch (deg_) { case 0: { const double c1 = table_[idx]*0.5*step_; if (dcdf <= c1) { assert(table_[idx] > 0.0); return xlo + dcdf/table_[idx]; } else { assert(table_[idx+1] > 0.0); return xlo + 0.5*step_ + (dcdf - c1)/table_[idx+1]; } } case 1: { const double a = (table_[idx+1] - table_[idx])/step_/2.0; if (a == 0.0) { assert(table_[idx] > 0.0); return xlo + dcdf/table_[idx]; } else { double x1, x2; const unsigned nroots = solveQuadratic( table_[idx]/a, -dcdf/a, &x1, &x2); if (nroots != 2U) throw std::runtime_error( "In npstat::Tabulated1D::unscaledQuantile: " "unexpected number of solutions"); if (fabs(x1 - 0.5*step_) < fabs(x2 - 0.5*step_)) return xlo + x1; else return xlo + x2; } } default: { double xhi = xlo + step_; const double eps = 2.0*DBL_EPSILON; while ((xhi - xlo)/(xhi + xlo) > eps) { const double med = (xhi + xlo)/2.0; if (unscaledCdf(med) >= r1) xhi = med; else xlo = med; } return (xhi + xlo)/2.0; } } } bool BinnedDensity1D::isEqual(const AbsDistribution1D& otherBase) const { const double eps = 1.0e-12; const BinnedDensity1D& r = static_cast(otherBase); if (!AbsScalableDistribution1D::isEqual(r)) return false; if (!(deg_ == r.deg_)) return false; const unsigned long n = table_.size(); if (!(n == r.table_.size())) return false; for (unsigned long i=0; i eps) return false; return true; } BinnedDensity1D::BinnedDensity1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale) { const unsigned npara = params.size(); if (npara) initialize(¶ms[0], npara, std::min(1U, npara-1U)); else initialize(static_cast(0), 0U, 0U); } void BinnedDensity1D::normalize() { cdf_.clear(); cdf_.reserve(len_); long double sum = 0.0L; switch (deg_) { case 0U: for (unsigned i=0; i(sum)); } break; case 1U: { double oldval = 0.0; double* data = &table_[0]; for (unsigned i=0; i(sum)); } sum += oldval*0.5; } break; default: assert(0); } const double norm = static_cast(sum); assert(norm > 0.0); const double integ = norm*step_; for (unsigned i=0; i 1.0) return 0.0; switch (deg_) { case 0: { unsigned idx = static_cast(x/step_); if (idx > len_ - 1) idx = len_ - 1; return table_[idx]; } case 1: { const double xs = x - step_/2.0; if (xs <= 0.0) return table_[0]; const unsigned idx = static_cast(xs/step_); if (idx > len_ - 2) return table_[len_ - 1]; const double dx = xs/step_ - idx; return interpolate_linear(dx, table_[idx], table_[idx + 1]); } default: assert(0); return 0.0; } } double BinnedDensity1D::unscaledCdf(double x) const { if (x <= 0.0) return 0.0; if (x >= 1.0) return 1.0; double v = 0.0; switch (deg_) { case 0: { unsigned idx = static_cast(x/step_); if (idx > len_ - 1) idx = len_ - 1; v = (idx ? cdf_[idx - 1] : 0.0) + table_[idx]*(x - idx*step_); } break; case 1: { const double xs = x - step_/2.0; if (xs <= 0.0) v = table_[0]*x; else { const unsigned idx = static_cast(xs/step_); if (idx > len_ - 2) v = 1.0 - table_[len_ - 1]*(1.0 - x); else { const double dx = xs - idx*step_; const double slope = (table_[idx+1] - table_[idx])/step_; v = cdf_[idx] + table_[idx]*dx + slope*dx*dx/2.0; } } } break; default: assert(0); break; } if (v < 0.0) v = 0.0; else if (v > 1.0) v = 1.0; return v; } double BinnedDensity1D::unscaledQuantile(const double r1) const { if (!(r1 >= 0.0 && r1 <= 1.0)) throw std::domain_error( "In npstat::BinnedDensity1D::unscaledQuantile: " "cdf argument outside of [0, 1] interval"); if (r1 == 0.0 && deg_) return 0.0; if (r1 == 1.0 && deg_) return 1.0; double v = 0.0; switch (deg_) { case 0: { if (r1 == 0) v = firstNonZeroBin_; else if (r1 == 1.0) v = lastNonZeroBin_ + 1U; else if (r1 <= cdf_[0]) v = r1/cdf_[0]; else { double rem; const unsigned bin = quantileBinFromCdf(&cdf_[0],len_,r1,&rem) + 1U; assert(bin < len_); v = bin + rem; } } break; case 1: { if (r1 <= cdf_[0]) v = 0.5*r1/cdf_[0]; else if (r1 >= cdf_[len_ - 1]) v = (len_-0.5+0.5*(r1-cdf_[len_-1])/(1.0-cdf_[len_-1])); else { const unsigned idx = std::lower_bound(cdf_.begin(),cdf_.end(),r1) - cdf_.begin() - 1U; assert(idx < len_ - 1); const double k = (table_[idx+1] - table_[idx])/step_; const double y = r1 - cdf_[idx]; double x; if (fabs(k) < 1.e-10*table_[idx]) x = y/table_[idx]; else { const double b = 2.0*table_[idx]/k; const double c = -2.0*y/k; double x1, x2; if (solveQuadratic(b, c, &x1, &x2)) { if (fabs(x1 - step_*0.5) < fabs(x2 - step_*0.5)) x = x1; else x = x2; } else { // This can happen due to various round-off problems. // Assume that the quadratic equation determinant // should have been 0 instead of negative. x = -b/2.0; } } if (x < 0.0) x = 0.0; else if (x > step_) x = step_; v = x/step_ + idx + 0.5; } } break; default: assert(0); return 0.0; } v *= step_; if (v < 0.0) v = 0.0; else if (v > 1.0) v = 1.0; return v; } bool BinnedDensity1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, deg_); gs::write_pod_vector(os, table_); return !os.fail(); } BinnedDensity1D* BinnedDensity1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current( gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { unsigned deg = 0; gs::read_pod(in, °); std::vector table; gs::read_pod_vector(in, &table); if (!in.fail() && table.size()) return new BinnedDensity1D(location, scale, &table[0], table.size(), deg); } distributionReadError(in, classname()); return 0; } bool StudentsT1D::write(std::ostream& os) const { AbsScalableDistribution1D::write(os); gs::write_pod(os, nDoF_); return !os.fail(); } StudentsT1D* StudentsT1D::read(const gs::ClassId& id, std::istream& in) { static const gs::ClassId current(gs::ClassId::makeId()); current.ensureSameId(id); double location, scale; if (AbsScalableDistribution1D::read(in, &location, &scale)) { double nDoF; gs::read_pod(in, &nDoF); if (!in.fail()) return new StudentsT1D(location, scale, nDoF); } distributionReadError(in, classname()); return 0; } bool StudentsT1D::isEqual(const AbsDistribution1D& otherBase) const { const StudentsT1D& r = static_cast(otherBase); return AbsScalableDistribution1D::isEqual(r) && nDoF_ == r.nDoF_; } StudentsT1D::StudentsT1D(const double location, const double scale, const double degreesOfFreedom) : AbsScalableDistribution1D(location, scale), nDoF_(degreesOfFreedom) { initialize(); } StudentsT1D::StudentsT1D(const double location, const double scale, const std::vector& params) : AbsScalableDistribution1D(location, scale), nDoF_(params[0]) { initialize(); } void StudentsT1D::initialize() { if (nDoF_ <= 0.0) throw std::invalid_argument( "In npstat::StudentsT1D::initialize: invalid number " "of degrees of freedom"); power_ = (nDoF_ + 1.0)/2.0; bignum_ = 0.0; normfactor_ = Gamma(power_)/Gamma(nDoF_/2.0)/sqrt(nDoF_)/SQRPI; } double StudentsT1D::unscaledDensity(const double x) const { return normfactor_*pow(1.0 + x*x/nDoF_, -power_); } double StudentsT1D::unscaledCdf(const double t) const { const double s = sqrt(t*t + nDoF_); return incompleteBeta(nDoF_/2.0, nDoF_/2.0, (t + s)/(2.0*s)); } double StudentsT1D::unscaledExceedance(const double x) const { return 1.0 - unscaledCdf(x); } double StudentsT1D::unscaledQuantile(const double r1) const { if (r1 == 0.5) return 0.0; const double c = inverseIncompleteBeta(nDoF_/2.0, nDoF_/2.0, r1); const double tmp = 2.0*c - 1.0; const double a = tmp*tmp; const double denom = 1.0 - a; if (denom > 0.0) { const double sqroot = sqrt(nDoF_*a/denom); return r1 > 0.5 ? sqroot : -sqroot; } else { if (bignum_ == 0.0) (const_cast(this))->bignum_ = effectiveSupport(); return r1 > 0.5 ? bignum_ : -bignum_; } } double StudentsT1D::effectiveSupport() const { // Figure out at which (positive) values of the argument // the density becomes effectively indistinguishable from 0 const double biglog = (log(normfactor_) - log(DBL_MIN) + power_*log(nDoF_))/(2.0*power_); if (biglog >= log(DBL_MAX)) return DBL_MAX; else return exp(biglog); } } Index: trunk/npstat/stat/EllipticalDistribution.cc =================================================================== --- trunk/npstat/stat/EllipticalDistribution.cc (revision 816) +++ trunk/npstat/stat/EllipticalDistribution.cc (revision 817) @@ -1,238 +1,238 @@ #include #include #include #include "geners/GenericIO.hh" #include "geners/vectorIO.hh" #include "npstat/nm/SpecialFunctions.hh" #include "npstat/rng/convertToSphericalRandom.hh" #include "npstat/stat/EllipticalDistribution.hh" #include "npstat/stat/Distributions1D.hh" #include "npstat/stat/distributionReadError.hh" static double calculateGNorm(const npstat::AbsDistribution1D& gDensity, const npstat::AbsDistribution1D& hDensity, const unsigned dim) { const double gmax = gDensity.quantile(1.0); const double hmax = hDensity.quantile(1.0); assert(gmax > 0.0 && hmax > 0.0); double rsq = 1.0; - while (gmax >= rsq || hmax >= rsq) + while (rsq >= gmax || rsq >= hmax) rsq /= 2.0; const double gd = gDensity.density(rsq); assert(gd > 0.0); const double hd = hDensity.density(rsq); assert(hd > 0.0); - const double dover2 = dim*0.5; + const double dover2 = dim/2.0; const double piu = M_PI*rsq; return hd*npstat::Gamma(dover2)*rsq/gd/pow(piu,dover2); } namespace npstat { EllipticalDistribution::EllipticalDistribution( const double* location, const unsigned i_dim, const Matrix& transformationMatrix, const AbsDistribution1D& gDensity, const AbsDistribution1D& hDensity) : AbsDistributionND(i_dim), mu_(location, location+i_dim), A_(transformationMatrix), g_(0), h_(0), det_(0.0), gNorm_(0.0), buf_(i_dim) { assert(location); assert(i_dim); if (!A_.isSquare()) throw std::invalid_argument( "In npstat::EllipticalDistribution constructor: " "transformation matrix is not square"); if (A_.nRows() != i_dim) throw std::invalid_argument( "In npstat::EllipticalDistribution constructor: " "incompatible dimensions of the shift and " "transformation matrix"); if (gDensity.quantile(0.0) != 0.0) throw std::invalid_argument( "In npstat::EllipticalDistribution constructor: " "inappropriate g distribution"); if (hDensity.quantile(0.0) != 0.0) throw std::invalid_argument( "In npstat::EllipticalDistribution constructor: " "inappropriate h distribution"); det_ = std::abs(A_.det()); assert(det_); const Matrix& covmat = A_.timesT(); InvCovmat_ = covmat.symPDInv(); g_ = gDensity.clone(); h_ = hDensity.clone(); gNorm_ = calculateGNorm(gDensity, hDensity, i_dim); } EllipticalDistribution::EllipticalDistribution( const EllipticalDistribution& r) : AbsDistributionND(r.dim_), mu_(r.mu_), A_(r.A_), InvCovmat_(r.InvCovmat_), g_(r.g_->clone()), h_(r.h_->clone()), det_(r.det_), gNorm_(r.gNorm_), buf_(r.dim_) { } EllipticalDistribution& EllipticalDistribution::operator=( const EllipticalDistribution& r) { if (&r != this) { AbsDistributionND::operator=(r); cleanup(); mu_ = r.mu_; A_ = r.A_; InvCovmat_ = r.InvCovmat_; g_ = r.g_->clone(); h_ = r.h_->clone(); det_ = r.det_; gNorm_ = r.gNorm_; buf_.resize(r.dim_); } return *this; } bool EllipticalDistribution::isEqual(const AbsDistributionND& other) const { const EllipticalDistribution& r = static_cast(other); return mu_ == r.mu_ && A_ == r.A_ && *g_ == *r.g_ && *h_ == *r.h_; } void EllipticalDistribution::cleanup() { A_.uninitialize(); InvCovmat_.uninitialize(); delete g_; g_ = 0; delete h_; h_ = 0; } double EllipticalDistribution::density( const double* x, const unsigned i_dim) const { assert(x); assert(i_dim == dim_); double* delta = &buf_[0]; const double* mu = &mu_[0]; for (unsigned i=0; idensity(chisq)*gNorm_/det_; } void EllipticalDistribution::unitMap( const double* rnd, const unsigned bufLen, double* x) const { if (bufLen) { if (bufLen < dim_) throw std::invalid_argument( "In npstat::EllipticalDistribution::unitMap: insufficient " "length of the random buffer"); if (bufLen > dim_) throw std::out_of_range( "In npstat::EllipticalDistribution::unitMap: " "input point dimensionality is out of range"); assert(rnd); assert(x); double *direction = &buf_[0]; const double rn = convertToSphericalRandom(rnd, dim_, direction); const double u = h_->quantile(rn); const double r = sqrt(u); for (unsigned i=0; i= dim_); Gauss1D normal(0.0, 1.0); unsigned callCount = 0U; double *direction = &buf_[0]; long double sumsql = 0.0L; for (unsigned i=0; iquantile(g()); ++callCount; const double r = sqrt(u); for (unsigned i=0; i()); current.ensureSameId(id); std::vector mu; gs::read_pod_vector(is, &mu); Matrix transform; gs::restore_item(is, &transform); CPP11_auto_ptr g = gs::read_obj(is); CPP11_auto_ptr h = gs::read_obj(is); if (is.fail()) { distributionReadError(is, classname()); return 0; } else return new EllipticalDistribution(&mu[0], mu.size(), transform, *g, *h); } } Index: trunk/npstat/stat/Distributions1D.hh =================================================================== --- trunk/npstat/stat/Distributions1D.hh (revision 816) +++ trunk/npstat/stat/Distributions1D.hh (revision 817) @@ -1,1067 +1,1068 @@ #ifndef NPSTAT_DISTRIBUTIONS1D_HH_ #define NPSTAT_DISTRIBUTIONS1D_HH_ /*! // \file Distributions1D.hh // // \brief A number of useful 1-d continuous statistical distributions // // Author: I. Volobouev // // November 2009 */ #include "npstat/stat/Distribution1DFactory.hh" namespace npstat { /** // The uniform distribution is defined here by a constant density // equal to 1 between 0 and 1 and equal to 0 everywhere else */ class Uniform1D : public AbsScalableDistribution1D { public: inline Uniform1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale) {} inline virtual Uniform1D* clone() const {return new Uniform1D(*this);} inline virtual ~Uniform1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Uniform1D";} static inline unsigned version() {return 1;} static Uniform1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; inline Uniform1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale) {} inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return 1.0 - unscaledCdf(x);} }; /** Isosceles triangle distribution: 1 - |x| supported on [-1, 1] */ class IsoscelesTriangle1D : public AbsScalableDistribution1D { public: inline IsoscelesTriangle1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale) {} inline virtual IsoscelesTriangle1D* clone() const {return new IsoscelesTriangle1D(*this);} inline virtual ~IsoscelesTriangle1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::IsoscelesTriangle1D";} static inline unsigned version() {return 1;} static IsoscelesTriangle1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; inline IsoscelesTriangle1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale) {} inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return 1.0 - unscaledCdf(x);} }; /** Exponential distribution. "scale" is the decay time. */ class Exponential1D : public AbsScalableDistribution1D { public: inline Exponential1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale) {} inline virtual Exponential1D* clone() const {return new Exponential1D(*this);} inline virtual ~Exponential1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Exponential1D";} static inline unsigned version() {return 1;} static Exponential1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; inline Exponential1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale) {} inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; }; /** Logistic distribution */ class Logistic1D : public AbsScalableDistribution1D { public: inline Logistic1D(const double location, const double scale) : AbsScalableDistribution1D(location, scale) {} inline virtual Logistic1D* clone() const {return new Logistic1D(*this);} inline virtual ~Logistic1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Logistic1D";} static inline unsigned version() {return 1;} static Logistic1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; inline Logistic1D(const double location, const double scale, const std::vector& /* params */) : AbsScalableDistribution1D(location, scale) {} inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; }; /** // A distribution whose density has a simple quadratic shape. // The support is from 0 to 1, and the coefficients "a" and "b" // are the coefficients for the Legendre polynomials of 1st // and 2nd degree translated to the support region. Note that // only those values of "a" and "b" that guarantee non-negativity // of the density are allowed, otherwise the code will generate // a run-time error. */ class Quadratic1D : public AbsScalableDistribution1D { public: Quadratic1D(double location, double scale, double a, double b); inline virtual Quadratic1D* clone() const {return new Quadratic1D(*this);} inline virtual ~Quadratic1D() {} inline double a() const {return a_;} inline double b() const {return b_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Quadratic1D";} static inline unsigned version() {return 3;} static Quadratic1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; Quadratic1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} void verifyNonNegative(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return 1.0 - unscaledCdf(x);} double a_; double b_; }; /** // A distribution whose density logarithm has a simple quadratic // shape. The support is from 0 to 1, and the coefficients "a" and "b" // are the coefficients for the Legendre polynomials of 1st and 2nd // degree translated to the support region. */ class LogQuadratic1D : public AbsScalableDistribution1D { public: LogQuadratic1D(double location, double scale, double a, double b); inline virtual LogQuadratic1D* clone() const {return new LogQuadratic1D(*this);} inline virtual ~LogQuadratic1D() {} inline double a() const {return a_;} inline double b() const {return b_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::LogQuadratic1D";} static inline unsigned version() {return 3;} static LogQuadratic1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; LogQuadratic1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} void normalize(); long double quadInteg(long double x) const; double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return 1.0 - unscaledCdf(x);} long double ref_; long double range_; double a_; double b_; double k_; double s_; double norm_; }; /** The Gaussian (or Normal) distribution */ class Gauss1D : public AbsScalableDistribution1D { public: Gauss1D(double location, double scale); inline virtual Gauss1D* clone() const {return new Gauss1D(*this);} inline virtual ~Gauss1D() {} // Higher quality generator than the one provided by // the quantile function virtual unsigned random(AbsRandomGenerator& g, double* generatedRandom) const; // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Gauss1D";} static inline unsigned version() {return 1;} static Gauss1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; Gauss1D(const double location, const double scale, const std::vector& params); inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; double xmin_; double xmax_; }; /** Gaussian distribution truncated at some number of sigmas */ class TruncatedGauss1D : public AbsScalableDistribution1D { public: TruncatedGauss1D(double location, double scale, double nsigma); inline virtual TruncatedGauss1D* clone() const {return new TruncatedGauss1D(*this);} inline virtual ~TruncatedGauss1D() {} inline double nsigma() const {return nsigma_;} // Higher quality generator than the one provided by // the quantile function virtual unsigned random(AbsRandomGenerator& g, double* generatedRandom) const; // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::TruncatedGauss1D";} static inline unsigned version() {return 1;} static TruncatedGauss1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; TruncatedGauss1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return unscaledCdf(-x);} double nsigma_; double norm_; double cdf0_; }; /** // Gaussian distribution on the [0, 1] interval, mirrored at the boundaries. // This is the Green's function of the diffusion equation on [0, 1]. The // interval can be shifted and scaled as for the uniform distribution. */ class MirroredGauss1D : public AbsScalableDistribution1D { public: MirroredGauss1D(double location, double scale, double meanOn0_1, double sigmaOn0_1); inline virtual MirroredGauss1D* clone() const {return new MirroredGauss1D(*this);} inline virtual ~MirroredGauss1D() {} inline double meanOn0_1() const {return mu0_;} inline double sigmaOn0_1() const {return sigma0_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::MirroredGauss1D";} static inline unsigned version() {return 1;} static MirroredGauss1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; MirroredGauss1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} void validate(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; long double ldCdf(double x) const; double mu0_; double sigma0_; }; /** // Bifurcated Gaussian distribution. Different sigmas // can be used on the left and on the right, with constructor // parameter "leftSigmaFraction" specifying the ratio of // the left sigma to the sum of sigmas (this ratio must be // between 0 and 1). Different truncations in terms of the // number of sigmas can be used as well. */ class BifurcatedGauss1D : public AbsScalableDistribution1D { public: BifurcatedGauss1D(double location, double scale, double leftSigmaFraction, double nSigmasLeft, double nSigmasRight); inline virtual BifurcatedGauss1D* clone() const {return new BifurcatedGauss1D(*this);} inline virtual ~BifurcatedGauss1D() {} inline double leftSigmaFraction() const {return leftSigma_/(leftSigma_ + rightSigma_);} inline double nSigmasLeft() const {return nSigmasLeft_;} inline double nSigmasRight() const {return nSigmasRight_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::BifurcatedGauss1D";} static inline unsigned version() {return 1;} static BifurcatedGauss1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: BifurcatedGauss1D(double location, double scale); friend class ScalableDistribution1DFactory; BifurcatedGauss1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 3;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; double leftSigma_; double rightSigma_; double nSigmasLeft_; double nSigmasRight_; double norm_; double leftCdfFrac_; double cdf0Left_; double cdf0Right_; }; /** Symmetric beta distribution */ class SymmetricBeta1D : public AbsScalableDistribution1D { public: SymmetricBeta1D(double location, double scale, double power); inline virtual SymmetricBeta1D* clone() const {return new SymmetricBeta1D(*this);} inline virtual ~SymmetricBeta1D() {} inline double power() const {return n_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::SymmetricBeta1D";} static inline unsigned version() {return 1;} static SymmetricBeta1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; SymmetricBeta1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return unscaledCdf(-x);} double calculateNorm(); double n_; double norm_; int intpow_; }; /** Beta1D density is proportional to x^(apha-1) * (1-x)^(beta-1) */ class Beta1D : public AbsScalableDistribution1D { public: Beta1D(double location, double scale, double alpha, double beta); inline virtual Beta1D* clone() const {return new Beta1D(*this);} inline virtual ~Beta1D() {} inline double alpha() const {return alpha_;} inline double beta() const {return beta_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Beta1D";} static inline unsigned version() {return 1;} static Beta1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; Beta1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 2;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; double calculateNorm() const; double alpha_; double beta_; double norm_; }; /** Shiftable gamma distribution */ class Gamma1D : public AbsScalableDistribution1D { public: Gamma1D(double location, double scale, double alpha); inline virtual Gamma1D* clone() const {return new Gamma1D(*this);} inline virtual ~Gamma1D() {} inline double alpha() const {return alpha_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Gamma1D";} static inline unsigned version() {return 1;} static Gamma1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; Gamma1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; void initialize(); double alpha_; double norm_; + double uplim_; }; /** // Pareto distribution. Location parameter is location of 0, scale // parameter is the distance between 0 and the start of the density // (like the normal Pareto distribution location parameter). */ class Pareto1D : public AbsScalableDistribution1D { public: Pareto1D(double location, double scale, double powerParameter); inline virtual Pareto1D* clone() const {return new Pareto1D(*this);} inline virtual ~Pareto1D() {} inline double powerParameter() const {return c_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Pareto1D";} static inline unsigned version() {return 1;} static Pareto1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; Pareto1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; double c_; double support_; }; /** // Uniform distribution with Pareto tail attached to the right, where // the support of the uniform would normally end. Location parameter // is location of 0, scale parameter is the width of the uniform part // (like the normal Pareto distribution location parameter). */ class UniPareto1D : public AbsScalableDistribution1D { public: UniPareto1D(double location, double scale, double powerParameter); inline virtual UniPareto1D* clone() const {return new UniPareto1D(*this);} inline virtual ~UniPareto1D() {} inline double powerParameter() const {return c_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::UniPareto1D";} static inline unsigned version() {return 1;} static UniPareto1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; UniPareto1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; double c_; double support_; double amplitude_; }; /** "Huber" distribution */ class Huber1D : public AbsScalableDistribution1D { public: Huber1D(double location, double scale, double tailWeight); inline virtual Huber1D* clone() const {return new Huber1D(*this);} inline virtual ~Huber1D() {} inline double tailWeight() const {return tailWeight_;} inline double tailStart() const {return a_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Huber1D";} static inline unsigned version() {return 1;} static Huber1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; Huber1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return unscaledCdf(-x);} double weight(double a) const; double tailWeight_; double a_; double normfactor_; double support_; double cdf0_; }; /** Cauchy (or Breit-Wigner) distribution */ class Cauchy1D : public AbsScalableDistribution1D { public: Cauchy1D(double location, double scale); inline virtual Cauchy1D* clone() const {return new Cauchy1D(*this);} inline virtual ~Cauchy1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Cauchy1D";} static inline unsigned version() {return 1;} static Cauchy1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; Cauchy1D(const double location, const double scale, const std::vector& params); inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return unscaledCdf(-x);} double support_; }; /** // Log-normal distribution represented by its mean, standard // deviation, and skewness. This representation is more useful // than other representations encountered in statistical literature. */ class LogNormal : public AbsScalableDistribution1D { public: LogNormal(double mean, double stdev, double skewness); inline virtual LogNormal* clone() const {return new LogNormal(*this);} inline virtual ~LogNormal() {} inline double skewness() const {return skew_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::LogNormal";} static inline unsigned version() {return 1;} static LogNormal* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; LogNormal(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; double skew_; double logw_; double s_; double xi_; double emgamovd_; }; /** // Moyal Distribution (originally derived by Moyal as an approximation // to the Landau distribution) */ class Moyal1D : public AbsScalableDistribution1D { public: Moyal1D(double location, double scale); inline virtual Moyal1D* clone() const {return new Moyal1D(*this);} inline virtual ~Moyal1D() {} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Moyal1D";} static inline unsigned version() {return 1;} static Moyal1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D& r) const {return AbsScalableDistribution1D::isEqual(r);} private: friend class ScalableDistribution1DFactory; Moyal1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 0;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; double xmax_; double xmin_; }; /** Student's t-distribution */ class StudentsT1D : public AbsScalableDistribution1D { public: StudentsT1D(double location, double scale, double nDegreesOfFreedom); inline virtual StudentsT1D* clone() const {return new StudentsT1D(*this);} inline virtual ~StudentsT1D() {} inline double nDegreesOfFreedom() const {return nDoF_;} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::StudentsT1D";} static inline unsigned version() {return 1;} static StudentsT1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; StudentsT1D(double location, double scale, const std::vector& params); inline static int nParameters() {return 1;} void initialize(); double effectiveSupport() const; double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledExceedance(double x) const; double unscaledQuantile(double x) const; double nDoF_; double normfactor_; double power_; double bignum_; }; /** Distribution defined by an interpolation table */ class Tabulated1D : public AbsScalableDistribution1D { public: // The "data" array gives (unnormalized) density values at // equidistant intervals. data[0] is density at 0.0, and // data[dataLen-1] is density at 1.0. If "dataLen" is less // than 2, uniform distribution will be created. Internally, // the data is kept in double precision. // // "interpolationDegree" must be less than 4 and less than "dataLen". // template Tabulated1D(double location, double scale, const Real* data, unsigned dataLen, unsigned interpolationDegree); inline Tabulated1D(const double location, const double scale, const std::vector& table, const unsigned interpolationDegree) : AbsScalableDistribution1D(location, scale) { const unsigned long sz = table.size(); initialize(sz ? &table[0] : (double*)0, sz, interpolationDegree); } inline virtual Tabulated1D* clone() const {return new Tabulated1D(*this);} inline virtual ~Tabulated1D() {} inline unsigned interpolationDegree() const {return deg_;} inline unsigned tableLength() const {return len_;} inline const double* tableData() const {return &table_[0];} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::Tabulated1D";} static inline unsigned version() {return 1;} static Tabulated1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; // The following constructor creates interpolator // of maximum degree possible Tabulated1D(double location, double scale, const std::vector& params); inline static int nParameters() {return -1;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; double unscaledExceedance(double x) const; template void initialize( const Real* data, unsigned dataLen, unsigned interpolationDegree); void normalize(); double interpolate(double x) const; double intervalInteg(unsigned intervalNumber) const; double interpIntegral(double x0, double x1) const; std::vector table_; std::vector cdf_; std::vector exceed_; double step_; unsigned len_; unsigned deg_; }; /** // Another interpolated distribution. For this one, we will assume // that the coordinates correspond to 1-d histogram bin centers. */ class BinnedDensity1D : public AbsScalableDistribution1D { public: // The "data" array gives density values at equidistant intervals. // data[0] is density at 0.5/dataLen, and data[dataLen-1] is density // at 1.0 - 0.5/dataLen. // // "interpolationDegree" must be less than 2. // template BinnedDensity1D(double location, double scale, const Real* data, unsigned dataLen, unsigned interpolationDegree); inline BinnedDensity1D(const double location, const double scale, const std::vector& table, const unsigned interpolationDegree) : AbsScalableDistribution1D(location, scale) { const unsigned long sz = table.size(); initialize(sz ? &table[0] : (double*)0, sz, interpolationDegree); } inline virtual BinnedDensity1D* clone() const {return new BinnedDensity1D(*this);} inline virtual ~BinnedDensity1D() {} inline unsigned interpolationDegree() const {return deg_;} inline unsigned tableLength() const {return len_;} inline const double* tableData() const {return &table_[0];} // Methods needed for I/O virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::BinnedDensity1D";} static inline unsigned version() {return 1;} static BinnedDensity1D* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistribution1D&) const; private: friend class ScalableDistribution1DFactory; // The following constructor creates interpolator // of maximum degree possible BinnedDensity1D(double location, double scale, const std::vector& params); inline static int nParameters() {return -1;} double unscaledDensity(double x) const; double unscaledCdf(double x) const; double unscaledQuantile(double x) const; inline double unscaledExceedance(const double x) const {return 1.0 - unscaledCdf(x);} template void initialize( const Real* data, unsigned dataLen, unsigned interpolationDegree); void normalize(); double interpolate(double x) const; std::vector table_; std::vector cdf_; double step_; unsigned len_; unsigned deg_; unsigned firstNonZeroBin_; unsigned lastNonZeroBin_; }; } #include "npstat/stat/Distributions1D.icc" #endif // NPSTAT_DISTRIBUTIONS1D_HH_ Index: trunk/npstat/stat/EllipticalDistribution.hh =================================================================== --- trunk/npstat/stat/EllipticalDistribution.hh (revision 816) +++ trunk/npstat/stat/EllipticalDistribution.hh (revision 817) @@ -1,96 +1,96 @@ #ifndef NPSTAT_ELLIPTICALDISTRIBUTION_HH_ #define NPSTAT_ELLIPTICALDISTRIBUTION_HH_ /*! // \file EllipticalDistribution.hh // // \brief Multivariate elliptical ditributions // // Author: I. Volobouev // // September 2022 */ +#include + #include "npstat/nm/Matrix.hh" #include "npstat/stat/AbsDistribution1D.hh" #include "npstat/stat/AbsDistributionND.hh" namespace npstat { class EllipticalDistribution : public AbsDistributionND { public: /** // The constructor arguments are as follows: // // location, dim The shift and the dimensionality of // the distribution. The array "location" // must have at least "dim" elements. // // transformationMatrix The square matrix for generating random // variables from this density according to // x = location + transformationMatrix*y, // where y is a spherically distributed // random variable. For multivariate normal, - // this is the square root of the covariance - // matrix. + // transformationMatrix is the square root + // of the covariance matrix. // - // gDensity The "generator" density for this - // distribution. The multivariate density - // is going to be proportional to - // gDensity.density(chi-square), where - // chi-square variable is constructed + // gDensity The "generator" distribution. The + // multivariate density is going to be + // proportional to gDensity.density(chi-square), + // where chi-square variable is constructed // as in the multivariate normal. // Must have gDensity.quantile(0.0) = 0.0. // - // hDensity Density of the r^2 variable used for - // generating random numbers from this density. - // Must have hDensity.quantile(0.0) = 0.0. + // hDensity Density of the r^2 variable. Must have + // hDensity.quantile(0.0) = 0.0. */ EllipticalDistribution(const double* location, unsigned dim, const Matrix& transformationMatrix, const AbsDistribution1D& gDensity, const AbsDistribution1D& hDensity); EllipticalDistribution(const EllipticalDistribution&); EllipticalDistribution& operator=(const EllipticalDistribution&); inline virtual EllipticalDistribution* clone() const {return new EllipticalDistribution(*this);} inline virtual ~EllipticalDistribution() {cleanup();} double density(const double* x, unsigned dim) const; void unitMap(const double* rnd, unsigned bufLen, double* x) const; inline bool mappedByQuantiles() const {return false;} unsigned random(AbsRandomGenerator& g, double* x, unsigned lenX) const; // Methods needed for I/O inline virtual gs::ClassId classId() const {return gs::ClassId(*this);} virtual bool write(std::ostream& os) const; static inline const char* classname() {return "npstat::EllipticalDistribution";} static inline unsigned version() {return 1;} static EllipticalDistribution* read(const gs::ClassId& id, std::istream& in); protected: virtual bool isEqual(const AbsDistributionND&) const; private: EllipticalDistribution(); void cleanup(); std::vector mu_; Matrix A_; Matrix InvCovmat_; AbsDistribution1D* g_; AbsDistribution1D* h_; double det_; double gNorm_; mutable std::vector buf_; }; } #endif // NPSTAT_ELLIPTICALDISTRIBUTION_HH_ Index: trunk/config.log =================================================================== --- trunk/config.log (revision 816) +++ trunk/config.log (revision 817) @@ -1,1181 +1,1177 @@ This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. It was created by npstat configure 5.7.0, which was generated by GNU Autoconf 2.71. Invocation command line was $ ./configure --disable-static --with-pic ## --------- ## ## Platform. ## ## --------- ## hostname = dawn uname -m = x86_64 uname -r = 5.15.0-47-generic uname -s = Linux uname -v = #51-Ubuntu SMP Thu Aug 11 07:51:15 UTC 2022 /usr/bin/uname -p = x86_64 /bin/uname -X = unknown /bin/arch = x86_64 /usr/bin/arch -k = unknown /usr/convex/getsysinfo = unknown /usr/bin/hostinfo = unknown /bin/machine = unknown /usr/bin/oslevel = unknown /bin/universe = unknown PATH: /home/igv/bin/ PATH: /home/igv/local/bin/ PATH: /usr/local/anaconda3/bin/ PATH: /usr/local/bin/ PATH: /usr/local/root/bin/ PATH: /usr/local/bin/ PATH: /bin/ PATH: /usr/bin/ PATH: /sbin/ PATH: /usr/sbin/ PATH: ./ ## ----------- ## ## Core tests. ## ## ----------- ## configure:2816: looking for aux files: compile ltmain.sh config.guess config.sub missing install-sh configure:2829: trying ./ configure:2858: ./compile found configure:2858: ./ltmain.sh found configure:2858: ./config.guess found configure:2858: ./config.sub found configure:2858: ./missing found configure:2840: ./install-sh found configure:2988: checking for a BSD-compatible install configure:3061: result: /bin/install -c configure:3072: checking whether build environment is sane configure:3127: result: yes configure:3286: checking for a race-free mkdir -p configure:3330: result: /bin/mkdir -p configure:3337: checking for gawk configure:3372: result: no configure:3337: checking for mawk configure:3358: found /bin/mawk configure:3369: result: mawk configure:3380: checking whether make sets $(MAKE) configure:3403: result: yes configure:3433: checking whether make supports nested variables configure:3451: result: yes configure:3650: checking for pkg-config configure:3673: found /bin/pkg-config configure:3685: result: /bin/pkg-config configure:3710: checking pkg-config is at least version 0.9.0 configure:3713: result: yes configure:3723: checking for fftw3 >= 3.1.2 geners >= 1.3.0 kstest >= 2.0.0 configure:3730: $PKG_CONFIG --exists --print-errors "fftw3 >= 3.1.2 geners >= 1.3.0 kstest >= 2.0.0" configure:3733: $? = 0 configure:3747: $PKG_CONFIG --exists --print-errors "fftw3 >= 3.1.2 geners >= 1.3.0 kstest >= 2.0.0" configure:3750: $? = 0 configure:3808: result: yes configure:3882: checking for g++ configure:3903: found /bin/g++ configure:3914: result: g++ configure:3941: checking for C++ compiler version configure:3950: g++ --version >&5 g++ (Ubuntu 11.2.0-19ubuntu1) 11.2.0 Copyright (C) 2021 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. configure:3961: $? = 0 configure:3950: g++ -v >&5 Using built-in specs. COLLECT_GCC=g++ COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa OFFLOAD_TARGET_DEFAULT=1 Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion='Ubuntu 11.2.0-19ubuntu1' --with-bugurl=file:///usr/share/doc/gcc-11/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-11 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --enable-cet --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-gcn/usr --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu --with-build-config=bootstrap-lto-lean --enable-link-serialization=2 Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 11.2.0 (Ubuntu 11.2.0-19ubuntu1) ... rest of stderr output deleted ... configure:3961: $? = 0 configure:3950: g++ -V >&5 g++: error: unrecognized command-line option '-V' g++: fatal error: no input files compilation terminated. configure:3961: $? = 1 configure:3950: g++ -qversion >&5 g++: error: unrecognized command-line option '-qversion'; did you mean '--version'? g++: fatal error: no input files compilation terminated. configure:3961: $? = 1 configure:3981: checking whether the C++ compiler works configure:4003: g++ -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:4007: $? = 0 configure:4057: result: yes configure:4060: checking for C++ compiler default output file name configure:4062: result: a.out configure:4068: checking for suffix of executables configure:4075: g++ -o conftest -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:4079: $? = 0 configure:4102: result: configure:4124: checking whether we are cross compiling configure:4132: g++ -o conftest -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:4136: $? = 0 configure:4143: ./conftest configure:4147: $? = 0 configure:4162: result: no configure:4167: checking for suffix of object files configure:4190: g++ -c -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:4194: $? = 0 configure:4216: result: o configure:4220: checking whether the compiler supports GNU C++ configure:4240: g++ -c -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:4240: $? = 0 configure:4250: result: yes configure:4261: checking whether g++ accepts -g configure:4282: g++ -c -g conftest.cpp >&5 configure:4282: $? = 0 configure:4326: result: yes configure:4346: checking for g++ option to enable C++11 features configure:4361: g++ -c -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 conftest.cpp: In function 'int main(int, char**)': conftest.cpp:134:8: error: unused variable 'a1' [-Werror=unused-variable] 134 | auto a1 = 6538; | ^~ conftest.cpp:141:16: error: unused variable 'a4' [-Werror=unused-variable] 141 | decltype(a2) a4 = 34895.034; | ^~ conftest.cpp:145:9: error: unused variable 'sa' [-Werror=unused-variable] 145 | short sa[cxx11test::get_val()] = { 0 }; | ^~ conftest.cpp:149:23: error: unused variable 'il' [-Werror=unused-variable] 149 | cxx11test::testinit il = { 4323, 435234.23544 }; | ^~ conftest.cpp:170:8: error: unused variable 'a' [-Werror=unused-variable] 170 | auto a = sum(1); | ^ conftest.cpp:171:8: error: unused variable 'b' [-Werror=unused-variable] 171 | auto b = sum(1, 2); | ^ conftest.cpp:172:8: error: unused variable 'c' [-Werror=unused-variable] 172 | auto c = sum(1.0, 2.0, 3.0); | ^ conftest.cpp:177:25: error: empty parentheses were disambiguated as a function declaration [-Werror=vexing-parse] 177 | cxx11test::delegate d2(); | ^~ conftest.cpp:177:25: note: remove parentheses to default-initialize a variable 177 | cxx11test::delegate d2(); | ^~ | -- conftest.cpp:177:25: note: or replace parentheses with braces to value-initialize a variable conftest.cpp:186:9: error: unused variable 'c' [-Werror=unused-variable] 186 | char *c = nullptr; | ^ conftest.cpp:194:15: error: unused variable 'utf8' [-Werror=unused-variable] 194 | char const *utf8 = u8"UTF-8 string \u2500"; | ^~~~ conftest.cpp:195:19: error: unused variable 'utf16' [-Werror=unused-variable] 195 | char16_t const *utf16 = u"UTF-8 string \u2500"; | ^~~~~ conftest.cpp:196:19: error: unused variable 'utf32' [-Werror=unused-variable] 196 | char32_t const *utf32 = U"UTF-32 string \u2500"; | ^~~~~ cc1plus: all warnings being treated as errors configure:4361: $? = 1 configure: failed program was: | /* confdefs.h */ | #define PACKAGE_NAME "npstat" | #define PACKAGE_TARNAME "npstat" | #define PACKAGE_VERSION "5.7.0" | #define PACKAGE_STRING "npstat 5.7.0" | #define PACKAGE_BUGREPORT "" | #define PACKAGE_URL "" | #define PACKAGE "npstat" | #define VERSION "5.7.0" | /* end confdefs.h. */ | | // Does the compiler advertise C++98 conformance? | #if !defined __cplusplus || __cplusplus < 199711L | # error "Compiler does not advertise C++98 conformance" | #endif | | // These inclusions are to reject old compilers that | // lack the unsuffixed header files. | #include | #include | | // and are *not* freestanding headers in C++98. | extern void assert (int); | namespace std { | extern int strcmp (const char *, const char *); | } | | // Namespaces, exceptions, and templates were all added after "C++ 2.0". | using std::exception; | using std::strcmp; | | namespace { | | void test_exception_syntax() | { | try { | throw "test"; | } catch (const char *s) { | // Extra parentheses suppress a warning when building autoconf itself, | // due to lint rules shared with more typical C programs. | assert (!(strcmp) (s, "test")); | } | } | | template struct test_template | { | T const val; | explicit test_template(T t) : val(t) {} | template T add(U u) { return static_cast(u) + val; } | }; | | } // anonymous namespace | | | // Does the compiler advertise C++ 2011 conformance? | #if !defined __cplusplus || __cplusplus < 201103L | # error "Compiler does not advertise C++11 conformance" | #endif | | namespace cxx11test | { | constexpr int get_val() { return 20; } | | struct testinit | { | int i; | double d; | }; | | class delegate | { | public: | delegate(int n) : n(n) {} | delegate(): delegate(2354) {} | | virtual int getval() { return this->n; }; | protected: | int n; | }; | | class overridden : public delegate | { | public: | overridden(int n): delegate(n) {} | virtual int getval() override final { return this->n * 2; } | }; | | class nocopy | { | public: | nocopy(int i): i(i) {} | nocopy() = default; | nocopy(const nocopy&) = delete; | nocopy & operator=(const nocopy&) = delete; | private: | int i; | }; | | // for testing lambda expressions | template Ret eval(Fn f, Ret v) | { | return f(v); | } | | // for testing variadic templates and trailing return types | template auto sum(V first) -> V | { | return first; | } | template auto sum(V first, Args... rest) -> V | { | return first + sum(rest...); | } | } | | | int | main (int argc, char **argv) | { | int ok = 0; | | assert (argc); | assert (! argv[0]); | { | test_exception_syntax (); | test_template tt (2.0); | assert (tt.add (4) == 6.0); | assert (true && !false); | } | | | { | // Test auto and decltype | auto a1 = 6538; | auto a2 = 48573953.4; | auto a3 = "String literal"; | | int total = 0; | for (auto i = a3; *i; ++i) { total += *i; } | | decltype(a2) a4 = 34895.034; | } | { | // Test constexpr | short sa[cxx11test::get_val()] = { 0 }; | } | { | // Test initializer lists | cxx11test::testinit il = { 4323, 435234.23544 }; | } | { | // Test range-based for | int array[] = {9, 7, 13, 15, 4, 18, 12, 10, 5, 3, | 14, 19, 17, 8, 6, 20, 16, 2, 11, 1}; | for (auto &x : array) { x += 23; } | } | { | // Test lambda expressions | using cxx11test::eval; | assert (eval ([](int x) { return x*2; }, 21) == 42); | double d = 2.0; | assert (eval ([&](double x) { return d += x; }, 3.0) == 5.0); | assert (d == 5.0); | assert (eval ([=](double x) mutable { return d += x; }, 4.0) == 9.0); | assert (d == 5.0); | } | { | // Test use of variadic templates | using cxx11test::sum; | auto a = sum(1); | auto b = sum(1, 2); | auto c = sum(1.0, 2.0, 3.0); | } | { | // Test constructor delegation | cxx11test::delegate d1; | cxx11test::delegate d2(); | cxx11test::delegate d3(45); | } | { | // Test override and final | cxx11test::overridden o1(55464); | } | { | // Test nullptr | char *c = nullptr; | } | { | // Test template brackets | test_template<::test_template> v(test_template(12)); | } | { | // Unicode literals | char const *utf8 = u8"UTF-8 string \u2500"; | char16_t const *utf16 = u"UTF-8 string \u2500"; | char32_t const *utf32 = U"UTF-32 string \u2500"; | } | | return ok; | } | configure:4379: result: none needed configure:4446: checking whether make supports the include directive configure:4461: make -f confmf.GNU && cat confinc.out this is the am__doit target configure:4464: $? = 0 configure:4483: result: yes (GNU style) configure:4509: checking dependency style of g++ configure:4621: result: gcc3 configure:4695: checking for g77 configure:4716: found /home/igv/bin/g77 configure:4727: result: g77 configure:4753: checking for Fortran 77 compiler version configure:4762: g77 --version >&5 GNU Fortran (Ubuntu 11.2.0-19ubuntu1) 11.2.0 Copyright (C) 2021 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. configure:4773: $? = 0 configure:4762: g77 -v >&5 Using built-in specs. COLLECT_GCC=g77 COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa OFFLOAD_TARGET_DEFAULT=1 Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion='Ubuntu 11.2.0-19ubuntu1' --with-bugurl=file:///usr/share/doc/gcc-11/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-11 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --enable-cet --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-gcn/usr --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu --with-build-config=bootstrap-lto-lean --enable-link-serialization=2 Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 11.2.0 (Ubuntu 11.2.0-19ubuntu1) ... rest of stderr output deleted ... configure:4773: $? = 0 configure:4762: g77 -V >&5 g77: error: unrecognized command-line option '-V' g77: fatal error: no input files compilation terminated. configure:4773: $? = 1 configure:4762: g77 -qversion >&5 g77: error: unrecognized command-line option '-qversion'; did you mean '--version'? g77: fatal error: no input files compilation terminated. configure:4773: $? = 1 configure:4782: checking whether the compiler supports GNU Fortran 77 configure:4796: g77 -c conftest.F >&5 configure:4796: $? = 0 configure:4806: result: yes configure:4814: checking whether g77 accepts -g configure:4826: g77 -c -g conftest.f >&5 configure:4826: $? = 0 configure:4835: result: yes configure:4870: checking build system type configure:4885: result: x86_64-pc-linux-gnu configure:4905: checking host system type configure:4919: result: x86_64-pc-linux-gnu configure:4944: checking how to get verbose linking output from g77 configure:4955: g77 -c -g -O2 conftest.f >&5 configure:4955: $? = 0 configure:4974: g77 -o conftest -g -O2 -v conftest.f Using built-in specs. Target: x86_64-linux-gnu Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 11.2.0 (Ubuntu 11.2.0-19ubuntu1) - /usr/lib/gcc/x86_64-linux-gnu/11/f951 conftest.f -ffixed-form -quiet -dumpbase conftest.f -dumpbase-ext .f -mtune=generic -march=x86-64 -g -O2 -version -fintrinsic-modules-path /usr/lib/gcc/x86_64-linux-gnu/11/finclude -fpre-include=/usr/include/finclude/math-vector-fortran.h -o /tmp/ccpJ01a8.s + /usr/lib/gcc/x86_64-linux-gnu/11/f951 conftest.f -ffixed-form -quiet -dumpbase conftest.f -dumpbase-ext .f -mtune=generic -march=x86-64 -g -O2 -version -fintrinsic-modules-path /usr/lib/gcc/x86_64-linux-gnu/11/finclude -fpre-include=/usr/include/finclude/math-vector-fortran.h -o /tmp/ccWFyWpt.s GNU Fortran (Ubuntu 11.2.0-19ubuntu1) version 11.2.0 (x86_64-linux-gnu) compiled by GNU C version 11.2.0, GMP version 6.2.1, MPFR version 4.1.0, MPC version 1.2.1, isl version isl-0.24-GMP GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072 GNU Fortran2008 (Ubuntu 11.2.0-19ubuntu1) version 11.2.0 (x86_64-linux-gnu) compiled by GNU C version 11.2.0, GMP version 6.2.1, MPFR version 4.1.0, MPC version 1.2.1, isl version isl-0.24-GMP GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072 - as -v --gdwarf-5 --64 -o /tmp/cc1E7ymw.o /tmp/ccpJ01a8.s + as -v --gdwarf-5 --64 -o /tmp/cccpWGsk.o /tmp/ccWFyWpt.s GNU assembler version 2.38 (x86_64-linux-gnu) using BFD version (GNU Binutils for Ubuntu) 2.38 Reading specs from /usr/lib/gcc/x86_64-linux-gnu/11/libgfortran.spec rename spec lib to liborig - /usr/lib/gcc/x86_64-linux-gnu/11/collect2 -plugin /usr/lib/gcc/x86_64-linux-gnu/11/liblto_plugin.so -plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper -plugin-opt=-fresolution=/tmp/ccZqVhst.res -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lquadmath -plugin-opt=-pass-through=-lm -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lc -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc --build-id --eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker /lib64/ld-linux-x86-64.so.2 -pie -z now -z relro -o conftest /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/Scrt1.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crti.o /usr/lib/gcc/x86_64-linux-gnu/11/crtbeginS.o -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. /tmp/cc1E7ymw.o -lgfortran -lm -lgcc_s -lgcc -lquadmath -lm -lgcc_s -lgcc -lc -lgcc_s -lgcc /usr/lib/gcc/x86_64-linux-gnu/11/crtendS.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crtn.o + /usr/lib/gcc/x86_64-linux-gnu/11/collect2 -plugin /usr/lib/gcc/x86_64-linux-gnu/11/liblto_plugin.so -plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper -plugin-opt=-fresolution=/tmp/ccGauNbE.res -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lquadmath -plugin-opt=-pass-through=-lm -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lc -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc --build-id --eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker /lib64/ld-linux-x86-64.so.2 -pie -z now -z relro -o conftest /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/Scrt1.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crti.o /usr/lib/gcc/x86_64-linux-gnu/11/crtbeginS.o -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. /tmp/cccpWGsk.o -lgfortran -lm -lgcc_s -lgcc -lquadmath -lm -lgcc_s -lgcc -lc -lgcc_s -lgcc /usr/lib/gcc/x86_64-linux-gnu/11/crtendS.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crtn.o configure:5057: result: -v configure:5059: checking for Fortran 77 libraries of g77 configure:5083: g77 -o conftest -g -O2 -v conftest.f Using built-in specs. Target: x86_64-linux-gnu Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 11.2.0 (Ubuntu 11.2.0-19ubuntu1) - /usr/lib/gcc/x86_64-linux-gnu/11/f951 conftest.f -ffixed-form -quiet -dumpbase conftest.f -dumpbase-ext .f -mtune=generic -march=x86-64 -g -O2 -version -fintrinsic-modules-path /usr/lib/gcc/x86_64-linux-gnu/11/finclude -fpre-include=/usr/include/finclude/math-vector-fortran.h -o /tmp/cc0cEw3q.s + /usr/lib/gcc/x86_64-linux-gnu/11/f951 conftest.f -ffixed-form -quiet -dumpbase conftest.f -dumpbase-ext .f -mtune=generic -march=x86-64 -g -O2 -version -fintrinsic-modules-path /usr/lib/gcc/x86_64-linux-gnu/11/finclude -fpre-include=/usr/include/finclude/math-vector-fortran.h -o /tmp/cczaUUNt.s GNU Fortran (Ubuntu 11.2.0-19ubuntu1) version 11.2.0 (x86_64-linux-gnu) compiled by GNU C version 11.2.0, GMP version 6.2.1, MPFR version 4.1.0, MPC version 1.2.1, isl version isl-0.24-GMP GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072 GNU Fortran2008 (Ubuntu 11.2.0-19ubuntu1) version 11.2.0 (x86_64-linux-gnu) compiled by GNU C version 11.2.0, GMP version 6.2.1, MPFR version 4.1.0, MPC version 1.2.1, isl version isl-0.24-GMP GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072 - as -v --gdwarf-5 --64 -o /tmp/ccI8PK9A.o /tmp/cc0cEw3q.s + as -v --gdwarf-5 --64 -o /tmp/ccT8pSJA.o /tmp/cczaUUNt.s GNU assembler version 2.38 (x86_64-linux-gnu) using BFD version (GNU Binutils for Ubuntu) 2.38 Reading specs from /usr/lib/gcc/x86_64-linux-gnu/11/libgfortran.spec rename spec lib to liborig - /usr/lib/gcc/x86_64-linux-gnu/11/collect2 -plugin /usr/lib/gcc/x86_64-linux-gnu/11/liblto_plugin.so -plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper -plugin-opt=-fresolution=/tmp/ccoRoUrh.res -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lquadmath -plugin-opt=-pass-through=-lm -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lc -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc --build-id --eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker /lib64/ld-linux-x86-64.so.2 -pie -z now -z relro -o conftest /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/Scrt1.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crti.o /usr/lib/gcc/x86_64-linux-gnu/11/crtbeginS.o -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. /tmp/ccI8PK9A.o -lgfortran -lm -lgcc_s -lgcc -lquadmath -lm -lgcc_s -lgcc -lc -lgcc_s -lgcc /usr/lib/gcc/x86_64-linux-gnu/11/crtendS.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crtn.o + /usr/lib/gcc/x86_64-linux-gnu/11/collect2 -plugin /usr/lib/gcc/x86_64-linux-gnu/11/liblto_plugin.so -plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper -plugin-opt=-fresolution=/tmp/cccp7i9n.res -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lquadmath -plugin-opt=-pass-through=-lm -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lc -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc --build-id --eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker /lib64/ld-linux-x86-64.so.2 -pie -z now -z relro -o conftest /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/Scrt1.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crti.o /usr/lib/gcc/x86_64-linux-gnu/11/crtbeginS.o -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. /tmp/ccT8pSJA.o -lgfortran -lm -lgcc_s -lgcc -lquadmath -lm -lgcc_s -lgcc -lc -lgcc_s -lgcc /usr/lib/gcc/x86_64-linux-gnu/11/crtendS.o /usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu/crtn.o configure:5299: result: -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. -lgfortran -lm -lquadmath configure:5362: checking how to print strings configure:5389: result: printf configure:5472: checking for gcc configure:5493: found /bin/gcc configure:5504: result: gcc configure:5857: checking for C compiler version configure:5866: gcc --version >&5 gcc (Ubuntu 11.2.0-19ubuntu1) 11.2.0 Copyright (C) 2021 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. configure:5877: $? = 0 configure:5866: gcc -v >&5 Using built-in specs. COLLECT_GCC=gcc COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/11/lto-wrapper OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa OFFLOAD_TARGET_DEFAULT=1 Target: x86_64-linux-gnu Configured with: ../src/configure -v --with-pkgversion='Ubuntu 11.2.0-19ubuntu1' --with-bugurl=file:///usr/share/doc/gcc-11/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-11 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --enable-cet --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-11-gBFGDP/gcc-11-11.2.0/debian/tmp-gcn/usr --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu --with-build-config=bootstrap-lto-lean --enable-link-serialization=2 Thread model: posix Supported LTO compression algorithms: zlib zstd gcc version 11.2.0 (Ubuntu 11.2.0-19ubuntu1) ... rest of stderr output deleted ... configure:5877: $? = 0 configure:5866: gcc -V >&5 gcc: error: unrecognized command-line option '-V' gcc: fatal error: no input files compilation terminated. configure:5877: $? = 1 configure:5866: gcc -qversion >&5 gcc: error: unrecognized command-line option '-qversion'; did you mean '--version'? gcc: fatal error: no input files compilation terminated. configure:5877: $? = 1 configure:5866: gcc -version >&5 gcc: error: unrecognized command-line option '-version' gcc: fatal error: no input files compilation terminated. configure:5877: $? = 1 configure:5881: checking whether the compiler supports GNU C configure:5901: gcc -c conftest.c >&5 configure:5901: $? = 0 configure:5911: result: yes configure:5922: checking whether gcc accepts -g configure:5943: gcc -c -g conftest.c >&5 configure:5943: $? = 0 configure:5987: result: yes configure:6007: checking for gcc option to enable C11 features configure:6022: gcc -c -g -O2 conftest.c >&5 configure:6022: $? = 0 configure:6040: result: none needed configure:6156: checking whether gcc understands -c and -o together configure:6179: gcc -c conftest.c -o conftest2.o configure:6182: $? = 0 configure:6179: gcc -c conftest.c -o conftest2.o configure:6182: $? = 0 configure:6194: result: yes configure:6213: checking dependency style of gcc configure:6325: result: gcc3 configure:6340: checking for a sed that does not truncate output configure:6410: result: /bin/sed configure:6428: checking for grep that handles long lines and -e configure:6492: result: /bin/grep configure:6497: checking for egrep configure:6565: result: /bin/grep -E configure:6570: checking for fgrep configure:6638: result: /bin/grep -F configure:6674: checking for ld used by gcc configure:6742: result: /bin/ld configure:6749: checking if the linker (/bin/ld) is GNU ld configure:6765: result: yes configure:6777: checking for BSD- or MS-compatible name lister (nm) configure:6832: result: /bin/nm -B configure:6972: checking the name lister (/bin/nm -B) interface configure:6980: gcc -c -g -O2 conftest.c >&5 configure:6983: /bin/nm -B "conftest.o" configure:6986: output 0000000000000000 B some_variable configure:6993: result: BSD nm configure:6996: checking whether ln -s works configure:7000: result: yes configure:7008: checking the maximum length of command line arguments configure:7140: result: 1572864 configure:7188: checking how to convert x86_64-pc-linux-gnu file names to x86_64-pc-linux-gnu format configure:7229: result: func_convert_file_noop configure:7236: checking how to convert x86_64-pc-linux-gnu file names to toolchain format configure:7257: result: func_convert_file_noop configure:7264: checking for /bin/ld option to reload object files configure:7272: result: -r configure:7351: checking for objdump configure:7372: found /bin/objdump configure:7383: result: objdump configure:7415: checking how to recognize dependent libraries configure:7616: result: pass_all configure:7706: checking for dlltool configure:7741: result: no configure:7771: checking how to associate runtime and link libraries configure:7799: result: printf %s\n configure:7865: checking for ar configure:7886: found /bin/ar configure:7897: result: ar configure:7934: checking for archiver @FILE support configure:7952: gcc -c -g -O2 conftest.c >&5 configure:7952: $? = 0 configure:7956: ar cr libconftest.a @conftest.lst >&5 configure:7959: $? = 0 configure:7964: ar cr libconftest.a @conftest.lst >&5 ar: conftest.o: No such file or directory configure:7967: $? = 1 configure:7979: result: @ configure:8042: checking for strip configure:8063: found /bin/strip configure:8074: result: strip configure:8151: checking for ranlib configure:8172: found /bin/ranlib configure:8183: result: ranlib configure:8285: checking command to parse /bin/nm -B output from gcc object configure:8439: gcc -c -g -O2 conftest.c >&5 configure:8442: $? = 0 configure:8446: /bin/nm -B conftest.o | sed -n -e 's/^.*[ ]\([ABCDGIRSTW][ABCDGIRSTW]*\)[ ][ ]*\([_A-Za-z][_A-Za-z0-9]*\)$/\1 \2 \2/p' | sed '/ __gnu_lto/d' > conftest.nm configure:8512: gcc -o conftest -g -O2 conftest.c conftstm.o >&5 configure:8515: $? = 0 configure:8553: result: ok configure:8600: checking for sysroot configure:8631: result: no configure:8638: checking for a working dd configure:8682: result: /bin/dd configure:8686: checking how to truncate binary pipes configure:8702: result: /bin/dd bs=4096 count=1 configure:8839: gcc -c -g -O2 conftest.c >&5 configure:8842: $? = 0 configure:9039: checking for mt configure:9060: found /bin/mt configure:9071: result: mt configure:9094: checking if mt is a manifest tool configure:9101: mt '-?' configure:9109: result: no configure:9839: checking for stdio.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for stdlib.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for string.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for inttypes.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for stdint.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for strings.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for sys/stat.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for sys/types.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9839: checking for unistd.h configure:9839: gcc -c -g -O2 conftest.c >&5 configure:9839: $? = 0 configure:9839: result: yes configure:9864: checking for dlfcn.h configure:9864: gcc -c -g -O2 conftest.c >&5 configure:9864: $? = 0 configure:9864: result: yes configure:10134: checking for objdir configure:10150: result: .libs configure:10414: checking if gcc supports -fno-rtti -fno-exceptions configure:10433: gcc -c -g -O2 -fno-rtti -fno-exceptions conftest.c >&5 cc1: warning: command-line option '-fno-rtti' is valid for C++/D/ObjC++ but not for C configure:10437: $? = 0 configure:10450: result: no configure:10814: checking for gcc option to produce PIC configure:10822: result: -fPIC -DPIC configure:10830: checking if gcc PIC flag -fPIC -DPIC works configure:10849: gcc -c -g -O2 -fPIC -DPIC -DPIC conftest.c >&5 configure:10853: $? = 0 configure:10866: result: yes configure:10895: checking if gcc static flag -static works configure:10924: result: yes configure:10939: checking if gcc supports -c -o file.o configure:10961: gcc -c -g -O2 -o out/conftest2.o conftest.c >&5 configure:10965: $? = 0 configure:10987: result: yes configure:10995: checking if gcc supports -c -o file.o configure:11043: result: yes configure:11076: checking whether the gcc linker (/bin/ld -m elf_x86_64) supports shared libraries configure:12346: result: yes configure:12383: checking whether -lc should be explicitly linked in configure:12392: gcc -c -g -O2 conftest.c >&5 configure:12395: $? = 0 configure:12410: gcc -shared -fPIC -DPIC conftest.o -v -Wl,-soname -Wl,conftest -o conftest 2\>\&1 \| /bin/grep -lc \>/dev/null 2\>\&1 configure:12413: $? = 0 configure:12427: result: no configure:12587: checking dynamic linker characteristics configure:13169: gcc -o conftest -g -O2 -Wl,-rpath -Wl,/foo conftest.c >&5 configure:13169: $? = 0 configure:13420: result: GNU/Linux ld.so configure:13542: checking how to hardcode library paths into programs configure:13567: result: immediate configure:14119: checking whether stripping libraries is possible configure:14124: result: yes configure:14159: checking if libtool supports shared libraries configure:14161: result: yes configure:14164: checking whether to build shared libraries configure:14189: result: yes configure:14192: checking whether to build static libraries configure:14196: result: no configure:14219: checking how to run the C++ preprocessor configure:14241: g++ -E conftest.cpp configure:14241: $? = 0 configure:14256: g++ -E conftest.cpp conftest.cpp:23:10: fatal error: ac_nonexistent.h: No such file or directory 23 | #include | ^~~~~~~~~~~~~~~~~~ compilation terminated. configure:14256: $? = 1 configure: failed program was: | /* confdefs.h */ | #define PACKAGE_NAME "npstat" | #define PACKAGE_TARNAME "npstat" | #define PACKAGE_VERSION "5.7.0" | #define PACKAGE_STRING "npstat 5.7.0" | #define PACKAGE_BUGREPORT "" | #define PACKAGE_URL "" | #define PACKAGE "npstat" | #define VERSION "5.7.0" | #define HAVE_STDIO_H 1 | #define HAVE_STDLIB_H 1 | #define HAVE_STRING_H 1 | #define HAVE_INTTYPES_H 1 | #define HAVE_STDINT_H 1 | #define HAVE_STRINGS_H 1 | #define HAVE_SYS_STAT_H 1 | #define HAVE_SYS_TYPES_H 1 | #define HAVE_UNISTD_H 1 | #define STDC_HEADERS 1 | #define HAVE_DLFCN_H 1 | #define LT_OBJDIR ".libs/" | /* end confdefs.h. */ | #include configure:14283: result: g++ -E configure:14297: g++ -E conftest.cpp configure:14297: $? = 0 configure:14312: g++ -E conftest.cpp conftest.cpp:23:10: fatal error: ac_nonexistent.h: No such file or directory 23 | #include | ^~~~~~~~~~~~~~~~~~ compilation terminated. configure:14312: $? = 1 configure: failed program was: | /* confdefs.h */ | #define PACKAGE_NAME "npstat" | #define PACKAGE_TARNAME "npstat" | #define PACKAGE_VERSION "5.7.0" | #define PACKAGE_STRING "npstat 5.7.0" | #define PACKAGE_BUGREPORT "" | #define PACKAGE_URL "" | #define PACKAGE "npstat" | #define VERSION "5.7.0" | #define HAVE_STDIO_H 1 | #define HAVE_STDLIB_H 1 | #define HAVE_STRING_H 1 | #define HAVE_INTTYPES_H 1 | #define HAVE_STDINT_H 1 | #define HAVE_STRINGS_H 1 | #define HAVE_SYS_STAT_H 1 | #define HAVE_SYS_TYPES_H 1 | #define HAVE_UNISTD_H 1 | #define STDC_HEADERS 1 | #define HAVE_DLFCN_H 1 | #define LT_OBJDIR ".libs/" | /* end confdefs.h. */ | #include configure:14477: checking for ld used by g++ configure:14545: result: /bin/ld -m elf_x86_64 configure:14552: checking if the linker (/bin/ld -m elf_x86_64) is GNU ld configure:14568: result: yes configure:14623: checking whether the g++ linker (/bin/ld -m elf_x86_64) supports shared libraries configure:15700: result: yes configure:15736: g++ -c -std=c++11 -O3 -Wall -W -Werror conftest.cpp >&5 configure:15739: $? = 0 configure:16220: checking for g++ option to produce PIC configure:16228: result: -fPIC -DPIC configure:16236: checking if g++ PIC flag -fPIC -DPIC works configure:16255: g++ -c -std=c++11 -O3 -Wall -W -Werror -fPIC -DPIC -DPIC conftest.cpp >&5 configure:16259: $? = 0 configure:16272: result: yes configure:16295: checking if g++ static flag -static works configure:16324: result: yes configure:16336: checking if g++ supports -c -o file.o configure:16358: g++ -c -std=c++11 -O3 -Wall -W -Werror -o out/conftest2.o conftest.cpp >&5 configure:16362: $? = 0 configure:16384: result: yes configure:16389: checking if g++ supports -c -o file.o configure:16437: result: yes configure:16467: checking whether the g++ linker (/bin/ld -m elf_x86_64) supports shared libraries configure:16510: result: yes configure:16652: checking dynamic linker characteristics configure:17412: result: GNU/Linux ld.so configure:17477: checking how to hardcode library paths into programs configure:17502: result: immediate configure:17643: checking if libtool supports shared libraries configure:17645: result: yes configure:17648: checking whether to build shared libraries configure:17672: result: yes configure:17675: checking whether to build static libraries configure:17679: result: no configure:18037: checking for g77 option to produce PIC configure:18045: result: -fPIC configure:18053: checking if g77 PIC flag -fPIC works configure:18072: g77 -c -g -O2 -fPIC conftest.f >&5 configure:18076: $? = 0 configure:18089: result: yes configure:18112: checking if g77 static flag -static works configure:18141: result: yes configure:18153: checking if g77 supports -c -o file.o configure:18175: g77 -c -g -O2 -o out/conftest2.o conftest.f >&5 configure:18179: $? = 0 configure:18201: result: yes configure:18206: checking if g77 supports -c -o file.o configure:18254: result: yes configure:18284: checking whether the g77 linker (/bin/ld -m elf_x86_64) supports shared libraries configure:19503: result: yes configure:19645: checking dynamic linker characteristics configure:20399: result: GNU/Linux ld.so configure:20464: checking how to hardcode library paths into programs configure:20489: result: immediate configure:20681: checking that generated files are newer than configure configure:20687: result: done configure:20714: creating ./config.status ## ---------------------- ## ## Running config.status. ## ## ---------------------- ## This file was extended by npstat config.status 5.7.0, which was generated by GNU Autoconf 2.71. Invocation command line was CONFIG_FILES = CONFIG_HEADERS = CONFIG_LINKS = CONFIG_COMMANDS = $ ./config.status on dawn config.status:1138: creating Makefile config.status:1138: creating npstat/nm/Makefile config.status:1138: creating npstat/rng/Makefile config.status:1138: creating npstat/stat/Makefile config.status:1138: creating npstat/wrap/Makefile config.status:1138: creating npstat/interfaces/Makefile config.status:1138: creating npstat/emsunfold/Makefile config.status:1138: creating npstat/Makefile config.status:1138: creating examples/C++/Makefile config.status:1138: creating npstat/swig/Makefile config.status:1138: creating npstat.pc config.status:1310: executing depfiles commands config.status:1387: cd npstat/nm && sed -e '/# am--include-marker/d' Makefile | make -f - am--depfiles -make: Nothing to be done for 'am--depfiles'. config.status:1392: $? = 0 config.status:1387: cd npstat/rng && sed -e '/# am--include-marker/d' Makefile | make -f - am--depfiles -make: Nothing to be done for 'am--depfiles'. config.status:1392: $? = 0 config.status:1387: cd npstat/stat && sed -e '/# am--include-marker/d' Makefile | make -f - am--depfiles -make: Nothing to be done for 'am--depfiles'. config.status:1392: $? = 0 config.status:1387: cd examples/C++ && sed -e '/# am--include-marker/d' Makefile | make -f - am--depfiles -make: Nothing to be done for 'am--depfiles'. config.status:1392: $? = 0 config.status:1387: cd npstat/swig && sed -e '/# am--include-marker/d' Makefile | make -f - am--depfiles make: Nothing to be done for 'am--depfiles'. config.status:1392: $? = 0 config.status:1310: executing libtool commands ## ---------------- ## ## Cache variables. ## ## ---------------- ## ac_cv_build=x86_64-pc-linux-gnu ac_cv_c_compiler_gnu=yes ac_cv_cxx_compiler_gnu=yes ac_cv_env_CCC_set= ac_cv_env_CCC_value= ac_cv_env_CC_set= ac_cv_env_CC_value= ac_cv_env_CFLAGS_set= ac_cv_env_CFLAGS_value= ac_cv_env_CPPFLAGS_set= ac_cv_env_CPPFLAGS_value= ac_cv_env_CXXCPP_set= ac_cv_env_CXXCPP_value= ac_cv_env_CXXFLAGS_set=set ac_cv_env_CXXFLAGS_value='-std=c++11 -O3 -Wall -W -Werror' ac_cv_env_CXX_set= ac_cv_env_CXX_value= ac_cv_env_DEPS_CFLAGS_set= ac_cv_env_DEPS_CFLAGS_value= ac_cv_env_DEPS_LIBS_set= ac_cv_env_DEPS_LIBS_value= ac_cv_env_F77_set= ac_cv_env_F77_value= ac_cv_env_FFLAGS_set= ac_cv_env_FFLAGS_value= ac_cv_env_LDFLAGS_set= ac_cv_env_LDFLAGS_value= ac_cv_env_LIBS_set= ac_cv_env_LIBS_value= ac_cv_env_LT_SYS_LIBRARY_PATH_set= ac_cv_env_LT_SYS_LIBRARY_PATH_value= ac_cv_env_PKG_CONFIG_LIBDIR_set= ac_cv_env_PKG_CONFIG_LIBDIR_value= ac_cv_env_PKG_CONFIG_PATH_set=set ac_cv_env_PKG_CONFIG_PATH_value=/usr/local/lib/pkgconfig ac_cv_env_PKG_CONFIG_set= ac_cv_env_PKG_CONFIG_value= ac_cv_env_build_alias_set= ac_cv_env_build_alias_value= ac_cv_env_host_alias_set= ac_cv_env_host_alias_value= ac_cv_env_target_alias_set= ac_cv_env_target_alias_value= ac_cv_f77_compiler_gnu=yes ac_cv_f77_libs=' -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. -lgfortran -lm -lquadmath' ac_cv_header_dlfcn_h=yes ac_cv_header_inttypes_h=yes ac_cv_header_stdint_h=yes ac_cv_header_stdio_h=yes ac_cv_header_stdlib_h=yes ac_cv_header_string_h=yes ac_cv_header_strings_h=yes ac_cv_header_sys_stat_h=yes ac_cv_header_sys_types_h=yes ac_cv_header_unistd_h=yes ac_cv_host=x86_64-pc-linux-gnu ac_cv_objext=o ac_cv_path_EGREP='/bin/grep -E' ac_cv_path_FGREP='/bin/grep -F' ac_cv_path_GREP=/bin/grep ac_cv_path_SED=/bin/sed ac_cv_path_ac_pt_PKG_CONFIG=/bin/pkg-config ac_cv_path_install='/bin/install -c' ac_cv_path_lt_DD=/bin/dd ac_cv_path_mkdir=/bin/mkdir ac_cv_prog_AWK=mawk ac_cv_prog_CXXCPP='g++ -E' ac_cv_prog_ac_ct_AR=ar ac_cv_prog_ac_ct_CC=gcc ac_cv_prog_ac_ct_CXX=g++ ac_cv_prog_ac_ct_F77=g77 ac_cv_prog_ac_ct_MANIFEST_TOOL=mt ac_cv_prog_ac_ct_OBJDUMP=objdump ac_cv_prog_ac_ct_RANLIB=ranlib ac_cv_prog_ac_ct_STRIP=strip ac_cv_prog_cc_c11= ac_cv_prog_cc_g=yes ac_cv_prog_cc_stdc= ac_cv_prog_cxx_11=no ac_cv_prog_cxx_g=yes ac_cv_prog_cxx_stdcxx= ac_cv_prog_f77_g=yes ac_cv_prog_f77_v=-v ac_cv_prog_make_make_set=yes am_cv_CC_dependencies_compiler_type=gcc3 am_cv_CXX_dependencies_compiler_type=gcc3 am_cv_make_support_nested_variables=yes am_cv_prog_cc_c_o=yes lt_cv_ar_at_file=@ lt_cv_archive_cmds_need_lc=no lt_cv_deplibs_check_method=pass_all lt_cv_file_magic_cmd='$MAGIC_CMD' lt_cv_file_magic_test_file= lt_cv_ld_reload_flag=-r lt_cv_nm_interface='BSD nm' lt_cv_objdir=.libs lt_cv_path_LD=/bin/ld lt_cv_path_LDCXX='/bin/ld -m elf_x86_64' lt_cv_path_NM='/bin/nm -B' lt_cv_path_mainfest_tool=no lt_cv_prog_compiler_c_o=yes lt_cv_prog_compiler_c_o_CXX=yes lt_cv_prog_compiler_c_o_F77=yes lt_cv_prog_compiler_pic='-fPIC -DPIC' lt_cv_prog_compiler_pic_CXX='-fPIC -DPIC' lt_cv_prog_compiler_pic_F77=-fPIC lt_cv_prog_compiler_pic_works=yes lt_cv_prog_compiler_pic_works_CXX=yes lt_cv_prog_compiler_pic_works_F77=yes lt_cv_prog_compiler_rtti_exceptions=no lt_cv_prog_compiler_static_works=yes lt_cv_prog_compiler_static_works_CXX=yes lt_cv_prog_compiler_static_works_F77=yes lt_cv_prog_gnu_ld=yes lt_cv_prog_gnu_ldcxx=yes lt_cv_sharedlib_from_linklib_cmd='printf %s\n' lt_cv_shlibpath_overrides_runpath=yes lt_cv_sys_global_symbol_pipe='sed -n -e '\''s/^.*[ ]\([ABCDGIRSTW][ABCDGIRSTW]*\)[ ][ ]*\([_A-Za-z][_A-Za-z0-9]*\)$/\1 \2 \2/p'\'' | sed '\''/ __gnu_lto/d'\''' lt_cv_sys_global_symbol_to_c_name_address='sed -n -e '\''s/^: \(.*\) .*$/ {"\1", (void *) 0},/p'\'' -e '\''s/^[ABCDGIRSTW][ABCDGIRSTW]* .* \(.*\)$/ {"\1", (void *) \&\1},/p'\''' lt_cv_sys_global_symbol_to_c_name_address_lib_prefix='sed -n -e '\''s/^: \(.*\) .*$/ {"\1", (void *) 0},/p'\'' -e '\''s/^[ABCDGIRSTW][ABCDGIRSTW]* .* \(lib.*\)$/ {"\1", (void *) \&\1},/p'\'' -e '\''s/^[ABCDGIRSTW][ABCDGIRSTW]* .* \(.*\)$/ {"lib\1", (void *) \&\1},/p'\''' lt_cv_sys_global_symbol_to_cdecl='sed -n -e '\''s/^T .* \(.*\)$/extern int \1();/p'\'' -e '\''s/^[ABCDGIRSTW][ABCDGIRSTW]* .* \(.*\)$/extern char \1;/p'\''' lt_cv_sys_global_symbol_to_import= lt_cv_sys_max_cmd_len=1572864 lt_cv_to_host_file_cmd=func_convert_file_noop lt_cv_to_tool_file_cmd=func_convert_file_noop lt_cv_truncate_bin='/bin/dd bs=4096 count=1' pkg_cv_DEPS_CFLAGS=-I/usr/local/include pkg_cv_DEPS_LIBS='-L/usr/local/lib -lfftw3 -lgeners -lkstest' ## ----------------- ## ## Output variables. ## ## ----------------- ## ACLOCAL='${SHELL} '\''/home/igv/Hepforge/npstat/trunk/missing'\'' aclocal-1.16' AMDEPBACKSLASH='\' AMDEP_FALSE='#' AMDEP_TRUE='' AMTAR='$${TAR-tar}' AM_BACKSLASH='\' AM_DEFAULT_V='$(AM_DEFAULT_VERBOSITY)' AM_DEFAULT_VERBOSITY='1' AM_V='$(V)' AR='ar' AUTOCONF='${SHELL} '\''/home/igv/Hepforge/npstat/trunk/missing'\'' autoconf' AUTOHEADER='${SHELL} '\''/home/igv/Hepforge/npstat/trunk/missing'\'' autoheader' AUTOMAKE='${SHELL} '\''/home/igv/Hepforge/npstat/trunk/missing'\'' automake-1.16' AWK='mawk' CC='gcc' CCDEPMODE='depmode=gcc3' CFLAGS='-g -O2' CPPFLAGS='' CSCOPE='cscope' CTAGS='ctags' CXX='g++' CXXCPP='g++ -E' CXXDEPMODE='depmode=gcc3' CXXFLAGS='-std=c++11 -O3 -Wall -W -Werror' CYGPATH_W='echo' DEFS='-DPACKAGE_NAME=\"npstat\" -DPACKAGE_TARNAME=\"npstat\" -DPACKAGE_VERSION=\"5.7.0\" -DPACKAGE_STRING=\"npstat\ 5.7.0\" -DPACKAGE_BUGREPORT=\"\" -DPACKAGE_URL=\"\" -DPACKAGE=\"npstat\" -DVERSION=\"5.7.0\" -DHAVE_STDIO_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_STDINT_H=1 -DHAVE_STRINGS_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_UNISTD_H=1 -DSTDC_HEADERS=1 -DHAVE_DLFCN_H=1 -DLT_OBJDIR=\".libs/\"' DEPDIR='.deps' DEPS_CFLAGS='-I/usr/local/include' DEPS_LIBS='-L/usr/local/lib -lfftw3 -lgeners -lkstest' DLLTOOL='false' DSYMUTIL='' DUMPBIN='' ECHO_C='' ECHO_N='-n' ECHO_T='' EGREP='/bin/grep -E' ETAGS='etags' EXEEXT='' F77='g77' FFLAGS='-g -O2' FGREP='/bin/grep -F' FLIBS=' -L/usr/lib/gcc/x86_64-linux-gnu/11 -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/11/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/11/../../.. -lgfortran -lm -lquadmath' GREP='/bin/grep' INSTALL_DATA='${INSTALL} -m 644' INSTALL_PROGRAM='${INSTALL}' INSTALL_SCRIPT='${INSTALL}' INSTALL_STRIP_PROGRAM='$(install_sh) -c -s' LD='/bin/ld -m elf_x86_64' LDFLAGS='' LIBOBJS='' LIBS='' LIBTOOL='$(SHELL) $(top_builddir)/libtool' LIPO='' LN_S='ln -s' LTLIBOBJS='' LT_SYS_LIBRARY_PATH='' MAKEINFO='${SHELL} '\''/home/igv/Hepforge/npstat/trunk/missing'\'' makeinfo' MANIFEST_TOOL=':' MKDIR_P='/bin/mkdir -p' NM='/bin/nm -B' NMEDIT='' OBJDUMP='objdump' OBJEXT='o' OTOOL64='' OTOOL='' PACKAGE='npstat' PACKAGE_BUGREPORT='' PACKAGE_NAME='npstat' PACKAGE_STRING='npstat 5.7.0' PACKAGE_TARNAME='npstat' PACKAGE_URL='' PACKAGE_VERSION='5.7.0' PATH_SEPARATOR=':' PKG_CONFIG='/bin/pkg-config' PKG_CONFIG_LIBDIR='' PKG_CONFIG_PATH='/usr/local/lib/pkgconfig' RANLIB='ranlib' SED='/bin/sed' SET_MAKE='' SHELL='/bin/bash' STRIP='strip' VERSION='5.7.0' ac_ct_AR='ar' ac_ct_CC='gcc' ac_ct_CXX='g++' ac_ct_DUMPBIN='' ac_ct_F77='g77' am__EXEEXT_FALSE='' am__EXEEXT_TRUE='#' am__fastdepCC_FALSE='#' am__fastdepCC_TRUE='' am__fastdepCXX_FALSE='#' am__fastdepCXX_TRUE='' am__include='include' am__isrc='' am__leading_dot='.' am__nodep='_no' am__quote='' am__tar='$${TAR-tar} chof - "$$tardir"' am__untar='$${TAR-tar} xf -' bindir='${exec_prefix}/bin' build='x86_64-pc-linux-gnu' build_alias='' build_cpu='x86_64' build_os='linux-gnu' build_vendor='pc' datadir='${datarootdir}' datarootdir='${prefix}/share' docdir='${datarootdir}/doc/${PACKAGE_TARNAME}' dvidir='${docdir}' exec_prefix='${prefix}' host='x86_64-pc-linux-gnu' host_alias='' host_cpu='x86_64' host_os='linux-gnu' host_vendor='pc' htmldir='${docdir}' includedir='${prefix}/include' infodir='${datarootdir}/info' install_sh='${SHELL} /home/igv/Hepforge/npstat/trunk/install-sh' libdir='${exec_prefix}/lib' libexecdir='${exec_prefix}/libexec' localedir='${datarootdir}/locale' localstatedir='${prefix}/var' mandir='${datarootdir}/man' mkdir_p='$(MKDIR_P)' oldincludedir='/usr/include' pdfdir='${docdir}' prefix='/usr/local' program_transform_name='s,x,x,' psdir='${docdir}' runstatedir='${localstatedir}/run' sbindir='${exec_prefix}/sbin' sharedstatedir='${prefix}/com' sysconfdir='${prefix}/etc' target_alias='' ## ----------- ## ## confdefs.h. ## ## ----------- ## /* confdefs.h */ #define PACKAGE_NAME "npstat" #define PACKAGE_TARNAME "npstat" #define PACKAGE_VERSION "5.7.0" #define PACKAGE_STRING "npstat 5.7.0" #define PACKAGE_BUGREPORT "" #define PACKAGE_URL "" #define PACKAGE "npstat" #define VERSION "5.7.0" #define HAVE_STDIO_H 1 #define HAVE_STDLIB_H 1 #define HAVE_STRING_H 1 #define HAVE_INTTYPES_H 1 #define HAVE_STDINT_H 1 #define HAVE_STRINGS_H 1 #define HAVE_SYS_STAT_H 1 #define HAVE_SYS_TYPES_H 1 #define HAVE_UNISTD_H 1 #define STDC_HEADERS 1 #define HAVE_DLFCN_H 1 #define LT_OBJDIR ".libs/" configure: exit 0