Compact monotone component

MParT/BasisEvaluator.h

@@ -234,7 +234,7 @@ class BasisEvaluator<BasisHomogeneity::OffdiagHomogeneous,
    * @param offdiag Evaluator for offdiagonal input elements
    * @param diag Evaluator for diagonal input element
    */
-  BasisEvaluator(unsigned int dim, OffdiagEvaluatorType const &offdiag = OffdiagEvaluatorType(),
+  BasisEvaluator(unsigned int dim = 0, OffdiagEvaluatorType const &offdiag = OffdiagEvaluatorType(),
                  DiagEvaluatorType const &diag = DiagEvaluatorType())
       : offdiag_(offdiag), diag_(diag), dim_(dim) {}
 

MParT/MapFactory.h

@@ -217,14 +217,15 @@ namespace mpart{
         */
         template<typename MemorySpace>
         struct CompFactoryImpl{
-            typedef std::tuple<BasisTypes, bool, PosFuncTypes, QuadTypes> OptionsKeyType;
+            typedef std::tuple<BasisTypes, bool, PosFuncTypes, QuadTypes, bool> OptionsKeyType;
             typedef std::function<std::shared_ptr<ConditionalMapBase<MemorySpace>>(FixedMultiIndexSet<MemorySpace> const&, MapOptions options)> FactoryFunctionType;
             typedef std::map<OptionsKeyType, FactoryFunctionType> FactoryMapType;
 
             static FactoryFunctionType GetFactoryFunction(MapOptions opts)
             {
                 bool isLinearized = (!isinf(opts.basisLB)) ||(!isinf(opts.basisUB));
-                OptionsKeyType optionsKey(opts.basisType, isLinearized, opts.posFuncType, opts.quadType);
+                bool isCompact = opts.isCompact;
+                OptionsKeyType optionsKey(opts.basisType, isLinearized, opts.posFuncType, opts.quadType, isCompact);
 
                 auto factoryMap = GetFactoryMap();
 

MParT/MapOptions.h

@@ -12,7 +12,8 @@ namespace mpart{
     {
         ProbabilistHermite,
         PhysicistHermite,
-        HermiteFunctions
+        HermiteFunctions,
+        Legendre
     };
 
     enum class PosFuncTypes
@@ -76,6 +77,8 @@ namespace mpart{
         double basisLB = -std::numeric_limits<double>::infinity();
         double basisUB =  std::numeric_limits<double>::infinity();
 
+        /** Whether the map should be compactly supported (cannot be linearized, only works with certain other map options */
+        bool isCompact = false;
 
         /** The type of positive bijector used inside the monotonicity-guaranteeing integral
             formulation.

MParT/MultivariateExpansionWorker.h

@@ -129,6 +129,27 @@ class MultivariateExpansionWorker
     */
     KOKKOS_INLINE_FUNCTION unsigned int InputSize() const {return multiSet_.Length();};
 
+    /**
+     @brief Fills cache to evaluate f(0)
+     @details
+     @param polyCache A pointer to the start of the cache.  This memory must be allocated before calling this function.
+     @param derivType
+
+     @see FillCache1, FillCache2
+     */
+    KOKKOS_FUNCTION void FillCache0(double*          polyCache,
+                                    DerivativeFlags::DerivativeType derivType) const
+    {
+        // Only allowed to evaluate, fail if you want derivatives; while possible this is a safety precaution
+        if(derivType == DerivativeFlags::None || derivType == DerivativeFlags::Parameters){
+            for(unsigned int d=0; d<dim_; ++d) {
+                basis1d_.EvaluateAll(d, &polyCache[startPos_(d)], maxDegrees_(d), 0.);
+            }
+        } else {
+            ProcAgnosticError<std::invalid_argument>("Cannot get derivatives for FillCache0");
+        }
+    }
+
     /**
      @brief Precomputes parts of the cache using all but the last component of the point, i.e., using only \f$x_1,x_2,\ldots,x_{d-1}\f$, not \f$x_d\f$.
      @details

MParT/OrthogonalPolynomial.h

@@ -1,5 +1,5 @@
-#ifndef ORTHOGONALPOLYNOMIAL_H
-#define ORTHOGONALPOLYNOMIAL_H
+#ifndef MPART_ORTHOGONALPOLYNOMIAL_H
+#define MPART_ORTHOGONALPOLYNOMIAL_H
 
 #include <Kokkos_Core.hpp>
 #include <cmath>
@@ -327,6 +327,21 @@ class PhysicistHermiteMixer{
 typedef OrthogonalPolynomial<PhysicistHermiteMixer> PhysicistHermite;
 
 
+class ShiftedLegendreMixer{ // Polynomials orthogonal on U[0,1]
+public:
+    KOKKOS_INLINE_FUNCTION double Normalization(unsigned int polyOrder) const { return 1./(2*polyOrder+1);}
+protected:
+    // \f[ p_{k}(x) = (a_k x + b_k) p_{k-1}(x) - c_k p_{k-2}(x) \f]
+    KOKKOS_INLINE_FUNCTION double ak(unsigned int k) const {return  2. * (2. * k - 1)/k;}
+    KOKKOS_INLINE_FUNCTION double bk(unsigned int k) const {return -(2. * k - 1)/k;}
+    KOKKOS_INLINE_FUNCTION double ck(unsigned int k) const {return (k-1.)/k;}
+    KOKKOS_INLINE_FUNCTION double phi0(double) const {return 1.0;}
+    KOKKOS_INLINE_FUNCTION double phi1(double x) const {return 2*x-1;}
+    KOKKOS_INLINE_FUNCTION double phi1_deriv(double) const{return 2.0;};
+};
+
+using ShiftedLegendre = OrthogonalPolynomial<ShiftedLegendreMixer>;
+
 } // namespace mpart
 
 #endif

MParT/UnivariateBases.h

@@ -0,0 +1,155 @@
+#ifndef MPART_UNIVARIATEBASES_H
+#define MPART_UNIVARIATEBASES_H
+
+#include <Kokkos_Core.hpp>
+#include <cmath>
+
+#include <iostream>
+
+#include "MParT/Utilities/MathFunctions.h"
+
+namespace mpart{
+
+/**
+ * @brief Generic class to represent functions
+ * \f[p_0(x)\equiv 1,p_k(x)=sin(2\pi k x)\f]
+ */
+class SineBasis
+{
+public:
+    static constexpr double PI2 = 2*M_PI;
+
+    /* Evaluates all polynomials up to a specified order. */
+    KOKKOS_FUNCTION void EvaluateAll(double*              output,
+                                     unsigned int         maxOrder,
+                                     double               x) const
+    {
+        output[0] = 1.;
+
+        double sin_x, cos_x, cos_kx;
+
+        if(maxOrder>0) {
+            sin_x = sin(PI2*x);
+            cos_x = cos(PI2*x);
+            cos_kx = cos_x;
+            output[1] = sin_x;
+        }
+
+        for(unsigned int order=2; order<=maxOrder; ++order) {
+            output[order] = output[order-1]*cos_x + cos_kx*sin_x;
+            cos_kx = cos_kx*cos_x - output[order-1]*sin_x;
+        }
+    }
+
+    /** Evaluates the derivative of every polynomial in this family up to degree maxOrder (inclusive).
+        The results are stored in the memory pointed to by the derivs pointer.
+    */
+    KOKKOS_FUNCTION void EvaluateDerivatives(double*      derivs,
+                             unsigned int maxOrder,
+                             double       x) const
+    {
+        derivs[0] = 0.;
+
+        double sin_x, cos_x, sin_kx;
+
+        if(maxOrder>0) {
+            sin_x = sin(PI2*x);
+            cos_x = cos(PI2*x);
+            sin_kx = sin_x;
+            derivs[1] = cos_x;
+        }
+
+        // d_x sin(2pi k x) = 2pi k cos(2pi k x)
+        for(unsigned int order=2; order<=maxOrder; ++order) {
+            derivs[order] = derivs[order-1]*cos_x - sin_kx*sin_x;
+            sin_kx = sin_kx*cos_x + derivs[order-1]*sin_x;
+            derivs[order-2] *= PI2*(order-2);
+        }
+        if(maxOrder>0) derivs[maxOrder-1] *= PI2*(maxOrder-1);
+        derivs[maxOrder] *= PI2*(maxOrder);
+    }
+
+    /** Evaluates the value and derivative of every polynomial in this family up to degree maxOrder (inclusive).
+        The results are stored in the memory pointed to by the derivs pointer.
+    */
+    KOKKOS_FUNCTION void EvaluateDerivatives(double*      vals,
+                           double*      derivs,
+                           unsigned int maxOrder,
+                           double       x) const
+    {
+        vals[0] = 1.;
+        derivs[0] = 0.;
+
+        double sin_x, cos_x;
+
+        if(maxOrder>0) {
+            sin_x = sin(PI2*x);
+            cos_x = cos(PI2*x);
+            vals[1] = sin_x;
+            derivs[1] = cos_x;
+        }
+
+        for(unsigned int order=2; order<=maxOrder; ++order) {
+            vals[order] = vals[order-1]*cos_x + derivs[order-1]*sin_x;
+            derivs[order] = derivs[order-1]*cos_x - vals[order-1]*sin_x;
+            derivs[order-2] *= PI2*(order-2);
+        }
+        if(maxOrder>0) derivs[maxOrder-1] *= PI2*(maxOrder-1);
+        derivs[maxOrder] *= PI2*(maxOrder);
+    }
+
+    KOKKOS_FUNCTION void EvaluateSecondDerivatives(double*      vals,
+                                   double*      derivs,
+                                   double*      secondDerivs,
+                                   unsigned int maxOrder,
+                                   double       x) const
+    {
+        vals[0] = 1.;
+        derivs[0] = 0.;
+        secondDerivs[0] = 0.;
+
+        double sin_x, cos_x;
+
+        if(maxOrder>0) {
+            sin_x = sin(PI2*x);
+            cos_x = cos(PI2*x);
+            vals[1] = sin_x;
+            derivs[1] = cos_x;
+            secondDerivs[1] = -PI2*PI2*vals[1];
+        }
+
+        for(unsigned int order=2; order<=maxOrder; ++order) {
+            vals[order] = vals[order-1]*cos_x + derivs[order-1]*sin_x;
+            derivs[order] = derivs[order-1]*cos_x - vals[order-1]*sin_x;
+            derivs[order-2] *= PI2*(order-2);
+            secondDerivs[order] = -(PI2*order)*(PI2*order)*vals[order];
+        }
+        if(maxOrder > 0) derivs[maxOrder-1] *= PI2*(maxOrder-1);
+        derivs[maxOrder] *= PI2*(maxOrder);
+    }
+
+
+
+    KOKKOS_FUNCTION double Evaluate(unsigned int const order,
+                    double const x) const
+    {
+        return order == 0 ? 1. : sin(PI2*order*x);
+    }
+
+    KOKKOS_FUNCTION double Derivative(unsigned int const order,
+                      double const x) const
+    {
+        return order == 0 ? 0. : PI2*order*cos(PI2*order*x);
+    }
+
+    KOKKOS_FUNCTION double SecondDerivative(unsigned int const order,
+                      double const x) const
+    {
+        return -(PI2*order)*(PI2*order)*sin(PI2*order*x);
+    }
+
+};
+
+} // namespace mpart
+
+#endif

MParT/Utilities/Serialization.h

@@ -11,8 +11,10 @@
 
 // A macro that can be used for registering the various MonotoneComponent classes with CEREAL
 // This macro is used in the MapFactoryImpl*.cpp files 
-#define REGISTER_MONO_COMP(BASIS_HOMOGENEITY, BASIS_TYPE, POS_TYPE, QUAD_TYPE, MEMORY_SPACE) \
-    CEREAL_REGISTER_TYPE(mpart::MonotoneComponent<mpart::MultivariateExpansionWorker<mpart::BasisEvaluator<BASIS_HOMOGENEITY, BASIS_TYPE>, MEMORY_SPACE>, mpart::POS_TYPE, mpart::QUAD_TYPE<MEMORY_SPACE>, MEMORY_SPACE>)
+#define REGISTER_HOMOGENEOUS_MONO_COMP(BASIS_TYPE, POS_TYPE, QUAD_TYPE, MEMORY_SPACE, IS_COMPACT) \
+    CEREAL_REGISTER_TYPE(mpart::MonotoneComponent<mpart::MultivariateExpansionWorker<mpart::BasisEvaluator<BasisHomogeneity::Homogeneous, BASIS_TYPE>, MEMORY_SPACE>, mpart::POS_TYPE, mpart::QUAD_TYPE<MEMORY_SPACE>, MEMORY_SPACE, IS_COMPACT>)
+#define REGISTER_OFFDIAGHOMOGENEOUS_MONO_COMP(BASIS_TYPE_1, BASIS_TYPE_2, RECTIFIER, POS_TYPE, QUAD_TYPE, MEMORY_SPACE, IS_COMPACT) \
+    CEREAL_REGISTER_TYPE(mpart::MonotoneComponent<mpart::MultivariateExpansionWorker<mpart::BasisEvaluator<BasisHomogeneity::OffdiagHomogeneous, Kokkos::pair<BASIS_TYPE_1, BASIS_TYPE_2>, RECTIFIER>, MEMORY_SPACE>, mpart::POS_TYPE, mpart::QUAD_TYPE<MEMORY_SPACE>, MEMORY_SPACE, IS_COMPACT>)
 
 
 namespace cereal {

src/CMakeLists.txt

@@ -54,6 +54,7 @@ target_sources(mpart
     MapFactoryImpl16.cpp
     MapFactoryImpl17.cpp
     MapFactoryImpl18.cpp
+    MapFactoryImpl19.cpp
 
     ${MPART_OPT_FILES}
     Initialization.cpp

src/MapFactoryImpl1.cpp

@@ -13,7 +13,7 @@
 using namespace mpart;
 
 
-template<typename MemorySpace, typename PosFuncType>
+template<typename MemorySpace, typename PosFuncType, bool isCompact>
 std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateComponentImpl_Phys_ACC(FixedMultiIndexSet<MemorySpace> const& mset, MapOptions opts)
 {
     BasisEvaluator<BasisHomogeneity::Homogeneous,PhysicistHermite> basis1d(opts.basisNorm);
@@ -24,27 +24,38 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateComponentImpl_Phys_ACC(Fi
     MultivariateExpansionWorker<BasisEvaluator<BasisHomogeneity::Homogeneous,PhysicistHermite>,MemorySpace> expansion(mset, basis1d);
     std::shared_ptr<ConditionalMapBase<MemorySpace>> output;
 
-    output = std::make_shared<MonotoneComponent<decltype(expansion), PosFuncType, decltype(quad), MemorySpace>>(expansion, quad, opts.contDeriv, opts.nugget);
+    output = std::make_shared<MonotoneComponent<decltype(expansion), PosFuncType, decltype(quad), MemorySpace, isCompact>>(expansion, quad, opts.contDeriv, opts.nugget);
 
     Kokkos::View<const double*, MemorySpace> coeffs = Kokkos::View<double*,MemorySpace>("Component Coefficients", mset.Size());
     output->SetCoeffs(coeffs);
     return output;
 }
 
-static auto reg_host_phys_acc_exp = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, Exp>));
-static auto reg_host_phys_acc_splus = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, SoftPlus>));
+static auto reg_host_phys_acc_exp = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis, false), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, Exp, false>));
+static auto reg_host_phys_acc_splus = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis, false), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, SoftPlus, false>));
 #if defined(MPART_ENABLE_GPU)
