diff --git a/.github/workflows/build-external-lib-tests.yml b/.github/workflows/build-external-lib-tests.yml
index ebe21dac..3e880ed4 100644
--- a/.github/workflows/build-external-lib-tests.yml
+++ b/.github/workflows/build-external-lib-tests.yml
@@ -59,7 +59,7 @@ jobs:
                 -DCMAKE_CXX_STANDARD=17    \
                 -DCMAKE_INSTALL_PREFIX=$GITHUB_WORKSPACE/KOKKOS_INSTALL/ \
                 ../
-          sudo make -j2 install
+          sudo make -j$(nproc) install
 
       - name: Install Conda Dependencies
         shell: bash -l {0}
@@ -74,7 +74,7 @@ jobs:
 
       - name: Build MParT
         shell: bash -l {0}
-        run: cd $GITHUB_WORKSPACE/mpart/build; make -j2
+        run: cd $GITHUB_WORKSPACE/mpart/build; make -j$(nproc)
 
       - name: Run Tests
         shell: bash -l {0}
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<MemorySpace>
     void LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                      StridedMatrix<double, MemorySpace>              output) override;
 
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> 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 ebac1ead..c1da168c 100644
--- a/MParT/ComposedMap.h
+++ b/MParT/ComposedMap.h
@@ -102,6 +102,8 @@ class ComposedMap : public ConditionalMapBase<MemorySpace>
     void GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
                       StridedMatrix<const double, MemorySpace> const& sens,
                       StridedMatrix<double, MemorySpace>              output) override;
+
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> Slice(int a, int b) override;
 private:
 
     unsigned int maxChecks_;
diff --git a/MParT/ConditionalMapBase.h b/MParT/ConditionalMapBase.h
index 3689bf3a..9fbab75b 100644
--- a/MParT/ConditionalMapBase.h
+++ b/MParT/ConditionalMapBase.h
@@ -149,6 +149,11 @@ namespace mpart {
 
         virtual void LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                  StridedMatrix<double, MemorySpace>              output) = 0;
+
+        virtual std::shared_ptr<ConditionalMapBase<MemorySpace>> Slice(int a, int b) {
+            throw std::runtime_error("Slice not implemented for this map type.");
+        }
+
 #if defined(MPART_HAS_CEREAL)
         template<typename Archive>
         void serialize(Archive& ar) {
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<MemorySpace>
     void LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                      StridedMatrix<double, MemorySpace>              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<ConditionalMapBase<MemorySpace>> Slice(int a, int b) override;
+
 }; // class IdentityMap
 
 }
diff --git a/MParT/MapFactory.h b/MParT/MapFactory.h
index 64001aef..d8e38248 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/SummarizedMap.h"
 //#include "MParT/DebugMap.h"
 #include "MParT/MultiIndices/FixedMultiIndexSet.h"
diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h
index 251deb5c..0812e5ae 100644
--- a/MParT/MonotoneComponent.h
+++ b/MParT/MonotoneComponent.h
@@ -1,6 +1,8 @@
 #ifndef MPART_MONOTONECOMPONENT_H
 #define MPART_MONOTONECOMPONENT_H
 
+#include <memory>
+
 #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 ExpansionType, class PosFuncType, class QuadratureType, typename MemorySpace>
-class MonotoneComponent : public ConditionalMapBase<MemorySpace>
+class MonotoneComponent : public ConditionalMapBase<MemorySpace>, public std::enable_shared_from_this<MonotoneComponent<ExpansionType, PosFuncType, QuadratureType, MemorySpace>>
 {
 
 public:
@@ -1078,6 +1080,11 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
         return 0.5*(xub+xlb);
     }
 
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> Slice(int a, int b) override {
+        assert(a == 0 && b == 1);
+        return MonotoneComponent<ExpansionType, PosFuncType, QuadratureType, MemorySpace>::shared_from_this();
+    }
+
     /** Give access to the underlying FixedMultiIndexSet
      * @return The FixedMultiIndexSet
      */
