Merge pull request MeasureTransport#386 from dannys4/newBindings

Bindings for multi-index methods, factory methods and bindings for sigmoids, new 1d map

.github/workflows/build-bindings.yml

@@ -61,15 +61,15 @@ jobs:
         id: cache-kokkos
         with:
           path: "${{ github.workspace }}/KOKKOS_INSTALL"
-          key: kokkos3.7.00
+          key: kokkos4.2.00
 
       - if: ${{steps.cache-kokkos.outputs.cache-hit != 'true'}}
         name: Checkout Kokkos
         uses: actions/checkout@v3
         with:
           repository: kokkos/kokkos
           path: kokkos
-          ref: '3.7.00'
+          ref: '4.2.00'
 
       - if: ${{steps.cache-kokkos.outputs.cache-hit != 'true'}}
         name: Install Kokkos

.github/workflows/build-external-lib-tests.yml

@@ -44,15 +44,15 @@ jobs:
         id: cache-kokkos
         with:
           path: "${{ github.workspace }}/KOKKOS_INSTALL"
-          key: kokkos3.7.00
+          key: kokkos4.2.00
 
       - if: ${{steps.cache-kokkos.outputs.cache-hit != 'true'}}
         name: Checkout Kokkos
         uses: actions/checkout@v3
         with:
           repository: kokkos/kokkos
           path: kokkos
-          ref: '3.7.00'
+          ref: '4.2.00'
 
       - if: ${{steps.cache-kokkos.outputs.cache-hit != 'true'}}
         name: Install Kokkos

MParT/BasisEvaluator.h

@@ -1,3 +1,6 @@
+#ifndef MPART_BASISEVALUATOR_H
+#define MPART_BASISEVALUATOR_H
+
 #include <Kokkos_Core.hpp>
 #include "MParT/PositiveBijectors.h"
 
@@ -327,4 +330,5 @@ class BasisEvaluator<BasisHomogeneity::Heterogeneous,
   // NOTE: shared_ptr is necessary to avoid type slicing
   std::vector<std::shared_ptr<CommonBasisEvaluatorType>> basis1d_;
 };
-}  // namespace mpart
+}  // namespace mpart
+#endif // MPART_BASISEVALUATOR_H

MParT/MapFactory.h

@@ -67,8 +67,7 @@ namespace mpart{
                                                                               std::shared_ptr<ConditionalMapBase<MemorySpace>> const &comp);
 
 
-                                                                                  /**
-            @brief Constructs a triangular map with MonotoneComponents for each block.  A total order multiindex
+        /** @brief Constructs a triangular map with MonotoneComponents for each block.  A total order multiindex
                    set is used to define the MonotoneComponent.
 
             @details For more control over the individual components, consider constructing components with
@@ -116,31 +115,80 @@ namespace mpart{
                                                                               std::shared_ptr<ConditionalMapBase<MemorySpace>> const &comp);
 
         /**
-         * @brief Create a MultivariateExpansion using a sigmoid object; INTERNAL USE ONLY
+         * @brief Create a RectifiedMultivariateExpansion using sigmoids
          * 
          * @details Using the parameters described in \c`opts`, this function creates
-         * a "multivariate expansion" object that uses a collection of sigmoid functions
+         * a "rectified multivariate expansion" object that uses a collection of sigmoid functions
          * to map the input space to the output space. The pertinent options are
          * - \c`opts.basisType` : The type of basis function to use in the expansion for the first \f$d-1\f$ inputs
          * - \c`opts.posFuncType` : The type of positive function to use in the rectified basis @ref MultivariateExpansionWorker
          * - \c`opts.edgeWidth` : The width of the "edge terms" on the last input @ref Sigmoid1d
          * - \c`opts.sigmoidType` : The type of sigmoid function to use in the expansion @ref Sigmoid1d
          * 
-         * By default, this constructs a total-order multi-index set. Note that, since the basis
-         * evaluator is rectified, this *will always* be a positive-valued function in the first \f$d-1\f$ inputs.
-         * As such, this factory function is largely for internal testing.
+         * By default, this constructs total-order multi-index sets. This is only creating a real-valued function.
+         * To create a map, use `TriangularMap` function.
          * 
          * @tparam MemorySpace 
-         * @param inputDim 
-         * @param centers 
+         * @param inputDim Dimension of the input space to component
+         * @param totalOrder Total order of your multi-index set. If 0, deduce from length of centers.
+         * @param centers Set of centers of the sigmoids. Should be length 2 + (1 + 2 + ... + p) where p is the default total order.
          * @param opts 
-         * @return std::shared_ptr<ParameterizedFunctionBase<MemorySpace>> 
+         * @return std::shared_ptr<ConditionalMapBase<MemorySpace>> 
          */
         template<typename MemorySpace>
