diff --git a/Makefile b/Makefile index 4abd5c6..f99c5e8 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) diff --git a/PhysTools/likelihood/likelihood.h b/PhysTools/likelihood/likelihood.h index fa64d25..2d0713b 100644 --- a/PhysTools/likelihood/likelihood.h +++ b/PhysTools/likelihood/likelihood.h @@ -2040,7 +2040,7 @@ namespace likelihood{ lastBestDiscreteIndex=dn; } } - //std::cout << " llh: " << bestLLH << std::endl; + std::cout << bestLLH << std::endl; return(bestLLH); } @@ -2338,32 +2338,6 @@ namespace likelihood{ }; - double CalcDeterminant(ublas::matrix m){ - /* - http://www.richelbilderbeek.nl/CppUblasMatrixExample7.htm - - This calculates the determinant for an arbitrary ublas matrix! Matrix `m` needs to be square - */ - 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); - - if (is_singular) return 0.0; - - double d = 1.0; - const std::size_t sz = pivots.size(); - for (std::size_t i=0; i != sz; ++i) - { - if (pivots(i) != i) - { - d *= -1.0; - } - d *= m(i,i); - } - return d; - } - template bool InvertMatrix (const ublas::matrix& input, ublas::matrix& inverse) { /* @@ -2422,7 +2396,31 @@ namespace likelihood{ } } - determinant = CalcDeterminant(*correlations); + /* + 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!"); @@ -2431,8 +2429,12 @@ namespace likelihood{ lnorm = log(sqrt(pow(boost::math::constants::one_div_two_pi(),static_cast(_size))/determinant)); } - template - DataType operator()(std::vector _x){ + + 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 */