Skip to content
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

Integration Nugget #353

Merged
merged 4 commits into from
Oct 20, 2023
Merged

Integration Nugget #353

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
05146b9
Added nugget to integrand so that map can have a guaranteed minimum s…
Oct 20, 2023
30a91ff
Updated binding code with new nugget option.
Oct 20, 2023
173c913
Fix to matlab bindings.
Oct 20, 2023
772de3d
More fixes to matlab bindings.
Oct 20, 2023
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
1 change: 1 addition & 0 deletions MParT/MapFactory.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ namespace mpart{
options.quadRelTol = 1e-6; // Optional. Default = 1e-6
options.quadMaxSub = 10; // Optional. Default = 30
options.contDeriv = true; // Optional. Default = true
options.nugget = 1e-4; // Optional. Default = 0.0

// Create a triangular map with these options
unsigned int inDim = 4;
Expand Down
9 changes: 7 additions & 2 deletions MParT/MapOptions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,11 +96,14 @@ namespace mpart{
/** If orthogonal polynomial basis functions should be normalized. */
bool basisNorm = true;

/** The minimum slope of the monotone component. This nugget is added to the g(df) integrand. Must be non-negative. */
double nugget = 0.0;

#if defined(MPART_HAS_CEREAL)
template<class Archive>
void serialize(Archive &archive)
{
archive( basisType, basisLB, basisUB, posFuncType, quadType, quadAbsTol, quadRelTol, quadMaxSub, quadMinSub, quadPts, contDeriv, basisNorm );
archive( basisType, basisLB, basisUB, posFuncType, quadType, quadAbsTol, quadRelTol, quadMaxSub, quadMinSub, quadPts, contDeriv, basisNorm, nugget);
}
#endif // MPART_HAS_CEREAL

Expand All @@ -118,6 +121,7 @@ namespace mpart{
ret &= (quadPts == opts2.quadPts);
ret &= (contDeriv == opts2.contDeriv);
ret &= (basisNorm == opts2.basisNorm);
ret &= (nugget == opts2.nugget);
return ret;
}

Expand All @@ -137,7 +141,8 @@ namespace mpart{
ss << "quadMaxSub = " << quadMaxSub << "\n";
ss << "quadMinSub = " << quadMinSub << "\n";
ss << "quadPts = " << quadPts << "\n";
ss << "contDeriv = " << (contDeriv ? "true" : "false");
ss << "contDeriv = " << (contDeriv ? "true" : "false") << "\n";
ss << "nugget = " << nugget << "\n";
return ss.str();
}
};
Expand Down
51 changes: 30 additions & 21 deletions MParT/MonotoneComponent.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,11 @@ namespace mpart{
@brief Defines a function \f$T:R^N\rightarrow R\f$ such that \f$\partial T / \partial x_N >0\f$ is strictly positive.

@details
The function \f$T\f$ is based on another (generally non-monotone) function \f$f : R^N\rightarrow R\f$ and a strictly positve function
\f$g : R\rightarrow R_{>0}\f$. Together, these functions define the monotone component $T$ through
The function \f$T\f$ is based on another (generally non-monotone) function \f$f : R^N\rightarrow R\f$, a strictly positve function
\f$g : R\rightarrow R_{>0}\f$, and a non-negative scalar \f$\epsilon\f$. Together, these functions define the monotone component $T$ through

\f[
T(x_1, x_2, ..., x_D) = f(x_1,x_2,..., x_{D-1}, 0) + \int_0^{x_D} g\left( \partial_D f(x_1,x_2,..., x_{D-1}, t) \right) dt
T(x_1, x_2, ..., x_D) = f(x_1,x_2,..., x_{D-1}, 0) + \int_0^{x_D} g\left( \partial_D f(x_1,x_2,..., x_{D-1}, t) + \epsilon \right) dt
\f]

@tparam ExpansionType A class defining the function \f$f\f$. It must satisfy the cached parameterization concept.
Expand All @@ -55,28 +55,33 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
@param expansion The expansion used to define the function \f$f\f$.
@param quad The quadrature rule used to approximate \f$\int_0^{x_D} g\left( \partial_D f(x_1,x_2,..., x_{D-1}, t) \right) dt\f$
@param useContDeriv A flag to specify whether the analytic derivative of \f$T(x_1, x_2, ..., x_D)\f$ should be used by default, or if the derivative of the discretized integral should be used. If "true", the analytic or "continuous" derivative will be used. If "false", the derivative of the numerically approximated integral will be used.
@param nugget The value of the non negative quantity \f$\epsilon\f$ defining the minimum slope of the transport map.
@verbatim embed:rst
See the :ref:`diag_deriv_section` mathematical background for more details.
@endverbatim
*/
MonotoneComponent(ExpansionType const& expansion,
QuadratureType const& quad,
bool useContDeriv=true) : ConditionalMapBase<MemorySpace>(expansion.InputSize(), 1, expansion.NumCoeffs()),
bool useContDeriv=true,
double nugget=0.0) : ConditionalMapBase<MemorySpace>(expansion.InputSize(), 1, expansion.NumCoeffs()),
expansion_(expansion),
quad_(quad),
dim_(expansion.InputSize()),
useContDeriv_(useContDeriv){};
useContDeriv_(useContDeriv),
nugget_(nugget){};



MonotoneComponent(ExpansionType const& expansion,
QuadratureType const& quad,
bool useContDeriv,
double nugget,
Kokkos::View<const double*, MemorySpace> coeffsIn) : ConditionalMapBase<MemorySpace>(expansion.InputSize(), 1, expansion.NumCoeffs(), coeffsIn),
expansion_(expansion),
quad_(quad),
dim_(expansion.InputSize()),
useContDeriv_(useContDeriv){};
useContDeriv_(useContDeriv),
nugget_(nugget){};

virtual std::shared_ptr<ParameterizedFunctionBase<MemorySpace>> GetBaseFunction() override{return std::make_shared<MultivariateExpansion<typename ExpansionType::BasisType, typename ExpansionType::KokkosSpace>>(1,expansion_);};

Expand Down Expand Up @@ -426,7 +431,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>

// Compute the inverse
Kokkos::View<double*,MemorySpace> workspace(team_member.thread_scratch(1), workspaceSize);
auto eval = SingleEvaluator<decltype(pt),decltype(coeffs)>(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_);
auto eval = SingleEvaluator<decltype(pt),decltype(coeffs)>(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_, nugget_);
output(ptInd) = RootFinding::InverseSingleBracket<MemorySpace>(ys(ptInd), eval, pt(pt.extent(0)-1), xtol, ytol, info);
}
};
Expand Down Expand Up @@ -614,7 +619,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None);