-    static auto reg_device_phys_acc_exp = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, Exp>));
-    static auto reg_device_phys_acc_splus = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, SoftPlus>));
+    static auto reg_device_phys_acc_exp = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis, false), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, Exp, false>));
+    static auto reg_device_phys_acc_splus = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis, false), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, SoftPlus, false>));
+#endif
+
+static auto reg_host_phys_acc_exp_compact = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis, true), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, Exp, true>));
+static auto reg_host_phys_acc_splus_compact = mpart::MapFactory::CompFactoryImpl<Kokkos::HostSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis, true), CreateComponentImpl_Phys_ACC<Kokkos::HostSpace, SoftPlus, true>));
+#if defined(MPART_ENABLE_GPU)
+    static auto reg_device_phys_acc_exp_compact = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::Exp, QuadTypes::AdaptiveClenshawCurtis, true), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, Exp, true>));
+    static auto reg_device_phys_acc_splus_compact = mpart::MapFactory::CompFactoryImpl<mpart::DeviceSpace>::GetFactoryMap()->insert(std::make_pair(std::make_tuple(BasisTypes::PhysicistHermite, false, PosFuncTypes::SoftPlus, QuadTypes::AdaptiveClenshawCurtis, true), CreateComponentImpl_Phys_ACC<mpart::DeviceSpace, SoftPlus, true>));
 #endif
 
 #if defined(MPART_HAS_CEREAL)
