diff --git a/Makefile b/Makefile index 4abd5c6..0203628 100644 --- a/Makefile +++ b/Makefile @@ -38,7 +38,7 @@ all: $(STAT_PRODUCT) $(DYN_PRODUCT) $(DYN_PRODUCT) : $(OBJS) @echo Linking $(DYN_PRODUCT) - @$(CXX) $(LDFLAGS) $(DYN_OPT) -shared $(OBJS) -o $(DYN_PRODUCT) + @$(CXX) $(DYN_OPT) -shared $(OBJS) -o $(DYN_PRODUCT) $(LDFLAGS) $(STAT_PRODUCT) : $(OBJS) @echo Linking $(STAT_PRODUCT) @$(AR) -rcs $(STAT_PRODUCT) $(OBJS) @@ -77,6 +77,7 @@ install : cp $(LBFGSB_HEADERS) $(PREFIX)/include/$(NAME)/lbfgsb/ mkdir -p $(PREFIX)/include/$(NAME)/likelihood cp $(LIKELIHOOD_HEADERS) $(PREFIX)/include/$(NAME)/likelihood/ + cp $(BASEDIR)/phystools.pc $(PREFIX)/lib/pkgconfig/ docs : doxygen diff --git a/PhysTools/autodiff.h b/PhysTools/autodiff.h index 7009c5b..7857559 100644 --- a/PhysTools/autodiff.h +++ b/PhysTools/autodiff.h @@ -793,8 +793,9 @@ namespace autodiff{ T te(e); result.v=pow(b.v,te); const unsigned int n=detail::dimensionExtractor::nVars(result); + T scale = std::pow(b.v, te - T(1)) * e; std::transform(result.g,result.g+n,result.g, - std::bind2nd(std::multiplies(),pow(b.v,te-T(1))*e)); + [scale](T x) {return x*scale;}); return(result); } @@ -807,8 +808,9 @@ namespace autodiff{ T tb(b); result.v=pow(tb,e.v); const unsigned int n=detail::dimensionExtractor::nVars(result); + T scale = result.v*std::log(tb); std::transform(result.g,result.g+n,result.g, - std::bind2nd(std::multiplies(),result.v*log(tb))); + [scale](T x) {return x*scale;}); return(result); } @@ -1074,7 +1076,7 @@ namespace autodiff{ result.v=-result.v; const unsigned int n=detail::dimensionExtractor::nVars(result); std::transform(result.g,result.g+n,result.g, - std::bind2nd(std::multiplies(),T(-1))); + [](T x) {return -x;}); } return(result); } @@ -1088,7 +1090,7 @@ namespace autodiff{ result.v=-result.v; const unsigned int n=detail::dimensionExtractor::nVars(result); std::transform(result.g,result.g+n,result.g, - std::bind2nd(std::multiplies(),T(-1))); + [](T x) {return -x;}); } return(result); } diff --git a/PhysTools/bin_types.h b/PhysTools/bin_types.h index c1b7a82..96732c3 100644 --- a/PhysTools/bin_types.h +++ b/PhysTools/bin_types.h @@ -4,6 +4,7 @@ #ifndef PHYSTOOLS_BIN_TYPES_H #define PHYSTOOLS_BIN_TYPES_H +#include #include "PhysTools/hdf5_serialization.h" #include "PhysTools/histogram.h" diff --git a/PhysTools/dynamic_histogram.h b/PhysTools/dynamic_histogram.h index 4e57d45..8a956ae 100644 --- a/PhysTools/dynamic_histogram.h +++ b/PhysTools/dynamic_histogram.h @@ -569,7 +569,7 @@ class histogram : public histogramBase : public histogramBase #include -#ifdef __APPLE__ +#if defined(__x86_64__) || defined(_M_X64) #include -#elif __linux__ - #include + #define ENABLE_FLOAT_EXCEPTIONS() \ + _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | \ + _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | \ + _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW) + #define DISABLE_FLOAT_EXCEPTIONS() \ + _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( \ + _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | \ + _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)) +#elif defined(__linux__) + #include + #define ENABLE_FLOAT_EXCEPTIONS() \ + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) + #define DISABLE_FLOAT_EXCEPTIONS() \ + fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) +#else + #define ENABLE_FLOAT_EXCEPTIONS() + #define DISABLE_FLOAT_EXCEPTIONS() #endif #include @@ -22,6 +37,13 @@ #include #include +#include +#include +#include +#include +#include +// #define BOOST_UBLAS_TYPE_CHECK 0 + #include "../lbfgsb/lbfgsb.h" #include "../histogram.h" @@ -35,6 +57,8 @@ namespace phys_tools{ ///Tools for performing binned maximum likelihood fits namespace likelihood{ + namespace ublas = boost::numeric::ublas; + //some useful type traits template struct remove_reference_wrapper{ using type=T; }; @@ -571,7 +595,10 @@ namespace likelihood{ // Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3 // Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8 - return x-pow(x,2.0)/2.0+pow(x,3.0)/3.0-pow(x,4.0)/4.0; + T x2 = x*x; + T x3 = x2*x; + T x4 = x3*x; + return x-x2/2.0+x3/3.0-x4/4.0; }; template @@ -593,7 +620,10 @@ namespace likelihood{ // Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3 // Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8 - return -x-pow(x,2.0)/2.0-pow(x,3.0)/3.0-pow(x,4.0)/4.0; + T x2 = x*x; + T x3 = x2*x; + T x4 = x3*x; + return -x-x2/2.0-x3/3.0-x4/4.0; }; template @@ -714,13 +744,7 @@ namespace likelihood{ struct thorstenLikelihood { template T operator()(double k, const std::vector& raw_w, int n_events) const { -#ifdef __APPLE__ - //_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW); - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW); -#elif __linux__ - //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); -#endif + ENABLE_FLOAT_EXCEPTIONS(); double log_percentage = log(0.99); @@ -836,22 +860,14 @@ namespace likelihood{ return accumulate(terms.begin(), terms.end()); -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)); -#elif __linux__ - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + DISABLE_FLOAT_EXCEPTIONS(); } }; /* struct thorstenLikelihood { template T operator()(double k, const std::vector& raw_w, int n_events) const { -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW); -#elif __linux__ - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + ENABLE_FLOAT_EXCEPTIONS(); unsigned int count = 1; std::vector kmc; @@ -927,20 +943,11 @@ namespace likelihood{ L += lf; //std::cout << "L += " << lf << std::endl; //std::cout << "L = " << L << std::endl; - -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)); -#elif __linux__ - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + DISABLE_FLOAT_EXCEPTIONS(); return L; } else { -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)); -#elif __linux__ - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + DISABLE_FLOAT_EXCEPTIONS(); return T(0); } } @@ -989,7 +996,7 @@ namespace likelihood{ } }; - struct bayesianSAYLikelihood { + struct SAYLikelihood { template T operator()(double k, T const & w_sum, T const & w2_sum) const { if(w_sum <= 0 || w2_sum < 0) { @@ -1011,7 +1018,7 @@ namespace likelihood{ } } - T alpha = pow(w_sum, T(2))/w2_sum; + T alpha = w_sum*w_sum/w2_sum + one; T beta = w_sum/w2_sum; T L = gammaPriorPoissonLikelihood()(k, alpha, beta); @@ -1019,15 +1026,26 @@ namespace likelihood{ } }; - struct frequentistSAYLikelihood { + struct SAYLikelihoodConstUncertaintyMod { + double sigma; + double sigma2; + SAYLikelihoodConstUncertaintyMod(double sigma):sigma(sigma),sigma2(sigma*sigma){} + void SetSigma(double sigma) { + this->sigma = sigma; + this->sigma2 = sigma*sigma; + } + double GetSigma() { + return this->sigma; + } template T operator()(double k, T const & w_sum, T const & w2_sum) const { - if(w_sum <= 0 || w2_sum < 0) { + T new_sigma2 = w2_sum + this->sigma2; + if(w_sum <= 0 || new_sigma2 < 0) { return(k==0?0:-std::numeric_limits::max()); } - if(w2_sum == 0) { - return poissonLikelihood()(k, w_sum, w2_sum); + if(new_sigma2 == 0) { + return poissonLikelihood()(k, w_sum, new_sigma2); } T one(1); @@ -1041,30 +1059,50 @@ namespace likelihood{ } } - const T & mu = w_sum; - T mu2 = pow(mu, T(2)); - const T & sigma2 = w2_sum; - - T beta = (mu + sqrt(mu2+sigma2*4.0))/(sigma2*2); - T alpha = (mu*sqrt(mu2+sigma2*4.0)/sigma2 + mu2/sigma2 + 2.0) / 2.0; + T alpha = w_sum*w_sum/new_sigma2 + one; + T beta = w_sum/new_sigma2; T L = gammaPriorPoissonLikelihood()(k, alpha, beta); return L; } }; - struct SAYLikelihood { - const bool is_frequentist; - SAYLikelihood():is_frequentist(false) {} - SAYLikelihood(bool is_frequentist):is_frequentist(is_frequentist) {} - bayesianSAYLikelihood bayesian; - frequentistSAYLikelihood frequentist; + struct SAYLikelihoodRelativeUncertaintyMod { + double sigma_over_mu; + SAYLikelihoodRelativeUncertaintyMod(double sigma_over_mu):sigma_over_mu(sigma_over_mu){} + void SetSigmaOverMu(double sigma_over_mu) { + this->sigma_over_mu = sigma_over_mu; + } + double GetSigmaOverMu() { + return this->sigma_over_mu; + } template T operator()(double k, T const & w_sum, T const & w2_sum) const { - if(is_frequentist) - return frequentist(k, w_sum, w2_sum); - else - return bayesian(k, w_sum, w2_sum); + T new_sigma2 = w2_sum + pow(this->sigma_over_mu*w_sum, 2); + if(w_sum <= 0 || new_sigma2 < 0) { + return(k==0?0:-std::numeric_limits::max()); + } + + if(new_sigma2 == 0) { + return poissonLikelihood()(k, w_sum, new_sigma2); + } + + T one(1); + T zero(0); + if(w_sum == zero) { + if(k == 0) { + return zero; + } + else { + return T(-std::numeric_limits::infinity()); + } + } + + T alpha = w_sum*w_sum/new_sigma2 + one; + T beta = w_sum/new_sigma2; + T L = gammaPriorPoissonLikelihood()(k, alpha, beta); + + return L; } }; @@ -1405,6 +1443,11 @@ namespace likelihood{ return(parameterSeeds); } + void setSeed(const std::vector& newSeed) { + assert(newSeed.size() == parameterSeeds.size()); + parameterSeeds = newSeed; + } + void setObservation(const HistogramsType& newObs){ observation=newObs; } @@ -1450,7 +1493,7 @@ namespace likelihood{ for(const RawEvent& e : ((entryStoringBin)*expIt)){ w = weighter(e); assert(w >= 0); - w2 = pow(w, DataType(2.0)); + w2 = w * w; assert(w2 >= 0); assert(e.num_events > 0); n_events += e.num_events; @@ -1531,7 +1574,7 @@ namespace likelihood{ for(const RawEvent& e : ((entryStoringBin)*it)) { w = weighter(e); assert(w >= 0); - w2 = pow(w, DataType(2.0)); + w2 = w*w; assert(w2 >= 0); assert(e.num_events > 0); n_events += e.num_events; @@ -1989,7 +2032,7 @@ namespace likelihood{ lastBestDiscreteIndex=dn; } } - //std::cout << " llh: " << bestLLH << std::endl; + std::cout << bestLLH << std::endl; return(bestLLH); } @@ -2172,6 +2215,38 @@ namespace likelihood{ } }; + struct PowerPrior{ + private: + double power, min, max, lnorm; + public: + PowerPrior(double power=0.0, double min=-std::numeric_limits::infinity(), + double max=std::numeric_limits::infinity()): + power(power), min(min), max(max) { + assert(min >= 0.0); + if(min == 0.0 and power <= 0.0) { + lnorm = 0.0; + } + else if(std::isfinite(max) && std::isfinite(min)) { + lnorm = log((std::pow(min, power+1.0) - std::pow(max, power+1.0))/(power+1.0)); + } + else if(std::isfinite(min) && power < 1.0) { + lnorm = log(std::pow(min, power+1.0)/(power+1.0)); + } + else { + lnorm = 0.0; + } + } + + template + DataType operator()(DataType x) const{ + if(xmax) + return(DataType(-std::numeric_limits::infinity())); + else if(power == 0.0) + return 0.0; + return power*log(x) - lnorm; + } + }; + struct GaussianPrior{ private: double mean; @@ -2222,7 +2297,7 @@ namespace likelihood{ Gaussian2DPrior(double mean0, double mean1, double stddev0, double stddev1, double correlation): mean0(mean0),mean1(mean1),stddev0(stddev0),stddev1(stddev1),correlation(correlation), lnorm(log(boost::math::constants::one_div_two_pi()/(stddev0*stddev1*sqrt(1.0-correlation*correlation)))), - prefactor(-1.0/(2.0*sqrt(1.0-correlation*correlation))){ + prefactor(-1.0/(2.0*(1.0-correlation*correlation))){ if(std::isinf(stddev0) || std::isinf(stddev0) || std::isnan(stddev0) || std::isnan(stddev1) || std::isnan(correlation)) { lnorm = 0.0; prefactor = 0.0; @@ -2254,6 +2329,124 @@ namespace likelihood{ } }; + + template + bool InvertMatrix (const ublas::matrix& input, ublas::matrix& inverse) { + /* + Using code snippet from https://stackoverflow.com/questions/17240624/matrix-inversion-in-boost + + Assigns inverse to provided matrix 'inverse' + returns True if successful, False if failed + */ + + typedef ublas::permutation_matrix pmatrix; + ublas::matrix A(input); + // create a permutation matrix for the LU-factorization + pmatrix pm(A.size1()); + + // perform LU-factorization + int res = ublas::lu_factorize(A,pm); + if( res != 0 ) + return false; + + // create identity matrix of "inverse" + inverse.assign(ublas::identity_matrix(A.size1())); + // backsubstitute to get the inverse + ublas::lu_substitute(A, pm, inverse); + return true; + } + + template + struct GaussianNDPrior{ + /* + This is a generalization of the `Gaussian2DPrior` class + */ + private: + // these are matrices rather than vectors, makes it easier! + std::shared_ptr> means = std::make_shared>( ublas::vector(_size)); + std::shared_ptr> stddevs = std::make_shared>( ublas::vector(_size)); + std::shared_ptr> correlations = std::make_shared>( ublas::matrix(_size, _size)); + std::shared_ptr> inverse = std::make_shared>( ublas::matrix(_size, _size)); + + double determinant; + double lnorm; + public: + + GaussianNDPrior(std::vector& _means, std::vector& _stddevs, std::vector>& _correlations){ + /* + We need to calculate the determinant of this 2D matrix, and we need to verify the lengths of these arrays are the right shape! + See: https://en.wikipedia.org/wiki/Covariance_matrix#Covariance_matrix_as_a_parameter_of_a_distribution + + We pre-calcualte the determinant and inverse since they're relevant later. Then we calculate the log-likelihood that a sample of values (x) is drawn from this distribution + */ + + for (unsigned int i = 0 ; i< _size; i++){ + (*stddevs)(i) = _stddevs[i]; + (*means)(i) = _means[i]; + for (unsigned int j=0; j<_size; j++){ + (*correlations)(i,j) = _correlations[i][j]; + } + } + + /* + http://www.richelbilderbeek.nl/CppUblasMatrixExample7.htm + + This calculates the determinant for an arbitrary ublas matrix! Matrix `m` needs to be square + */ + ublas::matrix m = *correlations; + assert(m.size1() == m.size2() && "Can only calculate the determinant of square matrices"); + ublas::permutation_matrix pivots(m.size1() ); + + const int is_singular = ublas::lu_factorize(m, pivots); + + double determinant = 1.0; + if (is_singular) determinant = 0.0; + else { + const std::size_t sz = pivots.size(); + for (std::size_t i=0; i != sz; ++i) + { + if (pivots(i) != i) + { + determinant *= -1.0; + } + determinant *= m(i,i); + } + } + + bool success = InvertMatrix(*correlations, *inverse); + if (!success){ + throw std::invalid_argument("Correlations matrix was un-invertible!"); + } + + lnorm = log(sqrt(pow(boost::math::constants::one_div_two_pi(),static_cast(_size))/determinant)); + } + + + template + DataType operator()(DataType const & arg0, ARGS const & ... args){ + + std::vector _x { arg0, args... }; + + /* + This calculates the log likelihood of measuring _x given this correlated prior + */ + if(_x.size()!= _size){ + throw std::invalid_argument("Vector should be same size as priors"); + } + + ublas::vector divvy(_size); + + for (size_t i=0;i<_size;i++){ + divvy(i) = ((*means)(i) - _x[i])/ (*stddevs)(i); + } + + auto first = ublas::prod(*inverse, divvy); + auto final = ublas::inner_prod(divvy, first); + return lnorm - 0.5*final; + + } + }; + struct ZigZagPrior{ private: double point; diff --git a/PhysTools/residuals.h b/PhysTools/residuals.h index b5ed16f..5c6413e 100644 --- a/PhysTools/residuals.h +++ b/PhysTools/residuals.h @@ -6,11 +6,28 @@ #include #include #include -#ifdef __APPLE__ + +#if defined(__x86_64__) || defined(_M_X64) #include -#elif __linux__ - #include + #define ENABLE_FLOAT_EXCEPTIONS() \ + _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | \ + _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | \ + _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW) + #define DISABLE_FLOAT_EXCEPTIONS() \ + _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( \ + _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | \ + _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)) +#elif defined(__linux__) + #include + #define ENABLE_FLOAT_EXCEPTIONS() \ + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) + #define DISABLE_FLOAT_EXCEPTIONS() \ + fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) +#else + #define ENABLE_FLOAT_EXCEPTIONS() + #define DISABLE_FLOAT_EXCEPTIONS() #endif + #include #include #include @@ -115,11 +132,7 @@ struct almost_equal > { template struct residual_computer { void operator()(std::vectorconst& z, std::vectorconst& n, std::vectorconst& s, std::vectorconst& m, std::vector& res) { -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW); -#elif __linux__ - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + ENABLE_FLOAT_EXCEPTIONS(); //std::cout << std::setprecision(16); std::cout << "residual_computer" << std::endl; // Useful things to precompute @@ -316,11 +329,7 @@ struct residual_computer { // Clean up the nasty stuff we defined #undef c -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)); -#elif __linux__ - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + DISABLE_FLOAT_EXCEPTIONS(); //return c; } }; @@ -342,7 +351,7 @@ template struct contour_integral_bignum { using big_type=boost::multiprecision::number >; big_type operator()(std::vectorconst& z, std::vectorconst& n, std::vectorconst& s, std::vectorconst& m) { - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + ENABLE_FLOAT_EXCEPTIONS(); const unsigned int M = m.size(); const unsigned int max_m = *std::max_element(m.begin(), m.end()); std::vector big_z(z.begin(), z.end()); @@ -361,7 +370,7 @@ struct contour_integral_bignum { // residual_sum += c(k,0); //} //#undef c - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + DISABLE_FLOAT_EXCEPTIONS(); return residual_sum; } }; @@ -383,7 +392,7 @@ struct contour_integral_bignum > { using result_type=phys_tools::autodiff::FD; using big_type=phys_tools::autodiff::FD > >; big_type operator()(std::vectorconst& z, std::vectorconst& n, std::vectorconst& s, std::vectorconst& m) { - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + ENABLE_FLOAT_EXCEPTIONS(); const unsigned int M = m.size(); const unsigned int max_m = *std::max_element(m.begin(), m.end()); std::vector big_z(z.begin(), z.end()); @@ -402,7 +411,7 @@ struct contour_integral_bignum > { // residual_sum += c(k,0); //} //#undef c - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + DISABLE_FLOAT_EXCEPTIONS(); return residual_sum; } }; @@ -473,7 +482,7 @@ template struct thorsten_fast { typedef typename make_big_type<16, T>::big_type big_type; big_type operator()(std::vectorconst& w, std::vectorconst& s, unsigned int D) { - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + ENABLE_FLOAT_EXCEPTIONS(); unsigned int M = s.size(); //const unsigned int max_m = *std::max_element(m.begin(), m.end()); @@ -545,7 +554,7 @@ struct thorsten_fast { std::cout << "result:" << result << std::endl; //T result = static_cast(big_result); - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); + DISABLE_FLOAT_EXCEPTIONS(); return result; } }; @@ -554,11 +563,7 @@ struct thorsten_fast { template struct contour_integral { T operator()(std::vectorconst& z, std::vectorconst& n, std::vectorconst& s, std::vectorconst& m) { -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW); -#elif __linux__ - feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + ENABLE_FLOAT_EXCEPTIONS(); const unsigned int M = m.size(); const unsigned int max_m = *std::max_element(m.begin(), m.end()); std::vector c; @@ -571,11 +576,7 @@ struct contour_integral { residual_sum += c(k,0); } #undef c -#ifdef __APPLE__ - _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~( _MM_MASK_DIV_ZERO | _MM_MASK_INVALID | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW)); -#elif __linux__ - fedisableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW); -#endif + DISABLE_FLOAT_EXCEPTIONS(); return residual_sum; } }; diff --git a/README.md b/README.md new file mode 100644 index 0000000..d93f0d8 --- /dev/null +++ b/README.md @@ -0,0 +1,25 @@ +# PhysTools +Bin-likelihood and histogram tools + +Author: C. Weaver + +Contributors: A. Schneider, C. Argüelles + +# Installation + +``` +./configure --prefix=$INSTALLATION_PATH +``` + +for most users `$INSTALLATION_PATH=/usr/local/` or their preferred installation directory. If the configuration use + +``` +./configure --help +``` + +to see the available flags to specify specific libraries. + + + + + diff --git a/configure b/configure index 4d73fa0..528620d 100755 --- a/configure +++ b/configure @@ -311,6 +311,22 @@ if [ ! -d ./build/ ]; then mkdir build; fi +echo "Generating pkg-config file..." +echo "prefix=$PREFIX" > phystools.pc + +echo ' +libdir=${prefix}/lib +includedir=${prefix}/include + +Name: PhysTools +Description: Likelihood tools. +URL: https://github.com/icecube/PhysTools +Version: 0.0.1 +Requires: gsl >= 1.15 hdf5 >= 1.8 +Libs: -L${libdir}/lib -lPhysTools +Cflags: -I${includedir}/PhysTools +' >> phystools.pc + echo "Generating config.mk..." echo "# Compiler CC=$CC @@ -325,7 +341,7 @@ PREFIX=$PREFIX # General Settings VERSION=${VERSION} -CXXFLAGS= -std=c++11 -O2 -g -fPIC +CXXFLAGS= -std=c++14 -O3 -g -fPIC # Libraries CXXFLAGS+=${BOOST_CFLAGS} ${HDF5_CFLAGS} diff --git a/test/gaussian_nd.output.txt b/test/gaussian_nd.output.txt new file mode 100644 index 0000000..e69de29 diff --git a/test/gaussian_nd.test.cpp b/test/gaussian_nd.test.cpp new file mode 100644 index 0000000..76b6929 --- /dev/null +++ b/test/gaussian_nd.test.cpp @@ -0,0 +1,51 @@ +#include +#include +#include +#include +#include + + +bool double_compare(double x, double y){ + double diff = std::abs(x-y); + double eps=(1e-6); + return (diff> corr{{1, alpha}, + {alpha, 1}}; + std::vector means{0.0, 0.0}; + std::vector sigmas{1.0, 1.0}; + + const size_t n_vals = static_cast(2); + + GaussianNDPrior<2> testprio(means, sigmas, corr); + Gaussian2DPrior otherprio(means[0], means[1], sigmas[0], sigmas[1], alpha); + + std::vector test_values = {0,0.5,1.0,1.5,2.0,2.5}; + + for (unsigned int i=0; i this_test = {test_values[i], test_values[j]}; + auto result = testprio(this_test); + auto other = otherprio(test_values[i], test_values[j]); + bool are_equal = double_compare(result, other); + if (!are_equal){ + throw std::runtime_error("Test failed!"); + } + + } + } + return 0; + +} \ No newline at end of file