Skip to content

Commit

Permalink
Merge pull request #339 from MeasureTransport/dgsharp/gracefulFail
Browse files Browse the repository at this point in the history 
Remove asserts
  • Loading branch information
@dannys4
dannys4 authored Jul 26, 2023
2 parents d666310 + 39b9aca commit 7356892
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 49 deletions.
101 changes: 67 additions & 34 deletions MParT/MonotoneComponent.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,15 +101,32 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
});
}

bool isGradFunctionInputValid(int sensRows, int sensCols, int ptsRows, int ptsCols, int outputRows, int outputCols, int expectedOutputRows) {
bool isSensRowsValid = sensRows==this->outputDim;
bool isInputColsValid = sensCols==ptsCols;
bool isPtsValid = ptsRows==this->inputDim;
bool isOutputRowsValid = outputRows==expectedOutputRows;
bool isOutputColsValid = outputCols==ptsCols;
return isSensRowsValid && isInputColsValid && isPtsValid && isOutputRowsValid && isOutputColsValid;
}

void checkGradFunctionInput(const std::string method, int sensRows, int sensCols, int ptsRows, int ptsCols, int outputRows, int outputCols, int expectedOutputRows) {
bool isInputValid = isGradFunctionInputValid(sensRows, sensCols, ptsRows, ptsCols, outputRows, outputCols, expectedOutputRows);
if(!isInputValid) {
std::stringstream ss;
ss << method << ": Invalid dimensions of input args." <<
"sens: (" << sensRows << "," << sensCols << "), expected: " << this->outputDim << ", " << ptsCols << "), "
<< "pts: (" << ptsRows << "," << ptsCols << "), expected: (" << this->inputDim << "," << ptsCols << "), "
<< "output: (" << outputRows << "," << outputCols << "), expected: (" << expectedOutputRows << "," << ptsCols << ")";
ProcAgnosticError<MemorySpace,std::invalid_argument>::error(ss.str().c_str());
}
}

void GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
StridedMatrix<const double, MemorySpace> const& sens,
StridedMatrix<double, MemorySpace> output) override
{
assert(sens.extent(0)==this->outputDim);
assert(sens.extent(1)==pts.extent(1));
assert(pts.extent(0)==this->inputDim);
assert(output.extent(0)==this->inputDim);
assert(output.extent(1)==pts.extent(1));
checkGradFunctionInput("Gradient", sens.extent(0), sens.extent(1), pts.extent(0), pts.extent(1), output.extent(0), output.extent(1), this->inputDim);

Kokkos::View<double*,MemorySpace> evals("Map output", pts.extent(1));

Expand All @@ -127,11 +144,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
StridedMatrix<const double, MemorySpace> const& sens,
StridedMatrix<double, MemorySpace> output) override
{
assert(sens.extent(0)==this->outputDim);
assert(sens.extent(1)==pts.extent(1));
assert(pts.extent(0)==this->inputDim);
assert(output.extent(0)==this->numCoeffs);
assert(output.extent(1)==pts.extent(1));
checkGradFunctionInput("CoeffGradImpl", sens.extent(0), sens.extent(1), pts.extent(0), pts.extent(1), output.extent(0), output.extent(1), this->numCoeffs);

Kokkos::View<double*,MemorySpace> evals("Map output", pts.extent(1));

Expand Down Expand Up @@ -222,8 +235,12 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
StridedVector<double,MemorySpace> output)
{
const unsigned int numPts = pts.extent(1);

assert(output.extent(0)==numPts);
if(output.extent(0)!=numPts) {
std::stringstream ss;
ss << "EvaluateImpl: output has incorrect number of columns. "
<< "Expected: " << pts.extent(1) << ", got " << output.extent(0);
ProcAgnosticError<MemorySpace,std::invalid_argument>::error(ss.str().c_str());
}

// Ask the expansion how much memory it would like for its one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -541,8 +558,6 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numPts = pts.extent(1);
const unsigned int numTerms = coeffs.extent(0);

assert(coeffs.extent(0)==numTerms);

Kokkos::View<double*,MemorySpace> output("ExpansionOutput", numPts);

// Ask the expansion how much memory it would like for it's one-point cache
Expand Down Expand Up @@ -600,6 +615,25 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
// DiscreteDerivative(pts2, coeffs2, evals, derivs);
// }

bool isJacobianInputValid(int jacRows, int jacCols, int evalRows, int expectJacRows, int expectJacCols, int expectEvalRows) {
bool isJacRowsCorrect = jacRows == expectJacRows;
bool isJacColsCorrect = jacCols == expectJacCols;
bool isEvalRowsCorrect = evalRows == expectEvalRows;
return isJacRowsCorrect && isJacColsCorrect && isEvalRowsCorrect;
}

void checkJacobianInput(const std::string method, int jacRows, int jacCols, int evalRows, int expectJacRows, int expectJacCols, int expectEvalRows) {
bool isInputValid = isJacobianInputValid(jacRows, jacCols, evalRows, expectJacRows, expectJacCols, expectEvalRows);
if(!isInputValid) {
std::stringstream ss;
ss << method << ": Incorrect input arg sizes. "
<< "jacobian: (" << jacRows << "," << jacCols << "), expected: (" << expectJacRows << "," << expectJacCols << "), ";
if(expectEvalRows > 0)
ss << "evaluations: (" << evalRows << "), expected: (" << expectEvalRows << ")";
ProcAgnosticError<MemorySpace,std::invalid_argument>::error(ss.str().c_str());
}
}

/** @brief Returns the gradient of the map with respect to the parameters \f$\mathbf{w}\f$ at multiple points.
@details
Expand All @@ -609,8 +643,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
@param[in] pts A \f$D\times N\f$ matrix containing the points \f$x^{(1)},\ldots,x^{(N)}\f$. Each column is a point.
@param[in] coeffs A vector of coefficients defining the function \f$f(\mathbf{x}; \mathbf{w})\f$.
@param[out] evaluations A vector containing the \f$N\f$ predictions \f$y_d^{(i)}\f$. The vector must be preallocated and have \f$N\f$ components when passed to this function. An assertion will be thrown in this vector is not the correct size.
@param[out] jacobian A matrix containing the \f$M\times N\f$ Jacobian matrix, where \f$M\f$ is the length of the parameter vector \f$\mathbf{w}\f$. This matrix must be sized correctly or an assertion will be thrown.
@param[out] evaluations A vector containing the \f$N\f$ predictions \f$y_d^{(i)}\f$. The vector must be preallocated and have \f$N\f$ components when passed to this function. An error will occur if this vector is not the correct size.
@param[out] jacobian A matrix containing the \f$M\times N\f$ Jacobian matrix, where \f$M\f$ is the length of the parameter vector \f$\mathbf{w}\f$. This matrix must be sized correctly or an error will occur.
@see CoeffGradient
*/
Expand All @@ -623,9 +657,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numPts = pts.extent(1);
const unsigned int numTerms = coeffs.extent(0);

assert(jacobian.extent(1)==numPts);
assert(jacobian.extent(0)==numTerms);
assert(evaluations.extent(0)==numPts);
checkJacobianInput("CoeffJacobian", jacobian.extent(0), jacobian.extent(1), evaluations.extent(0), numTerms, numPts, numPts);

// Ask the expansion how much memory it would like for it's one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -684,8 +716,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
@param[in] pts A \f$D\times N\f$ matrix containing the points \f$x^{(1)},\ldots,x^{(N)}\f$. Each column is a point.
@param[in] coeffs A vector of coefficients defining the function \f$f(\mathbf{x}; \mathbf{w})\f$.
@param[out] evaluations A vector containing the \f$N\f$ predictions \f$y_d^{(i)}\f$. The vector must be preallocated and have \f$N\f$ components when passed to this function. An assertion will be thrown in this vector is not the correct size.
@param[out] jacobian A matrix containing the \f$d\times N\f$ Jacobian matrix, where \f$d\f$ is the length of the input vector \f$\mathbf{x}\f$. This matrix must be sized correctly or an assertion will be thrown.
@param[out] evaluations A vector containing the \f$N\f$ predictions \f$y_d^{(i)}\f$. The vector must be preallocated and have \f$N\f$ components when passed to this function. An error will occur if this vector is not the correct size.
@param[out] jacobian A matrix containing the \f$d\times N\f$ Jacobian matrix, where \f$d\f$ is the length of the input vector \f$\mathbf{x}\f$. This matrix must be sized correctly or an error will occur.
@see CoeffGradient
*/
Expand All @@ -698,9 +730,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numPts = pts.extent(1);
const unsigned int numTerms = coeffs.extent(0);

assert(jacobian.extent(1)==numPts);
assert(jacobian.extent(0)==dim_);
assert(evaluations.extent(0)==numPts);
checkJacobianInput("InputJacobian",jacobian.extent(0), jacobian.extent(1), evaluations.extent(0), dim_, numPts, numPts);

// Ask the expansion how much memory it would like for it's one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -754,6 +784,10 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
Kokkos::parallel_for(policy, functor);
}