-REGISTER_MONO_COMP(BasisHomogeneity::Homogeneous, PhysicistHermite, Exp, AdaptiveClenshawCurtis, Kokkos::HostSpace)
-REGISTER_MONO_COMP(BasisHomogeneity::Homogeneous, PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, Kokkos::HostSpace)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, Exp, AdaptiveClenshawCurtis, Kokkos::HostSpace, false)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, Kokkos::HostSpace, false)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, Exp, AdaptiveClenshawCurtis, Kokkos::HostSpace, true)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, Kokkos::HostSpace, true)
 #if defined(MPART_ENABLE_GPU)
-REGISTER_MONO_COMP(BasisHomogeneity::Homogeneous, PhysicistHermite, Exp, AdaptiveClenshawCurtis, mpart::DeviceSpace)
-REGISTER_MONO_COMP(BasisHomogeneity::Homogeneous, PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, mpart::DeviceSpace)
-#endif 
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, Exp, AdaptiveClenshawCurtis, mpart::DeviceSpace, false)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, mpart::DeviceSpace, false)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, Exp, AdaptiveClenshawCurtis, mpart::DeviceSpace, true)
+REGISTER_HOMOGENEOUS_MONO_COMP(PhysicistHermite, SoftPlus, AdaptiveClenshawCurtis, mpart::DeviceSpace, true)
+#endif
 
 CEREAL_REGISTER_DYNAMIC_INIT(mpartInitMapFactory1)
 #endif