diff --git a/MParT/TriangularMap.h b/MParT/TriangularMap.h
index 086ff3f8..95a2b6dc 100644
--- a/MParT/TriangularMap.h
+++ b/MParT/TriangularMap.h
@@ -60,6 +60,18 @@ class TriangularMap : public ConditionalMapBase<MemorySpace>{
     void WrapCoeffs(Kokkos::View<double*, mpart::DeviceSpace> coeffs) override;
     #endif
 
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> 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.
+     *          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 ConditionalMapBase object containing the subsection of the map.
+     */
+    virtual std::shared_ptr<ConditionalMapBase<MemorySpace>> BlockSlice(int a, int b);
+
     virtual std::shared_ptr<ConditionalMapBase<MemorySpace>> GetComponent(unsigned int i){ return comps_.at(i);}
 
     /** @brief Computes the log determinant of the Jacobian matrix of this map.
diff --git a/bindings/julia/src/ConditionalMapBase.cpp b/bindings/julia/src/ConditionalMapBase.cpp
index c1f727e5..56ce235a 100644
--- a/bindings/julia/src/ConditionalMapBase.cpp
+++ b/bindings/julia/src/ConditionalMapBase.cpp
@@ -39,6 +39,7 @@ void mpart::binding::ConditionalMapBaseWrapper(jlcxx::Module &mod) {
             map.InverseImpl(JuliaToKokkos(x1), JuliaToKokkos(r), JuliaToKokkos(output));
             return output;
         })
+        .method("Slice", [](ConditionalMapBase<Kokkos::HostSpace>& map, int begin, int end){ return map.Slice(begin-1, end); })
         ;
     jlcxx::stl::apply_stl<ConditionalMapBase<Kokkos::HostSpace>*>(mod);
 }
\ No newline at end of file
diff --git a/bindings/julia/src/TriangularMap.cpp b/bindings/julia/src/TriangularMap.cpp
index 7ac54973..defb643d 100644
--- a/bindings/julia/src/TriangularMap.cpp
+++ b/bindings/julia/src/TriangularMap.cpp
@@ -14,6 +14,6 @@ void mpart::binding::TriangularMapWrapper(jlcxx::Module &mod) {
        .method("InverseInplace", [](TriangularMap<Kokkos::HostSpace> &map, jlcxx::ArrayRef<double,2> x, jlcxx::ArrayRef<double,2> r){
             map.InverseInplace(JuliaToKokkos(x), JuliaToKokkos(r));
        })
-       .method("GetComponent", &TriangularMap<Kokkos::HostSpace>::GetComponent)
+       .method("GetComponent", [](TriangularMap<Kokkos::HostSpace>& map, int i){ return map.GetComponent(i-1); })
     ;
 }
\ No newline at end of file
diff --git a/bindings/python/src/ConditionalMapBase.cpp b/bindings/python/src/ConditionalMapBase.cpp
index 22267d66..a44fc8d2 100644
--- a/bindings/python/src/ConditionalMapBase.cpp
+++ b/bindings/python/src/ConditionalMapBase.cpp
@@ -23,6 +23,12 @@ void mpart::binding::ConditionalMapBaseWrapper(py::module &m)
         .def("LogDeterminantCoeffGrad", static_cast<Eigen::RowMatrixXd (ConditionalMapBase<MemorySpace>::*)(Eigen::Ref<const Eigen::RowMatrixXd> const&)>(&ConditionalMapBase<MemorySpace>::LogDeterminantCoeffGrad))
         .def("LogDeterminantInputGrad", static_cast<Eigen::RowMatrixXd (ConditionalMapBase<MemorySpace>::*)(Eigen::Ref<const Eigen::RowMatrixXd> const&)>(&ConditionalMapBase<MemorySpace>::LogDeterminantInputGrad))
         .def("GetBaseFunction", &ConditionalMapBase<MemorySpace>::GetBaseFunction)
+        .def("__getitem__", [](ConditionalMapBase<MemorySpace> &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 8ab340c9..f3a3ceff 100644
--- a/bindings/python/src/TriangularMap.cpp
+++ b/bindings/python/src/TriangularMap.cpp
@@ -20,6 +20,8 @@ void mpart::binding::TriangularMapWrapper(py::module &m)
     py::class_<TriangularMap<MemorySpace>, ConditionalMapBase<MemorySpace>, std::shared_ptr<TriangularMap<MemorySpace>>>(m, tName.c_str())
         .def(py::init<std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>>, bool>(), py::arg("comps"), py::arg("moveCoeffs") = false)
         .def("GetComponent", &TriangularMap<MemorySpace>::GetComponent)
+        .def("Slice", &TriangularMap<MemorySpace>::Slice)
+        .def("__getitem__", &TriangularMap<MemorySpace>::GetComponent)
         ;
 
 }
diff --git a/src/AffineMap.cpp b/src/AffineMap.cpp
index 90b81e9d..d5427991 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<typename MemorySpace>
-AffineMap<MemorySpace>::AffineMap(StridedVector<double,MemorySpace> b) : ConditionalMapBase<MemorySpace>(b.size(),b.size(),0), 
-                                                                  b_("b",b.layout()), 
+AffineMap<MemorySpace>::AffineMap(StridedVector<double,MemorySpace> b) : ConditionalMapBase<MemorySpace>(b.size(),b.size(),0),
+                                                                  b_("b",b.layout()),
                                                                   logDet_(0.0)
 {
     Kokkos::deep_copy(b_, b);
@@ -30,7 +30,7 @@ AffineMap<MemorySpace>::AffineMap(StridedMatrix<double,MemorySpace> A,
                                   StridedVector<double,MemorySpace> b) : ConditionalMapBase<MemorySpace>(A.extent(1),A.extent(0),0),
                                                                   A_("A", A.layout()),
                                                                   b_("b",b.layout())
-{   
+{
     Kokkos::deep_copy(A_, A);
     Kokkos::deep_copy(b_, b);
     assert(A_.extent(0)<=A_.extent(1));
@@ -91,7 +91,7 @@ void AffineMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpace
     }else{
         outLeft = output;
     }
-    
+
     // Bias part
     if(b_.size()>0){
 
@@ -99,7 +99,7 @@ void AffineMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpace
 
         // After this for loop, we will have out = r - b
         Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(const int& j, const int& i) {
-            outLeft(i,j) = r(i,j) - b_(i); 
+            outLeft(i,j) = r(i,j) - b_(i);
         });
 
     }else{
@@ -115,7 +115,7 @@ void AffineMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpace
         // If the matrix is rectangular, treat it as the lower part of a block triangular matrix
         // The value of x1 contains the compute inverse for the first block
         if(nrows != ncols){
-            
+
             StridedMatrix<const double, MemorySpace> xsub = Kokkos::subview(x1, std::make_pair(0,ncols-nrows), Kokkos::ALL());
             StridedMatrix<const double, MemorySpace> Asub = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(0,ncols-nrows));
 
@@ -129,13 +129,13 @@ void AffineMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpace
         if(copyOut)
             Kokkos::deep_copy(output, outLeft);
     }
-    
+
 }
 
 
 
 template<typename MemorySpace>
-void AffineMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+void AffineMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                          StridedMatrix<double, MemorySpace>              output)
 {
     return;
@@ -151,7 +151,7 @@ void AffineMap<MemorySpace>::EvaluateImpl(StridedMatrix<const double, MemorySpac
     }else{
         Kokkos::deep_copy(output, pts);
     }
-    
+
     // Bias part
     if(b_.size()>0){
 
@@ -165,22 +165,22 @@ void AffineMap<MemorySpace>::EvaluateImpl(StridedMatrix<const double, MemorySpac
 }
 
 template<typename MemorySpace>
-void AffineMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void AffineMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                            StridedMatrix<const double, MemorySpace> const& sens,
                                            StridedMatrix<double, MemorySpace>              output)
 {
     return;
 }
 
-template<typename MemorySpace> 
-void AffineMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+template<typename MemorySpace>
+void AffineMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                          StridedMatrix<double, MemorySpace>              output)
 {
     return;
 }
 
 template<typename MemorySpace>
-void AffineMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void AffineMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                           StridedMatrix<const double, MemorySpace> const& sens,
                                           StridedMatrix<double, MemorySpace>              output)
 {
@@ -192,8 +192,16 @@ void AffineMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpac
     }
 }
 
+template<typename MemorySpace>
+std::shared_ptr<ConditionalMapBase<MemorySpace>> AffineMap<MemorySpace>::Slice(int a, int b) {
+    StridedMatrix<double, MemorySpace> A_new;
+    StridedVector<double, MemorySpace> b_new;
+    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<AffineMap<MemorySpace>>(A_new, b_new);
+}
 
 template class mpart::AffineMap<Kokkos::HostSpace>;
 #if defined(MPART_ENABLE_GPU)
     template class mpart::AffineMap<mpart::DeviceSpace>;
-#endif 
+#endif
diff --git a/src/ComposedMap.cpp b/src/ComposedMap.cpp
index c2eb0aca..381a1a98 100644
--- a/src/ComposedMap.cpp
+++ b/src/ComposedMap.cpp
@@ -8,7 +8,7 @@
 using namespace mpart;
 
 template<typename MemorySpace>
-ComposedMap<MemorySpace>::Checkpointer::Checkpointer(unsigned int maxSaves, 
+ComposedMap<MemorySpace>::Checkpointer::Checkpointer(unsigned int maxSaves,
                                                      StridedMatrix<const double, MemorySpace> initialPts,
                                                      std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>>& comps) : maxSaves_(maxSaves),
                                                                                                                             maps_(comps)
@@ -23,7 +23,7 @@ ComposedMap<MemorySpace>::Checkpointer::Checkpointer(unsigned int maxSaves,
 
     checkpointLayers_.push_back(0);
 
-    // Allocate the workspace 
+    // Allocate the workspace
     workspace1_ = Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>("xi", initialPts.extent(0), initialPts.extent(1));
     workspace2_ = Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>("xj", initialPts.extent(0), initialPts.extent(1));
 }
@@ -37,7 +37,7 @@ int ComposedMap<MemorySpace>::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<MemorySpace>::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<MemorySpace>::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<MemorySpace>::Checkpointer::GetNextCheckpoint(unsigned int layer
 
 template<typename MemorySpace>
 Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace> ComposedMap<MemorySpace>::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<double**, Kokkos::LayoutLeft, MemorySpace> ComposedMap<MemorySpace>
         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<layerInd; ++i){
-            
+
             if(i==nextCheckLayer){
-                
+
                 checkpoints_.push_back(Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>("xi", workspace1_.extent(0), workspace1_.extent(1)));
                 Kokkos::deep_copy(checkpoints_.back(),workspace1_);
 
@@ -133,7 +133,7 @@ Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace> ComposedMap<MemorySpace>
 
 template<typename MemorySpace>
 ComposedMap<MemorySpace>::ComposedMap(std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> const& maps, bool moveCoeffs,
-                                      int maxChecks) : ConditionalMapBase<MemorySpace>(maps.front()->inputDim, 
+                                      int maxChecks) : ConditionalMapBase<MemorySpace>(maps.front()->inputDim,
                                                                                        maps.front()->inputDim,
                                                                                        std::accumulate(maps.begin(), maps.end(), 0, [](size_t sum, std::shared_ptr<ConditionalMapBase<MemorySpace>> const& comp){ return sum + comp->numCoeffs; })),
                                                         maps_(maps),
@@ -159,7 +159,7 @@ ComposedMap<MemorySpace>::ComposedMap(std::vector<std::shared_ptr<ConditionalMap
 
         Kokkos::View<double*,MemorySpace> coeffs("coeffs", this->numCoeffs);
         unsigned int cumNumCoeffs = 0;
-        
+
         for(unsigned int i=0; i<maps_.size(); ++i){
 
             if(!maps_.at(i)->CheckCoefficients()){
@@ -241,14 +241,14 @@ void ComposedMap<MemorySpace>::LogDeterminantImpl(StridedMatrix<const double, Me
     if(maps_.size()==1)
         return;
 
-    // intermediate points and variable to hold logdet increments 
+    // intermediate points and variable to hold logdet increments
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intPts1("intermediate points 1", pts.extent(0), pts.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intPts2("intermediate points 2", pts.extent(0), pts.extent(1));
     Kokkos::deep_copy(intPts1, pts);
 
     Kokkos::View<double*, Kokkos::LayoutLeft, MemorySpace> compDetIncrement("Log Determinant", output.extent(0));
     for(unsigned int i=1; i<maps_.size(); ++i){
-        
+
         // Compute x_i = T_{i-1}(x_{i-1})
         maps_.at(i-1)->EvaluateImpl(intPts1, intPts2);
 
@@ -271,7 +271,7 @@ void ComposedMap<MemorySpace>::EvaluateImpl(StridedMatrix<const double, MemorySp
     // intermediate output
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intPts1("intermediate points 1", pts.extent(0), pts.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intPts2("intermediate points 2", pts.extent(0), pts.extent(1));
-    
+
     maps_.at(0)->EvaluateImpl(pts, intPts1);
 
     for(int j = 1; j<maps_.size(); j++){
@@ -283,7 +283,7 @@ void ComposedMap<MemorySpace>::EvaluateImpl(StridedMatrix<const double, MemorySp
 }
 
 template<typename MemorySpace>
-void ComposedMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void ComposedMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                             StridedMatrix<const double, MemorySpace> const& sens,
                                             StridedMatrix<double, MemorySpace>              output)
 {
@@ -296,14 +296,14 @@ void ComposedMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySp
 
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intSens1("intermediate Sens", sens.extent(0), sens.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  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);
@@ -320,7 +320,7 @@ void ComposedMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpa
                                            StridedMatrix<double, MemorySpace>              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<double**, Kokkos::LayoutLeft, MemorySpace> intR1("intermediate r1", r.extent(0), r.extent(1));
@@ -330,7 +330,7 @@ void ComposedMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpa
     Kokkos::deep_copy(intR1, r);
 
     for(int i = maps_.size() - 1; i>=0; --i){
-        
+
         maps_.at(i)->InverseImpl(x1, intR1, intR2); // <- Note the x1 doesn't matter here because the layer is square
         simple_swap<decltype(intR1)>(intR1,intR2);
     }
@@ -339,13 +339,13 @@ void ComposedMap<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpa
 }
 
 template<typename MemorySpace>
-void ComposedMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void ComposedMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                              StridedMatrix<const double, MemorySpace> const& sens,
                                              StridedMatrix<double, MemorySpace>              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<double**, Kokkos::LayoutLeft, MemorySpace>  intSens1("intermediate sens 1", sens.extent(0), sens.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intSens2("intermediate sens 2", sens.extent(0), sens.extent(1));
     Kokkos::deep_copy(intSens1, sens);
@@ -353,15 +353,15 @@ void ComposedMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemoryS
     Checkpointer checker(maxChecks_, pts, maps_);
 
     StridedMatrix<double, MemorySpace> 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());
 
 
@@ -381,43 +381,43 @@ void ComposedMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemoryS
 
 
 template<typename MemorySpace>
-void ComposedMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+void ComposedMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                            StridedMatrix<double, MemorySpace>              output)
 {
     StridedMatrix<double, MemorySpace> subOut;
 
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intSens1("intermediate Sens", pts.extent(0), pts.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  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<double**, Kokkos::LayoutLeft, MemorySpace> subOut2("temp", maps_.at(i)->numCoeffs, pts.extent(1));
         maps_.at(i)->CoeffGradImpl(input, intSens1, subOut2);
         subOut += subOut2;
@@ -428,31 +428,31 @@ void ComposedMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const d
             simple_swap<decltype(intSens1)>(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<typename MemorySpace>
-void ComposedMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+void ComposedMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                            StridedMatrix<double, MemorySpace>              output)
 {
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  intSens1("intermediate Sens", pts.extent(0), pts.extent(1));
     Kokkos::View<double**, Kokkos::LayoutLeft, MemorySpace>  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);
 
@@ -460,13 +460,27 @@ void ComposedMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const d
         maps_.at(i)->GradientImpl(input, intSens1, intSens2);
         simple_swap<decltype(intSens1)>(intSens1, intSens2);
 
-        maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2);  
-        intSens1 += intSens2;      
+        maps_.at(i)->LogDeterminantInputGradImpl(input, intSens2);
+        intSens1 += intSens2;
     }
 
     Kokkos::deep_copy(output, intSens1);
 }
 
+template<typename MemorySpace>
+std::shared_ptr<ConditionalMapBase<MemorySpace>> ComposedMap<MemorySpace>::Slice(int a, int b) {
+    assert(a < b && 0 <= a && b <= this->outputDim);
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> lastMap = maps_[maps_.size()-1];
+    if (lastMap->outputDim == b-a) {
+        return std::shared_ptr<ComposedMap<MemorySpace>>(this);
+    } else {
+        std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> newMaps;
+        newMaps = maps_;
+        newMaps[maps_.size()-1] = lastMap->Slice(a, b);
+        return std::make_shared<ComposedMap<MemorySpace>>(newMaps);
+    }
+}
+
 // Explicit template instantiation
 template class mpart::ComposedMap<Kokkos::HostSpace>;
 #if defined(MPART_ENABLE_GPU)
diff --git a/src/ConditionalMapBase.cpp b/src/ConditionalMapBase.cpp
index fcf0193a..ff09ddc0 100644
--- a/src/ConditionalMapBase.cpp
+++ b/src/ConditionalMapBase.cpp
@@ -40,27 +40,27 @@ Kokkos::View<double*, mpart::DeviceSpace> ConditionalMapBase<mpart::DeviceSpace>
 template<>
 template<>
 Kokkos::View<double*, Kokkos::HostSpace> ConditionalMapBase<mpart::DeviceSpace>::LogDeterminant(StridedMatrix<const double, Kokkos::HostSpace> const& pts)
-{   
+{
     return ToHost( this->LogDeterminant( ToDevice<mpart::DeviceSpace>(pts) ));
 }
 
 template<>
 template<>
 Kokkos::View<double*, mpart::DeviceSpace> ConditionalMapBase<Kokkos::HostSpace>::LogDeterminant(StridedMatrix<const double, mpart::DeviceSpace> const& pts)
-{   
+{
     return ToDevice<mpart::DeviceSpace>( this->LogDeterminant( ToHost(pts) ));
 }
 
 template<>
 Eigen::VectorXd ConditionalMapBase<mpart::DeviceSpace>::LogDeterminant(Eigen::Ref<const Eigen::RowMatrixXd> const& pts)
-{   
+{
     StridedMatrix<const double,Kokkos::HostSpace> pts_view = ConstRowMatToKokkos<double,Kokkos::HostSpace>(pts);
     Kokkos::View<double*, Kokkos::HostSpace> outView = ToHost( this->LogDeterminant( ToDevice<mpart::DeviceSpace>(pts_view) ));
     return KokkosToVec<double>(outView);
 }
 
 
-#endif 
+#endif
 
 template<>
 template<>
@@ -138,7 +138,7 @@ Eigen::RowMatrixXd ConditionalMapBase<mpart::DeviceSpace>::Inverse(Eigen::Ref<co
     return KokkosToMat( ToHost( Inverse(xView, rView) ));
 }
 
-#endif 
+#endif
 
 template<>
 template<>
@@ -177,10 +177,10 @@ template<>
 template<>
 StridedMatrix<double, Kokkos::HostSpace> ConditionalMapBase<mpart::DeviceSpace>::LogDeterminantCoeffGrad(StridedMatrix<const double, Kokkos::HostSpace> const& pts)
 {
-    // Copy the points to the device space 
+    // Copy the points to the device space
     StridedMatrix<const double, mpart::DeviceSpace> pts_device = ToDevice<mpart::DeviceSpace>(pts);
 
-    // Evaluate on the device space 
+    // Evaluate on the device space
     StridedMatrix<double, mpart::DeviceSpace> evals_device = this->LogDeterminantCoeffGrad(pts_device);
 
     // Copy back to the host space
@@ -191,10 +191,10 @@ template<>
 template<>
 StridedMatrix<double, mpart::DeviceSpace> ConditionalMapBase<Kokkos::HostSpace>::LogDeterminantCoeffGrad(StridedMatrix<const double, mpart::DeviceSpace> const& pts)
 {
-    // Copy the points to the host 
+    // Copy the points to the host
     StridedMatrix<const double, Kokkos::HostSpace> pts_host = ToHost(pts);
 
-    // Evaluate on the host 
+    // Evaluate on the host
     StridedMatrix<double, Kokkos::HostSpace> evals_host = this->LogDeterminantCoeffGrad(pts_host);
 
     // Copy back to the device
@@ -254,10 +254,10 @@ template<>
 template<>
 StridedMatrix<double, Kokkos::HostSpace> ConditionalMapBase<mpart::DeviceSpace>::LogDeterminantInputGrad(StridedMatrix<const double, Kokkos::HostSpace> const& pts)
 {
-    // Copy the points to the device space 
+    // Copy the points to the device space
     StridedMatrix<const double, mpart::DeviceSpace> pts_device = ToDevice<mpart::DeviceSpace,const double>(pts);
 
-    // Evaluate on the device space 
+    // Evaluate on the device space
     StridedMatrix<double, mpart::DeviceSpace> evals_device = this->LogDeterminantInputGrad(pts_device);
 
     // Copy back to the host space
@@ -268,10 +268,10 @@ template<>
 template<>
 StridedMatrix<double, mpart::DeviceSpace> ConditionalMapBase<Kokkos::HostSpace>::LogDeterminantInputGrad(StridedMatrix<const double, mpart::DeviceSpace> const& pts)
 {
-    // Copy the points to the host 
+    // Copy the points to the host
     StridedMatrix<const double, Kokkos::HostSpace> pts_host = ToHost(pts);
 
-    // Evaluate on the host 
+    // Evaluate on the host
     StridedMatrix<double, Kokkos::HostSpace> evals_host = this->LogDeterminantInputGrad(pts_host);
 
     // Copy back to the device
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<MemorySpace>::InverseImpl(StridedMatrix<const double, MemorySpa
 
 
 template<typename MemorySpace>
-void IdentityMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void IdentityMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                StridedMatrix<const double, MemorySpace> const& sens,
                                                StridedMatrix<double, MemorySpace>              output)
 {
@@ -57,7 +57,7 @@ void IdentityMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, MemoryS
 }
 
 template<typename MemorySpace>
-void IdentityMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+void IdentityMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
                             StridedMatrix<const double, MemorySpace> const& sens,
                             StridedMatrix<double, MemorySpace>              output)
 {
@@ -77,16 +77,16 @@ void IdentityMap<MemorySpace>::GradientImpl(StridedMatrix<const double, MemorySp
 }
 
 template<typename MemorySpace>
-void IdentityMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+void IdentityMap<MemorySpace>::LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                                              StridedMatrix<double, MemorySpace>              output)
 {
     assert(false);
 }
 
 template<typename MemorySpace>
-void IdentityMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+void IdentityMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                             StridedMatrix<double, MemorySpace>              output)
-{   
+{
     Kokkos::MDRangePolicy<Kokkos::Rank<2>, typename MemoryToExecution<MemorySpace>::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<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const d
 
 }
 
+template<typename MemorySpace>
+std::shared_ptr<ConditionalMapBase<MemorySpace>> IdentityMap<MemorySpace>::Slice(int a, int b) {
+    assert(a >= 0);
+    assert(b <= this->outputDim);
+    assert(a < b);
+    return std::make_shared<IdentityMap<MemorySpace>>(this->inputDim, b-a);
+}
+
 // Explicit template instantiation
 template class mpart::IdentityMap<Kokkos::HostSpace>;
 #if defined(KOKKOS_ENABLE_CUDA ) || defined(KOKKOS_ENABLE_SYCL)
diff --git a/src/MapFactory.cpp b/src/MapFactory.cpp
index 70df3ec0..47dc1032 100644
--- a/src/MapFactory.cpp
+++ b/src/MapFactory.cpp
@@ -100,9 +100,9 @@ std::shared_ptr<ConditionalMapBase<MemorySpace>> mpart::MapFactory::CreateSingle
 
 template<typename MemorySpace>
 std::shared_ptr<ConditionalMapBase<MemorySpace>> mpart::MapFactory::CreateTriangular(unsigned int inputDim,
-                                                                         unsigned int outputDim,
-                                                                         unsigned int totalOrder,
-                                                                         MapOptions options)
+                                                                                     unsigned int outputDim,
+                                                                                     unsigned int totalOrder,
+                                                                                     MapOptions options)
 {
 
     std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> comps(outputDim);
diff --git a/src/TriangularMap.cpp b/src/TriangularMap.cpp
index 5203719e..f3126db6 100644
--- a/src/TriangularMap.cpp
+++ b/src/TriangularMap.cpp
@@ -12,10 +12,10 @@ TriangularMap<MemorySpace>::TriangularMap(std::vector<std::shared_ptr<Conditiona
                         std::accumulate(components.begin(), components.end(), 0, [](size_t sum, std::shared_ptr<ConditionalMapBase<MemorySpace>> const& comp){ return sum + comp->numCoeffs; })),
                         comps_(components)
 {
-    
+
 
     // Check the sizes of all the inputs
-    
+
     for(unsigned int i=0; i<comps_.size(); ++i){
         if(comps_.at(i)->outputDim > comps_.at(i)->inputDim){
             std::stringstream msg;
@@ -42,7 +42,7 @@ TriangularMap<MemorySpace>::TriangularMap(std::vector<std::shared_ptr<Conditiona
 
         Kokkos::View<double*,MemorySpace> coeffs("coeffs", this->numCoeffs);
         unsigned int cumNumCoeffs = 0;
-        
+
         for(unsigned int i=0; i<comps_.size(); ++i){
 
             if(!comps_.at(i)->CheckCoefficients()){
@@ -281,7 +281,7 @@ void TriangularMap<MemorySpace>::CoeffGradImpl(StridedMatrix<const double, Memor
             subOut = Kokkos::subview(output, std::make_pair(startParamDim,int(startParamDim+comps_.at(i)->numCoeffs)), Kokkos::ALL());
             comps_.at(i)->CoeffGradImpl(subPts, subSens, subOut);
 
-            
+
             startParamDim += comps_.at(i)->numCoeffs;
         }
 
@@ -345,6 +345,63 @@ void TriangularMap<MemorySpace>::LogDeterminantInputGradImpl(StridedMatrix<const
     }
 }
 
+template<typename MemorySpace>
+std::shared_ptr<ConditionalMapBase<MemorySpace>> TriangularMap<MemorySpace>::Slice(int a, int b) {
+    std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> components;
+    // TODO: Handle empty case
+    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
+    //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; // 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)
+            break;
+        accum_b += this->comps_[k_b]->outputDim;
+        k_b++;
+    }
+
+    if(k_a == k_b){
+        components.push_back(this->comps_[k_a]->Slice(a-accum_a, b-accum_b));
+    } else {
+        components = std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>>(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<TriangularMap<MemorySpace>>(components);
+    output->SetCoeffs(Kokkos::View<double*,MemorySpace>("Component Coefficients", output->numCoeffs));
+    return output;
+}
+
+template<typename MemorySpace>
+std::shared_ptr<ConditionalMapBase<MemorySpace>> TriangularMap<MemorySpace>::BlockSlice(int a, int b) {
+    std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> components;
+    for(int k = a; k < b; k++){
+        components.push_back(this->comps_[k]);
+    }
+    return std::make_shared<TriangularMap<MemorySpace>>(components);
+}
+
 // Explicit template instantiation
 template class mpart::TriangularMap<Kokkos::HostSpace>;
 #if defined(MPART_ENABLE_GPU)
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<double*, Kokkos::HostSpace> 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<double**, Kokkos::HostSpace> 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<double**, Kokkos::HostSpace> 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<AffineMap<Kokkos::HostSpace>>(A);
@@ -139,7 +139,7 @@ TEST_CASE( "Testing Full AffineMap", "[FullAffineMap]" ) {
 
     Kokkos::View<double**, Kokkos::HostSpace> A("A", 2,2);
     Kokkos::View<double*, Kokkos::HostSpace> 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<double**, Kokkos::HostSpace> A("A", 3,3);
+    Kokkos::View<double*, Kokkos::HostSpace> 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<AffineMap<Kokkos::HostSpace>>(A,b);
+    auto slice = map->Slice(0,2);
+    REQUIRE(slice->inputDim==3);
+    REQUIRE(slice->outputDim==2);
+
+    unsigned int numPts = 3;
+    Kokkos::View<double**, Kokkos::HostSpace> pts("Point", A.extent(1), numPts);
+    for(unsigned int i=0; i<numPts; ++i){
+        for(unsigned int d=0; d<A.extent(0);++d)
+            pts(d,i) = double(i)/double(numPts-1);
+    }
+
+    // Check to make sure the log determinant is constant and what it should be
+    Kokkos::View<double*, Kokkos::HostSpace> logDet = slice->LogDeterminant(pts);
+    for(unsigned int i=0; i<numPts; ++i){
+        CHECK(logDet(i) == Approx(std::log(A(0,1)*A(1,2)-A(0,2)*A(1,1))).epsilon(1e-14));
+    }
+
+    // Test the forward evaluation
+    Kokkos::View<double**, Kokkos::HostSpace> evals = slice->Evaluate(pts);
+    REQUIRE(evals.extent(0)==2);
+    REQUIRE(evals.extent(1)==pts.extent(1));
+
+    Kokkos::View<double**, Kokkos::HostSpace> pts2 = slice->Inverse(pts, evals);
+    REQUIRE(pts2.extent(0)==2);
+    REQUIRE(pts2.extent(1)==pts.extent(1));
+
+    for(unsigned int i=0; i<numPts; ++i){
+        double trueOut1 = A(0,0)*pts(0,i) + A(0,1)*pts(1,i) + A(0,2)*pts(2,i);
+        double trueOut2 = A(1,0)*pts(0,i) + A(1,1)*pts(1,i) + A(1,2)*pts(2,i);
+        CHECK(evals(0,i)==Approx(trueOut1).epsilon(1e-14));
+        CHECK(evals(1,i)==Approx(trueOut2).epsilon(1e-14));
+        CHECK(pts2(0,i)==Approx(pts(1,i)).epsilon(1e-14));
+        CHECK(pts2(1,i)==Approx(pts(2,i)).epsilon(1e-14));
+    }
+}
+
 
 #if defined(MPART_ENABLE_GPU)
 
 TEST_CASE( "Testing Shift-only AffineMap on Device", "[DeviceShiftMap]" ) {
 
     Kokkos::View<double*, Kokkos::HostSpace> 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<AffineMap<mpart::DeviceSpace>>(db);
 
-    
+
     unsigned int numPts = 10;
     Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> hpts("Point", 2, numPts);
     for(unsigned int i=0; i<numPts; ++i){
@@ -230,7 +287,7 @@ TEST_CASE( "Testing Shift-only AffineMap on Device", "[DeviceShiftMap]" ) {
 TEST_CASE( "Testing Linear-only AffineMap on Device", "[DeviceLinearMap]" ) {
 
     Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> 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<numPts; ++i){
         double trueOut1 = hA(0,0)*hpts(0,i) + hA(0,1)*hpts(1,i);
         double trueOut2 = hA(1,0)*hpts(0,i) + hA(1,1)*hpts(1,i);
@@ -279,8 +336,8 @@ TEST_CASE( "Testing Linear-only AffineMap on Device", "[DeviceLinearMap]" ) {
 TEST_CASE( "Testing Rectangular AffineMap on Device", "[DeviceRectangularLinearMap]" ) {
 
     Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> 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<mpart::DeviceSpace>(hA);
     auto map = std::make_shared<AffineMap<mpart::DeviceSpace>>(dA);
 
@@ -338,7 +395,7 @@ TEST_CASE( "Testing Full AffineMap on Device", "[DeviceFullAffineMap]" ) {
 
     Kokkos::View<double**, Kokkos::LayoutLeft, Kokkos::HostSpace> hA("A", 2,2);
     Kokkos::View<double*, Kokkos::HostSpace> 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<mpart::DeviceSpace>(hA);
     auto db = ToDevice<mpart::DeviceSpace>(hb);
-    
+
     auto map = std::make_shared<AffineMap<mpart::DeviceSpace>>(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<numPts; ++i){
         double trueOut1 = hA(0,0)*hpts(0,i) + hA(0,1)*hpts(1,i) + hb(0);
         double trueOut2 = hA(1,0)*hpts(0,i) + hA(1,1)*hpts(1,i) + hb(1);
@@ -387,4 +444,4 @@ TEST_CASE( "Testing Full AffineMap on Device", "[DeviceFullAffineMap]" ) {
 }
 
 
-#endif 
+#endif
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<MemorySpace>{
 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 double, MemorySpace> const& pts,
-                              StridedMatrix<double, MemorySpace>              output) override{Kokkos::deep_copy(output,pts);};
+                              StridedMatrix<double, MemorySpace>              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 double, MemorySpace> const& pts,  
+    virtual void GradientImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                StridedMatrix<const double, MemorySpace> const& sens,
                                StridedMatrix<double, MemorySpace>              output) override
     {
-        assert(false);  
+        assert(false);
     }
 
     virtual void LogDeterminantImpl(StridedMatrix<const double, MemorySpace> const&,
@@ -32,25 +36,30 @@ class MyIdentityMap : public ConditionalMapBase<MemorySpace>{
                              StridedMatrix<const double, MemorySpace> const& r,
                              StridedMatrix<double, MemorySpace>              output) override{Kokkos::deep_copy(output,r);};
 
-    virtual void CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,  
+    virtual void CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                StridedMatrix<const double, MemorySpace> const& sens,
                                StridedMatrix<double, MemorySpace>              output) override
     {
-        assert(false);  
+        assert(false);
     }
 
 
-    virtual void LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+    virtual void LogDeterminantCoeffGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                              StridedMatrix<double, MemorySpace>              output) override
-    {   
+    {
         assert(false);
     }
 
-    virtual void LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts, 
+    virtual void LogDeterminantInputGradImpl(StridedMatrix<const double, MemorySpace> const& pts,
                                              StridedMatrix<double, MemorySpace>              output) override
-    {   
+    {
         assert(false);
     }
+
+    // Only creates a slice of the tail of the input
+    std::shared_ptr<ConditionalMapBase<MemorySpace>> Slice(int a, int b) override{
+        return std::make_shared<MyIdentityMap>(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<ConditionalMapBase<Kokkos::HostSpace>> 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<double**, Kokkos::HostSpace> pts("pts", dim, numPts);
+
+        for(unsigned int i=0; i<dim; ++i){
+            for(unsigned int j=0; j<numPts; ++j){
+                pts(i,j) = i;
+            }
+        }
+
+        Kokkos::View<double**, Kokkos::HostSpace> 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<sliceOutdim; ++i){
+            for(unsigned int j=0; j<numPts; ++j){
+                CHECK(output(i,j) == idx_start+i);
+            }
+        }
+    }
 }
\ No newline at end of file
diff --git a/tests/Test_TriangularMap.cpp b/tests/Test_TriangularMap.cpp
index 07b8ed1f..0b146ac2 100644
--- a/tests/Test_TriangularMap.cpp
+++ b/tests/Test_TriangularMap.cpp
@@ -12,7 +12,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
     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 with moveCoeffs=fa
         blocks.at(i) = MapFactory::CreateComponent<MemorySpace>(mset, options);
     }
 
-    std::shared_ptr<ConditionalMapBase<MemorySpace>> triMap = std::make_shared<TriangularMap<MemorySpace>>(blocks);
+    std::shared_ptr<TriangularMap<MemorySpace>> triMap = std::make_shared<TriangularMap<MemorySpace>>(blocks);
 
     CHECK(triMap->outputDim == numBlocks);
     CHECK(triMap->inputDim == numBlocks+extraInputs);
@@ -37,9 +37,9 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
     Kokkos::View<double*,Kokkos::HostSpace> coeffs("Coefficients", triMap->numCoeffs);
     for(unsigned int i=0; i<triMap->numCoeffs; ++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 with moveCoeffs=fa
 
     triMap->SetCoeffs(coeffs);
     auto out = triMap->Evaluate(in);
-    
+
     SECTION("Evaluation"){
 
         for(unsigned int i=0; i<numBlocks; ++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);
@@ -134,16 +134,16 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
             coeffs(i) -= fdstep;
         }
-        
+
     }
 
 
@@ -173,18 +173,18 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
                 in(i,ptInd) -= fdstep;
         }
-        
+
     }
 
     SECTION("LogDeterminantCoeffGrad"){
@@ -199,7 +199,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
         Kokkos::View<double**,Kokkos::HostSpace> detGrad = triMap->LogDeterminantCoeffGrad(in);
         REQUIRE(detGrad.extent(0)==triMap->numCoeffs);
         REQUIRE(detGrad.extent(1)==numSamps);
-        
+
         Kokkos::View<double*,Kokkos::HostSpace> logDet = triMap->LogDeterminant(in);
         Kokkos::View<double*,Kokkos::HostSpace> logDet2;
 
@@ -212,8 +212,8 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=fa
             logDet2 = triMap->LogDeterminant(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
-                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4)); 
-            
+                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4));
+
             coeffs(i) -= fdstep;
         }
 
@@ -226,7 +226,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
     MapOptions options;
     options.basisType = BasisTypes::ProbabilistHermite;
     options.basisNorm = false;
-    
+
     unsigned int numBlocks = 3;
     unsigned int maxDegree = 2;
     unsigned int extraInputs = 1;
@@ -238,7 +238,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
     std::vector<std::shared_ptr<ConditionalMapBase<MemorySpace>>> blocks(numBlocks);
     std::vector<Kokkos::View<double*,Kokkos::HostSpace>> coeffs_(numBlocks);
     for(unsigned int i=0;i<numBlocks;++i){
-        
+
         FixedMultiIndexSet<MemorySpace> mset(i+extraInputs+1,maxDegree);
         coeffSize += mset.Size();
 
@@ -248,9 +248,9 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
 
         for(unsigned int j=0; j<blocks.at(i)->numCoeffs; ++j)
             coeffs_.at(i)(j) = 0.1*(j+1);
-        
+
         blocks.at(i)->SetCoeffs(coeffs_.at(i));
-        
+
     }
     bool moveCoeffs=true;
     std::shared_ptr<ConditionalMapBase<MemorySpace>> triMap = std::make_shared<TriangularMap<MemorySpace>>(blocks, moveCoeffs);
@@ -269,7 +269,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
         for(unsigned int i=0; i<numBlocks; ++i){
             for(unsigned int j=0; j<blocks.at(i)->numCoeffs; ++j){
                 CHECK(blocks.at(i)->Coeffs()(j) == coeffs(cumCoeffInd)); // Values of coefficients should be equal
-                CHECK(blocks.at(i)->Coeffs()(j) == coeffs_.at(i)(j)); 
+                CHECK(blocks.at(i)->Coeffs()(j) == coeffs_.at(i)(j));
                 CHECK(&blocks.at(i)->Coeffs()(j) == &coeffs(cumCoeffInd)); // Memory location should also be the same (no copy)
                 cumCoeffInd++;
             }
@@ -286,13 +286,13 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
     }
 
 
-    
+
     auto out = triMap->Evaluate(in);
-    
+
     SECTION("Evaluation"){
 
         for(unsigned int i=0; i<numBlocks; ++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);
@@ -357,16 +357,16 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
             coeffs(i) -= fdstep;
         }
-        
+
     }
 
 
@@ -396,18 +396,18 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
                 in(i,ptInd) -= fdstep;
         }
-        
+
     }
 
     SECTION("LogDeterminantCoeffGrad"){
@@ -422,7 +422,7 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
         Kokkos::View<double**,Kokkos::HostSpace> detGrad = triMap->LogDeterminantCoeffGrad(in);
         REQUIRE(detGrad.extent(0)==triMap->numCoeffs);
         REQUIRE(detGrad.extent(1)==numSamps);
-        
+
         Kokkos::View<double*,Kokkos::HostSpace> logDet = triMap->LogDeterminant(in);
         Kokkos::View<double*,Kokkos::HostSpace> logDet2;
 
@@ -435,8 +435,8 @@ TEST_CASE( "Testing 3d triangular map from MonotoneComponents with moveCoeffs=tr
             logDet2 = triMap->LogDeterminant(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
-                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4)); 
-            
+                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4));
+
             coeffs(i) -= fdstep;
         }
 
@@ -477,9 +477,9 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
     Kokkos::View<double*,Kokkos::HostSpace> coeffs("Coefficients", triMap->numCoeffs);
     for(unsigned int i=0; i<triMap->numCoeffs; ++i)
         coeffs(i) = 0.1*(i+1);
-        
+
     SECTION("Coefficients"){
-        
+
         // Set the coefficients of the triangular map
         triMap->SetCoeffs(coeffs);
 
@@ -505,18 +505,18 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
 
     triMap->SetCoeffs(coeffs);
     auto out = triMap->Evaluate(in);
-    
+
     SECTION("Evaluation"){
 
         unsigned int start = 0;
 
         for(unsigned int b=0; b<numBlocks; ++b){
-            
+
             auto subIn = Kokkos::subview(in, std::make_pair(0, int(blocks.at(b)->inputDim)), Kokkos::ALL());
             auto subOut = Kokkos::subview(out, std::make_pair(int(start), int(start + blocks.at(b)->outputDim)), Kokkos::ALL());
 
             auto outBlock = blocks.at(b)->Evaluate(subIn);
-            
+
             REQUIRE(outBlock.extent(1)==subOut.extent(1));
             REQUIRE(outBlock.extent(0)==subOut.extent(0));
 
@@ -584,16 +584,16 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
             coeffs(i) -= fdstep;
         }
-        
+
     }
 
 
@@ -623,18 +623,18 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++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<numSamps; ++ptInd)
                 in(i,ptInd) -= fdstep;
         }
-        
+
     }
 
     SECTION("LogDeterminantCoeffGrad"){
@@ -649,7 +649,7 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
         Kokkos::View<double**,Kokkos::HostSpace> detGrad = triMap->LogDeterminantCoeffGrad(in);
         REQUIRE(detGrad.extent(0)==triMap->numCoeffs);
         REQUIRE(detGrad.extent(1)==numSamps);
-        
+
         Kokkos::View<double*,Kokkos::HostSpace> logDet = triMap->LogDeterminant(in);
         Kokkos::View<double*,Kokkos::HostSpace> logDet2;
 
@@ -662,8 +662,8 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             logDet2 = triMap->LogDeterminant(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
-                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).epsilon(1e-4)); 
-            
+                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).epsilon(1e-4));
+
             coeffs(i) -= fdstep;
         }
 
@@ -710,16 +710,16 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
     CHECK(triMap->numCoeffs == coeffSize);
 
 
-        
+
     SECTION("Coefficients"){
-        
+
 
         // Now make sure that the coefficients of each block were set
         unsigned int cumCoeffInd = 0;
         for(unsigned int i=0; i<numBlocks; ++i){
             for(unsigned int j=0; j<blocks.at(i)->numCoeffs; ++j){
                 CHECK(blocks.at(i)->Coeffs()(j) == triMap->Coeffs()(cumCoeffInd)); // Values of coefficients should be equal
-                CHECK(blocks.at(i)->Coeffs()(j) == coeffs_.at(i)(j)); 
+                CHECK(blocks.at(i)->Coeffs()(j) == coeffs_.at(i)(j));
                 CHECK(&blocks.at(i)->Coeffs()(j) == &triMap->Coeffs()(cumCoeffInd)); // Memory location should also be the same (no copy)
                 cumCoeffInd++;
             }
@@ -736,18 +736,18 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
     }
 
     auto out = triMap->Evaluate(in);
-    
+
     SECTION("Evaluation"){
 
         unsigned int start = 0;
 
         for(unsigned int b=0; b<numBlocks; ++b){
-            
+
             auto subIn = Kokkos::subview(in, std::make_pair(0, int(blocks.at(b)->inputDim)), Kokkos::ALL());
             auto subOut = Kokkos::subview(out, std::make_pair(int(start), int(start + blocks.at(b)->outputDim)), Kokkos::ALL());
 
             auto outBlock = blocks.at(b)->Evaluate(subIn);
-            
+
             REQUIRE(outBlock.extent(1)==subOut.extent(1));
             REQUIRE(outBlock.extent(0)==subOut.extent(0));
 
@@ -815,16 +815,16 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
             coeffs(i) -= fdstep;
         }
-        
+
     }
 
 
@@ -854,18 +854,18 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++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<numSamps; ++ptInd)
                 in(i,ptInd) -= fdstep;
         }
-        
+
     }
 
     SECTION("LogDeterminantCoeffGrad"){
@@ -880,7 +880,7 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
         Kokkos::View<double**,Kokkos::HostSpace> detGrad = triMap->LogDeterminantCoeffGrad(in);
         REQUIRE(detGrad.extent(0)==triMap->numCoeffs);
         REQUIRE(detGrad.extent(1)==numSamps);
-        
+
         Kokkos::View<double*,Kokkos::HostSpace> logDet = triMap->LogDeterminant(in);
         Kokkos::View<double*,Kokkos::HostSpace> logDet2;
 
@@ -893,16 +893,36 @@ TEST_CASE( "Testing TriangularMap made from smaller TriangularMaps with moveCoef
             logDet2 = triMap->LogDeterminant(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
-                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).epsilon(1e-4)); 
-            
+                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).epsilon(1e-4));
+
             coeffs(i) -= fdstep;
         }
 
     }
 
-}
+    SECTION("Slice"){
+        int sliceBegin = 1;
+        int sliceEnd = 5;
+        unsigned int extraInputs = 2;
+        auto slice = triMap->Slice(sliceBegin, sliceEnd);
+        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; i<sliceEnd; ++i){
+            for(unsigned int j=0; j<numSamps; ++j) {
+                CHECK( outSlice(i-sliceBegin,j) == Approx(outFull(i,j)).epsilon(1e-6));
+            }
+        }
 
+        // This is used for testing the gradient
+        auto tmapSlice = std::dynamic_pointer_cast<TriangularMap<Kokkos::HostSpace>>(slice);
+    }
 
+}
 
 TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[TriangularMap_CreateSingleEntryMap]" ) {
 
@@ -925,9 +945,9 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
     Kokkos::View<double*,Kokkos::HostSpace> coeffs("Coefficients", triMap->numCoeffs);
     for(unsigned int i=0; i<triMap->numCoeffs; ++i)
         coeffs(i) = 0.1*(i+1);
-        
+
     SECTION("Coefficients"){
-        
+
         // Set the coefficients of the triangular map
         triMap->SetCoeffs(coeffs);
 
@@ -937,7 +957,7 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
                 CHECK(comp->Coeffs()(i) == triMap->Coeffs()(i)); // Values of coefficients should be equal
                 CHECK(&comp->Coeffs()(i) == &triMap->Coeffs()(i)); // Memory location should also be the same (no copy)
         }
-        
+
     }
 
     unsigned int numSamps = 10;
@@ -950,7 +970,7 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
 
     triMap->SetCoeffs(coeffs);
     auto out = triMap->Evaluate(in);
-    
+
     SECTION("Evaluation"){
 
         unsigned int start = 0;
@@ -991,9 +1011,9 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
     }
 
     SECTION("LogDeterminant"){
-        
+
         auto inTopAndMid = Kokkos::subview(in, std::make_pair(0, int(activeInd)), Kokkos::ALL());
-        
+
         auto logDet = triMap->LogDeterminant(in);
         auto logDet_ = comp->LogDeterminant(inTopAndMid);
 
@@ -1030,16 +1050,16 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( coeffGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
             coeffs(i) -= fdstep;
         }
-        
+
     }
 
 
@@ -1069,18 +1089,18 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
             evals2 = triMap->Evaluate(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd){
-                
+
                 double fdDeriv = 0.0;
                 for(unsigned int j=0; j<triMap->outputDim; ++j)
                     fdDeriv += sens(j,ptInd) * (evals2(j,ptInd)-evals(j,ptInd))/fdstep;
 
-                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3)); 
+                CHECK( inputGrad(i,ptInd) == Approx(fdDeriv).margin(1e-3));
             }
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
                 in(i,ptInd) -= fdstep;
         }
-        
+
     }
 
     SECTION("LogDeterminantCoeffGrad"){
@@ -1095,7 +1115,7 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
         Kokkos::View<double**,Kokkos::HostSpace> detGrad = triMap->LogDeterminantCoeffGrad(in);
         REQUIRE(detGrad.extent(0)==triMap->numCoeffs);
         REQUIRE(detGrad.extent(1)==numSamps);
-        
+
         Kokkos::View<double*,Kokkos::HostSpace> logDet = triMap->LogDeterminant(in);
         Kokkos::View<double*,Kokkos::HostSpace> logDet2;
 
@@ -1108,8 +1128,8 @@ TEST_CASE( "Testing TriangularMap made using CreateSingleEntryMap", "[Triangular
             logDet2 = triMap->LogDeterminant(in);
 
             for(unsigned int ptInd=0; ptInd<numSamps; ++ptInd)
-                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4)); 
-            
+                CHECK( detGrad(i,ptInd) == Approx((logDet2(ptInd)-logDet(ptInd))/fdstep).margin(1e-4));
+
             coeffs(i) -= fdstep;
         }