-        std::shared_ptr<ParameterizedFunctionBase<MemorySpace>> CreateSigmoidExpansion(
-            unsigned int inputDim, StridedVector<double, MemorySpace> centers,
+        std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidComponent(
+            unsigned int inputDim, unsigned int totalOrder, StridedVector<const double, MemorySpace> centers,
             MapOptions opts);
 
+        template<typename MemorySpace, std::enable_if_t<std::is_same_v<MemorySpace, Kokkos::HostSpace>, bool> = true>
+        std::shared_ptr<ConditionalMapBase<Kokkos::HostSpace>> CreateSigmoidComponent(
+            unsigned int inputDim, unsigned int totalOrder, Eigen::Ref<const Eigen::RowVectorXd> centers,
+            MapOptions opts) {
+            StridedVector<const double, Kokkos::HostSpace> centersVec = ConstVecToKokkos<double, Kokkos::HostSpace>(centers);
+            return CreateSigmoidComponent<Kokkos::HostSpace>(inputDim, totalOrder, centersVec, opts);
+        }
+
+        template<typename MemorySpace>
+        std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidComponent(
+            FixedMultiIndexSet<MemorySpace> mset_offdiag, FixedMultiIndexSet<MemorySpace> mset_diag,
+            StridedVector<const double, MemorySpace> centers, MapOptions opts);
+
+        template<typename MemorySpace, std::enable_if_t<std::is_same_v<MemorySpace, Kokkos::HostSpace>, bool> = true>
+        std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidComponent(
+            FixedMultiIndexSet<MemorySpace> mset_offdiag, FixedMultiIndexSet<MemorySpace> mset_diag,
+            Eigen::Ref<const Eigen::RowVectorXd> centers, MapOptions opts) {
+            StridedVector<const double, Kokkos::HostSpace> centersVec = ConstVecToKokkos<double, Kokkos::HostSpace>(centers);
+            return CreateSigmoidComponent<Kokkos::HostSpace>(mset_offdiag, mset_diag, centersVec, opts);
+        }
+
+        template<typename MemorySpace>
+        std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidTriangular(
+            unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder,
+            std::vector<StridedVector<const double, MemorySpace>> const& centers, MapOptions opts
+        );
+
+        template<typename MemorySpace>
+        std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidTriangular(
+            unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder,
+            StridedMatrix<const double, MemorySpace> const& centers, MapOptions opts
+        ) {
+            std::vector<StridedVector<const double, MemorySpace>> centersVecs;
+            for(unsigned int i = 0; i < centers.extent(1); i++){
+                StridedVector<const double, MemorySpace> center_i = Kokkos::subview(centers, Kokkos::ALL(), i);
+                centersVecs.push_back(center_i);
+            }
+            return CreateSigmoidTriangular<MemorySpace>(inputDim, outputDim, totalOrder, centersVecs, opts);
+        }
+
+        template<typename MemorySpace, std::enable_if_t<std::is_same_v<MemorySpace, Kokkos::HostSpace>, bool> = true>
+        std::shared_ptr<ConditionalMapBase<Kokkos::HostSpace>> CreateSigmoidTriangular(
+            unsigned int inputDim, unsigned int outputDim, unsigned int totalOrder,
+            Eigen::Ref<const Eigen::RowMatrixXd> const& centers, MapOptions opts
+        ) {
+            StridedMatrix<const double, Kokkos::HostSpace> centersMat = ConstRowMatToKokkos<double,Kokkos::HostSpace>(centers);
+            return CreateSigmoidTriangular<Kokkos::HostSpace>(inputDim, outputDim, totalOrder, centersMat, opts);
+        }
+
         /** This struct is used to map the options to functions that can create a map component with types corresponding
             to the options.
         */

MParT/MultiIndices/MultiIndexLimiter.h

@@ -27,50 +27,40 @@ namespace MultiIndexLimiter{
 
     };
 
-    /** @class SeparableTotalOrderLimiter
-        @brief Same as TotalOrder without cross-terms
-        @details This limter only allows terms that satisfy
-        \f$\|\mathbf{j}\|_1\leq p_U\f$, where \f$\mathbf{j}\f$
-        is the multiindex, and \f$p_U\f$ is a nonnegative integer passed to the
-        constructor of this class.
+    /** @class Separable
+        @brief Restricts multi-indices to refuse cross-terms
+        @details This limiter only allows terms that satisfy
+        \f$\mathbf{j}_d = 0\f$ *or* \f$\mathbf{j}_d=\|\mathbf{j}\|_1\f$,
+        where \f$\mathbf{j}\in\mathbb{N}_0^d\f$ is the multiindex.
     */
-    class SeparableTotalOrder{
+    class Separable{
 
     public:
 
-        SeparableTotalOrder(unsigned int totalOrderIn) : totalOrder(totalOrderIn){};
+        Separable(){};
 
         bool operator()(MultiIndex const& multi){
           unsigned int sum = multi.Sum();
-          return (sum <= totalOrder) && (!multi.HasNonzeroEnd() || sum == multi.Get(multi.Length()-1));
+          return !multi.HasNonzeroEnd() || sum == multi.Get(multi.Length()-1);
         };
 
-    private:
-        const unsigned int totalOrder;
-
     };
 
-    /** @class NonzeroDiagTotalOrderLimiter
-        @brief Same as TotalOrder, except without any term that has nonzero diagonal entries
-        @details This limter only allows terms that satisfy
-        \f$\|\mathbf{j}\|_1\leq p_U\f$, where \f$\mathbf{j}\f$
-        is the multiindex, and \f$p_U\f$ is a nonnegative integer passed to the
-        constructor of this class.
+    /** @class NonzeroDiag
+        @brief Restricts acceptable to any term that has nonzero diagonal entries
+        @details This limiter only allows terms that satisfy
+        \f$\mathbf{j}_d \neq 0\f$ when the multi-index \f$\mathbf{j}\in\mathbb{N}_0^d\f$.
     */
-    class NonzeroDiagTotalOrderLimiter{
+    class NonzeroDiag{
 
     public:
 
-        NonzeroDiagTotalOrderLimiter(unsigned int totalOrderIn) : totalOrder(totalOrderIn){};
+        NonzeroDiag() {};
 
         bool operator()(MultiIndex const& multi){
-          unsigned int sum = multi.Sum();
-          return (sum <= totalOrder) && (multi.HasNonzeroEnd());
+          return (multi.HasNonzeroEnd());
         };
 
-    private:
-        const unsigned int totalOrder;
-
     };
 
 

MParT/RectifiedMultivariateExpansion.h

@@ -9,6 +9,7 @@
 #include "MParT/Utilities/KokkosHelpers.h"
 
 #include <algorithm>
+#include <numeric>
 
 namespace mpart{
 
@@ -31,7 +32,8 @@ namespace mpart{
             Rectifier>, MemorySpace
         >;
         using OffdiagWorker_T = MultivariateExpansionWorker<
-            BasisEvaluator<BasisHomogeneity::Homogeneous, OffdiagEval>
+            BasisEvaluator<BasisHomogeneity::Homogeneous, OffdiagEval>,
+            MemorySpace
         >;
 
 
@@ -42,11 +44,7 @@ namespace mpart{
                                     setSize_diag(worker_diag_.NumCoeffs()),
                                     worker_off(worker_off_),
                                     worker_diag(worker_diag_)
-        {
-            // TODO: Check that the inputs are compatible
-            // MVE_diag has no terms constant in last input
-            // Then, ensure all of these have the appropriate coefficient views
-        };
+        {};
 
         ~RectifiedMultivariateExpansion() = default;
 
@@ -179,8 +177,8 @@ namespace mpart{
             // Figure out how much memory we'll need in the cache
             unsigned int cacheSize = std::max(worker_diag.CacheSize(), worker_off.CacheSize());
             unsigned int maxParams = coeff_off.size() + coeff_diag.size();
-            Kokkos::pair<unsigned int, unsigned int> coeff_off_idx {0u,coeff_off.size()};
-            Kokkos::pair<unsigned int, unsigned int> coeff_diag_idx {coeff_off.size(), maxParams};
+            Kokkos::pair<unsigned int, unsigned int> coeff_off_idx {0u,(unsigned int)coeff_off.size()};
+            Kokkos::pair<unsigned int, unsigned int> coeff_diag_idx {(unsigned int) coeff_off.size(), maxParams};
 
             auto functor = KOKKOS_CLASS_LAMBDA (typename Kokkos::TeamPolicy<ExecutionSpace>::member_type team_member) {
 
@@ -424,6 +422,13 @@ namespace mpart{
             Kokkos::fence();
         }
 
+        std::vector<unsigned int> DiagonalCoeffIndices() const
+        {
+            std::vector<unsigned int> diagIndices(setSize_diag);
+            std::iota(diagIndices.begin(), diagIndices.end(), setSize_off);
+            return diagIndices;
+        }
+
     private:
         template<typename PointType, typename CoeffType>
         struct SingleWorkerEvaluator {

MParT/Sigmoid.h

@@ -106,7 +106,7 @@ class Sigmoid1d {
 	 * @param widths
 	 */
 	Sigmoid1d(Kokkos::View<double*, MemorySpace> centers,
-						Kokkos::View<double*, MemorySpace> widths)
+			  Kokkos::View<double*, MemorySpace> widths)
 			: centers_(centers), widths_(widths) {
 		Kokkos::View<double*, MemorySpace> weights ("Sigmoid weights", centers.extent(0));
 		Kokkos::parallel_for(centers.extent(0), KOKKOS_LAMBDA(unsigned int i) { weights(i) = 1.; });