Skip to content

Update likelihood file #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
alexwenym merged 3 commits into from
May 19, 2025
Merged

Update likelihood file #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1e87f1c
updated likelihood file
alexwenym Oct 10, 2024
83c4784
Fixed bug in compilation by putting LDFLAGS at the end
lhennig Dec 3, 2024
e8b2087
Merge pull request #1 from alexwenym/compile-in-docker
alexwenym Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
62 changes: 32 additions & 30 deletions PhysTools/likelihood/likelihood.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2040,7 +2040,7 @@ namespace likelihood{
lastBestDiscreteIndex=dn;
}
}
//std::cout << " llh: " << bestLLH << std::endl;
std::cout << bestLLH << std::endl;
return(bestLLH);
}

Expand Down Expand Up @@ -2338,32 +2338,6 @@ namespace likelihood{
};


double CalcDeterminant(ublas::matrix<double> 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<size_t> 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<class T>
bool InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
/*
Expand Down Expand Up @@ -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<double> m = *correlations;
assert(m.size1() == m.size2() && "Can only calculate the determinant of square matrices");
ublas::permutation_matrix<size_t> 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<double>(*correlations, *inverse);
if (!success){
throw std::invalid_argument("Correlations matrix was un-invertible!");
Expand All @@ -2431,8 +2429,12 @@ namespace likelihood{
lnorm = log(sqrt(pow(boost::math::constants::one_div_two_pi<double>(),static_cast<int>(_size))/determinant));
}

template <typename DataType>
DataType operator()(std::vector<DataType> _x){

template <typename DataType, typename ... ARGS>
DataType operator()(DataType const & arg0, ARGS const & ... args){

std::vector<DataType> _x { arg0, args... };

/*
This calculates the log likelihood of measuring _x given this correlated prior
*/
Expand Down