Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
5 changes: 3 additions & 2 deletions inc/PMNS_DensityMatrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -26,14 +26,15 @@ namespace OscProb {
using PMNS_Base::ProbMatrix;
virtual matrixD ProbMatrix(int nflvi, int nflvf);

/// Set the density matrix from an arbitrary state
virtual void SetInitialRho(matrixC rho_in);

protected:
// Resetting and propagating
/// Reset neutrino state to pure flavour flv
virtual void ResetToFlavour(int flv);
/// Set the density matrix from a pure state
virtual void SetPureState(vectorC nu_in);
/// Set the density matrix from an arbitrary state
virtual void SetInitialRho(matrixC rho_in);

/// Propagate neutrino through a single path
virtual void PropagatePath(NuPath p) = 0;
Expand Down
48 changes: 43 additions & 5 deletions inc/exceptions.h
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,10 @@ std::string format_args(const std::string& names, const Args&... args)
// Lambda to iterate and print each argument.
auto printer = [&](const auto& arg_value) {
if (i < arg_names.size()) {
if (i > 0) ss << ", ";
if (i > 0)
ss << ", ";
else
ss << "\n With: ";
ss << arg_names[i] << " = " << arg_value;
}
i++;
Expand All @@ -72,6 +75,18 @@ std::string format_args(const std::string& names, const Args&... args)
return ss.str();
}

//.............................................................................
///
/// @brief Helper function to handle empty arguments
///
template <typename... Args> auto format_args_wrapper(const Args&... args)
{
if constexpr (sizeof...(args) == 0) { return format_args(""); }
else {
return format_args(args...);
}
}

// Macro to get the most descriptive function name available
#if defined(__GNUC__) || defined(__clang__)
#define FUNCTION_NAME __PRETTY_FUNCTION__
Expand All @@ -82,17 +97,40 @@ std::string format_args(const std::string& names, const Args&... args)
#endif

// The core macro for throwing an exception with detailed info
#define THROW_ON_INVALID_ARG(condition, ...) \
#define THROW_ON_FAIL(exception_type, condition, message, ...) \
do { \
if (!(condition)) { \
std::stringstream ss; \
ss << "\n Condition '" << #condition << "' failed." \
<< "\n With: " << format_args(#__VA_ARGS__, __VA_ARGS__) \
ss << message << "\n Condition '" << #condition << "' failed." \
<< format_args_wrapper(#__VA_ARGS__, ##__VA_ARGS__) \
<< "\n In function: " << FUNCTION_NAME \
<< "\n At file: " << __FILE__ << "\n At line: " << __LINE__; \
throw std::invalid_argument(ss.str()); \
throw exception_type(ss.str()); \
} \
} \
while (0)

// Shorthand macro for invalid argument
#define THROW_ON_INVALID_ARG(condition, message, ...) \
THROW_ON_FAIL(std::invalid_argument, condition, message, ##__VA_ARGS__)

// Shorthand macro for logic error
#define THROW_ON_LOGIC_ERR(condition, message, ...) \
THROW_ON_FAIL(std::logic_error, condition, message, ##__VA_ARGS__)

#define TRY_TO_RETURN(function_call) \
([&]() -> decltype(function_call) { \
try { \
return function_call; \
} \
catch (const std::exception& e) { \
std::stringstream ss; \
ss << "An exception was thrown from " << #function_call \
<< "\n In function: " << FUNCTION_NAME \
<< "\n At file: " << __FILE__ << "\n At line: " << __LINE__ \
<< "\n Original error message: " << e.what(); \
throw std::runtime_error(ss.str()); \
} \
}())

#endif // EXCEPTIONS_H
10 changes: 5 additions & 5 deletions src/PMNS_Deco.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ PMNS_Deco::~PMNS_Deco() {}
///
void PMNS_Deco::SetGamma(int j, double val)
{
THROW_ON_INVALID_ARG(j == 2 || j == 3, j);
THROW_ON_INVALID_ARG(j == 2 || j == 3, "Gamma" << j << "1 not valid", j);

if (val < 0) {
cerr << "WARNING: Gamma_" << j << 1 << " must be positive. "
Expand Down Expand Up @@ -121,7 +121,8 @@ void PMNS_Deco::SetGamma32(double gamma32)

// Sanity check
double get_gamma32 = GetGamma(3, 2);
THROW_ON_INVALID_ARG(fabs(gamma32 - get_gamma32) < 1e-6, gamma32, gamma32);
THROW_ON_LOGIC_ERR(fabs(gamma32 - get_gamma32) < 1e-6, "Failed sanity check",
gamma32, get_gamma32);
}

//.............................................................................
Expand Down Expand Up @@ -182,9 +183,8 @@ double PMNS_Deco::GetGamma(int i, int j)
i = j;
j = temp;
}
THROW_ON_INVALID_ARG(i == 2 || i == 3, i);
THROW_ON_INVALID_ARG(j == 1 || j == 2, i);
THROW_ON_INVALID_ARG(i > j, i, j);
THROW_ON_INVALID_ARG(2 <= i && i <= 3 && 1 <= j && j <= 2 && i > j,
"Gamma" << i << j << " not valid.", i, j);

if (j == 1) { return fGamma[i - 1]; }
else {
Expand Down
31 changes: 18 additions & 13 deletions src/PMNS_DensityMatrix.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
// jcoelho\@apc.in2p3.fr
///////////////////////////////////////////////////////////////////////////////

#include <cassert>
//#include <cassert>

#include "MatrixDecomp/zheevh3.h"

Expand Down Expand Up @@ -86,7 +86,7 @@ void PMNS_DensityMatrix::RotateState(bool to_basis, matrixC U)
///
void PMNS_DensityMatrix::ResetToFlavour(int flv)
{
assert(flv >= 0 && flv < fNumNus);
THROW_ON_LOGIC_ERR(0 <= flv && flv < fNumNus, "flavour not valid,", flv);

PMNS_Base::ResetToFlavour(flv);

Expand Down Expand Up @@ -115,7 +115,7 @@ void PMNS_DensityMatrix::ResetToFlavour(int flv)
///
double PMNS_DensityMatrix::P(int flv)
{
assert(flv >= 0 && flv < fNumNus);
THROW_ON_LOGIC_ERR(0 <= flv && flv < fNumNus, "flavour not valid.", flv);
return abs(fRho[flv][flv]);
}

Expand All @@ -127,7 +127,8 @@ double PMNS_DensityMatrix::P(int flv)
///
void PMNS_DensityMatrix::SetPureState(vectorC nu_in)
{
assert(nu_in.size() == fNumNus);
THROW_ON_LOGIC_ERR(nu_in.size() == fNumNus, "Invalid state vector size.",
nu_in.size());

for (int i = 0; i < fNumNus; i++) {
for (int j = 0; j < fNumNus; j++) {
Expand All @@ -149,8 +150,9 @@ void PMNS_DensityMatrix::SetPureState(vectorC nu_in)
///
void PMNS_DensityMatrix::SetInitialRho(matrixC rho_in)
{
assert(rho_in.size() == fNumNus);
assert(rho_in[0].size() == fNumNus);
THROW_ON_LOGIC_ERR(rho_in.size() == fNumNus && rho_in[0].size() == fNumNus,
"Invalid density matrix size", rho_in.size(),
rho_in[0].size());

double eps = 1e-12;

Expand All @@ -160,13 +162,14 @@ void PMNS_DensityMatrix::SetInitialRho(matrixC rho_in)
trace += abs(rho_in[i][i]);
for (int j = i; j < fNumNus; j++) {
// Ensure rho_in is hermitian
assert(abs(rho_in[i][j] - conj(rho_in[j][i])) < eps);
THROW_ON_LOGIC_ERR(abs(rho_in[i][j] - conj(rho_in[j][i])) < eps,
"Matrix not hermitian", rho_in[i][j], rho_in[j][i]);
rho_array[i][j] = rho_in[i][j];
}
}

// Ensure rho_in has trace 1
assert(abs(trace - 1) < eps);
THROW_ON_LOGIC_ERR(abs(trace - 1) < eps, "Total probability not 1", trace);

double fEvalGLoBES[3];
complexD fEvecGLoBES[3][3];
Expand All @@ -175,7 +178,9 @@ void PMNS_DensityMatrix::SetInitialRho(matrixC rho_in)
zheevh3(rho_array, fEvecGLoBES, fEvalGLoBES);

// Ensure rho_in is positive definite
for (int i = 0; i < fNumNus; i++) assert(fEvalGLoBES[i] >= -eps);
for (int i = 0; i < fNumNus; i++)
THROW_ON_LOGIC_ERR(fEvalGLoBES[i] >= -eps, "Matrix not positive definite",
fEvalGLoBES[i], i);

fRho = rho_in;
}
Expand All @@ -196,10 +201,10 @@ void PMNS_DensityMatrix::SetInitialRho(matrixC rho_in)
///
matrixD PMNS_DensityMatrix::ProbMatrix(int nflvi, int nflvf)
{
THROW_ON_INVALID_ARG(nflvi >= 0, nflvi);
THROW_ON_INVALID_ARG(nflvi <= fNumNus, nflvi, fNumNus);
THROW_ON_INVALID_ARG(nflvf >= 0, nflvf);
THROW_ON_INVALID_ARG(nflvf <= fNumNus, nflvf, fNumNus);
THROW_ON_INVALID_ARG(nflvi >= 0, "", nflvi);
THROW_ON_INVALID_ARG(nflvi <= fNumNus, "", nflvi, fNumNus);
THROW_ON_INVALID_ARG(nflvf >= 0, "", nflvf);
THROW_ON_INVALID_ARG(nflvf <= fNumNus, "", nflvf, fNumNus);

// Output probabilities
matrixD probs(nflvi, vectorD(nflvf));
Expand Down
52 changes: 22 additions & 30 deletions src/PMNS_OQS.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -57,10 +57,7 @@ void PMNS_OQS::SetPower(int n) { fPower = n; }
///
void PMNS_OQS::SetDecoElement(int i, double val)
{
THROW_ON_INVALID_ARG(i > 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);

fBuiltDissipator *= (fa[i] == abs(val));
fBuiltDissipator *= (TRY_TO_RETURN(GetDecoElement(i)) == abs(val));
fa[i] = abs(val);
}

Expand All @@ -74,14 +71,8 @@ void PMNS_OQS::SetDecoElement(int i, double val)
///
void PMNS_OQS::SetDecoAngle(int i, int j, double th)
{
THROW_ON_INVALID_ARG(i > 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
THROW_ON_INVALID_ARG(j > 0, i);
THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
THROW_ON_INVALID_ARG(i != j, i, j);

double val = cos(th);
fBuiltDissipator *= (fcos[i][j] == val);
fBuiltDissipator *= (TRY_TO_RETURN(GetDecoAngle(i, j)) == val);
fcos[i][j] = val;
fcos[j][i] = val;
}
Expand All @@ -104,8 +95,10 @@ int PMNS_OQS::GetPower() { return fPower; }
///
double PMNS_OQS::GetDecoElement(int i)
{
THROW_ON_INVALID_ARG(i > 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
THROW_ON_INVALID_ARG(0 < i && i < SU3_DIM,
"a_" << i << " is not a valid element. "
<< "i must be in range [1, 8]",
i, SU3_DIM);

return fa[i];
}
Expand All @@ -121,11 +114,10 @@ double PMNS_OQS::GetDecoElement(int i)
///
double PMNS_OQS::GetDecoAngle(int i, int j)
{
THROW_ON_INVALID_ARG(i > 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
THROW_ON_INVALID_ARG(j > 0, i);
THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
THROW_ON_INVALID_ARG(i != j, i, j);
THROW_ON_INVALID_ARG(0 < i && i < SU3_DIM && 0 < j && j < SU3_DIM && i != j,
"th_" << i << j << " is not a valid angle. "
<< "i,j must be in range [1, 8] with i<j",
i, j, SU3_DIM);

return acos(fcos[i][j]);
}
Expand All @@ -141,10 +133,10 @@ double PMNS_OQS::GetDecoAngle(int i, int j)
///
double PMNS_OQS::GetHGM(int i, int j)
{
THROW_ON_INVALID_ARG(i >= 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
THROW_ON_INVALID_ARG(j >= 0, i);
THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
THROW_ON_INVALID_ARG(0 <= i && i < SU3_DIM && 0 <= j && j < SU3_DIM,
"HGM_" << i << j << " is not a valid element. "
<< "i,j must be in range [0, 8]",
i, j, SU3_DIM);

return fHGM[i][j];
}
Expand All @@ -160,10 +152,10 @@ double PMNS_OQS::GetHGM(int i, int j)
///
double PMNS_OQS::GetDissipatorElement(int i, int j)
{
THROW_ON_INVALID_ARG(i >= 0, i);
THROW_ON_INVALID_ARG(i < SU3_DIM, i, SU3_DIM);
THROW_ON_INVALID_ARG(j >= 0, i);
THROW_ON_INVALID_ARG(j < SU3_DIM, j, SU3_DIM);
THROW_ON_INVALID_ARG(0 <= i && i < SU3_DIM && 0 <= j && j < SU3_DIM,
"D_" << i << j << " is not a valid element. "
<< "i,j must be in range [0, 8]",
i, j, SU3_DIM);

return fD[i][j];
}
Expand Down Expand Up @@ -513,10 +505,10 @@ void PMNS_OQS::PropagatePath(NuPath p)
///
matrixD PMNS_OQS::ProbMatrix(int nflvi, int nflvf)
{
THROW_ON_INVALID_ARG(nflvi >= 0, nflvi);
THROW_ON_INVALID_ARG(nflvi <= fNumNus, nflvi, fNumNus);
THROW_ON_INVALID_ARG(nflvf >= 0, nflvf);
THROW_ON_INVALID_ARG(nflvf <= fNumNus, nflvf, fNumNus);
THROW_ON_INVALID_ARG(1 <= nflvi <= fNumNus,
"Invalid number of initial flavours", nflvi);
THROW_ON_INVALID_ARG(1 <= nflvi <= fNumNus,
"Invalid number of final flavours", nflvi);

// Output probabilities
matrixD probs(nflvi, vectorD(nflvf));
Expand Down