// Create the integrand g( \partial_D f(x_1,...,x_{D-1},t))
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt), decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Diagonal);
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt), decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Diagonal, nugget_);

// Compute \int_0^x g( \partial_D f(x_1,...,x_{D-1},t)) dt
quad_.Integrate(workspace.data(), integrand, 0, 1, both.data());
Expand Down Expand Up @@ -713,7 +718,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None);

// Create the integrand g( \partial_D f(x_1,...,x_{D-1},t))
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt),decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Parameters);
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt),decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Parameters, nugget_);

// Compute \int_0^x g( \partial_D f(x_1,...,x_{D-1},t)) dt as well as the gradient of this term wrt the coefficients of f
quad_.Integrate(workspace.data(), integrand, 0, 1, integral.data());
Expand Down Expand Up @@ -786,7 +791,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
expansion_.FillCache1(cache.data(), pt, DerivativeFlags::Input);

// Create the integrand g( \partial_D f(x_1,...,x_{D-1},t))
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt),decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Input);
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt),decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Input, nugget_);

// Compute \int_0^x g( \partial_D f(x_1,...,x_{D-1},t)) dt as well as the gradient of this term wrt the map input
quad_.Integrate(workspace.data(), integrand, 0, 1, integral.data());
Expand Down Expand Up @@ -957,7 +962,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>

// Create the integrand g( \partial_D f(x_1,...,x_{D-1},t))
Kokkos::View<double*,MemorySpace> integrandWork(team_member.thread_scratch(1), numTerms);
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt), decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Mixed, integrandWork);
MonotoneIntegrand<ExpansionType, PosFuncType, decltype(pt), decltype(coeffs), MemorySpace> integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Mixed, nugget_, integrandWork);

// Compute \int_0^x g( \partial_D f(x_1,...,x_{D-1},t)) dt as well as the gradient of this term wrt the coefficients of f
quad_.Integrate(workspace.data(), integrand, 0, 1, integral.data());
Expand Down Expand Up @@ -991,7 +996,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
double xd,
CoeffsType const& coeffs,
QuadratureType const& quad,
ExpansionType const& expansion)
ExpansionType const& expansion,
double nugget=0.0)
{
double output = 0.0;
// Compute the integral \int_0^1 g( \partial_D f(x_1,...,x_{D-1},t*x_d)) dt
Expand All @@ -1000,8 +1006,8 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
pt,
xd,
coeffs,
DerivativeFlags::None);

DerivativeFlags::None,
nugget);
quad.Integrate(workspace, integrand, 0, 1, &output);

expansion.FillCache2(cache, pt, 0.0, DerivativeFlags::None);
Expand All @@ -1023,7 +1029,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
void save( Archive & ar ) const
{
ar( cereal::base_class<ConditionalMapBase<MemorySpace>>( this ));
ar( expansion_, quad_, useContDeriv_);
ar( expansion_, quad_, useContDeriv_, nugget_);
ar( this->savedCoeffs );
}

Expand All @@ -1033,15 +1039,16 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
ExpansionType expansion;
QuadratureType quad;
bool useContDeriv;
ar(expansion, quad, useContDeriv);
double nugget;
ar(expansion, quad, useContDeriv, nugget);

Kokkos::View<double*, MemorySpace> coeffs;
ar( coeffs );

if(coeffs.size() == expansion.NumCoeffs()){
construct( expansion, quad, useContDeriv, coeffs);
construct( expansion, quad, useContDeriv, nugget, coeffs);
}else{
construct( expansion, quad, useContDeriv);
construct( expansion, quad, useContDeriv, nugget);
}
}

Expand All @@ -1053,6 +1060,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
QuadratureType quad_;
unsigned int dim_;
bool useContDeriv_;
double nugget_;


template<typename PointType, typename CoeffType>
Expand All @@ -1063,11 +1071,12 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
CoeffType coeffs;
QuadratureType quad;
ExpansionType expansion;
double nugget;

SingleEvaluator(double* workspace_, double* cache_, PointType pt_, CoeffType coeffs_, QuadratureType quad_, ExpansionType expansion_):
workspace(workspace_), cache(cache_), pt(pt_), coeffs(coeffs_), quad(quad_), expansion(expansion_) {};
SingleEvaluator(double* workspace_, double* cache_, PointType pt_, CoeffType coeffs_, QuadratureType quad_, ExpansionType expansion_, double nugget_):
workspace(workspace_), cache(cache_), pt(pt_), coeffs(coeffs_), quad(quad_), expansion(expansion_), nugget(nugget_) {};
double operator()(double x) {
return EvaluateSingle(cache, workspace, pt, x, coeffs, quad, expansion);
return EvaluateSingle(cache, workspace, pt, x, coeffs, quad, expansion, nugget);
}
};
};
Expand Down
Loading