From f763802af41ea9dbfa8d64cc1b7bacedb95947a8 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Thu, 29 Sep 2022 20:55:19 -0400 Subject: [PATCH 01/16] Create C++ for Slice of triangular map --- MParT/TriangularMap.h | 32 +++++++++++++++++------- tests/Test_TriangularMap.cpp | 48 +++++++++++++++++++++++++----------- 2 files changed, 56 insertions(+), 24 deletions(-) diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index d2ce80a2..74d7735e 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -29,7 +29,7 @@ is positive definite. */ template class TriangularMap : public ConditionalMapBase{ - + public: /** @brief Construct a block triangular map from a collection of other ConditionalMapBase objects. @@ -55,8 +55,22 @@ class TriangularMap : public ConditionalMapBase{ #if defined(MPART_ENABLE_GPU) virtual void SetCoeffs(Kokkos::View coeffs) override; virtual void WrapCoeffs(Kokkos::View coeffs) override; - #endif - + #endif + + /** @brief Returns a subsection of the map + * @details This function returns a subsection (or "slice") of the map defined by some contiguous subset of the components. + * @param a The index of the first component to take + * @param b The index of the last component to take + * @return A shared pointer to a new TriangularMap object containing the subsection of the map. + */ + std::shared_ptr> Slice(int a, int b) const { + std::vector>> components; + for(int i=a; i<=b; i++){ + components.push_back(this->comps_[i]); + } + return std::make_shared>(components); + } + virtual std::shared_ptr> GetComponent(unsigned int i){ return comps_.at(i);} /** @brief Computes the log determinant of the Jacobian matrix of this map. @@ -80,10 +94,10 @@ class TriangularMap : public ConditionalMapBase{ void EvaluateImpl(StridedMatrix const& pts, StridedMatrix output) override; - virtual void GradientImpl(StridedMatrix const& pts, + 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 @@ -105,15 +119,15 @@ class TriangularMap : public ConditionalMapBase{ StridedMatrix const& r); - virtual void CoeffGradImpl(StridedMatrix const& pts, + virtual void CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override; - virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, + virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) override; - - virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, + + virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) override; private: diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index 7ec02599..1abca323 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -12,7 +12,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ MapOptions options; options.basisType = BasisTypes::ProbabilistHermite; options.basisNorm = false; - + unsigned int numBlocks = 3; unsigned int maxDegree = 2; unsigned int extraInputs = 1; @@ -27,7 +27,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ blocks.at(i) = MapFactory::CreateComponent(mset, options); } - std::shared_ptr> triMap = std::make_shared>(blocks); + std::shared_ptr> triMap = std::make_shared>(blocks); CHECK(triMap->outputDim == numBlocks); CHECK(triMap->inputDim == numBlocks+extraInputs); @@ -37,9 +37,9 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ Kokkos::View coeffs("Coefficients", triMap->numCoeffs); for(unsigned int i=0; inumCoeffs; ++i) coeffs(i) = 0.1*(i+1); - + SECTION("Coefficients"){ - + // Set the coefficients of the triangular map triMap->SetCoeffs(coeffs); @@ -65,11 +65,11 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ triMap->SetCoeffs(coeffs); auto out = triMap->Evaluate(in); - + SECTION("Evaluation"){ for(unsigned int i=0; iEvaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); REQUIRE(outBlock.extent(1)==numSamps); @@ -134,16 +134,16 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ evals2 = triMap->Evaluate(in); for(unsigned int ptInd=0; ptIndoutputDim; ++j) fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep; - CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).epsilon(1e-3)); + CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).epsilon(1e-3)); } coeffs(i) -= fdstep; } - + } @@ -173,18 +173,18 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ evals2 = triMap->Evaluate(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)); + CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).epsilon(1e-3)); } for(unsigned int ptInd=0; ptInd detGrad = triMap->LogDeterminantCoeffGrad(in); REQUIRE(detGrad.extent(0)==triMap->numCoeffs); REQUIRE(detGrad.extent(1)==numSamps); - + Kokkos::View logDet = triMap->LogDeterminant(in); Kokkos::View logDet2; @@ -212,11 +212,29 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ logDet2 = triMap->LogDeterminant(in); for(unsigned int ptInd=0; ptIndSlice(sliceBegin, sliceEnd); + REQUIRE(slice->inputDim == triMap->inputDim); + REQUIRE(slice->outputDim == (sliceEnd-sliceBegin+1)); + // Just test the evaluation + for(unsigned int i=sliceBegin; i<=sliceEnd; ++i){ + + auto outBlock = blocks.at(i)->Evaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); + + REQUIRE(outBlock.extent(1)==numSamps); + REQUIRE(outBlock.extent(0)==1); + for(unsigned int j=0; j Date: Thu, 29 Sep 2022 21:10:25 -0400 Subject: [PATCH 02/16] Add Python & Julia Wrappers --- bindings/julia/src/TriangularMap.cpp | 1 + bindings/python/src/TriangularMap.cpp | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/bindings/julia/src/TriangularMap.cpp b/bindings/julia/src/TriangularMap.cpp index 4ba390c7..d2e4f815 100644 --- a/bindings/julia/src/TriangularMap.cpp +++ b/bindings/julia/src/TriangularMap.cpp @@ -15,5 +15,6 @@ void mpart::binding::TriangularMapWrapper(jlcxx::Module &mod) { map.InverseInplace(JuliaToKokkos(x), JuliaToKokkos(r)); }) .method("GetComponent", &TriangularMap::GetComponent) + .method("Slice", &TriangularMap::Slice) ; } \ No newline at end of file diff --git a/bindings/python/src/TriangularMap.cpp b/bindings/python/src/TriangularMap.cpp index 5dcc5312..1b3c9d22 100644 --- a/bindings/python/src/TriangularMap.cpp +++ b/bindings/python/src/TriangularMap.cpp @@ -20,6 +20,12 @@ void mpart::binding::TriangularMapWrapper(py::module &m) py::class_, ConditionalMapBase, std::shared_ptr>>(m, tName.c_str()) .def(py::init>>>()) .def("GetComponent", &TriangularMap::GetComponent) + .def("Slice", &TriangularMap::Slice) + .def("__getitem__", &TriangularMap::GetComponent) + .def("__getitem__", [](const TriangularMap &m, const py::slice &s) { + if(s.stride() != 1) throw std::runtime_error("TriangularMap slice stride must be 1"); + return m.Slice(s.start(), s.start() + s.size()); + }) ; } From 3408dc0857c9bddfcda99a819770b6e8cc4c280d Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Thu, 29 Sep 2022 21:52:46 -0400 Subject: [PATCH 03/16] Fix slicing --- MParT/MapFactory.h | 11 ++++++----- MParT/TriangularMap.h | 6 +++--- bindings/julia/src/TriangularMap.cpp | 4 ++-- bindings/python/src/TriangularMap.cpp | 6 ++++-- src/MapFactory.cpp | 12 ++++++------ tests/Test_TriangularMap.cpp | 6 +++--- 6 files changed, 24 insertions(+), 21 deletions(-) diff --git a/MParT/MapFactory.h b/MParT/MapFactory.h index 752fd60c..e5355892 100644 --- a/MParT/MapFactory.h +++ b/MParT/MapFactory.h @@ -4,6 +4,7 @@ #include "MParT/MapOptions.h" #include "MParT/ConditionalMapBase.h" +#include "MParT/TriangularMap.h" #include "MParT/MultiIndices/FixedMultiIndexSet.h" #include @@ -34,7 +35,7 @@ namespace mpart{ auto map = MapFactory::CreateTriangular(inDim, outDim, totalOrder, options); @endcode - + */ namespace MapFactory{ @@ -61,7 +62,7 @@ namespace mpart{ @param options Additional options that will be passed on to CreateComponent to construct each MonotoneComponent. */ template - std::shared_ptr> CreateTriangular(unsigned int inputDim, + std::shared_ptr> CreateTriangular(unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder, MapOptions options = MapOptions()); @@ -79,7 +80,7 @@ namespace mpart{ - /** This struct is used to map the options to functions that can create a map component with types corresponding + /** This struct is used to map the options to functions that can create a map component with types corresponding to the options. */ template @@ -98,14 +99,14 @@ namespace mpart{ auto iter = factoryMap->find(optionsKey); if(iter == factoryMap->end()) throw std::runtime_error("Could not find registered factory method for given MapOptions."); - + return iter->second; } static std::shared_ptr GetFactoryMap() { static std::shared_ptr map; - if( !map ) + if( !map ) map = std::make_shared(); return map; } diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index 74d7735e..11c3e748 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -58,14 +58,14 @@ class TriangularMap : public ConditionalMapBase{ #endif /** @brief Returns a subsection of the map - * @details This function returns a subsection (or "slice") of the map defined by some contiguous subset of the components. + * @details This function returns a subsection (or "slice") of the map defined by some contiguous subset of the components, python-style indices. * @param a The index of the first component to take - * @param b The index of the last component to take + * @param b The index AFTER the last component to take * @return A shared pointer to a new TriangularMap object containing the subsection of the map. */ std::shared_ptr> Slice(int a, int b) const { std::vector>> components; - for(int i=a; i<=b; i++){ + for(int i=a; icomps_[i]); } return std::make_shared>(components); diff --git a/bindings/julia/src/TriangularMap.cpp b/bindings/julia/src/TriangularMap.cpp index d2e4f815..8545906a 100644 --- a/bindings/julia/src/TriangularMap.cpp +++ b/bindings/julia/src/TriangularMap.cpp @@ -14,7 +14,7 @@ void mpart::binding::TriangularMapWrapper(jlcxx::Module &mod) { .method("InverseInplace", [](TriangularMap &map, jlcxx::ArrayRef x, jlcxx::ArrayRef r){ map.InverseInplace(JuliaToKokkos(x), JuliaToKokkos(r)); }) - .method("GetComponent", &TriangularMap::GetComponent) - .method("Slice", &TriangularMap::Slice) + .method("GetComponent", [](TriangularMap& map, int i){ return map.GetComponent(i-1); }) + .method("Slice", [](TriangularMap& map, int begin, int end){ return map.Slice(begin-1, end); }) ; } \ No newline at end of file diff --git a/bindings/python/src/TriangularMap.cpp b/bindings/python/src/TriangularMap.cpp index 1b3c9d22..5c9b45e7 100644 --- a/bindings/python/src/TriangularMap.cpp +++ b/bindings/python/src/TriangularMap.cpp @@ -23,8 +23,10 @@ void mpart::binding::TriangularMapWrapper(py::module &m) .def("Slice", &TriangularMap::Slice) .def("__getitem__", &TriangularMap::GetComponent) .def("__getitem__", [](const TriangularMap &m, const py::slice &s) { - if(s.stride() != 1) throw std::runtime_error("TriangularMap slice stride must be 1"); - return m.Slice(s.start(), s.start() + s.size()); + size_t start = 0, stop = 0, step = 0, slicelength = 0; + s.compute(m.outputDim, &start, &stop, &step, &slicelength); + if(step != 1) throw std::runtime_error("TriangularMap slice stride must be 1"); + return m.Slice(start, stop); }) ; diff --git a/src/MapFactory.cpp b/src/MapFactory.cpp index 3f89b482..d35f30fc 100644 --- a/src/MapFactory.cpp +++ b/src/MapFactory.cpp @@ -15,12 +15,12 @@ using namespace mpart; template std::shared_ptr> mpart::MapFactory::CreateComponent(FixedMultiIndexSet const& mset, MapOptions opts) -{ +{ return CompFactoryImpl::GetFactoryFunction(opts)(mset,opts); } template -std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int inputDim, +std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder, MapOptions options) @@ -41,14 +41,14 @@ std::shared_ptr> mpart::MapFactory::CreateTriang template -std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int outputDim, +std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int outputDim, FixedMultiIndexSet const& mset, MapOptions opts) { std::shared_ptr> output; if(opts.basisType==BasisTypes::ProbabilistHermite){ - + if(isinf(opts.basisLB) && isinf(opts.basisUB)){ ProbabilistHermite basis1d(opts.basisNorm); output = std::make_shared>(outputDim, mset, basis1d); @@ -91,9 +91,9 @@ std::shared_ptr> mpart::MapFactory::Creat template std::shared_ptr> mpart::MapFactory::CreateComponent(FixedMultiIndexSet const&, MapOptions); template std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int, FixedMultiIndexSet const&, MapOptions); -template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); +template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); #if defined(MPART_ENABLE_GPU) template std::shared_ptr> mpart::MapFactory::CreateComponent(FixedMultiIndexSet const&, MapOptions); template std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int, FixedMultiIndexSet const&, MapOptions); - template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); + template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); #endif \ No newline at end of file diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index 1abca323..4aeaffb4 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -221,12 +221,12 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ SECTION("Slice"){ int sliceBegin = 1; - int sliceEnd = 2; + int sliceEnd = 3; auto slice = triMap->Slice(sliceBegin, sliceEnd); REQUIRE(slice->inputDim == triMap->inputDim); - REQUIRE(slice->outputDim == (sliceEnd-sliceBegin+1)); + REQUIRE(slice->outputDim == sliceEnd-sliceBegin); // Just test the evaluation - for(unsigned int i=sliceBegin; i<=sliceEnd; ++i){ + for(unsigned int i=sliceBegin; iEvaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); From fc90fc73639bf56947c539b2fac0a01ff1ea464e Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Thu, 29 Sep 2022 23:24:42 -0400 Subject: [PATCH 04/16] Fix wording --- bindings/python/src/TriangularMap.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bindings/python/src/TriangularMap.cpp b/bindings/python/src/TriangularMap.cpp index 5c9b45e7..94daa23e 100644 --- a/bindings/python/src/TriangularMap.cpp +++ b/bindings/python/src/TriangularMap.cpp @@ -25,7 +25,7 @@ void mpart::binding::TriangularMapWrapper(py::module &m) .def("__getitem__", [](const TriangularMap &m, const py::slice &s) { size_t start = 0, stop = 0, step = 0, slicelength = 0; s.compute(m.outputDim, &start, &stop, &step, &slicelength); - if(step != 1) throw std::runtime_error("TriangularMap slice stride must be 1"); + if(step != 1) throw std::runtime_error("TriangularMap slice step must be 1"); return m.Slice(start, stop); }) ; From c0b14f3eecb4fc800dfc122541de8c39b4b8ce19 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Tue, 18 Oct 2022 17:23:57 -0400 Subject: [PATCH 05/16] Work on Slicing --- MParT/TriangularMap.h | 61 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 56 insertions(+), 5 deletions(-) diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index 11c3e748..598b987d 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -57,16 +57,67 @@ class TriangularMap : public ConditionalMapBase{ virtual void WrapCoeffs(Kokkos::View coeffs) override; #endif - /** @brief Returns a subsection of the map + /** @brief Returns a subsection of the map according to index of outputs * @details This function returns a subsection (or "slice") of the map defined by some contiguous subset of the components, python-style indices. - * @param a The index of the first component to take - * @param b The index AFTER the last component to take + * If the \f$j\f$th component has \f$M_j\f$ outputs, and \f$S_n=\sum_{j=1}^nM_j\f$, then the slice \f$(a,b)\f$ will find the indices + * \f$k_a,k_b\f$ such that \f$S_{k_a} \leq a < S_{k_a+1}\f$ and, similarly, \f$S_{k_b}\leq b < S_{k_b+1}\f$. Then, we recursively slice + * the components \f$k_a,k_b\f$ such that, setting \f$T_{k_a}^\prime := T_{k_a}(x)[a-S_{k_a}:S_{k_a+1}]\f$ and \f$T_{k_b}^\prime := T_{k_b}(x)[0:b-S_{k_b}]\f$, + * the returned map is block triangular and defined by T_{a:b}(x) = [T_{k_a}^\prime(x), T_{k_a+1}(x), \ldots, T_{k_b-1}(x), T_{k_b}^\prime(x)]. Further, the + * returned map will have the same number of inputs as the original map and the number of outputs will be \f$b-a\f$. + * @param a The first output index to take + * @param b The index AFTER the last output index to take * @return A shared pointer to a new TriangularMap object containing the subsection of the map. */ std::shared_ptr> Slice(int a, int b) const { std::vector>> components; - for(int i=a; icomps_[i]); + // TODO: Handle empty case + if( a < 0 || a >= b || b > this->GetNumOutputs() ){ + throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= GetNumOutputs() must be satisfied."); + } + int accum_a = 0; + int k_a = 0; + //TODO: Check that this is correct + while(accum_a < a){ + k_a++; + accum_a += this->components[k_a]->GetNumOutputs(); + } + if(a != accum_a){ + accum_a -= this->components[k_a]->GetNumOutputs(); + } + int accum_b = 0; + int k_b = k_a; + //TODO: Check that this is correct + while(accum_b < b){ + k_b++; + accum_b += this->components[k_b]->GetNumOutputs(); + } + if(b != accum_b){ + accum_b -= this->components[k_b]->GetNumOutputs(); + } + if(k_a == k_b){ + components.push_back(this->components[k_a]->Slice(a-accum_a, b-accum_b)); + } else { + components.push_back(this->components[k_a]->Slice(a-accum_a, this->components[k_a]->GetNumOutputs())); + for(int k = k_a+1; k < k_b; k++){ + components.push_back(this->components[k]); + } + components.push_back(this->components[k_b]->Slice(0, b-accum_b)); + } + return std::make_shared>(components); + } + + /** @brief Returns a subsection of the map according to index of components + * @details This function returns a subsection (or "slice") of a block triangular map by simply just taking a subset of the blocks. + * The returned map will have the same number of inputs as the original map, and this is equivalent to {@link #Slice Slice} when + * the number of outputs of each component is 1. + * @param a The first component block to take + * @param b The index AFTER the last component block to take + * @return A shared pointer to a new TriangularMap object containing the subsection of the map. + */ + std::shared_ptr> BlockSlice(int a, int b) const { + std::vector>> components; + for(int k = a; k < b; k++){ + components.push_back(this->components[k]); } return std::make_shared>(components); } From b9d6fcf137c89b2abb066417cdc36ad2f89f3a0e Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Wed, 19 Oct 2022 09:30:17 -0400 Subject: [PATCH 06/16] Creating slice within ConditionalMapBase --- MParT/ConditionalMapBase.h | 35 ++++++++++++++++++----------------- MParT/TriangularMap.h | 13 +------------ 2 files changed, 19 insertions(+), 29 deletions(-) diff --git a/MParT/ConditionalMapBase.h b/MParT/ConditionalMapBase.h index 173d9edb..bcb91df8 100644 --- a/MParT/ConditionalMapBase.h +++ b/MParT/ConditionalMapBase.h @@ -18,7 +18,7 @@ namespace mpart { Let \f$r\in\mathbb{R}^M\f$ denote the map output, \f$r=T(x_2; x_1)\f$. The conditional maps implemented by children of this class guarantee that for fixed \f$x_1\f$, the conditional mapping from \f$x_1 \rightarrow r\f$ is invertible and the Jacobian matrix with respect to \f$x_2\f$, \f$\nabla_{x_2} T\f$, is positive definite. - + @see mpart::MapFactory */ template @@ -28,7 +28,7 @@ namespace mpart { public: /** - @brief Construct a new Conditional Map Base object. + @brief Construct a new Conditional Map Base object. @details Typically the ConditionalMapBase class is not constructed directly. The preferred way of creating a ConditionalMapBase object is through factory methods in the mpart::MapFactory namespace. @param inDim The dimension \f$N\f$ of the input to this map. @param outDim The dimension \f$M\f$ of the output from this map. @@ -37,8 +37,8 @@ namespace mpart { ConditionalMapBase(unsigned int inDim, unsigned int outDim, unsigned int nCoeffs) : ParameterizedFunctionBase(inDim, outDim, nCoeffs){}; virtual ~ConditionalMapBase() = default; - - /** For Monotone parameterizations that are based on a non-monotone base function, this function will return the base function. If the monotone parameterization is + + /** For Monotone parameterizations that are based on a non-monotone base function, this function will return the base function. If the monotone parameterization is not constructed from a non-monotone base, then this function will return a nullptr. */ virtual std::shared_ptr> GetBaseFunction(){return nullptr;}; @@ -73,14 +73,14 @@ namespace mpart { /** @brief Computes the inverse of the map. - @details The ConditionalMapBase class defines invertible functions \f$T:\mathbb{R}^{d_{in}}\rightarrow\mathbb{R}^{d_{out}}\f$ + @details The ConditionalMapBase class defines invertible functions \f$T:\mathbb{R}^{d_{in}}\rightarrow\mathbb{R}^{d_{out}}\f$ where the input dimension \f$d_{in}\geq d_{out}\f$. An input \f$x\in\mathbb{R}^{d_{in}}\f$ can then be split - into two parts \f$x=[x_1,x_2]\f$, where \f$x_1\in\mathbb{R}^{d_{in}-d_{out}}\f$ represents the "extra" inputs + into two parts \f$x=[x_1,x_2]\f$, where \f$x_1\in\mathbb{R}^{d_{in}-d_{out}}\f$ represents the "extra" inputs and \f$x_2\in\mathbb{R}^{d_{out}}\f$ has the same size as the map output. This function computes the second block \f$x_2\f$ of the input given the first block \f$x_1\f$ and a vector \f$r=T(x_1,x_2)\f$. - + @param x1 A \f$d_{in}-d_{out}\times N\f$ or \f$d_{in}\times N\f$ matrix containing \f$N\f$ values of the first input block: - \f$\{x_1^{(1)},\ldots x_1^{(N)}\}\f$. If \f$d_{in}\f$ rows are given, then the last \f$d_{out}\f$ rows may be used as an initial + \f$\{x_1^{(1)},\ldots x_1^{(N)}\}\f$. If \f$d_{in}\f$ rows are given, then the last \f$d_{out}\f$ rows may be used as an initial guess for \f$x_2\f$ by certain solvers. @param r A \f$d_{out}\times N\f$ matrix containing \f$N\f$ values of the map output: \f$\{r^{(1)},\ldots, r^{(N)}\}\f$. @return A \f$d_{out} \times N\f$ matrix containing the computed values of \f$\{x_2^{(1)},\ldots,x_2^{(N)}\}\f$. @@ -92,13 +92,13 @@ namespace mpart { /** Inverse function with conversion between general view type and const strided matrix. */ template StridedMatrix Inverse(ViewType1 x, ViewType2 r){ - StridedMatrix newx(x); - StridedMatrix newr(r); + StridedMatrix newx(x); + StridedMatrix newr(r); return this->Inverse(newx,newr); } /** Inverse function with conversion between eigen matrix and Kokkos view. */ - Eigen::RowMatrixXd Inverse(Eigen::Ref const& x1, + Eigen::RowMatrixXd Inverse(Eigen::Ref const& x1, Eigen::Ref const& r); @@ -109,10 +109,10 @@ namespace mpart { /** @brief Computes the gradient of the log determinant with respect to the map coefficients. @details For 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 + this function computes \f[ \nabla_w \det{\nabla_x T(x_i; w)}, - \f] + \f] at multiple points \f$x_i\f$. @param pts A collection of points where we want to evaluate the gradient. Each column corresponds to a point. @return A matrix containing the coefficient gradient at each input point. The \f$i^{th}\f$ column contains \f$\nabla_w \det{\nabla_x T(x_i; w)}\f$. @@ -123,14 +123,14 @@ namespace mpart { /** Include conversion between general view type and Strided matrix. */ template StridedMatrix LogDeterminantCoeffGrad(ViewType pts){ - StridedMatrix newpts(pts); + StridedMatrix newpts(pts); return this->LogDeterminantCoeffGrad(newpts); } /** Evaluation with additional conversion of Eigen matrix to Kokkos unmanaged view. */ Eigen::RowMatrixXd LogDeterminantCoeffGrad(Eigen::Ref const& pts); - virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, + virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) = 0; @@ -140,16 +140,17 @@ namespace mpart { /** Include conversion between general view type and Strided matrix. */ template StridedMatrix LogDeterminantInputGrad(ViewType pts){ - StridedMatrix newpts(pts); + StridedMatrix newpts(pts); return this->LogDeterminantInputGrad(newpts); } /** Evaluation with additional conversion of Eigen matrix to Kokkos unmanaged view. */ Eigen::RowMatrixXd LogDeterminantInputGrad(Eigen::Ref const& pts); - virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, + virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) = 0; + virtual void Slice(int a, int b) = 0; }; // class ConditionalMapBase } diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index 598b987d..73089896 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -57,18 +57,7 @@ class TriangularMap : public ConditionalMapBase{ virtual void WrapCoeffs(Kokkos::View coeffs) override; #endif - /** @brief Returns a subsection of the map according to index of outputs - * @details This function returns a subsection (or "slice") of the map defined by some contiguous subset of the components, python-style indices. - * If the \f$j\f$th component has \f$M_j\f$ outputs, and \f$S_n=\sum_{j=1}^nM_j\f$, then the slice \f$(a,b)\f$ will find the indices - * \f$k_a,k_b\f$ such that \f$S_{k_a} \leq a < S_{k_a+1}\f$ and, similarly, \f$S_{k_b}\leq b < S_{k_b+1}\f$. Then, we recursively slice - * the components \f$k_a,k_b\f$ such that, setting \f$T_{k_a}^\prime := T_{k_a}(x)[a-S_{k_a}:S_{k_a+1}]\f$ and \f$T_{k_b}^\prime := T_{k_b}(x)[0:b-S_{k_b}]\f$, - * the returned map is block triangular and defined by T_{a:b}(x) = [T_{k_a}^\prime(x), T_{k_a+1}(x), \ldots, T_{k_b-1}(x), T_{k_b}^\prime(x)]. Further, the - * returned map will have the same number of inputs as the original map and the number of outputs will be \f$b-a\f$. - * @param a The first output index to take - * @param b The index AFTER the last output index to take - * @return A shared pointer to a new TriangularMap object containing the subsection of the map. - */ - std::shared_ptr> Slice(int a, int b) const { + virtual std::shared_ptr> Slice(int a, int b) const override { std::vector>> components; // TODO: Handle empty case if( a < 0 || a >= b || b > this->GetNumOutputs() ){ From fd5dd9101e3f99705155c0f17b077579472bc762 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Wed, 19 Oct 2022 11:39:11 -0400 Subject: [PATCH 07/16] Add Slicing for a lot of classes --- MParT/AffineMap.h | 2 + MParT/ComposedMap.h | 2 + MParT/IdentityMap.h | 4 + MParT/MonotoneComponent.h | 5 ++ MParT/TriangularMap.h | 42 +--------- src/AffineMap.cpp | 38 +++++---- src/ComposedMap.cpp | 132 +++++++++++++++++------------- src/IdentityMap.cpp | 18 ++-- src/TriangularMap.cpp | 53 ++++++++++-- tests/Test_ConditionalMapBase.cpp | 71 +++++++++++++--- tests/Test_TriangularMap.cpp | 3 +- 11 files changed, 234 insertions(+), 136 deletions(-) diff --git a/MParT/AffineMap.h b/MParT/AffineMap.h index e69bd131..42fcd10e 100644 --- a/MParT/AffineMap.h +++ b/MParT/AffineMap.h @@ -59,6 +59,8 @@ class AffineMap : public ConditionalMapBase void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) override; + std::shared_ptr> Slice(int a, int b) override; + /** Computes an LU factorization of the matrix A_ */ void Factorize(); diff --git a/MParT/ComposedMap.h b/MParT/ComposedMap.h index 9fb60432..c70e97a6 100644 --- a/MParT/ComposedMap.h +++ b/MParT/ComposedMap.h @@ -100,6 +100,8 @@ class ComposedMap : public ConditionalMapBase void GradientImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override; + + std::shared_ptr> Slice(int a, int b) override; private: unsigned int maxChecks_; diff --git a/MParT/IdentityMap.h b/MParT/IdentityMap.h index 0a9d9fad..bafea491 100644 --- a/MParT/IdentityMap.h +++ b/MParT/IdentityMap.h @@ -57,6 +57,10 @@ class IdentityMap : public ConditionalMapBase void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) override; + // This only slices the tail of the map due to how IdentityMap is defined. + // i.e. this is invariant on identical shifting of a,b. + std::shared_ptr> Slice(int a, int b) override; + }; // class IdentityMap } diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h index c4891072..e30c9c2e 100644 --- a/MParT/MonotoneComponent.h +++ b/MParT/MonotoneComponent.h @@ -1077,6 +1077,11 @@ class MonotoneComponent : public ConditionalMapBase assert(it> Slice(int a, int b) override { + assert(a == 0 && b == 1); + return std::shared_ptr>(this); + } private: ExpansionType expansion_; QuadratureType quad_; diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index 0893db30..8c7a4288 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -57,43 +57,7 @@ class TriangularMap : public ConditionalMapBase{ void WrapCoeffs(Kokkos::View coeffs) override; #endif - virtual std::shared_ptr> Slice(int a, int b) const override { - std::vector>> components; - // TODO: Handle empty case - if( a < 0 || a >= b || b > this->GetNumOutputs() ){ - throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= GetNumOutputs() must be satisfied."); - } - int accum_a = 0; - int k_a = 0; - //TODO: Check that this is correct - while(accum_a < a){ - k_a++; - accum_a += this->components[k_a]->GetNumOutputs(); - } - if(a != accum_a){ - accum_a -= this->components[k_a]->GetNumOutputs(); - } - int accum_b = 0; - int k_b = k_a; - //TODO: Check that this is correct - while(accum_b < b){ - k_b++; - accum_b += this->components[k_b]->GetNumOutputs(); - } - if(b != accum_b){ - accum_b -= this->components[k_b]->GetNumOutputs(); - } - if(k_a == k_b){ - components.push_back(this->components[k_a]->Slice(a-accum_a, b-accum_b)); - } else { - components.push_back(this->components[k_a]->Slice(a-accum_a, this->components[k_a]->GetNumOutputs())); - for(int k = k_a+1; k < k_b; k++){ - components.push_back(this->components[k]); - } - components.push_back(this->components[k_b]->Slice(0, b-accum_b)); - } - return std::make_shared>(components); - } + std::shared_ptr> Slice(int a, int b) override; /** @brief Returns a subsection of the map according to index of components * @details This function returns a subsection (or "slice") of a block triangular map by simply just taking a subset of the blocks. @@ -103,10 +67,10 @@ class TriangularMap : public ConditionalMapBase{ * @param b The index AFTER the last component block to take * @return A shared pointer to a new TriangularMap object containing the subsection of the map. */ - std::shared_ptr> BlockSlice(int a, int b) const { + std::shared_ptr> BlockSlice(int a, int b) { std::vector>> components; for(int k = a; k < b; k++){ - components.push_back(this->components[k]); + components.push_back(this->comps_[k]); } return std::make_shared>(components); } diff --git a/src/AffineMap.cpp b/src/AffineMap.cpp index 4c76e79f..72b774ac 100644 --- a/src/AffineMap.cpp +++ b/src/AffineMap.cpp @@ -3,13 +3,13 @@ #include "MParT/Utilities/ArrayConversions.h" #include "MParT/Initialization.h" - + using namespace mpart; template -AffineMap::AffineMap(StridedVector b) : ConditionalMapBase(b.size(),b.size(),0), - b_(b), +AffineMap::AffineMap(StridedVector b) : ConditionalMapBase(b.size(),b.size(),0), + b_(b), logDet_(0.0) {} @@ -27,7 +27,7 @@ AffineMap::AffineMap(StridedMatrix A, StridedVector b) : ConditionalMapBase(A.extent(1),A.extent(0),0), A_(A), b_(b) -{ +{ assert(A_.extent(0)<=A_.extent(1)); assert(A_.extent(0) == b_.extent(0)); Factorize(); @@ -86,7 +86,7 @@ void AffineMap::InverseImpl(StridedMatrix0){ @@ -94,7 +94,7 @@ void AffineMap::InverseImpl(StridedMatrix::InverseImpl(StridedMatrix xsub = Kokkos::subview(x1, std::make_pair(0,ncols-nrows), Kokkos::ALL()); StridedMatrix Asub = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(0,ncols-nrows)); @@ -124,13 +124,13 @@ void AffineMap::InverseImpl(StridedMatrix -void AffineMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, +void AffineMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) { return; @@ -146,7 +146,7 @@ void AffineMap::EvaluateImpl(StridedMatrix0){ @@ -160,22 +160,22 @@ void AffineMap::EvaluateImpl(StridedMatrix -void AffineMap::CoeffGradImpl(StridedMatrix const& pts, +void AffineMap::CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { return; } -template -void AffineMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, +template +void AffineMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) { return; } template -void AffineMap::GradientImpl(StridedMatrix const& pts, +void AffineMap::GradientImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { @@ -187,8 +187,16 @@ void AffineMap::GradientImpl(StridedMatrix +std::shared_ptr> AffineMap::Slice(int a, int b) { + StridedMatrix A_new; + StridedVector b_new; + if(A_.extent(0) > 0) A_new = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(a,b)); + if(b_.size() > 0) b_new = Kokkos::subview(b_, std::make_pair(a,b)); + return std::make_shared>(A_new, b_new); +} template class mpart::AffineMap; #if defined(MPART_ENABLE_GPU) template class mpart::AffineMap; -#endif +#endif diff --git a/src/ComposedMap.cpp b/src/ComposedMap.cpp index 42e2ded1..f0bd2760 100644 --- a/src/ComposedMap.cpp +++ b/src/ComposedMap.cpp @@ -8,7 +8,7 @@ using namespace mpart; template -ComposedMap::Checkpointer::Checkpointer(unsigned int maxSaves, +ComposedMap::Checkpointer::Checkpointer(unsigned int maxSaves, StridedMatrix initialPts, std::vector>>& comps) : maxSaves_(maxSaves), maps_(comps) @@ -23,7 +23,7 @@ ComposedMap::Checkpointer::Checkpointer(unsigned int maxSaves, checkpointLayers_.push_back(0); - // Allocate the workspace + // Allocate the workspace workspace1_ = Kokkos::View("xi", initialPts.extent(0), initialPts.extent(1)); workspace2_ = Kokkos::View("xj", initialPts.extent(0), initialPts.extent(1)); } @@ -37,7 +37,7 @@ int ComposedMap::Checkpointer::GetNextCheckpoint(unsigned int layer int currLayer = checkpointLayers_.back(); int oldLayer = currLayer; - + //ds = snaps - *check; int availSaves = int(maxSaves_) - checkpoints_.size(); if(availSaves==0){ @@ -46,12 +46,12 @@ int ComposedMap::Checkpointer::GetNextCheckpoint(unsigned int layer int reps = 0; int range = 1; - - while(range < layerInd+1 - currLayer) { + + while(range < layerInd+1 - currLayer) { reps += 1; - range = range*(reps + availSaves)/reps; + range = range*(reps + availSaves)/reps; } - + int bino1, bino2, bino3, bino4, bino5; bino1 = range*reps/(availSaves+reps); bino2 = (availSaves > 1) ? bino1*availSaves/(availSaves+reps-1) : 1; @@ -71,15 +71,15 @@ int ComposedMap::Checkpointer::GetNextCheckpoint(unsigned int layer if(layerInd-currLayer <= bino1 + bino3){ currLayer = currLayer + bino4; }else{ - if(layerInd-currLayer >= range - bino5){ - currLayer = currLayer + bino1; + if(layerInd-currLayer >= range - bino5){ + currLayer = currLayer + bino1; }else{ currLayer = layerInd-bino2-bino3; } } - if (currLayer == oldLayer) - currLayer = oldLayer + 1; + if (currLayer == oldLayer) + currLayer = oldLayer + 1; return currLayer; } @@ -87,7 +87,7 @@ int ComposedMap::Checkpointer::GetNextCheckpoint(unsigned int layer template Kokkos::View ComposedMap::Checkpointer::GetLayerInput(unsigned int layerInd) -{ +{ // Check if the layer asked for is the last checkpoint if(layerInd==checkpointLayers_.back()){ return checkpoints_.back(); @@ -99,21 +99,21 @@ Kokkos::View ComposedMap return GetLayerInput(layerInd); }else{ - - + + // Figure out how many more checkpoints we can make ("available saves") - - // Figure out the index of the next checkpoint + + // Figure out the index of the next checkpoint int nextCheckLayer = GetNextCheckpoint(layerInd); - + // Compute the input by reevaluating from the last checkpoint int i = checkpointLayers_.back(); maps_.at(i)->EvaluateImpl(checkpoints_.back(), workspace1_); ++i; for(; i("xi", workspace1_.extent(0), workspace1_.extent(1))); Kokkos::deep_copy(checkpoints_.back(),workspace1_); @@ -133,7 +133,7 @@ Kokkos::View ComposedMap template ComposedMap::ComposedMap(std::vector>> const& maps, - int maxChecks) : ConditionalMapBase(maps.front()->inputDim, + int maxChecks) : ConditionalMapBase(maps.front()->inputDim, maps.front()->inputDim, std::accumulate(maps.begin(), maps.end(), 0, [](size_t sum, std::shared_ptr> const& comp){ return sum + comp->numCoeffs; })), maps_(maps), @@ -202,14 +202,14 @@ void ComposedMap::LogDeterminantImpl(StridedMatrix intPts1("intermediate points 1", pts.extent(0), pts.extent(1)); Kokkos::View intPts2("intermediate points 2", pts.extent(0), pts.extent(1)); Kokkos::deep_copy(intPts1, pts); Kokkos::View compDetIncrement("Log Determinant", output.extent(0)); for(unsigned int i=1; iEvaluateImpl(intPts1, intPts2); @@ -232,7 +232,7 @@ void ComposedMap::EvaluateImpl(StridedMatrix intPts1("intermediate points 1", pts.extent(0), pts.extent(1)); Kokkos::View intPts2("intermediate points 2", pts.extent(0), pts.extent(1)); - + maps_.at(0)->EvaluateImpl(pts, intPts1); for(int j = 1; j::EvaluateImpl(StridedMatrix -void ComposedMap::GradientImpl(StridedMatrix const& pts, +void ComposedMap::GradientImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { @@ -257,14 +257,14 @@ void ComposedMap::GradientImpl(StridedMatrix intSens1("intermediate Sens", sens.extent(0), sens.extent(1)); Kokkos::View intSens2("intermediate Sens", sens.extent(0), sens.extent(1)); - Kokkos::deep_copy(intSens1,sens); - + Kokkos::deep_copy(intSens1,sens); + Checkpointer checker(maxChecks_, pts, maps_); for(int i = maps_.size() - 1; i>=0; --i){ - + // x* = T_{i-1} o ... o T_1(x) - auto input = checker.GetLayerInput(i); + auto input = checker.GetLayerInput(i); //s_{i-1}^T = s_{i}^T J_i(x*) maps_.at(i)->GradientImpl(input, intSens1, intSens2); @@ -281,7 +281,7 @@ void ComposedMap::InverseImpl(StridedMatrix output) { // We assume each component of the composed map is square, so the x1 input is not used. - // We can just loop backward through the layers + // We can just loop backward through the layers // intermediate vals Kokkos::View intR1("intermediate r1", r.extent(0), r.extent(1)); @@ -291,7 +291,7 @@ void ComposedMap::InverseImpl(StridedMatrix=0; --i){ - + maps_.at(i)->InverseImpl(x1, intR1, intR2); // <- Note the x1 doesn't matter here because the layer is square simple_swap(intR1,intR2); } @@ -300,13 +300,13 @@ void ComposedMap::InverseImpl(StridedMatrix -void ComposedMap::CoeffGradImpl(StridedMatrix const& pts, +void ComposedMap::CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { // T(x) = T_{n} o ... o T_1(x), T_i coeffs w_i - // s^T \nabla_{w_i} T(x) = s^T J_n * J_{n-1} *...*J_{i-1}* \nabla_w_i T_i + // s^T \nabla_{w_i} T(x) = s^T J_n * J_{n-1} *...*J_{i-1}* \nabla_w_i T_i Kokkos::View intSens1("intermediate sens 1", sens.extent(0), sens.extent(1)); Kokkos::View intSens2("intermediate sens 2", sens.extent(0), sens.extent(1)); Kokkos::deep_copy(intSens1, sens); @@ -314,15 +314,15 @@ void ComposedMap::CoeffGradImpl(StridedMatrix subOut; - int endParamDim = this->numCoeffs; + int endParamDim = this->numCoeffs; for(int i = maps_.size() - 1; i>=0; --i){ - // intPts1 = T_{i-1} o ... o T_1(x) + // intPts1 = T_{i-1} o ... o T_1(x) auto input = checker.GetLayerInput(i); // finish g_i = s^T \nabla_{w_i} T(x) - subOut = Kokkos::subview(output, - std::make_pair(int(endParamDim-maps_.at(i)->numCoeffs), endParamDim), + subOut = Kokkos::subview(output, + std::make_pair(int(endParamDim-maps_.at(i)->numCoeffs), endParamDim), Kokkos::ALL()); @@ -342,43 +342,43 @@ void ComposedMap::CoeffGradImpl(StridedMatrix -void ComposedMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, +void ComposedMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) { StridedMatrix subOut; Kokkos::View intSens1("intermediate Sens", pts.extent(0), pts.extent(1)); Kokkos::View intSens2("intermediate Sens", pts.extent(0), pts.extent(1)); - + // Get the gradient of the log determinant contribution from the last component Checkpointer checker(maxChecks_, pts, maps_); - auto input = checker.GetLayerInput(maps_.size()-1); - + auto input = checker.GetLayerInput(maps_.size()-1); + int endParamDim = this->numCoeffs; - subOut = Kokkos::subview(output, - std::make_pair(int(endParamDim-maps_.back()->numCoeffs), endParamDim), + subOut = Kokkos::subview(output, + std::make_pair(int(endParamDim-maps_.back()->numCoeffs), endParamDim), Kokkos::ALL()); maps_.back()->LogDeterminantCoeffGradImpl(input, subOut); // Get the sensitivity of this log determinant term wrt changes in the input maps_.back()->LogDeterminantInputGradImpl(input, intSens1); - - endParamDim -= maps_.back()->numCoeffs; + + endParamDim -= maps_.back()->numCoeffs; for(int i = maps_.size() - 2; i>=0; --i){ - + // Compute input to this layer input = checker.GetLayerInput(i); // Gradient for direct contribution of these parameters on the log determinant - subOut = Kokkos::subview(output, - std::make_pair(int(endParamDim-maps_.at(i)->numCoeffs), endParamDim), + subOut = Kokkos::subview(output, + std::make_pair(int(endParamDim-maps_.at(i)->numCoeffs), endParamDim), Kokkos::ALL()); - + maps_.at(i)->LogDeterminantCoeffGradImpl(input, subOut); - // Gradient of later log determinant terms on the coefficients of this layer + // Gradient of later log determinant terms on the coefficients of this layer Kokkos::View subOut2("temp", maps_.at(i)->numCoeffs, pts.extent(1)); maps_.at(i)->CoeffGradImpl(input, intSens1, subOut2); subOut += subOut2; @@ -389,31 +389,31 @@ void ComposedMap::LogDeterminantCoeffGradImpl(StridedMatrix(intSens1, intSens2); // Add sensitivity of log determinant to input - maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2); - - intSens1 += intSens2; + maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2); + + intSens1 += intSens2; } - endParamDim -= maps_.at(i)->numCoeffs; + endParamDim -= maps_.at(i)->numCoeffs; } } template -void ComposedMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, +void ComposedMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) { Kokkos::View intSens1("intermediate Sens", pts.extent(0), pts.extent(1)); Kokkos::View intSens2("intermediate Sens", pts.extent(0), pts.extent(1)); - + // Get the gradient of the log determinant contribution from the last component Checkpointer checker(maxChecks_, pts, maps_); auto input = checker.GetLayerInput(maps_.size()-1); - + maps_.back()->LogDeterminantInputGradImpl(input, intSens1); for(int i = maps_.size() - 2; i>=0; --i){ - + // reset intPts1 to initial pts input = checker.GetLayerInput(i); @@ -421,13 +421,27 @@ void ComposedMap::LogDeterminantInputGradImpl(StridedMatrixGradientImpl(input, intSens1, intSens2); simple_swap(intSens1, intSens2); - maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2); - intSens1 += intSens2; + maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2); + intSens1 += intSens2; } Kokkos::deep_copy(output, intSens1); } +template +std::shared_ptr> ComposedMap::Slice(int a, int b) { + assert(a < b && 0 <= a && b <= this->outputDim); + std::shared_ptr> lastMap = maps_[maps_.size()-1]; + if (lastMap->outputDim == b-a) { + return std::shared_ptr>(this); + } else { + std::vector>> newMaps; + newMaps = maps_; + newMaps[maps_.size()-1] = lastMap->Slice(a, b); + return std::make_shared>(newMaps); + } +} + // Explicit template instantiation template class mpart::ComposedMap; #if defined(MPART_ENABLE_GPU) diff --git a/src/IdentityMap.cpp b/src/IdentityMap.cpp index 01c6297c..538351ed 100644 --- a/src/IdentityMap.cpp +++ b/src/IdentityMap.cpp @@ -49,7 +49,7 @@ void IdentityMap::InverseImpl(StridedMatrix -void IdentityMap::CoeffGradImpl(StridedMatrix const& pts, +void IdentityMap::CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { @@ -57,7 +57,7 @@ void IdentityMap::CoeffGradImpl(StridedMatrix -void IdentityMap::GradientImpl(StridedMatrix const& pts, +void IdentityMap::GradientImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) { @@ -77,16 +77,16 @@ void IdentityMap::GradientImpl(StridedMatrix -void IdentityMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, +void IdentityMap::LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) { assert(false); } template -void IdentityMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, +void IdentityMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) -{ +{ Kokkos::MDRangePolicy, typename MemoryToExecution::Space> zeroPolicy({0, 0}, {output.extent(0), output.extent(1)}); Kokkos::parallel_for(zeroPolicy, KOKKOS_LAMBDA(const int& i, const int& j) { output(i,j) = 0.0; @@ -94,6 +94,14 @@ void IdentityMap::LogDeterminantInputGradImpl(StridedMatrix +std::shared_ptr> IdentityMap::Slice(int a, int b) { + assert(a >= 0); + assert(b <= this->outputDim); + assert(a < b); + return std::make_shared>(this->inputDim, b-a); +} + // Explicit template instantiation template class mpart::IdentityMap; #if defined(KOKKOS_ENABLE_CUDA ) || defined(KOKKOS_ENABLE_SYCL) diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index 0dc1b896..c83e49be 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -106,7 +106,7 @@ void TriangularMap::WrapCoeffs(Kokkos::View void TriangularMap::LogDeterminantImpl(StridedMatrix const& pts, StridedVector output) -{ +{ // Evaluate the log determinant for the first component StridedMatrix subPts = Kokkos::subview(pts, std::make_pair(0,int(comps_.at(0)->inputDim)), Kokkos::ALL()); comps_.at(0)->LogDeterminantImpl(subPts, output); @@ -200,8 +200,8 @@ void TriangularMap::GradientImpl(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){ @@ -209,7 +209,7 @@ void TriangularMap::GradientImpl(StridedMatrix void TriangularMap::LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) { - // Initialize the output to zero + // Initialize the output to zero Kokkos::MDRangePolicy, typename MemoryToExecution::Space> zeroPolicy({0, 0}, {output.extent(0), output.extent(1)}); Kokkos::parallel_for(zeroPolicy, KOKKOS_LAMBDA(const int& i, const int& j) { output(i,j) = 0.0; @@ -289,10 +289,10 @@ void TriangularMap::LogDeterminantInputGradImpl(StridedMatrix subPts; StridedMatrix subOut; - int numPts = pts.extent(1); + int numPts = pts.extent(1); Kokkos::View compGrad("Component Gradient", this->inputDim, numPts); Kokkos::View subGrad; - + for(unsigned int i=0; iinputDim; subPts = Kokkos::subview(pts, std::make_pair(0,compDim), Kokkos::ALL()); @@ -309,6 +309,45 @@ void TriangularMap::LogDeterminantInputGradImpl(StridedMatrix +std::shared_ptr> TriangularMap::Slice(int a, int b) { + std::vector>> components; + // TODO: Handle empty case + if( a < 0 || a >= b || b > this->outputDim ){ + throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); + } + int accum_a = 0; + int k_a = 0; + //TODO: Check that this is correct + while(accum_a < a){ + k_a++; + accum_a += this->comps_[k_a]->outputDim; + } + if(a != accum_a){ + accum_a -= this->comps_[k_a]->outputDim; + } + int accum_b = 0; + int k_b = k_a; + //TODO: Check that this is correct + while(accum_b < b){ + k_b++; + accum_b += this->comps_[k_b]->outputDim; + } + if(b != accum_b){ + accum_b -= this->comps_[k_b]->outputDim; + } + if(k_a == k_b){ + components.push_back(this->comps_[k_a]->Slice(a-accum_a, b-accum_b)); + } else { + components.push_back(this->comps_[k_a]->Slice(a-accum_a, this->comps_[k_a]->outputDim)); + for(int k = k_a+1; k < k_b; k++){ + components.push_back(this->comps_[k]); + } + components.push_back(this->comps_[k_b]->Slice(0, b-accum_b)); + } + return std::make_shared>(components); + } + // Explicit template instantiation template class mpart::TriangularMap; #if defined(MPART_ENABLE_GPU) diff --git a/tests/Test_ConditionalMapBase.cpp b/tests/Test_ConditionalMapBase.cpp index f51fb546..bbb10de9 100644 --- a/tests/Test_ConditionalMapBase.cpp +++ b/tests/Test_ConditionalMapBase.cpp @@ -8,18 +8,22 @@ using MemorySpace = Kokkos::HostSpace; class MyIdentityMap : public ConditionalMapBase{ public: - MyIdentityMap(unsigned int dim, unsigned int numCoeffs) : ConditionalMapBase(dim,dim,numCoeffs){}; + MyIdentityMap(unsigned int dim, unsigned int numCoeffs, unsigned int dimOut = 0) : ConditionalMapBase(dim,!dimOut? dim : dimOut,numCoeffs){}; virtual ~MyIdentityMap() = default; virtual void EvaluateImpl(StridedMatrix const& pts, - StridedMatrix output) override{Kokkos::deep_copy(output,pts);}; + StridedMatrix output) override{ + int a = this->inputDim - this->outputDim; + auto pts_view = Kokkos::subview(pts, std::make_pair(a, int(this->inputDim)), Kokkos::ALL()); + Kokkos::deep_copy(output,pts_view); + }; - virtual void GradientImpl(StridedMatrix const& pts, + virtual void GradientImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override { - assert(false); + assert(false); } virtual void LogDeterminantImpl(StridedMatrix const&, @@ -32,25 +36,30 @@ class MyIdentityMap : public ConditionalMapBase{ StridedMatrix const& r, StridedMatrix output) override{Kokkos::deep_copy(output,r);}; - virtual void CoeffGradImpl(StridedMatrix const& pts, + virtual void CoeffGradImpl(StridedMatrix const& pts, StridedMatrix const& sens, StridedMatrix output) override { - assert(false); + assert(false); } - virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, + virtual void LogDeterminantCoeffGradImpl(StridedMatrix const& pts, StridedMatrix output) override - { + { assert(false); } - virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, + virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) override - { + { assert(false); } + + // Only creates a slice of the tail of the input + std::shared_ptr> Slice(int a, int b) override{ + return std::make_shared(this->inputDim, 0, b-a); + } }; @@ -241,5 +250,47 @@ TEST_CASE( "Testing inverse evaluation of an identity conditional map", "[Condit } } } +} + +TEST_CASE( "Testing slicing evaluation of an identity conditional map", "[ConditionalMapBaseSlice]" ) { + unsigned int dim = 7; + unsigned int numPts = 5; + MyIdentityMap map(dim,0); + + int a = 2; + int b = 5; + std::shared_ptr> mapSlice = map.Slice(a,b); + + int sliceOutdim = b-a; + + CHECK(map.inputDim == dim); + CHECK(map.outputDim == dim); + + CHECK(mapSlice->inputDim == dim); + CHECK(mapSlice->outputDim == sliceOutdim); + + SECTION("Using Kokkos"){ + + Kokkos::View pts("pts", dim, numPts); + + for(unsigned int i=0; i output = mapSlice->Evaluate(pts); + + REQUIRE(output.extent(0)==sliceOutdim); + REQUIRE(output.extent(1)==numPts); + + int idx_start = map.inputDim - (b-a); + + for(unsigned int i=0; i Date: Wed, 19 Oct 2022 12:28:29 -0400 Subject: [PATCH 08/16] affix kokkos to 3.6.01 --- .github/workflows/build-external-lib-tests.yml | 3 ++- src/TriangularMap.cpp | 4 ++-- tests/Test_TriangularMap.cpp | 14 ++++++-------- 3 files changed, 10 insertions(+), 11 deletions(-) diff --git a/.github/workflows/build-external-lib-tests.yml b/.github/workflows/build-external-lib-tests.yml index 5f0e66a0..14d89439 100644 --- a/.github/workflows/build-external-lib-tests.yml +++ b/.github/workflows/build-external-lib-tests.yml @@ -3,7 +3,7 @@ name: external-lib-tests on: push: branches: - - dannys4/issue50 + - main pull_request: {} jobs: @@ -27,6 +27,7 @@ jobs: with: repository: kokkos/kokkos path: kokkos + ref: '3.6.01' - name: Install Kokkos run: | diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index c83e49be..8f5a387d 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -319,7 +319,7 @@ std::shared_ptr> TriangularMap::Sli int accum_a = 0; int k_a = 0; //TODO: Check that this is correct - while(accum_a < a){ + while(accum_a < a && k_a < this->comps_.size()-1){ k_a++; accum_a += this->comps_[k_a]->outputDim; } @@ -329,7 +329,7 @@ std::shared_ptr> TriangularMap::Sli int accum_b = 0; int k_b = k_a; //TODO: Check that this is correct - while(accum_b < b){ + while(accum_b < b && k_b < this->comps_.size()-1){ k_b++; accum_b += this->comps_[k_b]->outputDim; } diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index 5852369a..121e5643 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -226,16 +226,14 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ auto slice = triMap->Slice(sliceBegin, sliceEnd); REQUIRE(slice->inputDim == triMap->inputDim); REQUIRE(slice->outputDim == sliceEnd-sliceBegin); + + auto outSlice = slice->Evaluate(in); // Just test the evaluation for(unsigned int i=sliceBegin; iEvaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); - - REQUIRE(outBlock.extent(1)==numSamps); - REQUIRE(outBlock.extent(0)==1); - for(unsigned int j=0; j Date: Wed, 19 Oct 2022 12:41:29 -0400 Subject: [PATCH 09/16] Fix Bindings for Slice --- .github/workflows/build-external-lib-tests.yml | 4 ++-- bindings/julia/src/ConditionalMapBase.cpp | 1 + bindings/julia/src/TriangularMap.cpp | 1 - bindings/python/src/ConditionalMapBase.cpp | 6 ++++++ bindings/python/src/TriangularMap.cpp | 6 ------ 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/.github/workflows/build-external-lib-tests.yml b/.github/workflows/build-external-lib-tests.yml index 14d89439..3b1e9ba7 100644 --- a/.github/workflows/build-external-lib-tests.yml +++ b/.github/workflows/build-external-lib-tests.yml @@ -39,7 +39,7 @@ jobs: -DKokkos_CXX_STANDARD=17 \ -DCMAKE_INSTALL_PREFIX=$GITHUB_WORKSPACE/KOKKOS_INSTALL/ \ ../ - sudo make install + sudo make -j$(nproc) install - name: Install Eigen, Catch2, Pybind11 shell: bash -l {0} @@ -54,7 +54,7 @@ jobs: - name: Build MParT shell: bash -l {0} - run: cd $GITHUB_WORKSPACE/mpart/build; make + run: cd $GITHUB_WORKSPACE/mpart/build; make -j$(nproc) - name: Run Tests shell: bash -l {0} diff --git a/bindings/julia/src/ConditionalMapBase.cpp b/bindings/julia/src/ConditionalMapBase.cpp index 99926541..f4dd67a5 100644 --- a/bindings/julia/src/ConditionalMapBase.cpp +++ b/bindings/julia/src/ConditionalMapBase.cpp @@ -32,6 +32,7 @@ void mpart::binding::ConditionalMapBaseWrapper(jlcxx::Module &mod) { map.InverseImpl(JuliaToKokkos(x1), JuliaToKokkos(r), JuliaToKokkos(output)); return output; }) + .method("Slice", [](ConditionalMapBase& map, int begin, int end){ return map.Slice(begin-1, end); }) ; jlcxx::stl::apply_stl*>(mod); } \ No newline at end of file diff --git a/bindings/julia/src/TriangularMap.cpp b/bindings/julia/src/TriangularMap.cpp index d4449753..defb643d 100644 --- a/bindings/julia/src/TriangularMap.cpp +++ b/bindings/julia/src/TriangularMap.cpp @@ -15,6 +15,5 @@ void mpart::binding::TriangularMapWrapper(jlcxx::Module &mod) { map.InverseInplace(JuliaToKokkos(x), JuliaToKokkos(r)); }) .method("GetComponent", [](TriangularMap& map, int i){ return map.GetComponent(i-1); }) - .method("Slice", [](TriangularMap& map, int begin, int end){ return map.Slice(begin-1, end); }) ; } \ No newline at end of file diff --git a/bindings/python/src/ConditionalMapBase.cpp b/bindings/python/src/ConditionalMapBase.cpp index 5cba8a09..5404178f 100644 --- a/bindings/python/src/ConditionalMapBase.cpp +++ b/bindings/python/src/ConditionalMapBase.cpp @@ -22,6 +22,12 @@ void mpart::binding::ConditionalMapBaseWrapper(py::module &m) .def("Inverse", static_cast::*)(Eigen::Ref const&, Eigen::Ref const&)>(&ConditionalMapBase::Inverse)) .def("LogDeterminantCoeffGrad", static_cast::*)(Eigen::Ref const&)>(&ConditionalMapBase::LogDeterminantCoeffGrad)) .def("GetBaseFunction", &ConditionalMapBase::GetBaseFunction) + .def("__getitem__", [](ConditionalMapBase &m, const py::slice &s) { + size_t start = 0, stop = 0, step = 0, slicelength = 0; + s.compute(m.outputDim, &start, &stop, &step, &slicelength); + if(step != 1) throw std::runtime_error("Map slice step size must be 1"); + return m.Slice(start, stop); + }) ; } diff --git a/bindings/python/src/TriangularMap.cpp b/bindings/python/src/TriangularMap.cpp index 94daa23e..379078b6 100644 --- a/bindings/python/src/TriangularMap.cpp +++ b/bindings/python/src/TriangularMap.cpp @@ -22,12 +22,6 @@ void mpart::binding::TriangularMapWrapper(py::module &m) .def("GetComponent", &TriangularMap::GetComponent) .def("Slice", &TriangularMap::Slice) .def("__getitem__", &TriangularMap::GetComponent) - .def("__getitem__", [](const TriangularMap &m, const py::slice &s) { - size_t start = 0, stop = 0, step = 0, slicelength = 0; - s.compute(m.outputDim, &start, &stop, &step, &slicelength); - if(step != 1) throw std::runtime_error("TriangularMap slice step must be 1"); - return m.Slice(start, stop); - }) ; } From cde0fe6f7e7efd672c8ad7734f51bb6b6e17d525 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Wed, 19 Oct 2022 15:32:43 -0400 Subject: [PATCH 10/16] Working on slice --- MParT/MapFactory.h | 2 +- MParT/TriangularMap.h | 10 +---- src/MapFactory.cpp | 12 +++--- src/TriangularMap.cpp | 81 ++++++++++++++++++++---------------- tests/Test_TriangularMap.cpp | 16 ++++--- 5 files changed, 66 insertions(+), 55 deletions(-) diff --git a/MParT/MapFactory.h b/MParT/MapFactory.h index e5355892..fbe46aeb 100644 --- a/MParT/MapFactory.h +++ b/MParT/MapFactory.h @@ -62,7 +62,7 @@ namespace mpart{ @param options Additional options that will be passed on to CreateComponent to construct each MonotoneComponent. */ template - std::shared_ptr> CreateTriangular(unsigned int inputDim, + std::shared_ptr> CreateTriangular(unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder, MapOptions options = MapOptions()); diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h index 8c7a4288..075eda2c 100644 --- a/MParT/TriangularMap.h +++ b/MParT/TriangularMap.h @@ -65,15 +65,9 @@ class TriangularMap : public ConditionalMapBase{ * the number of outputs of each component is 1. * @param a The first component block to take * @param b The index AFTER the last component block to take - * @return A shared pointer to a new TriangularMap object containing the subsection of the map. + * @return A shared pointer to a new ConditionalMapBase object containing the subsection of the map. */ - std::shared_ptr> BlockSlice(int a, int b) { - std::vector>> components; - for(int k = a; k < b; k++){ - components.push_back(this->comps_[k]); - } - return std::make_shared>(components); - } + virtual std::shared_ptr> BlockSlice(int a, int b); virtual std::shared_ptr> GetComponent(unsigned int i){ return comps_.at(i);} diff --git a/src/MapFactory.cpp b/src/MapFactory.cpp index d35f30fc..3c923fd4 100644 --- a/src/MapFactory.cpp +++ b/src/MapFactory.cpp @@ -20,10 +20,10 @@ std::shared_ptr> mpart::MapFactory::CreateCompon } template -std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int inputDim, - unsigned int outputDim, - unsigned int totalOrder, - MapOptions options) +std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int inputDim, + unsigned int outputDim, + unsigned int totalOrder, + MapOptions options) { std::vector>> comps(outputDim); @@ -91,9 +91,9 @@ std::shared_ptr> mpart::MapFactory::Creat template std::shared_ptr> mpart::MapFactory::CreateComponent(FixedMultiIndexSet const&, MapOptions); template std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int, FixedMultiIndexSet const&, MapOptions); -template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); +template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); #if defined(MPART_ENABLE_GPU) template std::shared_ptr> mpart::MapFactory::CreateComponent(FixedMultiIndexSet const&, MapOptions); template std::shared_ptr> mpart::MapFactory::CreateExpansion(unsigned int, FixedMultiIndexSet const&, MapOptions); - template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); + template std::shared_ptr> mpart::MapFactory::CreateTriangular(unsigned int, unsigned int, unsigned int, MapOptions); #endif \ No newline at end of file diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index 8f5a387d..92c7cea5 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -311,43 +311,54 @@ void TriangularMap::LogDeterminantInputGradImpl(StridedMatrix std::shared_ptr> TriangularMap::Slice(int a, int b) { - std::vector>> components; - // TODO: Handle empty case - if( a < 0 || a >= b || b > this->outputDim ){ - throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); - } - int accum_a = 0; - int k_a = 0; - //TODO: Check that this is correct - while(accum_a < a && k_a < this->comps_.size()-1){ - k_a++; - accum_a += this->comps_[k_a]->outputDim; - } - if(a != accum_a){ - accum_a -= this->comps_[k_a]->outputDim; - } - int accum_b = 0; - int k_b = k_a; - //TODO: Check that this is correct - while(accum_b < b && k_b < this->comps_.size()-1){ - k_b++; - accum_b += this->comps_[k_b]->outputDim; - } - if(b != accum_b){ - accum_b -= this->comps_[k_b]->outputDim; - } - if(k_a == k_b){ - components.push_back(this->comps_[k_a]->Slice(a-accum_a, b-accum_b)); - } else { - components.push_back(this->comps_[k_a]->Slice(a-accum_a, this->comps_[k_a]->outputDim)); - for(int k = k_a+1; k < k_b; k++){ - components.push_back(this->comps_[k]); - } - components.push_back(this->comps_[k_b]->Slice(0, b-accum_b)); - } - return std::make_shared>(components); + std::vector>> components; + // TODO: Handle empty case + if( a < 0 || a >= b || b > this->outputDim ) + throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); + + int accum_a = 0; + int k_a = 0; + //TODO: Check that this is correct + while(k_a < this->comps_.size()){ + if(accum_a + this->comps_[k_a]->outputDim > a) + break; + accum_a += this->comps_[k_a]->outputDim; + k_a++; + } + + int accum_b = accum_a; + int k_b = k_a; + //TODO: Check that this is correct + while(k_b < this->comps_.size()){ + if(accum_b + this->comps_[k_b]->outputDim >= b) + break; + accum_b += this->comps_[k_b]->outputDim; + k_b++; } + std::cerr << "a = " << a << ", accum_a = " << accum_a << ", k_a = " << k_a << std::endl; + std::cerr << "b = " << b << ", accum_b = " << accum_b << ", k_b = " << k_b << std::endl; + if(k_a == k_b){ + components.push_back(this->comps_[k_a]->Slice(a-accum_a, b-accum_b)); + } else { + components = std::vector>>(this->comps_.begin()+k_a, this->comps_.begin()+k_b+1); + components[0] = this->comps_[k_a]->Slice(a-accum_a, this->comps_[k_a]->outputDim); + components[components.size()-1] = this->comps_[k_b]->Slice(0, b-accum_b); + } + auto output = std::make_shared>(components); + output->SetCoeffs(Kokkos::View("Component Coefficients", output->numCoeffs)); + return output; +} + +template +std::shared_ptr> TriangularMap::BlockSlice(int a, int b) { + std::vector>> components; + for(int k = a; k < b; k++){ + components.push_back(this->comps_[k]); + } + return std::make_shared>(components); +} + // Explicit template instantiation template class mpart::TriangularMap; #if defined(MPART_ENABLE_GPU) diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index 121e5643..51153ae5 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -219,7 +219,6 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ } - /* SECTION("Slice"){ int sliceBegin = 1; int sliceEnd = 3; @@ -228,12 +227,19 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents", "[TriangularMap_ REQUIRE(slice->outputDim == sliceEnd-sliceBegin); auto outSlice = slice->Evaluate(in); + // Just test the evaluation for(unsigned int i=sliceBegin; iEvaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); + + REQUIRE(outBlock.extent(1)==numSamps); + REQUIRE(outBlock.extent(0)==1); + for(unsigned int j=0; j>(slice); + } } \ No newline at end of file From aa41b54483a1b6450ac0e890c6b9f39116a4912c Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Wed, 26 Oct 2022 09:43:30 -0400 Subject: [PATCH 11/16] CMB --- MParT/ConditionalMapBase.h | 4 +++- src/ConditionalMapBase.cpp | 31 ++++++++++++++++++------------- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/MParT/ConditionalMapBase.h b/MParT/ConditionalMapBase.h index 00d192e5..e98021a8 100644 --- a/MParT/ConditionalMapBase.h +++ b/MParT/ConditionalMapBase.h @@ -150,7 +150,9 @@ namespace mpart { virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) = 0; - virtual std::shared_ptr> Slice(int a, int b) = 0; + std::shared_ptr> Slice(int a, int b); + + virtual ConditionalMapBase SliceImpl(int a, int b); }; // class ConditionalMapBase } diff --git a/src/ConditionalMapBase.cpp b/src/ConditionalMapBase.cpp index fcf0193a..1704a0c7 100644 --- a/src/ConditionalMapBase.cpp +++ b/src/ConditionalMapBase.cpp @@ -40,27 +40,27 @@ Kokkos::View ConditionalMapBase template<> template<> Kokkos::View ConditionalMapBase::LogDeterminant(StridedMatrix const& pts) -{ +{ return ToHost( this->LogDeterminant( ToDevice(pts) )); } template<> template<> Kokkos::View ConditionalMapBase::LogDeterminant(StridedMatrix const& pts) -{ +{ return ToDevice( this->LogDeterminant( ToHost(pts) )); } template<> Eigen::VectorXd ConditionalMapBase::LogDeterminant(Eigen::Ref const& pts) -{ +{ StridedMatrix pts_view = ConstRowMatToKokkos(pts); Kokkos::View outView = ToHost( this->LogDeterminant( ToDevice(pts_view) )); return KokkosToVec(outView); } -#endif +#endif template<> template<> @@ -138,7 +138,7 @@ Eigen::RowMatrixXd ConditionalMapBase::Inverse(Eigen::Ref template<> @@ -150,6 +150,11 @@ StridedMatrix ConditionalMapBase:: return output; } +template +std::shared_ptr> Slice(int a, int b) { + return std::shared_ptr>(this->SliceImpl(a,b)); +} + template<> Eigen::RowMatrixXd ConditionalMapBase::LogDeterminantCoeffGrad(Eigen::Ref const& pts) { @@ -177,10 +182,10 @@ template<> template<> StridedMatrix ConditionalMapBase::LogDeterminantCoeffGrad(StridedMatrix const& pts) { - // Copy the points to the device space + // Copy the points to the device space StridedMatrix pts_device = ToDevice(pts); - // Evaluate on the device space + // Evaluate on the device space StridedMatrix evals_device = this->LogDeterminantCoeffGrad(pts_device); // Copy back to the host space @@ -191,10 +196,10 @@ template<> template<> StridedMatrix ConditionalMapBase::LogDeterminantCoeffGrad(StridedMatrix const& pts) { - // Copy the points to the host + // Copy the points to the host StridedMatrix pts_host = ToHost(pts); - // Evaluate on the host + // Evaluate on the host StridedMatrix evals_host = this->LogDeterminantCoeffGrad(pts_host); // Copy back to the device @@ -254,10 +259,10 @@ template<> template<> StridedMatrix ConditionalMapBase::LogDeterminantInputGrad(StridedMatrix const& pts) { - // Copy the points to the device space + // Copy the points to the device space StridedMatrix pts_device = ToDevice(pts); - // Evaluate on the device space + // Evaluate on the device space StridedMatrix evals_device = this->LogDeterminantInputGrad(pts_device); // Copy back to the host space @@ -268,10 +273,10 @@ template<> template<> StridedMatrix ConditionalMapBase::LogDeterminantInputGrad(StridedMatrix const& pts) { - // Copy the points to the host + // Copy the points to the host StridedMatrix pts_host = ToHost(pts); - // Evaluate on the host + // Evaluate on the host StridedMatrix evals_host = this->LogDeterminantInputGrad(pts_host); // Copy back to the device From 449668f7a2defddd3dccf4408a5b9306aed542ad Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Sun, 8 Jan 2023 14:22:46 -0500 Subject: [PATCH 12/16] Changing slice --- MParT/ConditionalMapBase.h | 3 +-- src/ConditionalMapBase.cpp | 5 ----- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/MParT/ConditionalMapBase.h b/MParT/ConditionalMapBase.h index 35f932c1..29acf142 100644 --- a/MParT/ConditionalMapBase.h +++ b/MParT/ConditionalMapBase.h @@ -150,9 +150,8 @@ namespace mpart { virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) = 0; - std::shared_ptr> Slice(int a, int b); + virtual std::shared_ptr> Slice(int a, int b) = 0; - virtual ConditionalMapBase SliceImpl(int a, int b); #if defined(MPART_HAS_CEREAL) template void serialize(Archive& ar) { diff --git a/src/ConditionalMapBase.cpp b/src/ConditionalMapBase.cpp index 1704a0c7..ff09ddc0 100644 --- a/src/ConditionalMapBase.cpp +++ b/src/ConditionalMapBase.cpp @@ -150,11 +150,6 @@ StridedMatrix ConditionalMapBase:: return output; } -template -std::shared_ptr> Slice(int a, int b) { - return std::shared_ptr>(this->SliceImpl(a,b)); -} - template<> Eigen::RowMatrixXd ConditionalMapBase::LogDeterminantCoeffGrad(Eigen::Ref const& pts) { From 608cc66eb38054f348106d351633f7c6dab6de14 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Mon, 9 Jan 2023 10:02:59 -0500 Subject: [PATCH 13/16] Start enable shared from this --- MParT/ConditionalMapBase.h | 4 +++- MParT/MonotoneComponent.h | 6 ++++-- tests/Test_TriangularMap.cpp | 1 + 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/MParT/ConditionalMapBase.h b/MParT/ConditionalMapBase.h index 29acf142..9fbab75b 100644 --- a/MParT/ConditionalMapBase.h +++ b/MParT/ConditionalMapBase.h @@ -150,7 +150,9 @@ namespace mpart { virtual void LogDeterminantInputGradImpl(StridedMatrix const& pts, StridedMatrix output) = 0; - virtual std::shared_ptr> Slice(int a, int b) = 0; + virtual std::shared_ptr> Slice(int a, int b) { + throw std::runtime_error("Slice not implemented for this map type."); + } #if defined(MPART_HAS_CEREAL) template diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h index 2f609b91..58457d46 100644 --- a/MParT/MonotoneComponent.h +++ b/MParT/MonotoneComponent.h @@ -1,6 +1,8 @@ #ifndef MPART_MONOTONECOMPONENT_H #define MPART_MONOTONECOMPONENT_H +#include + #include "MParT/MultiIndices/FixedMultiIndexSet.h" #include "MParT/MultiIndices/MultiIndexSet.h" @@ -37,7 +39,7 @@ T(x_1, x_2, ..., x_D) = f(x_1,x_2,..., x_{D-1}, 0) + \int_0^{x_D} g\left( \part @tparam QuadratureType A class defining the integration scheme used to approximate \f$\int_0^{x_N} g\left( \frac{\partial f}{\partial x_d}(x_1,x_2,..., x_{N-1}, t) \right) dt\f$. The type must have a function `Integrate(f,lb,ub)` that accepts a functor `f`, a double lower bound `lb`, a double upper bound `ub`, and returns a double with an estimate of the integral. The MParT::AdaptiveSimpson and MParT::RecursiveQuadrature classes provide this interface. */ template -class MonotoneComponent : public ConditionalMapBase +class MonotoneComponent : public ConditionalMapBase, public std::enable_shared_from_this> { public: @@ -1080,7 +1082,7 @@ class MonotoneComponent : public ConditionalMapBase std::shared_ptr> Slice(int a, int b) override { assert(a == 0 && b == 1); - return std::shared_ptr>(this); + return shared_from_this(); } /** Give access to the underlying FixedMultiIndexSet diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index bf3333fa..9164f118 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -903,6 +903,7 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef SECTION("Slice"){ int sliceBegin = 1; int sliceEnd = 3; + unsigned int extraInputs = 2; auto slice = triMap->Slice(sliceBegin, sliceEnd); REQUIRE(slice->inputDim == triMap->inputDim); REQUIRE(slice->outputDim == sliceEnd-sliceBegin); From 6202974d8ef145f7bccf71d8c9dc8f8c48235469 Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Mon, 9 Jan 2023 10:26:19 -0500 Subject: [PATCH 14/16] Working Triangular Test --- MParT/MonotoneComponent.h | 2 +- src/TriangularMap.cpp | 18 +++++++++--------- tests/Test_TriangularMap.cpp | 14 ++++++-------- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h index 58457d46..0812e5ae 100644 --- a/MParT/MonotoneComponent.h +++ b/MParT/MonotoneComponent.h @@ -1082,7 +1082,7 @@ class MonotoneComponent : public ConditionalMapBase, public std::en std::shared_ptr> Slice(int a, int b) override { assert(a == 0 && b == 1); - return shared_from_this(); + return MonotoneComponent::shared_from_this(); } /** Give access to the underlying FixedMultiIndexSet diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index 318f1b25..2e892667 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -12,10 +12,10 @@ TriangularMap::TriangularMap(std::vector> const& comp){ return sum + comp->numCoeffs; })), comps_(components) { - + // Check the sizes of all the inputs - + for(unsigned int i=0; ioutputDim > comps_.at(i)->inputDim){ std::stringstream msg; @@ -42,7 +42,7 @@ TriangularMap::TriangularMap(std::vector coeffs("coeffs", this->numCoeffs); unsigned int cumNumCoeffs = 0; - + for(unsigned int i=0; iCheckCoefficients()){ @@ -281,7 +281,7 @@ void TriangularMap::CoeffGradImpl(StridedMatrixnumCoeffs)), Kokkos::ALL()); comps_.at(i)->CoeffGradImpl(subPts, subSens, subOut); - + startParamDim += comps_.at(i)->numCoeffs; } @@ -350,10 +350,10 @@ std::shared_ptr> TriangularMap::Sli std::vector>> components; // TODO: Handle empty case if( a < 0 || a >= b || b > this->outputDim ) - throw std::runtime_error("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); + throw std::invalid_argument("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); - int accum_a = 0; - int k_a = 0; + int accum_a = 0; // Accumulated output dimension before a + int k_a = 0; // Index of component containing a //TODO: Check that this is correct while(k_a < this->comps_.size()){ if(accum_a + this->comps_[k_a]->outputDim > a) @@ -362,8 +362,8 @@ std::shared_ptr> TriangularMap::Sli k_a++; } - int accum_b = accum_a; - int k_b = k_a; + int accum_b = accum_a; // Accumulated output dimension before b + int k_b = k_a; // Index of component containing b //TODO: Check that this is correct while(k_b < this->comps_.size()){ if(accum_b + this->comps_[k_b]->outputDim >= b) diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp index 9164f118..0b146ac2 100644 --- a/tests/Test_TriangularMap.cpp +++ b/tests/Test_TriangularMap.cpp @@ -902,22 +902,20 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef SECTION("Slice"){ int sliceBegin = 1; - int sliceEnd = 3; + int sliceEnd = 5; unsigned int extraInputs = 2; auto slice = triMap->Slice(sliceBegin, sliceEnd); - REQUIRE(slice->inputDim == triMap->inputDim); + REQUIRE(slice->inputDim == sliceEnd); REQUIRE(slice->outputDim == sliceEnd-sliceBegin); auto outSlice = slice->Evaluate(in); + auto outFull = triMap->Evaluate(in); // Just test the evaluation for(unsigned int i=sliceBegin; iEvaluate(Kokkos::subview(in, std::make_pair(0,int(i+1+extraInputs)), Kokkos::ALL())); - - REQUIRE(outBlock.extent(1)==numSamps); - REQUIRE(outBlock.extent(0)==1); - for(unsigned int j=0; j Date: Mon, 9 Jan 2023 10:31:01 -0500 Subject: [PATCH 15/16] Add crust cases --- src/TriangularMap.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp index 2e892667..f3126db6 100644 --- a/src/TriangularMap.cpp +++ b/src/TriangularMap.cpp @@ -349,8 +349,17 @@ template std::shared_ptr> TriangularMap::Slice(int a, int b) { std::vector>> components; // TODO: Handle empty case - if( a < 0 || a >= b || b > this->outputDim ) + if( a < 0 || a >= b || b > this->outputDim ) { throw std::invalid_argument("TriangularMap::Slice: 0 <= a < b <= outputDim must be satisfied."); + } + // Special cases if the slice is at the end of the map + if( b <= this->comps_[0]->outputDim) { + return this->comps_[0]->Slice(a, b); + } + if( a >= this->outputDim - this->comps_[this->comps_.size()-1]->outputDim) { + unsigned int rest_of_output = this->outputDim - this->comps_[this->comps_.size()-1]->outputDim; + return this->comps_[this->comps_.size()-1]->Slice(a - rest_of_output, b - rest_of_output); + } int accum_a = 0; // Accumulated output dimension before a int k_a = 0; // Index of component containing a @@ -372,8 +381,6 @@ std::shared_ptr> TriangularMap::Sli k_b++; } - std::cerr << "a = " << a << ", accum_a = " << accum_a << ", k_a = " << k_a << std::endl; - std::cerr << "b = " << b << ", accum_b = " << accum_b << ", k_b = " << k_b << std::endl; if(k_a == k_b){ components.push_back(this->comps_[k_a]->Slice(a-accum_a, b-accum_b)); } else { From 7c4c91b02fc9fbcad3406fd357a3132af7b68a5f Mon Sep 17 00:00:00 2001 From: Danny Sharp Date: Mon, 9 Jan 2023 11:35:08 -0500 Subject: [PATCH 16/16] Work on affinemap issues --- src/AffineMap.cpp | 2 +- tests/Test_AffineMap.cpp | 91 ++++++++++++++++++++++++++++++++-------- 2 files changed, 75 insertions(+), 18 deletions(-) diff --git a/src/AffineMap.cpp b/src/AffineMap.cpp index bca7fb07..d5427991 100644 --- a/src/AffineMap.cpp +++ b/src/AffineMap.cpp @@ -196,7 +196,7 @@ template std::shared_ptr> AffineMap::Slice(int a, int b) { StridedMatrix A_new; StridedVector b_new; - if(A_.extent(0) > 0) A_new = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(a,b)); + if(A_.extent(0) > 0) A_new = Kokkos::subview(A_, std::make_pair(a,b), Kokkos::ALL()); if(b_.size() > 0) b_new = Kokkos::subview(b_, std::make_pair(a,b)); return std::make_shared>(A_new, b_new); } diff --git a/tests/Test_AffineMap.cpp b/tests/Test_AffineMap.cpp index b6dbffa5..554af1c1 100644 --- a/tests/Test_AffineMap.cpp +++ b/tests/Test_AffineMap.cpp @@ -9,7 +9,7 @@ using namespace Catch; TEST_CASE( "Testing Shift-only AffineMap", "[ShiftMap]" ) { Kokkos::View b("b", 2); - + b(0) = 1.0; b(1) = 2.0; @@ -46,7 +46,7 @@ TEST_CASE( "Testing Shift-only AffineMap", "[ShiftMap]" ) { TEST_CASE( "Testing Linear-only AffineMap", "[LinearMap]" ) { Kokkos::View A("A", 2,2); - + A(0,0) = 2.0; A(1,0) = 1.0; A(0,1) = 1.0; @@ -85,8 +85,8 @@ TEST_CASE( "Testing Linear-only AffineMap", "[LinearMap]" ) { TEST_CASE( "Testing Rectangular AffineMap", "[RectangularLinearMap]" ) { Kokkos::View A("A", 2,3); - - /* + + /* A = [[2.0, 3.0, 1.0], [1.0, 1.0, 4.0]] */ @@ -96,7 +96,7 @@ TEST_CASE( "Testing Rectangular AffineMap", "[RectangularLinearMap]" ) { A(1,1) = 1.0; A(0,2) = 1.0; A(1,2) = 4.0; - + auto map = std::make_shared>(A); @@ -139,7 +139,7 @@ TEST_CASE( "Testing Full AffineMap", "[FullAffineMap]" ) { Kokkos::View A("A", 2,2); Kokkos::View b("b", 2); - + A(0,0) = 2.0; A(1,0) = 1.0; A(0,1) = 1.0; @@ -178,13 +178,70 @@ TEST_CASE( "Testing Full AffineMap", "[FullAffineMap]" ) { } +TEST_CASE( "Testing AffineMap::Slice", "[AffineMapSlice]" ) { + Kokkos::View A("A", 3,3); + Kokkos::View b("b", 3); + std::fill(A.data(), A.data()+A.extent(0)*A.extent(1), 0.0); + + A(0,0) = 2.0; + A(1,0) = 1.0; + A(0,1) = 1.0; + A(1,1) = 4.0; + A(0,2) = 1.0; + A(1,2) = 4.0; + + A(2,0) = 2.0; + A(2,1) = 1.0; + A(2,2) = 3.0; + + b(0) = 1.0; + b(1) = 2.0; + b(2) = 3.0; + + auto map = std::make_shared>(A,b); + auto slice = map->Slice(0,2); + REQUIRE(slice->inputDim==3); + REQUIRE(slice->outputDim==2); + + unsigned int numPts = 3; + Kokkos::View pts("Point", A.extent(1), numPts); + for(unsigned int i=0; i logDet = slice->LogDeterminant(pts); + for(unsigned int i=0; i evals = slice->Evaluate(pts); + REQUIRE(evals.extent(0)==2); + REQUIRE(evals.extent(1)==pts.extent(1)); + + Kokkos::View pts2 = slice->Inverse(pts, evals); + REQUIRE(pts2.extent(0)==2); + REQUIRE(pts2.extent(1)==pts.extent(1)); + + for(unsigned int i=0; i hb("b", 2); - + hb(0) = 1.0; hb(1) = 2.0; @@ -193,7 +250,7 @@ TEST_CASE( "Testing Shift-only AffineMap on Device", "[DeviceShiftMap]" ) { auto map = std::make_shared>(db); - + unsigned int numPts = 10; Kokkos::View hpts("Point", 2, numPts); for(unsigned int i=0; i hA("A", 2,2); - + hA(0,0) = 2.0; hA(1,0) = 1.0; hA(0,1) = 1.0; @@ -263,7 +320,7 @@ TEST_CASE( "Testing Linear-only AffineMap on Device", "[DeviceLinearMap]" ) { auto dpts2 = map->Inverse(dpts, devals); auto hpts2 = ToHost(dpts2); - + for(unsigned int i=0; i hA("A", 2,3); - - /* + + /* A = [[2.0, 3.0, 1.0], [1.0, 1.0, 4.0]] */ @@ -290,7 +347,7 @@ TEST_CASE( "Testing Rectangular AffineMap on Device", "[DeviceRectangularLinearM hA(1,1) = 1.0; hA(0,2) = 1.0; hA(1,2) = 4.0; - + auto dA = ToDevice(hA); auto map = std::make_shared>(dA); @@ -338,7 +395,7 @@ TEST_CASE( "Testing Full AffineMap on Device", "[DeviceFullAffineMap]" ) { Kokkos::View hA("A", 2,2); Kokkos::View hb("b", 2); - + hA(0,0) = 2.0; hA(1,0) = 1.0; hA(0,1) = 1.0; @@ -349,7 +406,7 @@ TEST_CASE( "Testing Full AffineMap on Device", "[DeviceFullAffineMap]" ) { auto dA = ToDevice(hA); auto db = ToDevice(hb); - + auto map = std::make_shared>(dA,db); unsigned int numPts = 10; @@ -374,7 +431,7 @@ TEST_CASE( "Testing Full AffineMap on Device", "[DeviceFullAffineMap]" ) { auto dpts2 = map->Inverse(dpts, devals); auto hpts2 = ToHost(dpts2); - + for(unsigned int i=0; i