diff --git a/CMakeLists.txt b/CMakeLists.txt index 28f1785e..f7edc515 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.13) -project(MParT VERSION 0.0.1) +project(MParT VERSION 2.0.2) message(STATUS "Will install MParT to ${CMAKE_INSTALL_PREFIX}") @@ -53,6 +53,7 @@ option(BUILD_SHARED_LIBS "Build using shared libraries" ON) # ############################################################# # RPATH settings +# See https://gist.github.com/kprussing/db21614ca5b51cedff07dfb70059f280 for scikit-build example # use, i.e. don't skip the full RPATH for the build tree set(CMAKE_SKIP_BUILD_RPATH FALSE) @@ -70,9 +71,9 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) # the RPATH to be used when installing, but only if it's not a system directory list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir) -if("${isSystemDir}" STREQUAL "-1") - set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") -endif("${isSystemDir}" STREQUAL "-1") +# if("${isSystemDir}" STREQUAL "-1") +# set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") +# endif("${isSystemDir}" STREQUAL "-1") # ############################################################# # Dependencies diff --git a/MParT/AffineFunction.h b/MParT/AffineFunction.h index a8da5f31..37096a82 100644 --- a/MParT/AffineFunction.h +++ b/MParT/AffineFunction.h @@ -51,9 +51,9 @@ class AffineFunction : public ParameterizedFunctionBase protected: - - StridedMatrix A_; - StridedVector b_; + + Kokkos::View A_; + Kokkos::View b_; int ldA; diff --git a/MParT/AffineMap.h b/MParT/AffineMap.h index e69bd131..0def7313 100644 --- a/MParT/AffineMap.h +++ b/MParT/AffineMap.h @@ -13,6 +13,7 @@ namespace mpart{ /** @brief Defines transformations of the form \f$Ax+b\f$ for an invertible matrix \f$A\f$ and vector offset \f$b\f$. + * Makes a deep copy of any views passed to the constructor. */ template class AffineMap : public ConditionalMapBase @@ -64,8 +65,8 @@ class AffineMap : public ConditionalMapBase protected: - StridedMatrix A_; - StridedVector b_; + Kokkos::View A_; + Kokkos::View b_; mpart::PartialPivLU luSolver_; double logDet_; diff --git a/MParT/MapObjective.h b/MParT/MapObjective.h index 012a5158..86ce36e3 100644 --- a/MParT/MapObjective.h +++ b/MParT/MapObjective.h @@ -63,8 +63,9 @@ class MapObjective { */ double operator()(unsigned int n, const double* x, double* grad, std::shared_ptr> map); - unsigned int Dim(){return train_.extent(0);} - unsigned int NumSamples(){return train_.extent(1);} + unsigned int InputDim() const {return train_.extent(0);} + virtual unsigned int MapOutputDim() const {return train_.extent(0);} + unsigned int NumSamples() const {return train_.extent(1);} /** * @brief Shortcut to calculate the error of the map on the training dataset @@ -176,7 +177,7 @@ class KLObjective: public MapObjective { double ObjectivePlusCoeffGradImpl(StridedMatrix data, StridedVector grad, std::shared_ptr> map) const override; double ObjectiveImpl(StridedMatrix data, std::shared_ptr> map) const override; void CoeffGradImpl(StridedMatrix data, StridedVector grad, std::shared_ptr> map) const override; - + unsigned int MapOutputDim() const override {return density_->Dim();} private: /** * @brief Density \f$\mu\f$ to calculate the KL with respect to (i.e. \f$D(\cdot||\mu)\f$ ) diff --git a/bindings/matlab/tests/TrainMapAdaptiveTest.m b/bindings/matlab/tests/TrainMapAdaptiveTest.m new file mode 100644 index 00000000..8b916c84 --- /dev/null +++ b/bindings/matlab/tests/TrainMapAdaptiveTest.m @@ -0,0 +1,44 @@ +classdef TrainMapTest < matlab.unittest.TestCase + % Copyright 2014 - 2016 The MathWorks, Inc. + + methods ( Test ) + function train(testCase) + % Create data + dim = 2; + N=20000; + N_test = N/5; + data = randn(dim+1,N); + target = [data(1,:);data(2,:);data(3,:) + data(2,:).^2]; + test = target(:,1:N_test); + train = target(:,N_test+1:end); + + % Create objective and map + obj1 = GaussianKLObjective(train, test, 1); + obj2 = GaussianKLObjective(train, test); + map_options = ATMOptions(); + map_options.maxDegrees = MultiIndex([3,5]) + msets2 = [CreateTotalOrder(1,1), CreateTotalOrder(2,1)] + msets1 = [CreateTotalOrder(2,1)] + + % Print test error before + map1 = TrainMapAdaptive(msets1, obj1, map_options); + map2 = TrainMapAdaptive(msets2, obj2, map_options); + + % Evaluate test samples after training + KS_stat1 = KSStatistic(map1,test); + KS_stat2 = KSStatistic(map2,test); + testCase.verifyTrue(KS_stat1 < 0.1); + testCase.verifyTrue(KS_stat2 < 0.1); + end + end +end + + +% Perform Kolmogorov-Smirnov test +function KS_stat = KSStatistic(map,samples) + pullback_evals = Evaluate(map,samples); + sorted_samples = sort(pullback_evals(:)); + samps_cdf = (1 + erf(sorted_samples/sqrt(2)))/2; + samps_ecdf = (1:numel(sorted_samples))'/numel(sorted_samples); + KS_stat = max(abs(samps_ecdf - samps_cdf)); +end \ No newline at end of file diff --git a/bindings/python/CMakeLists.txt b/bindings/python/CMakeLists.txt index e2d04495..0fff3980 100644 --- a/bindings/python/CMakeLists.txt +++ b/bindings/python/CMakeLists.txt @@ -38,4 +38,14 @@ target_link_libraries(pympart PRIVATE mpart Kokkos::kokkos Eigen3::Eigen ${EXT_L # Add an installation target for the python bindings install(TARGETS pympart DESTINATION "${PYTHON_INSTALL_PREFIX}") -install(DIRECTORY package/ DESTINATION "${PYTHON_INSTALL_PREFIX}") \ No newline at end of file +install(DIRECTORY package/ DESTINATION "${PYTHON_INSTALL_PREFIX}") + +# See https://gist.github.com/kprussing/db21614ca5b51cedff07dfb70059f280 +#set(lib_path "${PYTHON_PREFIX}/${CMAKE_INSTALL_LIBDIR}") +#message(STATUS "LIB_PATH = ${lib_path}") +#list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${lib_path}" is_system) +if(SKBUILD) + if (SKBUILD_LIB_RPATH) + set_target_properties(pympart PROPERTIES INSTALL_RPATH "${SKBUILD_LIB_RPATH}") + endif() +endif() \ No newline at end of file diff --git a/bindings/python/src/AffineMap.cpp b/bindings/python/src/AffineMap.cpp index 994b7ea2..94a6bcaf 100644 --- a/bindings/python/src/AffineMap.cpp +++ b/bindings/python/src/AffineMap.cpp @@ -13,15 +13,15 @@ using namespace mpart::binding; void mpart::binding::AffineMapWrapperHost(py::module &m) { py::class_, ConditionalMapBase, std::shared_ptr>>(m, "AffineMap") - .def(py::init( [](Eigen::Ref const& b) + .def(py::init( [](py::EigenDRef const& b) { return new AffineMap(VecToKokkos(b)); })) - .def(py::init( [](Eigen::Ref const& A, Eigen::Ref const& b) + .def(py::init( [](py::EigenDRef const& A, py::EigenDRef const& b) { return new AffineMap(MatToKokkos(A), VecToKokkos(b)); })) - .def(py::init( [](Eigen::Ref const& A) + .def(py::init( [](py::EigenDRef const& A) { return new AffineMap(MatToKokkos(A)); })); @@ -31,15 +31,15 @@ void mpart::binding::AffineMapWrapperHost(py::module &m) void mpart::binding::AffineFunctionWrapperHost(py::module &m) { py::class_, ParameterizedFunctionBase, std::shared_ptr>>(m, "AffineFunction") - .def(py::init( [](Eigen::Ref const& b) + .def(py::init( [](py::EigenDRef const& b) { return new AffineFunction(VecToKokkos(b)); })) - .def(py::init( [](Eigen::Ref const& A, Eigen::Ref const& b) + .def(py::init( [](py::EigenDRef const& A, py::EigenDRef const& b) { return new AffineFunction(MatToKokkos(A), VecToKokkos(b)); })) - .def(py::init( [](Eigen::Ref const& A) + .def(py::init( [](py::EigenDRef const& A) { return new AffineFunction(MatToKokkos(A)); })); diff --git a/pyproject.toml b/pyproject.toml index f8a6d9a6..679ba838 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ license={file="LICENSE.txt"} readme="README.md" requires-python = ">=3.7" description="A Monotone Parameterization Toolkit" -version="1.4.4" +version="2.0.2" keywords=["Measure Transport", "Monotone", "Transport Map", "Isotonic Regression", "Triangular", "Knothe-Rosenblatt"] [project.urls] diff --git a/setup.py b/setup.py index b81a2121..0c55b3d1 100644 --- a/setup.py +++ b/setup.py @@ -1,9 +1,41 @@ from skbuild import setup +import os, sys, site +import warnings + +def get_install_locations(): + """Return the installation directory, or '' if no directory could be found + + Adapted from stack overflow post https://stackoverflow.com/a/36205159 + """ + + if '--user' in sys.argv: + paths = (site.getusersitepackages(),) + else: + py_version = f'{sys.version_info[0]}.{sys.version_info[1]}' + paths = [ + sys.prefix + f'/lib/python{py_version}/dist-packages/', + sys.prefix + f'/lib/python{py_version}/site-packages/', + sys.prefix + f'/local/lib/python{py_version}/dist-packages/', + sys.prefix + f'/local/lib/python{py_version}/site-packages/', + f'/Library/Frameworks/Python.framework/Versions/{py_version}/lib/python{py_version}/site-packages/' + ] + + for path in paths: + if os.path.exists(path): + parts = path.split('/') + lib_indices = [index for index, item in enumerate(parts) if item == 'lib'] + return path, '/'.join(parts[0:(lib_indices[-1]+1)]) + + return '' + +site_folder, lib_folder = get_install_locations() + + setup( packages=['mpart'], package_dir={'mpart': 'bindings/python/package'}, package_data={'mpart':['**/*pympart*']}, include_package_data=True, - cmake_args=['-DKokkos_ENABLE_THREADS:BOOL=ON', '-DPYTHON_INSTALL_SUFFIX=bindings/python/package/', '-DMPART_JULIA:BOOL=OFF', '-DMPART_MATLAB:BOOL=OFF', '-DMPART_BUILD_TESTS:BOOL=OFF', '-DMPART_PYTHON:BOOL=ON', '-DPYTHON_INSTALL_PREFIX='] + cmake_args=['-DKokkos_ENABLE_THREADS:BOOL=ON', f'-DSKBUILD_LIB_RPATH={lib_folder}', f'-DSKBUILD_SITE_PATH={site_folder}', '-DPYTHON_INSTALL_SUFFIX=bindings/python/package/', '-DMPART_JULIA:BOOL=OFF', '-DMPART_MATLAB:BOOL=OFF', '-DMPART_BUILD_TESTS:BOOL=OFF', '-DMPART_PYTHON:BOOL=ON', '-DPYTHON_INSTALL_PREFIX='] ) \ No newline at end of file diff --git a/src/AffineFunction.cpp b/src/AffineFunction.cpp index 2f3bce01..1969ecd7 100644 --- a/src/AffineFunction.cpp +++ b/src/AffineFunction.cpp @@ -11,7 +11,7 @@ using namespace mpart; template AffineFunction::AffineFunction(StridedVector b) : ParameterizedFunctionBase(b.size(),b.size(),0), - b_("b",b.layout()) + b_("b",b.extent(0)) { Kokkos::deep_copy(b_, b); @@ -19,7 +19,7 @@ AffineFunction::AffineFunction(StridedVector b) template AffineFunction::AffineFunction(StridedMatrix A) : ParameterizedFunctionBase(A.extent(1),A.extent(0),0), - A_("A", A.layout()) + A_("A", A.extent(0), A.extent(1)) { Kokkos::deep_copy(A_, A); assert(A_.extent(0)<=A_.extent(1)); @@ -29,8 +29,8 @@ AffineFunction::AffineFunction(StridedMatrix A) template AffineFunction::AffineFunction(StridedMatrix A, StridedVector b) : ParameterizedFunctionBase(A.extent(1),A.extent(0),0), - A_("A", A.layout()), - b_("b",b.layout()) + A_("A", A.extent(0), A.extent(1)), + b_("b", b.extent(0)) { Kokkos::deep_copy(A_, A); Kokkos::deep_copy(b_, b); @@ -43,19 +43,26 @@ template void AffineFunction::EvaluateImpl(StridedMatrix const& pts, StridedMatrix output) { + unsigned int numPts = pts.extent(1); + Kokkos::MDRangePolicy, typename MemoryToExecution::Space> policy({{0, 0}}, {{numPts, this->outputDim}}); + // Linear part if(A_.extent(0)>0){ + + // Initialize output to zeros + Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(const int& j, const int& i) { + output(i,j) = 0.0; + }); + dgemm(1.0, A_, pts, 0.0, output); + }else{ Kokkos::deep_copy(output, pts); } - + // Bias part if(b_.size()>0){ - unsigned int numPts = pts.extent(1); - Kokkos::MDRangePolicy, typename MemoryToExecution::Space> policy({0, 0}, {numPts, this->outputDim}); - Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(const int& j, const int& i) { output(i,j) += b_(i); }); diff --git a/src/AffineMap.cpp b/src/AffineMap.cpp index 0c2d3fb0..a09e3d76 100644 --- a/src/AffineMap.cpp +++ b/src/AffineMap.cpp @@ -9,7 +9,7 @@ using namespace mpart; template AffineMap::AffineMap(StridedVector b) : ConditionalMapBase(b.size(),b.size(),0), - b_("b",b.layout()), + b_("b", b.extent(0)), logDet_(0.0) { Kokkos::deep_copy(b_, b); @@ -17,7 +17,7 @@ AffineMap::AffineMap(StridedVector b) : Conditi template AffineMap::AffineMap(StridedMatrix A) : ConditionalMapBase(A.extent(1),A.extent(0),0), - A_("A", A.layout()) + A_("A", A.extent(0), A.extent(1)) { Kokkos::deep_copy(A_, A); assert(A_.extent(0)<=A_.extent(1)); @@ -28,8 +28,8 @@ AffineMap::AffineMap(StridedMatrix A) : Conditi template AffineMap::AffineMap(StridedMatrix A, StridedVector b) : ConditionalMapBase(A.extent(1),A.extent(0),0), - A_("A", A.layout()), - b_("b",b.layout()) + A_("A", A.extent(0), A.extent(1)), + b_("b", b.extent(0)) { Kokkos::deep_copy(A_, A); Kokkos::deep_copy(b_, b); @@ -42,13 +42,6 @@ AffineMap::AffineMap(StridedMatrix A, template void AffineMap::Factorize(){ - // If A_ is not column major, create a column major version - if(A_.stride_0()!=1){ - Kokkos::View anew("A_", A_.extent(0), A_.extent(1)); - Kokkos::deep_copy(anew, A_); - A_ = anew; - } - if(A_.extent(0)!=A_.extent(1)){ StridedMatrix Asub = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(A_.extent(1)-A_.extent(0),A_.extent(1))); luSolver_.compute(Asub); @@ -145,19 +138,25 @@ template void AffineMap::EvaluateImpl(StridedMatrix const& pts, StridedMatrix output) { + unsigned int numPts = pts.extent(1); + Kokkos::MDRangePolicy, typename MemoryToExecution::Space> policy({{0, 0}}, {{numPts, this->outputDim}}); + // Linear part if(A_.extent(0)>0){ + + // Initialize output to zeros + Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(const int& j, const int& i) { + output(i,j) = 0.0; + }); + dgemm(1.0, A_, pts, 0.0, output); + }else{ Kokkos::deep_copy(output, pts); } - + // Bias part if(b_.size()>0){ - - unsigned int numPts = pts.extent(1); - Kokkos::MDRangePolicy, typename MemoryToExecution::Space> policy({{0, 0}}, {{numPts, this->outputDim}}); - Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(const int& j, const int& i) { output(i,j) += b_(i); }); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1b04ef93..adda6569 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -50,6 +50,9 @@ target_sources(mpart MapFactoryImpl13.cpp MapFactoryImpl14.cpp MapFactoryImpl15.cpp + MapFactoryImpl16.cpp + MapFactoryImpl17.cpp + MapFactoryImpl18.cpp ${MPART_OPT_FILES} Initialization.cpp diff --git a/src/MapFactoryImpl16.cpp b/src/MapFactoryImpl16.cpp new file mode 100644 index 00000000..6a195b25 --- /dev/null +++ b/src/MapFactoryImpl16.cpp @@ -0,0 +1,35 @@ +#include "MParT/MapFactory.h" + +#include "MParT/MonotoneComponent.h" +#include "MParT/TriangularMap.h" +#include "MParT/Quadrature.h" + +#include "MParT/HermiteFunction.h" +#include "MParT/MultivariateExpansionWorker.h" +#include "MParT/PositiveBijectors.h" + +#include "MParT/LinearizedBasis.h" + +using namespace mpart; + +template +std::shared_ptr> CreateComponentImpl_LinHF_AS(FixedMultiIndexSet const& mset, MapOptions opts) +{ + LinearizedBasis basis1d(HermiteFunction(), opts.basisLB, opts.basisUB); + AdaptiveSimpson quad(opts.quadMaxSub, 1, nullptr, opts.quadAbsTol, opts.quadRelTol, QuadError::First, opts.quadMinSub); + + MultivariateExpansionWorker expansion(mset, basis1d); + std::shared_ptr> output; + + output = std::make_shared>(expansion, quad, opts.contDeriv); + + output->SetCoeffs(Kokkos::View("Component Coefficients", mset.Size())); + return output; +} + +static auto reg_host_linhf_as_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::AdaptiveSimpson), CreateComponentImpl_LinHF_AS)); +static auto reg_host_linhf_as_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveSimpson), CreateComponentImpl_LinHF_AS)); +#if defined(MPART_ENABLE_GPU) + static auto reg_device_linhf_as_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::AdaptiveSimpson), CreateComponentImpl_LinHF_AS)); + static auto reg_device_linhf_as_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveSimpson), CreateComponentImpl_LinHF_AS)); +#endif \ No newline at end of file diff --git a/src/MapFactoryImpl17.cpp b/src/MapFactoryImpl17.cpp new file mode 100644 index 00000000..2207b57f --- /dev/null +++ b/src/MapFactoryImpl17.cpp @@ -0,0 +1,35 @@ +#include "MParT/MapFactory.h" + +#include "MParT/MonotoneComponent.h" +#include "MParT/TriangularMap.h" +#include "MParT/Quadrature.h" + +#include "MParT/HermiteFunction.h" +#include "MParT/MultivariateExpansionWorker.h" +#include "MParT/PositiveBijectors.h" + +#include "MParT/LinearizedBasis.h" + +using namespace mpart; + +template +std::shared_ptr> CreateComponentImpl_LinHF_CC(FixedMultiIndexSet const& mset, MapOptions opts) +{ + LinearizedBasis basis1d(HermiteFunction(), opts.basisLB, opts.basisUB); + ClenshawCurtisQuadrature quad(opts.quadPts, 1); + + MultivariateExpansionWorker expansion(mset, basis1d); + std::shared_ptr> output; + + output = std::make_shared>(expansion, quad, opts.contDeriv); + + output->SetCoeffs(Kokkos::View("Component Coefficients", mset.Size())); + return output; +} + +static auto reg_host_linhf_cc_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::ClenshawCurtis), CreateComponentImpl_LinHF_CC)); +static auto reg_host_linhf_cc_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::ClenshawCurtis), CreateComponentImpl_LinHF_CC)); +#if defined(MPART_ENABLE_GPU) + static auto reg_device_linhf_cc_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::ClenshawCurtis), CreateComponentImpl_LinHF_CC)); + static auto reg_device_linhf_cc_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::ClenshawCurtis), CreateComponentImpl_LinHF_CC)); +#endif \ No newline at end of file diff --git a/src/MapFactoryImpl18.cpp b/src/MapFactoryImpl18.cpp new file mode 100644 index 00000000..34b5959f --- /dev/null +++ b/src/MapFactoryImpl18.cpp @@ -0,0 +1,36 @@ +#include "MParT/MapFactory.h" + +#include "MParT/MonotoneComponent.h" +#include "MParT/TriangularMap.h" +#include "MParT/Quadrature.h" + +#include "MParT/HermiteFunction.h" +#include "MParT/MultivariateExpansionWorker.h" +#include "MParT/PositiveBijectors.h" + +#include "MParT/LinearizedBasis.h" + +using namespace mpart; + +template +std::shared_ptr> CreateComponentImpl_LinHF_ACC(FixedMultiIndexSet const& mset, MapOptions opts) +{ + LinearizedBasis basis1d(HermiteFunction(), opts.basisLB, opts.basisUB); + unsigned int level = std::log2(opts.quadPts-2); + AdaptiveClenshawCurtis quad(level, opts.quadMaxSub, 1, nullptr, opts.quadAbsTol, opts.quadRelTol, QuadError::First, opts.quadMinSub); + + MultivariateExpansionWorker expansion(mset, basis1d); + std::shared_ptr> output; + + output = std::make_shared>(expansion, quad, opts.contDeriv); + + output->SetCoeffs(Kokkos::View("Component Coefficients", mset.Size())); + return output; +} + +static auto reg_host_linhf_acc_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_LinHF_ACC)); +static auto reg_host_linhf_acc_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_LinHF_ACC)); +#if defined(MPART_ENABLE_GPU) + static auto reg_device_linhf_acc_exp = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_LinHF_ACC)); + static auto reg_device_linhf_acc_splus = mpart::MapFactory::CompFactoryImpl::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::HermiteFunctions, true, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_LinHF_ACC)); +#endif \ No newline at end of file diff --git a/src/MapObjective.cpp b/src/MapObjective.cpp index f4e3ce95..d11123d9 100644 --- a/src/MapObjective.cpp +++ b/src/MapObjective.cpp @@ -61,7 +61,6 @@ double KLObjective::ObjectivePlusCoeffGradImpl(StridedMatrix rc(densityGradX, -1.0/((double) N_samps)); Kokkos::parallel_reduce(N_samps, rc, &grad(0)); - Kokkos::fence("End of MapObjective"); return sumDensity/N_samps; } diff --git a/src/MultiIndices/MultiIndexSet.cpp b/src/MultiIndices/MultiIndexSet.cpp index c178c94b..184b6e08 100644 --- a/src/MultiIndices/MultiIndexSet.cpp +++ b/src/MultiIndices/MultiIndexSet.cpp @@ -658,7 +658,7 @@ std::vector MultiIndexSet::Expand(unsigned int activeIndex) // loop through the forward neighbors of this index std::set tempSet = outEdges.at(globalIndex); for(int neighbor : tempSet) - { + { if(IsAdmissible(neighbor)&&(!IsActive(neighbor))){ Activate(neighbor); newIndices.push_back(global2active.at(neighbor)); @@ -672,13 +672,31 @@ std::vector MultiIndexSet::Expand(unsigned int activeIndex) std::vector MultiIndexSet::Expand() { - std::vector frontier = Frontier(); - std::vector newInds, allNewInds; - for(auto& ind : frontier){ - newInds = Expand(ind); - allNewInds.insert(allNewInds.end(), newInds.begin(), newInds.end()); - } - return allNewInds; + // Figure out which terms in the frontier should be activated + std::vector frontier = Frontier(); + std::vector toExpand; // Global indices to activate + for(auto& ind : frontier){ + + unsigned int globalIndex = active2global.at(ind); + + // loop through the forward neighbors of this index + std::set tempSet = outEdges.at(globalIndex); + for(int neighbor : tempSet) + { + if(IsAdmissible(neighbor)&&(!IsActive(neighbor))){ + toExpand.push_back(neighbor); + } + } + } + + // Activate those terms + std::vector allNewInds; + for(auto& neighbor : toExpand){ + Activate(neighbor); + allNewInds.push_back(global2active.at(neighbor)); + } + + return allNewInds; } std::vector MultiIndexSet::ForciblyExpand(unsigned int const activeIndex) diff --git a/src/TrainMapAdaptive.cpp b/src/TrainMapAdaptive.cpp index 9e2b9d68..7091db57 100644 --- a/src/TrainMapAdaptive.cpp +++ b/src/TrainMapAdaptive.cpp @@ -46,8 +46,8 @@ std::shared_ptr> mpart::TrainMapAdaptive(s ATMOptions options) { // Dimensions - unsigned int inputDim = objective->Dim(); - unsigned int outputDim = mset0.size(); + unsigned int inputDim = objective->InputDim(); + unsigned int outputDim = objective->MapOutputDim(); std::vector mset_sizes (outputDim); unsigned int currSz = 0; @@ -60,6 +60,13 @@ std::shared_ptr> mpart::TrainMapAdaptive(s std::vector mset_tmp {}; std::vector mset_best {}; + if(mset0.size() != outputDim) { + std::stringstream ss; + ss << "AdaptiveTransportMap: Number of MultiIndexSets must match dimension of Objective.\n"; + ss << "Expected Length " << inputDim << ", got " << mset0.size() << "."; + throw std::invalid_argument(ss.str()); + } + // Ensure the given vector of multiindexsets is valid unsigned int currMsetDim = mset0[0].Length(); for(unsigned int j = 0; j < outputDim; j++) { diff --git a/tests/MultiIndices/Test_MultiIndexSet.cpp b/tests/MultiIndices/Test_MultiIndexSet.cpp index 9f5fbe97..417e6adb 100644 --- a/tests/MultiIndices/Test_MultiIndexSet.cpp +++ b/tests/MultiIndices/Test_MultiIndexSet.cpp @@ -356,6 +356,16 @@ TEST_CASE("Testing the MultiIndexSet class", "[MultiIndexSet]" ) { // Check the result of IsAdmissable(). REQUIRE( indexFamily.IsAdmissible(multi)); + + // Check to make sure https://github.com/MeasureTransport/MParT/issues/308 is resolved + indexFamily = MultiIndexSet(2); + indexFamily += MultiIndex{0,0}; + indexFamily += MultiIndex{1,0}; + indexFamily.Expand(); + + REQUIRE( indexFamily.IsActive(MultiIndex{2,0}) ); + REQUIRE( indexFamily.IsActive(MultiIndex{0,1}) ); + REQUIRE(! indexFamily.IsActive(MultiIndex{1,1}) ); } /* diff --git a/tests/Test_TrainMapAdaptive.cpp b/tests/Test_TrainMapAdaptive.cpp index f09c7874..d9d061e3 100644 --- a/tests/Test_TrainMapAdaptive.cpp +++ b/tests/Test_TrainMapAdaptive.cpp @@ -140,4 +140,35 @@ TEST_CASE("Adaptive Transport Map","[ATM]") { StridedMatrix pullback_test = atm->Evaluate(testSamples); TestStandardNormalSamples(pullback_test); } + SECTION("TraditionalBananaOneComp") { + Kokkos::View targetSamples("targetSamples", 2, numPts); + Kokkos::parallel_for("Intializing targetSamples", numPts, KOKKOS_LAMBDA(const unsigned int i){ + targetSamples(0,i) = samples(0,i); + targetSamples(1,i) = samples(1,i) + samples(0,i)*samples(0,i); + }); + NormalizeSamples(targetSamples); + + StridedMatrix testSamples = Kokkos::subview(targetSamples, Kokkos::ALL, Kokkos::pair(0, testPts)); + StridedMatrix trainSamples = Kokkos::subview(targetSamples, Kokkos::ALL, Kokkos::pair(testPts, numPts)); + auto objective = ObjectiveFactory::CreateGaussianKLObjective(trainSamples,testSamples,1); + + std::vector mset0 {MultiIndexSet::CreateTotalOrder(2,0)}; + MultiIndexSet correctMset = (MultiIndexSet::CreateTotalOrder(2,0) + MultiIndex{0,1}) + MultiIndex{2,0}; + + ATMOptions opts; + opts.maxSize = 15; // Algorithm must add 3 correct terms in 6 iterations + opts.basisLB = -3.; + opts.basisUB = 3.; + opts.maxDegrees = MultiIndex{1000,4}; // Limit the second input to have cubic complexity or less + + std::shared_ptr> atm = TrainMapAdaptive(mset0, objective, opts); + MultiIndexSet finalMset = mset0[0]; + CHECK((finalMset + correctMset).Size() == finalMset.Size()); + std::vector bounded = finalMset.FilterBounded(opts.maxDegrees); + bool checkBound = false; + for(auto b1 : bounded) checkBound |= b1; + CHECK(!checkBound); + StridedMatrix pullback_test = atm->Evaluate(testSamples); + TestStandardNormalSamples(pullback_test); + } } \ No newline at end of file