void checkMixedJacobianInput(const std::string method, int jacRows, int jacCols, int expectJacRows, int expectJacCols) {
checkJacobianInput(method, jacRows, jacCols, 0, expectJacRows, expectJacCols, 0);
}

template<typename ExecutionSpace=typename MemoryToExecution<MemorySpace>::Space>
void ContinuousMixedJacobian(StridedMatrix<const double, MemorySpace> const& pts,
StridedVector<const double, MemorySpace> const& coeffs,
Expand All @@ -763,8 +797,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numTerms = coeffs.extent(0);
const unsigned int dim = pts.extent(0);

assert(jacobian.extent(1)==numPts);
assert(jacobian.extent(0)==numTerms);
checkMixedJacobianInput("ContinuousMixedJacobian", jacobian.extent(0), jacobian.extent(1), numTerms, numPts);

// Ask the expansion how much memory it would like for it's one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -815,9 +848,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numPts = pts.extent(1);
const unsigned int dim = pts.extent(0);


assert(jacobian.extent(1)==numPts);
assert(jacobian.extent(0)==dim);
checkMixedJacobianInput("ContinuousMixedInputJacobian", jacobian.extent(0), jacobian.extent(1), dim, numPts);

// Ask the expansion how much memory it would like for it's one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -869,8 +900,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
const unsigned int numPts = pts.extent(1);
const unsigned int numTerms = coeffs.extent(0);

assert(jacobian.extent(1)==numPts);
assert(jacobian.extent(0)==numTerms);
checkMixedJacobianInput("DiscreteMixedJacobian", jacobian.extent(0), jacobian.extent(1), numTerms, numPts);

// Ask the expansion how much memory it would like for it's one-point cache
const unsigned int cacheSize = expansion_.CacheSize();
Expand Down Expand Up @@ -1013,7 +1043,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
break;
}
}
assert(i<maxIts);
if(i>maxIts)
ProcAgnosticError<MemorySpace, std::runtime_error>::error("InverseSingleBracket: lower bound iterations exceed maxIts");

// We have a lower bound...
}else{
Expand All @@ -1030,7 +1061,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
break;
}
}
assert(i<10000);
if(i>maxIts)
ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: upper bound calculation exceeds maxIts");
}

assert(ylb<yub);
Expand Down Expand Up @@ -1074,7 +1106,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
break;
};

assert(it<maxIts);
if(it>maxIts)
ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
return 0.5*(xub+xlb);
}

Expand Down
13 changes: 8 additions & 5 deletions MParT/MonotoneIntegrand.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define MPART_MONOTONEINTEGRAND_H

#include "MParT/DerivativeFlags.h"
#include "MParT/Utilities/Miscellaneous.h"

#include <Kokkos_Core.hpp>

