Extended interfaces to sigmoids to leverage more flexible definitions…

… of 'order'.

MParT/MapOptions.h

@@ -4,6 +4,7 @@
 #include <string>
 #include <sstream>
 #include <limits>
+#include "MParT/Sigmoid.h"
 
 namespace mpart{
 
@@ -65,9 +66,12 @@ namespace mpart{
         /** The type of edge terms to use with sigmoids */
         EdgeTypes edgeType = EdgeTypes::SoftPlus;
 
-        /** The "shape" of the edge terms in a sigmoid expansion (larger is "sharper" edge terms)*/
+        /** The "shape" of the edge terms in a sigmoid basis (larger is "sharper" edge terms)*/
         double edgeShape = 1.5;
 
+        /** How arrays of centers should be interpreted in sigmoid basis functions. */
+        SigmoidSumSizeType sigmoidBasisSumType;
+
         /** Linearization bounds for the 1d basis function. The basis function is linearized outside [lb,ub] */
         double basisLB = -std::numeric_limits<double>::infinity();
         double basisUB =  std::numeric_limits<double>::infinity();
@@ -132,6 +136,8 @@ namespace mpart{
             ret &= (basisLB     == opts2.basisLB);
             ret &= (basisUB     == opts2.basisUB);
             ret &= (edgeType    == opts2.edgeType);
+            ret &= (edgeShape   == opts2.edgeShape);
+            ret &= (sigmoidBasisSumType == opts2.sigmoidBasisSumType);
             ret &= (sigmoidType == opts2.sigmoidType);
             ret &= (posFuncType == opts2.posFuncType);
             ret &= (quadType    == opts2.quadType);
@@ -153,6 +159,7 @@ namespace mpart{
             ss << "basisUB = " << basisUB << "\n";
             ss << "edgeType = " << etypes[static_cast<unsigned int>(edgeType)] << "\n";
             ss << "sigmoidType = " << stypes[static_cast<unsigned int>(sigmoidType)] << "\n";
+            ss << "sigmoidBasisSumType = " << sbstypes[static_cast<unsigned int>(sigmoidBasisSumType)] << "\n";
             ss << "basisNorm = " << (basisNorm ? "true" : "false") << "\n";
             ss << "posFuncType = " << pftypes[static_cast<unsigned int>(posFuncType)] << "\n";
             ss << "quadType = " << qtypes[static_cast<unsigned int>(quadType)] << "\n";
@@ -171,6 +178,7 @@ namespace mpart{
         inline static const std::string qtypes[3] = {"ClenshawCurtis", "AdaptiveSimpson", "AdaptiveClenshawCurtis"};
         inline static const std::string etypes[1] = {"SoftPlus"};
         inline static const std::string stypes[1] = {"Logistic"};
+        inline static const std::string sbstypes[2] = {"Linear", "Constant"};
     };
 };
 

bindings/python/CMakeLists.txt

@@ -13,6 +13,7 @@ set(PYTHON_BINDING_SOURCES
   # src/DebugMap.cpp
   src/MapFactory.cpp
   src/IdentityMap.cpp
+  src/Sigmoid.cpp
 
   ../common/src/CommonUtilities.cpp
 )

bindings/python/include/CommonPybindUtilities.h

@@ -43,6 +43,8 @@ void SummarizedMapWrapper(pybind11::module &m);
 template<typename MemorySpace>
 void IdentityMapWrapper(pybind11::module &m);
 
+void Sigmoid1dWrapper(pybind11::module &m);
+
 // template<typename MemorySpace>
 // void DebugMapWrapper(pybind11::module &m);
 

bindings/python/src/MapOptions.cpp

@@ -38,6 +38,11 @@ void mpart::binding::MapOptionsWrapper(py::module &m)
     py::enum_<SigmoidTypes>(m, "SigmoidTypes")
     .value("Logistic",SigmoidTypes::Logistic);
 
+    // SigmoidSumTypes
+    py::enum_<SigmoidSumSizeType>(m, "SigmoidSumSizeType")
+    .value("Linear",SigmoidSumSizeType::Linear)
+    .value("Constant",SigmoidSumSizeType::Constant);
+
     // EdgeTypes
     py::enum_<EdgeTypes>(m, "EdgeTypes")
     .value("SoftPlus",EdgeTypes::SoftPlus);
@@ -56,6 +61,7 @@ void mpart::binding::MapOptionsWrapper(py::module &m)
     .def_readwrite("edgeType", &MapOptions::edgeType)
     .def_readwrite("edgeShape", &MapOptions::edgeShape)
     .def_readwrite("sigmoidType", &MapOptions::sigmoidType)
+    .def_readwrite("sigmoidBasisSumType", &MapOptions::sigmoidBasisSumType)
     .def_readwrite("quadType", &MapOptions::quadType)
     .def_readwrite("quadAbsTol", &MapOptions::quadAbsTol)
     .def_readwrite("quadRelTol", &MapOptions::quadRelTol)

bindings/python/src/Sigmoid.cpp

@@ -0,0 +1,52 @@
+#include "CommonPybindUtilities.h"
+#include "MParT/Sigmoid.h"
+#include "MParT/Utilities/ArrayConversions.h"
+#include <Eigen/Core>
+#include <pybind11/stl.h>
+#include <pybind11/eigen.h>
+
+#include <Kokkos_Core.hpp>
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+using namespace mpart::binding;
+
+void mpart::binding::Sigmoid1dWrapper(py::module &m)
+{
+    typedef Sigmoid1d<Kokkos::HostSpace, mpart::SigmoidTypeSpace::Logistic, mpart::SoftPlus> SigmoidType;
+
+    py::class_<SigmoidType, std::shared_ptr<SigmoidType>>(m, "Sigmoid1d")
+        .def(py::init([](Eigen::Ref<Eigen::VectorXd> centers, 
+                         Eigen::Ref<Eigen::VectorXd> widths,
+                         Eigen::Ref<Eigen::VectorXd> weights, 
+                         Eigen::Ref<Eigen::Matrix<unsigned int, Eigen::Dynamic, 1>> startInds) {
+                return std::make_shared<SigmoidType>(mpart::VecToKokkos<double>(centers),
+                                                     mpart::VecToKokkos<double>(widths),
+                                                     mpart::VecToKokkos<double>(weights),
+                                                     mpart::VecToKokkos<unsigned int>(startInds));
+            }))
+        .def("EvaluateAll", [](const SigmoidType &sig,
+                               Eigen::Ref<Eigen::VectorXd> output, 
+                               double x){sig.EvaluateAll(output.data(),output.size()-1,x);})
+        .def("EvaluateDerivatives", [](const SigmoidType &sig,
+                               Eigen::Ref<Eigen::VectorXd> output, 
+                               Eigen::Ref<Eigen::VectorXd> derivs,
+                               double x){
+                                if(output.size()!=derivs.size()){
+                                    throw std::invalid_argument("Vector dimensions do not match.");
+                                }else{
+                                   sig.EvaluateDerivatives(output.data(),derivs.data(),output.size()-1,x);
+                                }})
+        .def("EvaluateSecondDerivatives", [](const SigmoidType &sig,
+                               Eigen::Ref<Eigen::VectorXd> output, 
+                               Eigen::Ref<Eigen::VectorXd> derivs,
+                               Eigen::Ref<Eigen::VectorXd> derivs2,
+                               double x){
+                                if((output.size()!=derivs.size())||(derivs.size()!=derivs2.size())){
+                                    throw std::invalid_argument("Vector dimensions do not match.");
+                                }else{
+                                   sig.EvaluateSecondDerivatives(output.data(),derivs.data(),derivs2.data(),output.size()-1,x);
+                                }})
+        ;
+
+}

bindings/python/src/Wrapper.cpp

@@ -21,7 +21,8 @@ PYBIND11_MODULE(pympart, m) {
     // DebugMapWrapper<Kokkos::HostSpace>(m);
     MapFactoryWrapper<Kokkos::HostSpace>(m);
     MapObjectiveWrapper<Kokkos::HostSpace>(m);
-
+    Sigmoid1dWrapper(m);
+
 #if defined(MPART_HAS_NLOPT)
     TrainOptionsWrapper(m);
     ATMOptionsWrapper(m);

src/MapFactory.cpp

@@ -174,63 +174,98 @@ std::shared_ptr<ParameterizedFunctionBase<MemorySpace>> mpart::MapFactory::Creat
     return nullptr;
 }
 
-
+/** Constructs a sigmoid basis from specified centers.  Chooses the widths of the 
+ *  sigmoid functions using a rule of thumb based on the distance between neighboring
+ *  centers.
+*/
 template <typename MemorySpace, typename SigmoidType, typename EdgeType>
-Sigmoid1d<MemorySpace, SigmoidType, EdgeType> CreateSigmoid(
-    unsigned int inputDim, StridedVector<double, MemorySpace> centers, double edgeShape) {
-    Kokkos::View<double*, MemorySpace> widths("Sigmoid Widths", centers.size());
-    Kokkos::View<double*, MemorySpace> weights("Sigmoid Weights", centers.size());
-    int sigmoid_count = centers.size() - 2;
-    double max_order_double = Sigmoid1d<MemorySpace, SigmoidType, EdgeType>::LengthToOrder(sigmoid_count);
-    int max_order = int(max_order_double);
-    if(max_order_double  != double(max_order)) {
-        std::stringstream ss;
-        ss << "Incorrect length of centers, " << centers.size() << ".\n";
-        ss << "Length should be of form 2+(1+2+3+...+n) for some order n";
-        ProcAgnosticError<std::invalid_argument>(
-            ss.str().c_str());
-    }
-    unsigned int N = max_order+2; 
-    Kokkos::RangePolicy<typename MemoryToExecution<MemorySpace>::Space> policy{0, N};
-    Kokkos::parallel_for(policy, KOKKOS_LAMBDA(unsigned int i) {
-        if (i == max_order) {
-            widths(0) = edgeShape;
-            weights(0) = 1./edgeShape;
-            return;
-        } else if (i == max_order+1) {
-            widths(1) = edgeShape;
-            weights(1) = 1./edgeShape;
-            return;
+Sigmoid1d<MemorySpace, SigmoidType, EdgeType> CreateSigmoid(unsigned int inputDim, 
+                                                            StridedVector<double, MemorySpace> centers, 
+                                                            double edgeShape,
+                                                            SigmoidSumSizeType sumType) 
+{
+    //Kokkos::View<double*, MemorySpace> widths("Sigmoid Widths", centers.size());
+    //Kokkos::View<double*, MemorySpace> weights("Sigmoid Weights", centers.size());
+
+
+    int numSigmoids = Sigmoid1d<MemorySpace, SigmoidType, EdgeType>::ComputeNumSigmoids(centers.size(), sumType);
+
+    // Fill vectors on host and then copy to device if necessary
+    Kokkos::View<double*, Kokkos::HostSpace> hWidths("Widths", centers.extent(0));
+    Kokkos::View<double*, Kokkos::HostSpace> hWeights("Weights", centers.extent(0));
+    Kokkos::View<double*, Kokkos::HostSpace> hCenters("Centers", centers.extent(0));
+    Kokkos::deep_copy(hCenters, centers);
+    Kokkos::deep_copy(hWeights, 1.0);
+
+    // Set widths and weights for edge terms
+    hWidths(0) = edgeShape;
+    hWeights(0) = 1.0/edgeShape;
+    hWidths(1) = edgeShape;
+    hWeights(1) = 1.0/edgeShape;
+
+    // Now set sigmoid widths/weights
+    if(sumType==SigmoidSumSizeType::Linear){
+
+        for(int i=0; i<numSigmoids; i++){
+            int start_idx = 2+(i*(i+1))/2;
+
+            for(unsigned int j = 0; j < i; j++) {
+                double prev_center, next_center;
+                if(j == 0 || i == 0) {// Use center for left edge term
+                    prev_center = centers(0);
+                } else {
+                    prev_center = centers(start_idx+j-1);
+                }
+                if(j == i-1 || i == 0) { // Use center for right edge term
+                    next_center = centers(1);
+                } else {
+                    next_center = centers(start_idx+j);
+                }
+                hWidths(start_idx + j) = 2/(next_center - prev_center);
+                hWeights(start_idx + j) = 1.;
+            }
         }
-        int start_idx = 2+(i*(i+1))/2;
-
-        for(unsigned int j = 0; j < i; j++) {
-            double prev_center, next_center;
-            if(j == 0 || i == 0) {// Use center for left edge term
-                prev_center = centers(0);
-            } else {
-                prev_center = centers(start_idx+j-1);
+    }else{
+        if(numSigmoids>0){
+            // Make sure the non-edge centers are sorted
+            std::sort(&hCenters(2), &hCenters(hCenters.extent(0)-1)+1);
+
+            // Set first sigmoid
+            if(numSigmoids>1){
+                hWidths(2) = 4.0/(hCenters(3)-hCenters(2)); // Second sigmoid center minus first
+            }else{
+                hWidths(2) = 2.0/std::abs(hCenters(1)-hCenters(0)); // Distance between edge centers
+            }
+
+            // Set "interior" sigmoids
+            for(int i=1; i<numSigmoids-1; i++){
+                hWidths[2+i] = 2.0/(hCenters[3+i]-hCenters[1+i]);
             }
-            if(j == i-1 || i == 0) { // Use center for right edge term
-                next_center = centers(1);
-            } else {
-                next_center = centers(start_idx+j);
+
+            // If there is more than one sigmoid, set the last one
+            if(numSigmoids>1){
+                hWidths(2+numSigmoids-1) = 4.0/(hCenters(2+numSigmoids-1)-hCenters(2+numSigmoids-2));
             }
-            widths(start_idx + j) = 2/(next_center - prev_center);
-            weights(start_idx + j) = 1.;
         }
-    });
-    Sigmoid1d<MemorySpace, SigmoidType, EdgeType> sig {centers, widths, weights};
-    return sig;
+    }
+
+    Kokkos::View<double*, MemorySpace> widths = Kokkos::create_mirror_view_and_copy(MemorySpace(), hWidths);
+    Kokkos::View<double*, MemorySpace> weights = Kokkos::create_mirror_view_and_copy(MemorySpace(), hWeights);
+    Kokkos::View<double*, MemorySpace> dCenters = Kokkos::create_mirror_view_and_copy(MemorySpace(), hCenters);
+
+    Sigmoid1d<MemorySpace, SigmoidType, EdgeType> sig {dCenters, widths, weights, sumType};
+    return sig;  
 }
 
 template<typename MemorySpace, typename OffdiagEval, typename Rectifier, typename SigmoidType, typename EdgeType>
 using SigmoidBasisEval = BasisEvaluator<BasisHomogeneity::OffdiagHomogeneous, Kokkos::pair<OffdiagEval, Sigmoid1d<MemorySpace, SigmoidType, EdgeType>>>;
 
 template <typename MemorySpace, typename OffdiagEval, typename Rectifier, typename SigmoidType, typename EdgeType>
-std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(
-    FixedMultiIndexSet<MemorySpace> mset_offdiag, FixedMultiIndexSet<MemorySpace> mset_diag,
-    StridedVector<double, MemorySpace> centers, double edgeWidth)
+std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(FixedMultiIndexSet<MemorySpace> mset_offdiag, 
+                                                                                FixedMultiIndexSet<MemorySpace> mset_diag,
+                                                                                StridedVector<double, MemorySpace> centers, 
+                                                                                double edgeWidth,
+                                                                                SigmoidSumSizeType sumType)
 {
     unsigned int inputDim = mset_diag.Length();
     if(inputDim != 1 && inputDim != mset_offdiag.Length() + 1) {
@@ -243,7 +278,7 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(
     using Sigmoid_T = Sigmoid1d<MemorySpace, SigmoidType, EdgeType>;
     using DiagBasisEval_T = BasisEvaluator<BasisHomogeneity::OffdiagHomogeneous, Kokkos::pair<OffdiagEval, Sigmoid_T>, Rectifier>;
     using OffdiagBasisEval_T = BasisEvaluator<BasisHomogeneity::Homogeneous, OffdiagEval>;
-    auto sigmoid = CreateSigmoid<MemorySpace, SigmoidType, EdgeType>(inputDim, centers, edgeWidth);
+    auto sigmoid = CreateSigmoid<MemorySpace, SigmoidType, EdgeType>(inputDim, centers, edgeWidth, sumType);
     if(inputDim == 1) {
         unsigned int maxOrder = mset_diag.Size() - 1;
         auto output = std::make_shared<UnivariateExpansion<MemorySpace, Sigmoid_T>>(maxOrder, sigmoid);
@@ -262,24 +297,41 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(
 }
 
 template <typename MemorySpace, typename OffdiagEval, typename Rectifier, typename SigmoidType, typename EdgeType>
-std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(
-    unsigned int inputDim, unsigned int totalOrder, StridedVector<double, MemorySpace> centers, double edgeWidth)
+std::shared_ptr<ConditionalMapBase<MemorySpace>> CreateSigmoidExpansionTemplate(unsigned int inputDim, 
+                                                                                unsigned int totalOrder, 
+                                                                                StridedVector<double, MemorySpace> centers, 
+                                                                                double edgeWidth,
+                                                                                SigmoidSumSizeType sumType)
 {
     using Sigmoid_T = Sigmoid1d<MemorySpace, SigmoidType, EdgeType>;
-    if(totalOrder == 0)
-        totalOrder = Sigmoid_T::LengthToOrder(centers.size()-2)+3;
+    int numSigmoids = Sigmoid_T::ComputeNumSigmoids(centers.size(), sumType);
+
     if(inputDim == 1) {
-        FixedMultiIndexSet<MemorySpace> mset {1, totalOrder};
+        FixedMultiIndexSet<MemorySpace> mset(1, numSigmoids+3);
         FixedMultiIndexSet<MemorySpace> mset_offdiag {1,0};
         return CreateSigmoidExpansionTemplate<MemorySpace, OffdiagEval, Rectifier, SigmoidType, EdgeType>(
-            mset_offdiag, mset, centers, edgeWidth);
+            mset_offdiag, mset, centers, edgeWidth, sumType);
     }
-    MultiIndexSet mset = MultiIndexSet::CreateTotalOrder(inputDim, totalOrder, MultiIndexLimiter::NonzeroDiag());
+
+    // Create a multiindex containing the tensor product of a total order set in the first n-1 inputs and all sigmoid terms
+    // Start with a total order set having all zeros on the diagonal
+    MultiIndexSet mset0 = MultiIndexSet::CreateTotalOrder(inputDim, totalOrder, MultiIndexLimiter::Dimension(0, inputDim-1));
+    // Now construct the tensor-product set
+    MultiIndexSet mset(inputDim);
+    for(int baseInd=0; baseInd<mset0.Size(); ++baseInd){
+        MultiIndex m = mset0.at(baseInd);
+        for(int sigInd=0; sigInd<numSigmoids+3; ++sigInd){
+            m.Set(inputDim-1, sigInd);
+            mset += m;
+        }
+    }
+
     FixedMultiIndexSet<MemorySpace> fmset_diag_d = mset.Fix(true).ToDevice<MemorySpace>();
     FixedMultiIndexSet<Kokkos::HostSpace> fmset_offdiag {inputDim-1, totalOrder};
     FixedMultiIndexSet<MemorySpace> fmset_offdiag_d = fmset_offdiag.ToDevice<MemorySpace>();
+
     return CreateSigmoidExpansionTemplate<MemorySpace, OffdiagEval, Rectifier, SigmoidType, EdgeType>(
-        fmset_offdiag_d, fmset_diag_d, centers, edgeWidth);
+        fmset_offdiag_d, fmset_diag_d, centers, edgeWidth, sumType);
 }
 
 template<typename MemorySpace>
@@ -319,9 +371,9 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> MapFactory::CreateSigmoidCompon
     Kokkos::deep_copy(centers_copy, centers);
     // Dispatch to the correct sigmoid expansion template
     if(opts.posFuncType == PosFuncTypes::Exp) {
-        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, Exp, SigmoidTypeSpace::Logistic, SoftPlus>(inputDim, totalOrder, centers_copy, opts.edgeShape);
+        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, Exp, SigmoidTypeSpace::Logistic, SoftPlus>(inputDim, totalOrder, centers_copy, opts.edgeShape, opts.sigmoidBasisSumType);
     } else if(opts.posFuncType == PosFuncTypes::SoftPlus) {
-        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, SoftPlus, SigmoidTypeSpace::Logistic, SoftPlus>(inputDim, totalOrder, centers_copy, opts.edgeShape);
+        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, SoftPlus, SigmoidTypeSpace::Logistic, SoftPlus>(inputDim, totalOrder, centers_copy, opts.edgeShape, opts.sigmoidBasisSumType);
     }
     else {
         return nullptr;
@@ -337,9 +389,9 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> MapFactory::CreateSigmoidCompon
     Kokkos::deep_copy(centers_copy, centers);
     // Dispatch to the correct sigmoid expansion template
     if(opts.posFuncType == PosFuncTypes::Exp) {
-        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, Exp, SigmoidTypeSpace::Logistic, SoftPlus>(mset_offdiag, mset_diag, centers_copy, opts.edgeShape);
+        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, Exp, SigmoidTypeSpace::Logistic, SoftPlus>(mset_offdiag, mset_diag, centers_copy, opts.edgeShape, opts.sigmoidBasisSumType);
     } else if(opts.posFuncType == PosFuncTypes::SoftPlus) {
-        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, SoftPlus, SigmoidTypeSpace::Logistic, SoftPlus>(mset_offdiag, mset_diag, centers_copy, opts.edgeShape);
+        return CreateSigmoidExpansionTemplate<MemorySpace, HermiteFunction, SoftPlus, SigmoidTypeSpace::Logistic, SoftPlus>(mset_offdiag, mset_diag, centers_copy, opts.edgeShape, opts.sigmoidBasisSumType);
     }
     else {
         return nullptr;