diff --git a/MParT/DerivativeFlags.h b/MParT/DerivativeFlags.h index be96113b..9d7c58da 100644 --- a/MParT/DerivativeFlags.h +++ b/MParT/DerivativeFlags.h @@ -9,7 +9,9 @@ namespace DerivativeFlags{ Parameters, //<- Deriv wrt coeffs Diagonal, //<- first derivative wrt diagonal Diagonal2, //<- second derivative wrt diagonal - Mixed //<- gradient wrt coeffs of first derivative wrt x_d + Mixed, //<- gradient wrt coeffs of first derivative wrt x_d + Input, //<- gradient wrt map input + MixedInput //<- gradient of diagonal wrt map input }; } diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h index 889241d1..476562c2 100644 --- a/MParT/MonotoneComponent.h +++ b/MParT/MonotoneComponent.h @@ -54,13 +54,13 @@ class MonotoneComponent : public ConditionalMapBase MonotoneComponent(ExpansionType const& expansion, QuadratureType const& quad, bool useContDeriv=true) : ConditionalMapBase(expansion.InputSize(), 1, expansion.NumCoeffs()), - _expansion(expansion), - _quad(quad), - _dim(expansion.InputSize()), - _useContDeriv(useContDeriv){}; + expansion_(expansion), + quad_(quad), + dim_(expansion.InputSize()), + useContDeriv_(useContDeriv){}; - virtual std::shared_ptr> GetBaseFunction() override{return std::make_shared>(1,_expansion);}; + virtual std::shared_ptr> GetBaseFunction() override{return std::make_shared>(1,expansion_);}; /** Override the ConditionalMapBase Evaluate function. */ virtual void EvaluateImpl(StridedMatrix const& pts, @@ -83,7 +83,7 @@ class MonotoneComponent : public ConditionalMapBase StridedVector output) override { // First, get the diagonal derivative - if(_useContDeriv){ + if(useContDeriv_){ ContinuousDerivative(pts, this->savedCoeffs, output); }else{ Kokkos::View evals("Evaluations", pts.extent(1)); @@ -101,6 +101,28 @@ class MonotoneComponent : public ConditionalMapBase }); } + virtual void GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix 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)); + + Kokkos::View evals("Map output", pts.extent(1)); + + InputJacobian(pts, this->savedCoeffs, evals, output); + + // Scale each column by the sensitivity + auto policy = Kokkos::RangePolicy::Space>(0,pts.extent(1)); + Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA (unsigned int ptInd) { + for(unsigned int i=0; iinputDim; ++i) + output(i,ptInd) *= sens(0,ptInd); + }); + } + virtual void CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override @@ -130,7 +152,7 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View derivs("Diagonal Derivative", pts.extent(1)); // First, get the diagonal derivative - if(_useContDeriv){ + if(useContDeriv_){ ContinuousMixedJacobian(pts,this->savedCoeffs, output); ContinuousDerivative(pts, this->savedCoeffs, derivs); }else{ @@ -178,9 +200,9 @@ class MonotoneComponent : public ConditionalMapBase assert(output.extent(0)==numPts); // Ask the expansion how much memory it would like for its one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); - _quad.SetDim(1); - const unsigned int workspaceSize = _quad.WorkspaceSize(); + const unsigned int cacheSize = expansion_.CacheSize(); + quad_.SetDim(1); + const unsigned int workspaceSize = quad_.WorkspaceSize(); auto functor = KOKKOS_CLASS_LAMBDA (typename Kokkos::TeamPolicy::member_type team_member) { @@ -195,8 +217,8 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View workspace(team_member.thread_scratch(1), workspaceSize); // Fill in entries in the cache that are independent of x_d. By passing DerivativeFlags::None, we are telling the expansion that no derivatives with wrt x_1,...x_{d-1} will be needed. - _expansion.FillCache1(cache.data(), pt, DerivativeFlags::None); - output(ptInd) = EvaluateSingle(cache.data(), workspace.data(), pt, pt(_dim-1), coeffs, _quad, _expansion); + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None); + output(ptInd) = EvaluateSingle(cache.data(), workspace.data(), pt, pt(dim_-1), coeffs, quad_, expansion_); } }; @@ -309,9 +331,9 @@ class MonotoneComponent : public ConditionalMapBase } // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); - _quad.SetDim(1); - const unsigned int workspaceSize = _quad.WorkspaceSize(); + const unsigned int cacheSize = expansion_.CacheSize(); + quad_.SetDim(1); + const unsigned int workspaceSize = quad_.WorkspaceSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations unsigned int cacheBytes = Kokkos::View::shmem_size(cacheSize+workspaceSize); @@ -330,11 +352,11 @@ class MonotoneComponent : public ConditionalMapBase // Fill in the cache with everything that doesn't depend on x_d Kokkos::View cache(team_member.thread_scratch(1), cacheSize); - _expansion.FillCache1(cache.data(), pt, DerivativeFlags::None); + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None); // Compute the inverse Kokkos::View workspace(team_member.thread_scratch(1), workspaceSize); - output(ptInd) = InverseSingleBracket(workspace.data(), cache.data(), pt, ys(ptInd), coeffs, xtol, ytol, _quad, _expansion); + output(ptInd) = InverseSingleBracket(workspace.data(), cache.data(), pt, ys(ptInd), coeffs, xtol, ytol, quad_, expansion_); } }; @@ -395,7 +417,7 @@ class MonotoneComponent : public ConditionalMapBase const unsigned int dim = pts.extent(0); // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); + const unsigned int cacheSize = expansion_.CacheSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations auto cacheBytes = Kokkos::View::shmem_size(cacheSize); @@ -413,13 +435,13 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View cache(team_member.thread_scratch(1), cacheSize); // Precompute anything that does not depend on x_d. The DerivativeFlags::None arguments specifies that we won't want to derivative wrt to x_i for i Kokkos::View output("ExpansionOutput", numPts); // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); + const unsigned int cacheSize = expansion_.CacheSize(); - _quad.SetDim(2); - const unsigned int workspaceSize = _quad.WorkspaceSize(); + quad_.SetDim(2); + const unsigned int workspaceSize = quad_.WorkspaceSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations auto cacheBytes = Kokkos::View::shmem_size(cacheSize+workspaceSize+2); @@ -520,19 +542,19 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View both(team_member.thread_scratch(1), 2); // Fill in the cache with anything that doesn't depend on x_d - _expansion.FillCache1(cache.data(), pt, DerivativeFlags::None); + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None); // Create the integrand g( \partial_D f(x_1,...,x_{D-1},t)) - MonotoneIntegrand integrand(cache.data(), _expansion, pt, coeffs, DerivativeFlags::Diagonal); + MonotoneIntegrand integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Diagonal); // Compute \int_0^x g( \partial_D f(x_1,...,x_{D-1},t)) dt - _quad.Integrate(workspace.data(), integrand, 0, 1, both.data()); + quad_.Integrate(workspace.data(), integrand, 0, 1, both.data()); evals(ptInd) = both(0); derivs(ptInd) = both(1); // Add f(x_1,x_2,...,x_{d-1},0) to the evaluation output - _expansion.FillCache2(cache.data(), pt, 0.0, DerivativeFlags::None); - evals(ptInd) += _expansion.Evaluate(cache.data(), coeffs); + expansion_.FillCache2(cache.data(), pt, 0.0, DerivativeFlags::None); + evals(ptInd) += expansion_.Evaluate(cache.data(), coeffs); } }; @@ -580,9 +602,9 @@ class MonotoneComponent : public ConditionalMapBase assert(evaluations.extent(0)==numPts); // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); - _quad.SetDim(numTerms+1); - const unsigned int workspaceSize = _quad.WorkspaceSize(); + const unsigned int cacheSize = expansion_.CacheSize(); + quad_.SetDim(numTerms+1); + const unsigned int workspaceSize = quad_.WorkspaceSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations auto cacheBytes = Kokkos::View::shmem_size(cacheSize+workspaceSize+numTerms+1); @@ -602,18 +624,18 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View integral(team_member.thread_scratch(1), numTerms+1); // Fill in the cache with anything that doesn't depend on x_d - _expansion.FillCache1(cache.data(), pt, DerivativeFlags::None); + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None); // Create the integrand g( \partial_D f(x_1,...,x_{D-1},t)) - MonotoneIntegrand integrand(cache.data(), _expansion, pt, coeffs, DerivativeFlags::Parameters); + MonotoneIntegrand integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Parameters); // 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()); + quad_.Integrate(workspace.data(), integrand, 0, 1, integral.data()); evaluations(ptInd) = integral(0); - _expansion.FillCache2(cache.data(), pt, 0.0, DerivativeFlags::None); - evaluations(ptInd) += _expansion.CoeffDerivative(cache.data(), coeffs, jacView); + expansion_.FillCache2(cache.data(), pt, 0.0, DerivativeFlags::None); + evaluations(ptInd) += expansion_.CoeffDerivative(cache.data(), coeffs, jacView); // Add the Integral to the coefficient gradient for(unsigned int termInd=0; termInd Kokkos::parallel_for(policy, functor); } + /** @brief Returns the gradient of the map with respect to the input \f$x_{1:d}\f$ at multiple points. + + @details + Consider \f$N\f$ points \f$\{\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(N)}\}\f$ and let + \f$y_d^{(i)} = T_d(\mathbf{x}^{(i)}; \mathbf{w})\f$. This function computes \f$\nabla_{\mathbf{x}} y_d^{(i)}\f$ + for each output \f$y_d^{(i)}\f$. + + @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. + + @see CoeffGradient + */ + template::Space> + void InputJacobian(StridedMatrix const& pts, + StridedVector const& coeffs, + StridedVector evaluations, + StridedMatrix jacobian) + { + 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); + + // Ask the expansion how much memory it would like for it's one-point cache + const unsigned int cacheSize = expansion_.CacheSize(); + quad_.SetDim(dim_+1); + const unsigned int workspaceSize = quad_.WorkspaceSize(); + + // Create a policy with enough scratch memory to cache the polynomial evaluations + auto cacheBytes = Kokkos::View::shmem_size(cacheSize+workspaceSize+dim_+1); + + auto functor = KOKKOS_CLASS_LAMBDA (typename Kokkos::TeamPolicy::member_type team_member) { + + unsigned int ptInd = team_member.league_rank () * team_member.team_size () + team_member.team_rank (); + + if(ptInd cache(team_member.thread_scratch(1), cacheSize); + Kokkos::View workspace(team_member.thread_scratch(1), workspaceSize); + Kokkos::View integral(team_member.thread_scratch(1), dim_+1); + + // Fill in the cache with anything that doesn't depend on x_d + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::Input); + + // Create the integrand g( \partial_D f(x_1,...,x_{D-1},t)) + MonotoneIntegrand integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Input); + + // 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()); + + evaluations(ptInd) = integral(0); + + expansion_.FillCache2(cache.data(), pt, 0.0, DerivativeFlags::Input); + evaluations(ptInd) += expansion_.InputDerivative(cache.data(), coeffs, jacView); + + // Add the Integral to the coefficient gradient + for(unsigned int d=0; d(numPts, cacheBytes, functor); + Kokkos::parallel_for(policy, functor); + } + template::Space> void ContinuousMixedJacobian(StridedMatrix const& pts, StridedVector const& coeffs, @@ -640,7 +741,7 @@ class MonotoneComponent : public ConditionalMapBase assert(jacobian.extent(0)==numTerms); // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); + const unsigned int cacheSize = expansion_.CacheSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations auto cacheBytes = Kokkos::View::shmem_size(cacheSize); @@ -659,13 +760,13 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View cache(team_member.thread_scratch(1), cacheSize); // Precompute anything that does not depend on x_d. The DerivativeFlags::None arguments specifies that we won't want to derivative wrt to x_i for i assert(jacobian.extent(0)==numTerms); // Ask the expansion how much memory it would like for it's one-point cache - const unsigned int cacheSize = _expansion.CacheSize(); - _quad.SetDim(numTerms+1); - const unsigned int workspaceSize = _quad.WorkspaceSize(); + const unsigned int cacheSize = expansion_.CacheSize(); + quad_.SetDim(numTerms+1); + const unsigned int workspaceSize = quad_.WorkspaceSize(); // Create a policy with enough scratch memory to cache the polynomial evaluations auto cacheBytes = Kokkos::View::shmem_size(cacheSize+workspaceSize+2*numTerms+1); @@ -714,14 +815,14 @@ class MonotoneComponent : public ConditionalMapBase Kokkos::View integral(team_member.thread_scratch(1), numTerms+1); // Fill in the cache with anything that doesn't depend on x_d - _expansion.FillCache1(cache.data(), pt, DerivativeFlags::None); + expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None); // Create the integrand g( \partial_D f(x_1,...,x_{D-1},t)) Kokkos::View integrandWork(team_member.thread_scratch(1), numTerms); - MonotoneIntegrand integrand(cache.data(), _expansion, pt, coeffs, DerivativeFlags::Mixed, integrandWork); + MonotoneIntegrand integrand(cache.data(), expansion_, pt, coeffs, DerivativeFlags::Mixed, 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()); + quad_.Integrate(workspace.data(), integrand, 0, 1, integral.data()); // Add the Integral to the coefficient gradient for(unsigned int termInd=0; termInd with faster convergence. The worst case performance of ITP, in terms of required evaluations of \f$T\f$, is the same as bisection. - @param cache Memory set up by Kokkos for this evaluation. The `_expansion.FillCache1` function must be + @param cache Memory set up by Kokkos for this evaluation. The `expansion_.FillCache1` function must be called before calling this function. @param pt Array of length \f$D\f$ containing the fixed values of \f$x_{1:D-1}\f$ and an initial guess for \f$x_D\f$ in the last component. @param coeffs An array of coefficients for the expansion. @@ -897,10 +998,10 @@ class MonotoneComponent : public ConditionalMapBase return 0.5*(xub+xlb); } private: - ExpansionType _expansion; - QuadratureType _quad; - const unsigned int _dim; - bool _useContDeriv; + ExpansionType expansion_; + QuadratureType quad_; + const unsigned int dim_; + bool useContDeriv_; }; } // namespace mpart diff --git a/MParT/MonotoneIntegrand.h b/MParT/MonotoneIntegrand.h index cddbf949..9699ae32 100644 --- a/MParT/MonotoneIntegrand.h +++ b/MParT/MonotoneIntegrand.h @@ -42,7 +42,7 @@ class MonotoneIntegrand{ /** @param cache A pointer to memory storing evaluations of \phi_{i,p}(x_i) for each i. These terms need - to be evaluated outside this class (e.g., using `_expansion.FillCache1` for \f$i\in\{0,\ldots,D-1\}\f$. + to be evaluated outside this class (e.g., using `expansion_.FillCache1` for \f$i\in\{0,\ldots,D-1\}\f$. @param expansion @param xd @param coeffs @@ -57,44 +57,45 @@ class MonotoneIntegrand{ } KOKKOS_INLINE_FUNCTION MonotoneIntegrand(double* cache, - ExpansionType const& expansion, - PointType const& pt, - double xd, - CoeffsType const& coeffs, - DerivativeFlags::DerivativeType derivType) : _dim(pt.extent(0)), - _cache(cache), - _expansion(expansion), - _pt(pt), - _xd(xd), - _coeffs(coeffs), - _derivType(derivType) + ExpansionType const& expansion, + PointType const& pt, + double xd, + CoeffsType const& coeffs, + DerivativeFlags::DerivativeType derivType) : dim_(pt.extent(0)), + cache_(cache), + expansion_(expansion), + pt_(pt), + xd_(xd), + coeffs_(coeffs), + derivType_(derivType) { assert(derivType!=DerivativeFlags::Mixed); + assert(derivType!=DerivativeFlags::MixedInput); } KOKKOS_INLINE_FUNCTION MonotoneIntegrand(double* cache, - ExpansionType const& expansion, - PointType const& pt, - CoeffsType const& coeffs, - DerivativeFlags::DerivativeType derivType, - Kokkos::View workspace) : MonotoneIntegrand(cache, expansion, pt, pt(pt.extent(0)-1), coeffs, derivType, workspace) + ExpansionType const& expansion, + PointType const& pt, + CoeffsType const& coeffs, + DerivativeFlags::DerivativeType derivType, + Kokkos::View workspace) : MonotoneIntegrand(cache, expansion, pt, pt(pt.extent(0)-1), coeffs, derivType, workspace) { } KOKKOS_INLINE_FUNCTION MonotoneIntegrand(double* cache, - ExpansionType const& expansion, - PointType const& pt, - double xd, - CoeffsType const& coeffs, - DerivativeFlags::DerivativeType derivType, - Kokkos::View workspace) : _dim(pt.extent(0)), - _cache(cache), - _expansion(expansion), - _pt(pt), - _xd(xd), - _coeffs(coeffs), - _derivType(derivType), - _workspace(workspace) + ExpansionType const& expansion, + PointType const& pt, + double xd, + CoeffsType const& coeffs, + DerivativeFlags::DerivativeType derivType, + Kokkos::View workspace) : dim_(pt.extent(0)), + cache_(cache), + expansion_(expansion), + pt_(pt), + xd_(xd), + coeffs_(coeffs), + derivType_(derivType), + _workspace(workspace) { if(derivType==DerivativeFlags::Mixed) assert(workspace.extent(0)>=coeffs.extent(0)); @@ -109,58 +110,69 @@ class MonotoneIntegrand{ */ KOKKOS_INLINE_FUNCTION void operator()(double t, double* output) const { - const unsigned int numTerms = _expansion.NumCoeffs(); + const unsigned int numTerms = expansion_.NumCoeffs(); + const unsigned int dim = pt_.size(); unsigned int numOutputs = 1; - if(_derivType==DerivativeFlags::Diagonal) + if(derivType_==DerivativeFlags::Diagonal) numOutputs++; - if((_derivType==DerivativeFlags::Parameters) || (_derivType==DerivativeFlags::Mixed)) + if((derivType_==DerivativeFlags::Parameters) || (derivType_==DerivativeFlags::Mixed)) numOutputs += numTerms; + if((derivType_==DerivativeFlags::Input) || (derivType_==DerivativeFlags::MixedInput)) + numOutputs += dim; // Finish filling in the cache at the quadrature point (FillCache1 is called outside this class) - if((_derivType==DerivativeFlags::Diagonal)||(_derivType==DerivativeFlags::Mixed)){ - _expansion.FillCache2(_cache, _pt, t*_xd, DerivativeFlags::Diagonal2); + if((derivType_==DerivativeFlags::Diagonal)||(derivType_==DerivativeFlags::Mixed)||(derivType_==DerivativeFlags::Input)){ + expansion_.FillCache2(cache_, pt_, t*xd_, DerivativeFlags::Diagonal2); }else{ - _expansion.FillCache2(_cache, _pt, t*_xd, DerivativeFlags::Diagonal); + expansion_.FillCache2(cache_, pt_, t*xd_, DerivativeFlags::Diagonal); } - // Use the cache to evaluate \partial_d f and, optionally, the gradient of \partial_d f wrt the coefficients. + // Use the cache to evaluate \partial_d f and, optionally, the gradient of \partial_d f wrt the coefficients or input. double df = 0; - if(_derivType==DerivativeFlags::Parameters){ + if(derivType_==DerivativeFlags::Parameters){ Kokkos::View> gradSeg(&output[1], numTerms); - df = _expansion.MixedDerivative(_cache, _coeffs, 1, gradSeg); + df = expansion_.MixedCoeffDerivative(cache_, coeffs_, 1, gradSeg); - double scale = _xd*PosFuncType::Derivative(df); + double scale = xd_*PosFuncType::Derivative(df); for(unsigned int i=0; i> gradSeg(&output[1], numTerms); - df = _expansion.MixedDerivative(_cache, _coeffs, 1, gradSeg); + df = expansion_.MixedCoeffDerivative(cache_, coeffs_, 1, gradSeg); - scale = _xd*t*df2*PosFuncType::SecondDerivative(df) + dgdf; + scale = xd_*t*df2*PosFuncType::SecondDerivative(df) + dgdf; for(unsigned int i=0; i> gradSeg(&output[1], dim); + df = expansion_.MixedInputDerivative(cache_, coeffs_, gradSeg); + + double scale = xd_*PosFuncType::Derivative(df); + for(unsigned int i=0; i _workspace; }; // class MonotoneIntegrand diff --git a/MParT/MultivariateExpansion.h b/MParT/MultivariateExpansion.h index 578ee199..4cc126dd 100644 --- a/MParT/MultivariateExpansion.h +++ b/MParT/MultivariateExpansion.h @@ -93,7 +93,67 @@ namespace mpart{ Kokkos::fence(); } - void CoeffGradImpl(StridedMatrix const& pts, + virtual void GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix output) override + { + using ExecutionSpace = typename MemoryToExecution::Space; + + const unsigned int numPts = pts.extent(1); + const unsigned int inDim = pts.extent(0); + + // Figure out how much memory we'll need in the cache + unsigned int cacheSize = worker.CacheSize(); + + auto functor = KOKKOS_CLASS_LAMBDA (typename Kokkos::TeamPolicy::member_type team_member) { + + unsigned int ptInd = team_member.league_rank () * team_member.team_size () + team_member.team_rank (); + + if(ptInd cache(team_member.thread_scratch(1), cacheSize); + Kokkos::View grad(team_member.thread_scratch(1), inDim); + + // Fill in entries in the cache that are independent of x_d. By passing DerivativeFlags::None, we are telling the expansion that no derivatives with wrt x_1,...x_{d-1} will be needed. + worker.FillCache1(cache.data(), pt, DerivativeFlags::Input); + worker.FillCache2(cache.data(), pt, pt(pt.size()-1), DerivativeFlags::Input); + + unsigned int coeffStartInd = 0; + for(unsigned int i=0; ioutputDim; ++d){ + + // Extract the coefficients for this output dimension + auto coeffs = Kokkos::subview(this->savedCoeffs, std::make_pair(coeffStartInd, coeffStartInd+worker.NumCoeffs())); + + // Evaluate the expansion + worker.InputDerivative(cache.data(), coeffs, grad); + + for(unsigned int i=0; i::shmem_size(cacheSize + inDim); + + // Paralel loop over each point computing T(x_1,...,x_D) for that point + auto policy = GetCachedRangePolicy(numPts, cacheBytes, functor); + Kokkos::parallel_for(policy, functor); + + Kokkos::fence(); + } + + + virtual void CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override { diff --git a/MParT/MultivariateExpansionWorker.h b/MParT/MultivariateExpansionWorker.h index 9d3b82d2..de67b584 100644 --- a/MParT/MultivariateExpansionWorker.h +++ b/MParT/MultivariateExpansionWorker.h @@ -16,14 +16,16 @@ namespace mpart{ template struct MultivariateExpansionMaxDegreeFunctor { - MultivariateExpansionMaxDegreeFunctor(unsigned int dim, Kokkos::View startPos, Kokkos::View maxDegrees) : dim(dim), startPos(startPos), maxDegrees(maxDegrees) {}; + MultivariateExpansionMaxDegreeFunctor(unsigned int dim, Kokkos::View startPos, Kokkos::View maxDegreesIn) : dim(dim), startPos(startPos), maxDegrees(maxDegreesIn) + { + }; KOKKOS_FUNCTION void operator()(const unsigned int i, unsigned int& update, const bool final) const{ if(final) startPos(i) = update; - if(i::Space>(0,dim_+3), MultivariateExpansionMaxDegreeFunctor(dim_,startPos_, maxDegrees_)); + Kokkos::parallel_scan(Kokkos::RangePolicy::Space>(0,2*dim_+2), MultivariateExpansionMaxDegreeFunctor(dim_,startPos_, maxDegrees_)); // Compute the cache size and store in a double on the host Kokkos::View dCacheSize("Temporary cache size",1); @@ -134,11 +136,18 @@ class MultivariateExpansionWorker template KOKKOS_FUNCTION void FillCache1(double* polyCache, PointType const& pt, - DerivativeFlags::DerivativeType) const + DerivativeFlags::DerivativeType derivType) const { + // Fill in first derivative information if needed + if((derivType == DerivativeFlags::Input)||(derivType==DerivativeFlags::MixedInput)){ + for(unsigned int d=0; d + KOKKOS_FUNCTION double InputDerivative(const double* polyCache, CoeffVecType const& coeffs, GradVecType& grad) const + { + const unsigned int numTerms = multiSet_.Size(); + double f = 0.0; + + for(int wrt=-1; wrt=0){ + grad(wrt) = 0; + } + + for(unsigned int termInd=0; termInd + KOKKOS_FUNCTION double MixedInputDerivative(const double* polyCache, CoeffVecType const& coeffs, GradVecType& grad) const + { + const unsigned int numTerms = multiSet_.Size(); + double df = 0.0; + unsigned int posInd; + + for(int wrt=-1; wrt=0){ + grad(wrt) = 0; + } + + for(unsigned int termInd=0; termInd - KOKKOS_FUNCTION double MixedDerivative(const double* cache, CoeffVecType const& coeffs, unsigned int derivOrder, GradVecType& grad) const + KOKKOS_FUNCTION double MixedCoeffDerivative(const double* cache, CoeffVecType const& coeffs, unsigned int derivOrder, GradVecType& grad) const { const unsigned int numTerms = multiSet_.Size(); @@ -288,7 +398,7 @@ class MultivariateExpansionWorker double df=0; - const unsigned int posIndex = dim_+derivOrder-1; + const unsigned int posIndex = 2*dim_+derivOrder-2; // Compute coeff * polyval for each term for(unsigned int termInd=0; termInd output) = 0; + /** Evaluate the gradient of the function with conversion between default view layout and const strided matrix. */ + template + StridedMatrix Gradient(ViewType1 pts, ViewType2 sens){ + StridedMatrix newpts(pts); + StridedMatrix newSens(sens); + return this->Gradient(newpts, newSens); + } + + /** Evaluate the gradient of the function with conversion from Eigen to Kokkos (and possibly copy to/from device.) */ + Eigen::RowMatrixXd Gradient(Eigen::Ref const& pts, + Eigen::Ref const& sens); + + /** @brief Evaluate the gradient of the function at multiple points. + @details For input points \f$x^{(i)}\f$ and sensitivity vectors \f$s^{(i)}\f$, this function computes + \f[ + g^{(i)} = \left[s^{(i)}\right]^T\nabla_x T(x^{(i)}; w), + \f] + where \f$\nabla_x T\f$ is the Jacobian of the function \f$T\f$ with respect to the input \f$x\f$. Note that this function can + be used to evaluate one step of the chain rule (e.g., one backpropagation step). Given a scalar-valued function \f$g\f$, the gradient of + \f$g(T(x))\f$ with respect to \f$x\f$ is given by \f$\left(\nabla g\right)^T \left(\nabla_x T\right)\f$. Passing \f$\nabla g\f$ as the + sensitivity input to this function then allows you to compute this product and thus the gradient of the composed function \f$g(T(x))\f$. + + @param pts A \f$d_{in}\times N\f$ matrix containining \f$N\f$ points in \f$\mathbb{R}^{d_{in}}\f$ where this function be evaluated. Each column is a point. + @param sens A \f$d_{out}\times N\f$ matrix containing \f$N\f$ sensitivity vectors in \f$\mathbb{R}^{d_{out}}\f$. + @return A \f$d_{in}\times N\f$ matrix containing the gradient of vectors. Each column corresponds to the gradient at a particular point and sensitivity. + */ + template + StridedMatrix Gradient(StridedMatrix const& pts, + StridedMatrix const& sens); + + + virtual void GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix output) = 0; + + /** @brief Computes the gradient of the map output with respect to the map coefficients. @details Consider a map \f$T(x; w) : \mathbb{R}^N \rightarrow \mathbb{R}^M\f$ parameterized by coefficients \f$w\in\mathbb{R}^K\f$. This function computes diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index efc8b35f..97751424 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -78,7 +78,10 @@ class TriangularMap : public ConditionalMapBase void EvaluateImpl(StridedMatrix const& pts, StridedMatrix output) override; - + virtual void GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix output) override; + /** @brief Evaluates the map inverse. @details To understand this function, consider splitting the map input \f$x_{1:N}\f$ into two parts so that \f$x_{1:N} = [x_{1:N-M},x_{N-M+1:M}]\f$. Note that the diff --git a/bindings/julia/src/ParameterizedFunctionBase.cpp b/bindings/julia/src/ParameterizedFunctionBase.cpp index 4c58ae5b..02a84448 100644 --- a/bindings/julia/src/ParameterizedFunctionBase.cpp +++ b/bindings/julia/src/ParameterizedFunctionBase.cpp @@ -25,5 +25,12 @@ void mpart::binding::ParameterizedFunctionBaseWrapper(jlcxx::Module &mod) { pfb.CoeffGradImpl(JuliaToKokkos(pts), JuliaToKokkos(sens), JuliaToKokkos(output)); return output; }) + .method("Gradient", [](ParameterizedFunctionBase &pfb, jlcxx::ArrayRef pts, jlcxx::ArrayRef sens) { + unsigned int numPts = size(pts,1); + unsigned int dim = size(pts,0); + jlcxx::ArrayRef output = jlMalloc(dim, numPts); + pfb.GradientImpl(JuliaToKokkos(pts), JuliaToKokkos(sens), JuliaToKokkos(output)); + return output; + }) ; } \ No newline at end of file diff --git a/bindings/matlab/mat/ConditionalMap.m b/bindings/matlab/mat/ConditionalMap.m index 328bbd38..069af50d 100644 --- a/bindings/matlab/mat/ConditionalMap.m +++ b/bindings/matlab/mat/ConditionalMap.m @@ -126,6 +126,11 @@ function SetCoeffs(this,coeffs) MParT_('ConditionalMap_CoeffGrad',this.id_,pts,sens,result); end + function result = Gradient(this,pts,sens) + result = zeros(size(pts,1), size(pts,2)); + MParT_('ConditionalMap_Gradient',this.id_,pts,sens,result); + end + function result = LogDeterminantCoeffGrad(this,pts) result = zeros(this.numCoeffs, size(pts,2)); MParT_('ConditionalMap_LogDeterminantCoeffGrad',this.id_,pts,result); diff --git a/bindings/matlab/mat/ParameterizedFunction.m b/bindings/matlab/mat/ParameterizedFunction.m index 8ca53d88..c51c6a89 100644 --- a/bindings/matlab/mat/ParameterizedFunction.m +++ b/bindings/matlab/mat/ParameterizedFunction.m @@ -76,6 +76,11 @@ function SetCoeffs(this,coeffs) MParT_('ParameterizedFunction_CoeffGrad',this.id_,pts,sens,result); end + function result = Gradient(this,pts,sens) + result = zeros(size(pts,1), size(pts,2)); + MParT_('ParameterizedFunction_Gradient',this.id_,pts,sens,result); + end + function result = get_id(this) result = this.id_; end diff --git a/bindings/matlab/src/ConditionalMap_mex.cpp b/bindings/matlab/src/ConditionalMap_mex.cpp index 291b47f4..98db0627 100644 --- a/bindings/matlab/src/ConditionalMap_mex.cpp +++ b/bindings/matlab/src/ConditionalMap_mex.cpp @@ -268,6 +268,21 @@ MEX_DEFINE(ConditionalMap_CoeffGrad) (int nlhs, mxArray* plhs[], condMap.map_ptr->CoeffGradImpl(pts,sens,out); } +MEX_DEFINE(ConditionalMap_Gradient) (int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) { + + InputArguments input(nrhs, prhs, 4); + OutputArguments output(nlhs, plhs, 0); + + const ConditionalMapMex& condMap = Session::getConst(input.get(0)); + + auto pts = MexToKokkos2d(prhs[1]); + auto sens = MexToKokkos2d(prhs[2]); + auto out = MexToKokkos2d(prhs[3]); + + condMap.map_ptr->GradientImpl(pts,sens,out); +} + MEX_DEFINE(ConditionalMap_LogDeterminantCoeffGrad) (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { InputArguments input(nrhs, prhs, 3); diff --git a/bindings/matlab/src/ParameterizedFunctionBase_mex.cpp b/bindings/matlab/src/ParameterizedFunctionBase_mex.cpp index ce422d95..0f28e0c5 100644 --- a/bindings/matlab/src/ParameterizedFunctionBase_mex.cpp +++ b/bindings/matlab/src/ParameterizedFunctionBase_mex.cpp @@ -138,4 +138,19 @@ MEX_DEFINE(ParameterizedFunction_CoeffGrad) (int nlhs, mxArray* plhs[], parFunc.fun_ptr->CoeffGradImpl(pts,sens,out); } +MEX_DEFINE(ParameterizedFunction_Gradient) (int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) { + + InputArguments input(nrhs, prhs, 4); + OutputArguments output(nlhs, plhs, 0); + + const ParameterizedFunctionMex& parFunc = Session::getConst(input.get(0)); + + auto pts = MexToKokkos2d(prhs[1]); + auto sens = MexToKokkos2d(prhs[2]); + auto out = MexToKokkos2d(prhs[3]); + + parFunc.fun_ptr->GradientImpl(pts,sens,out); +} + } // namespace \ No newline at end of file diff --git a/bindings/python/src/MultiIndex.cpp b/bindings/python/src/MultiIndex.cpp index 1fb554b5..79774875 100644 --- a/bindings/python/src/MultiIndex.cpp +++ b/bindings/python/src/MultiIndex.cpp @@ -20,6 +20,61 @@ using namespace mpart::binding; void mpart::binding::MultiIndexWrapper(py::module &m) { + // MultiIndexSetLimiters + //TotalOrder + py::class_>(m, "TotalOrder") + .def(py::init()) + .def("__call__", &MultiIndexLimiter::TotalOrder::operator()) + ; + + + //Dimension + py::class_>(m, "Dimension") + .def(py::init()) + .def("__call__", &MultiIndexLimiter::Dimension::operator()) + ; + + + //Anisotropic + py::class_>(m, "Anisotropic") + .def(py::init const&, double>()) + .def("__call__", &MultiIndexLimiter::Anisotropic::operator()) + ; + + + //MaxDegree + py::class_>(m, "MaxDegree") + .def(py::init()) + .def("__call__", &MultiIndexLimiter::MaxDegree::operator()) + ; + + + //None + py::class_>(m, "NoneLim") + .def(py::init<>()) + .def("__call__", &MultiIndexLimiter::None::operator()) + ; + + + //And + py::class_>(m, "And") + .def(py::init,std::function>()) + .def("__call__", &MultiIndexLimiter::And::operator()) + ; + + //Or + py::class_>(m, "Or") + .def(py::init,std::function>()) + .def("__call__", &MultiIndexLimiter::Or::operator()) + ; + + //Xor + py::class_>(m, "Xor") + .def(py::init,std::function>()) + .def("__call__", &MultiIndexLimiter::Xor::operator()) + ; + + // MultiIndex py::class_>(m, "MultiIndex") .def(py::init<>()) @@ -56,8 +111,8 @@ void mpart::binding::MultiIndexWrapper(py::module &m) .def("at", &MultiIndexSet::at) .def("Size", &MultiIndexSet::Size) - .def("CreateTotalOrder", &MultiIndexSet::CreateTotalOrder) - .def("CreateTensorProduct", &MultiIndexSet::CreateTensorProduct) + .def_static("CreateTotalOrder", &MultiIndexSet::CreateTotalOrder, py::arg("length"), py::arg("maxOrder"), py::arg("limiter")=MultiIndexLimiter::None()) + .def_static("CreateTensorProduct", &MultiIndexSet::CreateTensorProduct, py::arg("length"), py::arg("maxOrder"), py::arg("limiter")=MultiIndexLimiter::None()) .def("union", &MultiIndexSet::Union) .def("SetLimiter",&MultiIndexSet::SetLimiter) @@ -83,60 +138,7 @@ void mpart::binding::MultiIndexWrapper(py::module &m) ; - // MultiIndexSetLimiters - //TotalOrder - py::class_>(m, "TotalOrder") - .def(py::init()) - .def("__call__", &MultiIndexLimiter::TotalOrder::operator()) - ; - - - //Dimension - py::class_>(m, "Dimension") - .def(py::init()) - .def("__call__", &MultiIndexLimiter::Dimension::operator()) - ; - - - //Anisotropic - py::class_>(m, "Anisotropic") - .def(py::init const&, double>()) - .def("__call__", &MultiIndexLimiter::Anisotropic::operator()) - ; - - - //MaxDegree - py::class_>(m, "MaxDegree") - .def(py::init()) - .def("__call__", &MultiIndexLimiter::MaxDegree::operator()) - ; - - - //None - py::class_>(m, "NoneLim") - .def(py::init<>()) - .def("__call__", &MultiIndexLimiter::None::operator()) - ; - - - //And - py::class_>(m, "And") - .def(py::init,std::function>()) - .def("__call__", &MultiIndexLimiter::And::operator()) - ; - - //Or - py::class_>(m, "Or") - .def(py::init,std::function>()) - .def("__call__", &MultiIndexLimiter::Or::operator()) - ; - - //Xor - py::class_>(m, "Xor") - .def(py::init,std::function>()) - .def("__call__", &MultiIndexLimiter::Xor::operator()) - ; - + //========================================================================================================== //FixedMultiIndexSet diff --git a/bindings/python/src/ParameterizedFunctionBase.cpp b/bindings/python/src/ParameterizedFunctionBase.cpp index 062e7b65..a71b77f6 100644 --- a/bindings/python/src/ParameterizedFunctionBase.cpp +++ b/bindings/python/src/ParameterizedFunctionBase.cpp @@ -17,6 +17,7 @@ void mpart::binding::ParameterizedFunctionBaseWrapper(py::mod .def("CoeffMap", &ParameterizedFunctionBase::CoeffMap) .def("SetCoeffs", py::overload_cast>(&ParameterizedFunctionBase::SetCoeffs)) .def("Evaluate", static_cast::*)(Eigen::Ref const&)>(&ParameterizedFunctionBase::Evaluate)) + .def("Gradient", static_cast::*)(Eigen::Ref const&, Eigen::Ref const&)>(&ParameterizedFunctionBase::Gradient)) .def("CoeffGrad",static_cast::*)(Eigen::Ref const&, Eigen::Ref const&)>(&ParameterizedFunctionBase::CoeffGrad)) .def_readonly("numCoeffs", &ParameterizedFunctionBase::numCoeffs) .def_readonly("inputDim", &ParameterizedFunctionBase::inputDim) @@ -37,6 +38,7 @@ void mpart::binding::ParameterizedFunctionBaseWrapper(py::mo }) .def("SetCoeffs", py::overload_cast>(&ParameterizedFunctionBase::SetCoeffs)) .def("Evaluate", static_cast::*)(Eigen::Ref const&)>(&ParameterizedFunctionBase::Evaluate)) + .def("Gradient", static_cast::*)(Eigen::Ref const&, Eigen::Ref const&)>(&ParameterizedFunctionBase::Gradient)) .def("CoeffGrad",static_cast::*)(Eigen::Ref const&, Eigen::Ref const&)>(&ParameterizedFunctionBase::CoeffGrad)) .def_readonly("numCoeffs", &ParameterizedFunctionBase::numCoeffs) .def_readonly("inputDim", &ParameterizedFunctionBase::inputDim) diff --git a/src/ParameterizedFunctionBase.cpp b/src/ParameterizedFunctionBase.cpp index c80a1c18..6cc09e91 100644 --- a/src/ParameterizedFunctionBase.cpp +++ b/src/ParameterizedFunctionBase.cpp @@ -99,6 +99,91 @@ Eigen::RowMatrixXd ParameterizedFunctionBase::Evaluate(Eigen #endif + + +template<> +template<> +StridedMatrix ParameterizedFunctionBase::Gradient(StridedMatrix const& pts, StridedMatrix const& sens) +{ + CheckCoefficients("Gradient"); + + Kokkos::View output("Gradients", inputDim, pts.extent(1)); + GradientImpl(pts, sens, output); + return output; +} + +template<> +Eigen::RowMatrixXd ParameterizedFunctionBase::Gradient(Eigen::Ref const& pts, Eigen::Ref const& sens) +{ + CheckCoefficients("Gradient"); + + Eigen::RowMatrixXd output(inputDim, pts.cols()); + StridedMatrix ptsView = ConstRowMatToKokkos(pts); + StridedMatrix sensView = ConstRowMatToKokkos(sens); + StridedMatrix outView = MatToKokkos(output); + GradientImpl(ptsView, sensView, outView); + return output; +} + + +#if defined(MPART_ENABLE_GPU) +template<> +template<> +StridedMatrix ParameterizedFunctionBase::Gradient(StridedMatrix const& pts, StridedMatrix const& sens) +{ + CheckCoefficients("Gradient"); + + Kokkos::View output("Map Evaluations", outputDim, pts.extent(1)); + GradientImpl(pts, sens, output); + return output; +} + +template<> +template<> +StridedMatrix ParameterizedFunctionBase::Gradient(StridedMatrix const& pts, StridedMatrix const& sens) +{ + // Copy the points to the device space + StridedMatrix pts_device = ToDevice(pts); + StridedMatrix sens_device = ToDevice(sens); + // Evaluate on the device space + StridedMatrix evals_device = this->Gradient(pts_device, sens_device); + + // Copy back to the host space + return ToHost(evals_device); +} + +template<> +template<> +StridedMatrix ParameterizedFunctionBase::Gradient(StridedMatrix const& pts, StridedMatrix const& sens) +{ + // Copy the points to host + StridedMatrix pts_host = ToHost(pts); + StridedMatrix sens_host = ToHost(sens); + + // Evaluate on the host + StridedMatrix evals_host = this->Gradient(pts_host, sens_host); + + // Copy back to the device + return ToDevice(evals_host); +} + + +template<> +Eigen::RowMatrixXd ParameterizedFunctionBase::Gradient(Eigen::Ref const& pts, Eigen::Ref const& sens) +{ + CheckCoefficients("Evaluate"); + + Eigen::RowMatrixXd output(outputDim, pts.cols()); + StridedMatrix ptsView = ToDevice( ConstRowMatToKokkos(pts)); + StridedMatrix sensView = ToDevice( ConstRowMatToKokkos(sens)); + + return KokkosToMat( ToHost(this->Gradient(ptsView, sensView))); +} + +#endif + + + template void ParameterizedFunctionBase::SetCoeffs(Kokkos::View coeffs){ diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index fa98da64..d3d601b3 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -157,6 +157,45 @@ void TriangularMap::InverseInplace(StridedMatrix +void TriangularMap::GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix output) +{ + // Evaluate the output for each component + StridedMatrix subPts; + StridedMatrix subSens; + + + + Kokkos::RangePolicy::Space> policy(0,pts.extent(1)); + unsigned int dim = pts.extent(0); + Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int& ptInd){ + for(unsigned int d=0; dinputDim)), Kokkos::ALL()); + subSens = Kokkos::subview(sens, std::make_pair(startOutDim,int(startOutDim+comps_.at(i)->outputDim)), Kokkos::ALL()); + + Kokkos::View subOut("Component Jacobian", comps_.at(i)->inputDim, pts.extent(1)); + comps_.at(i)->GradientImpl(subPts, subSens, subOut); + + dim = comps_.at(i)->inputDim; + Kokkos::parallel_for(policy, KOKKOS_LAMBDA(const int& ptInd){ + for(unsigned int d=0; doutputDim; + } +} + template void TriangularMap::CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, diff --git a/tests/Test_ConditionalMapBase.cpp b/tests/Test_ConditionalMapBase.cpp index 3f3efb68..1b3e8b1e 100644 --- a/tests/Test_ConditionalMapBase.cpp +++ b/tests/Test_ConditionalMapBase.cpp @@ -15,6 +15,13 @@ class MyIdentityMap : public ConditionalMapBase{ virtual void EvaluateImpl(StridedMatrix const& pts, StridedMatrix output) override{Kokkos::deep_copy(output,pts);}; + virtual void GradientImpl(StridedMatrix const& pts, + StridedMatrix const& sens, + StridedMatrix output) override + { + assert(false); + } + virtual void LogDeterminantImpl(StridedMatrix const&, StridedVector output) override{ for(unsigned int i=0; i expansion(mset); - // // Make room for the cache + // Make room for the cache std::vector cache(expansion.CacheSize()); REQUIRE(cache.size() == 6); @@ -129,7 +129,7 @@ TEST_CASE( "MonotoneIntegrand2d", "[MonotoneIntegrand2d]") { pt(i) = 0.25; // Set up the first part of the cache - expansion.FillCache1(&cache[0], pt, DerivativeFlags::None); + expansion.FillCache1(&cache[0], pt, DerivativeFlags::MixedInput); SECTION("Integrand Only") { @@ -168,6 +168,35 @@ TEST_CASE( "MonotoneIntegrand2d", "[MonotoneIntegrand2d]") { } } + SECTION("Integrand Input Gradient") { + double fdStep = 1e-5; + + MonotoneIntegrand, Exp, Kokkos::View,Kokkos::View> integrand(&cache[0], expansion, pt, coeffs, DerivativeFlags::Input); + + // Evaluate the expansion + double df, d2f; + Kokkos::View fval("Integrand", dim+1); + Kokkos::View coeffGrad("Coefficient Gradient", dim); + + for(double t : std::vector{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, Exp, Kokkos::View,Kokkos::View> integrand(&cache[0], expansion, pt, coeffs, DerivativeFlags::Parameters); @@ -183,7 +212,7 @@ TEST_CASE( "MonotoneIntegrand2d", "[MonotoneIntegrand2d]") { expansion.FillCache2(&cache[0], pt, t*xd, DerivativeFlags::Diagonal); df = expansion.DiagonalDerivative(&cache[0], coeffs, 1); - expansion.MixedDerivative(&cache[0], coeffs, 1, coeffGrad); + expansion.MixedCoeffDerivative(&cache[0], coeffs, 1, coeffGrad); integrand(t, fval.data()); CHECK(fval(0) == Approx(xd*exp(df)).epsilon(testTol)); @@ -213,8 +242,8 @@ TEST_CASE( "MonotoneIntegrand2d", "[MonotoneIntegrand2d]") { expansion.FillCache2(&cache[0], pt, t*xd, DerivativeFlags::Diagonal2); df = expansion.DiagonalDerivative(&cache[0], coeffs, 1); d2f = expansion.DiagonalDerivative(&cache[0], coeffs, 2); - expansion.MixedDerivative(&cache[0], coeffs, 1, coeffGrad1); - double d2f2 = expansion.MixedDerivative(&cache[0], coeffs, 2, coeffGrad2); + expansion.MixedCoeffDerivative(&cache[0], coeffs, 1, coeffGrad1); + double d2f2 = expansion.MixedCoeffDerivative(&cache[0], coeffs, 2, coeffGrad2); CHECK(d2f == Approx(d2f2).epsilon(1e-15)); @@ -571,6 +600,35 @@ TEST_CASE( "Testing monotone component derivative", "[MonotoneComponentDerivativ } } + + SECTION("Input Jacobian"){ + + const double fdStep = 1e-4; + + Kokkos::View evals("Evaluations", numPts); + Kokkos::View evals2("Evaluations 2", numPts); + + Kokkos::View jac("Jacobian", dim, numPts); + + comp.InputJacobian(evalPts, coeffs, evals, jac); + + Kokkos::View evalPts2("Points2", evalPts.extent(0), evalPts.extent(1)); + Kokkos::deep_copy(evalPts2, evalPts); + + for(unsigned int j=0; j sens("Sensitivity", outDim, numPts); + Kokkos::View grads; + + for(unsigned int d=0; d expansion(mset); unsigned int cacheSize = expansion.CacheSize(); - CHECK(cacheSize == (maxDegree+1)*(dim+2)); + CHECK(cacheSize == (maxDegree+1)*(2*dim+1)); // Allocate some memory for the cache std::vector cache(cacheSize); @@ -90,7 +90,7 @@ TEST_CASE( "Testing multivariate expansion worker", "[MultivariateExpansionWorke CHECK( gradEig.dot(stepDir) == Approx((f2-f)/fdStep).epsilon(1e-4)); // Mixed first derivatives - df2 = expansion.MixedDerivative(&cache[0], coeffs, 1, grad); + df2 = expansion.MixedCoeffDerivative(&cache[0], coeffs, 1, grad); CHECK(df2==Approx(df).epsilon(1e-15)); df2 = expansion.DiagonalDerivative(&cache[0], coeffs2, 1); @@ -98,11 +98,55 @@ TEST_CASE( "Testing multivariate expansion worker", "[MultivariateExpansionWorke // Mixed second derivatives (grad of d2f wrt coeffs) - double d2f2 = expansion.MixedDerivative(&cache[0], coeffs, 2, grad); + double d2f2 = expansion.MixedCoeffDerivative(&cache[0], coeffs, 2, grad); CHECK(d2f2==Approx(d2f).epsilon(1e-15)); d2f2 = expansion.DiagonalDerivative(&cache[0], coeffs2, 2); CHECK( gradEig.dot(stepDir) == Approx((d2f2-d2f)/fdStep).epsilon(1e-4)); + + SECTION("Input Derivatives"){ + // Check input derivatives + expansion.FillCache1(&cache[0], pt, DerivativeFlags::Input); + expansion.FillCache2(&cache[0], pt, pt(dim-1), DerivativeFlags::Input); + + Kokkos::View inGrad("Input Gradient", dim); + double eval = expansion.Evaluate(&cache[0], coeffs); + double eval2 = expansion.InputDerivative(&cache[0], coeffs, inGrad); + CHECK(eval2 == Approx(eval).epsilon(1e-13)); + + for(unsigned int wrt=0; wrt inGrad("Input Gradient", dim); + double df = expansion.DiagonalDerivative(&cache[0], coeffs, 1); + double df2 = expansion.MixedInputDerivative(&cache[0], coeffs, inGrad); + CHECK(df2 == Approx(df).epsilon(1e-13)); + + for(unsigned int wrt=0; wrt sens("Sensitivities", triMap->outputDim, numSamps); + for(unsigned int j=0; joutputDim; ++i){ + sens(i,j) = 1.0 + 0.1*i + j; + } + } + + Kokkos::View evals = triMap->Evaluate(in); + Kokkos::View evals2; + + Kokkos::View inputGrad = triMap->Gradient(in, sens); + + REQUIRE(inputGrad.extent(0)==triMap->inputDim); + REQUIRE(inputGrad.extent(1)==numSamps); + + // Compare with finite differences + double fdstep = 1e-5; + for(unsigned int i=0; iinputDim; ++i){ + for(unsigned int ptInd=0; ptIndEvaluate(in); + + for(unsigned int ptInd=0; ptIndoutputDim; ++j) + fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep; + + CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).epsilon(1e-3)); + } + + for(unsigned int ptInd=0; ptInd sens("Sensitivities", triMap->outputDim, numSamps);