Expand Down Expand Up @@ -157,7 +158,7 @@ class MonotoneIntegrand{
gradSeg(i) = scale*gradSeg(i) + _workspace(i);

}else if(derivType_==DerivativeFlags::Input){

Kokkos::View<double*,MemorySpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>> gradSeg(&output[1], dim);
df = expansion_.MixedInputDerivative(cache_, coeffs_, gradSeg);

Expand All @@ -172,17 +173,17 @@ class MonotoneIntegrand{
// First output is always the integrand itself
double gf = PosFuncType::Evaluate(df);
output[0] = xd_*gf;

// Check for infs or nans
if(std::isinf(gf)){
printf("\nERROR: In MonotoneIntegrand, value of g(df(...)) is inf. The value of df(...) is %0.4f, and the value of f(df(...)) is %0.4f.\n\n", df, gf);
}else if(std::isnan(gf)){
printf("\nERROR: In MonotoneIntegrand, A nan was encountered in value of g(df(...)).\n\n");
}else if(failOnNaN && std::isnan(gf)){
ProcAgnosticError<MemorySpace,std::runtime_error>::error("MonotoneIntegrand: nan was encountered in value of g(df(...)). Use MonotoneIntegrand::setFailOnNaN for enabling NaN propagation.");
}

// Compute the derivative with respect to x_d
if((derivType_==DerivativeFlags::Diagonal) || (derivType_==DerivativeFlags::Input)){

unsigned int ind = (derivType_==DerivativeFlags::Diagonal) ? 1 : dim;
// Compute \partial^2_d f
output[ind] = expansion_.DiagonalDerivative(cache_, coeffs_, 2);
Expand All @@ -193,6 +194,7 @@ class MonotoneIntegrand{
}
}

void setFailOnNaN(bool shouldFailOnNaN) {failOnNaN = shouldFailOnNaN;}
private:

const unsigned int dim_;
Expand All @@ -203,6 +205,7 @@ class MonotoneIntegrand{
CoeffsType const& coeffs_;
DerivativeFlags::DerivativeType derivType_;
Kokkos::View<double*,MemorySpace> _workspace;
bool failOnNaN = true;

}; // class MonotoneIntegrand

Expand Down
14 changes: 14 additions & 0 deletions MParT/Utilities/Miscellaneous.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,20 @@ namespace mpart{
std::string GetOption(std::unordered_map<std::string,std::string> const& map,
std::string const& key,
std::string const& defaultValue);

template<typename MemorySpace, typename ErrorType>
struct ProcAgnosticError {
static void error(const char*) {
assert(0);
}
};

template<typename ErrorType>
struct ProcAgnosticError<Kokkos::HostSpace, ErrorType> {
static void error(const char* message) {
throw ErrorType(message);
}
};
}

#endif
28 changes: 18 additions & 10 deletions tests/Test_MonotoneComponent.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,20 +179,20 @@ TEST_CASE( "MonotoneIntegrand2d", "[MonotoneIntegrand2d]") {
Kokkos::View<double*,HostSpace> coeffGrad("Coefficient Gradient", dim);

for(double t : std::vector<double>{0.0, 0.5, -0.5, 1.0}){


integrand(t, fval.data());

// Compute what it should be
double xd = pt(dim-1);
expansion.FillCache2(&cache[0], pt, t*xd, DerivativeFlags::Diagonal);

df = expansion.MixedInputDerivative(&cache[0], coeffs, coeffGrad);
REQUIRE(fval(0) == Approx(std::abs(xd)*exp(df)).epsilon(testTol));

// Check the derivative against finite differences
for(unsigned int wrt=0; wrt<dim-1; ++wrt){
CHECK(fval(wrt+1) == Approx(xd * std::exp(df)*coeffGrad(wrt)).epsilon(1e-4));
CHECK(fval(wrt+1) == Approx(xd * std::exp(df)*coeffGrad(wrt)).epsilon(1e-4));
}
}
}
Expand Down Expand Up @@ -629,8 +629,16 @@ TEST_CASE( "Testing monotone component derivative", "[MonotoneComponentDerivativ
}

}
}

SECTION("GradientImpl") {

Kokkos::View<double**, HostSpace> evals("Evaluations", 1, numPts);

Kokkos::View<double**, HostSpace> sens("Jacobian", dim+1, numPts);
REQUIRE_THROWS_AS(comp.GradientImpl(evalPts, sens, evals), std::invalid_argument);

}
}


TEST_CASE( "Least squares test", "[MonotoneComponentRegression]" ) {
Expand Down Expand Up @@ -691,7 +699,7 @@ TEST_CASE("Testing MonotoneComponent CoeffGrad and LogDeterminantCoeffGrad", "[M
{
//const double testTol = 1e-4;
unsigned int dim = 2;

// Create points evently spaced on [lb,ub]
unsigned int numPts = 20;
//double lb = -0.5;
Expand All @@ -717,11 +725,11 @@ TEST_CASE("Testing MonotoneComponent CoeffGrad and LogDeterminantCoeffGrad", "[M
Kokkos::View<double*, HostSpace> coeffs("Expansion coefficients", mset.Size());
for(unsigned int i=0; i<coeffs.extent(0); ++i)
coeffs(i) = 0.1*std::cos( 0.01*i );

comp.SetCoeffs(coeffs);

SECTION("CoeffGrad"){

Kokkos::View<double**, HostSpace> sens("Sensitivity", 1, numPts);
for(unsigned int i=0; i<numPts; ++i)
sens(0,i) = 0.25*(i+1);
Expand Down Expand Up @@ -760,7 +768,7 @@ TEST_CASE("Testing MonotoneComponent CoeffGrad and LogDeterminantCoeffGrad", "[M
logDets2 = comp.LogDeterminant(evalPts);
for(unsigned int ptInd=0; ptInd<numPts; ++ptInd)
CHECK( grads(i,ptInd) == Approx((logDets2(ptInd)-logDets(ptInd))/fdstep).epsilon(1e-5));

coeffs(i) -= fdstep;
}

Expand Down

0 comments on commit 7356892

Please sign in to comment.