diff --git a/source/mir/ndslice/algorithm.d b/source/mir/ndslice/algorithm.d
index 9f2236a4..21239b9d 100644
--- a/source/mir/ndslice/algorithm.d
+++ b/source/mir/ndslice/algorithm.d
@@ -52,17 +52,25 @@ private void checkShapesMatch(
     }
 }
 
-
-private auto ref frontOf(alias slice)() { return slice.front; };
+template frontOf(size_t N)
+{
+    static if (N == 0)
+        enum frontOf = "";
+    else
+    {
+        enum i = N - 1;
+        enum frontOf = frontOf!i ~ "slices[" ~ i.stringof ~ "].front, ";
+    }
+}
 
 S reduceImpl(alias fun, S, Slices...)(S seed, Slices slices)
 {
     do
     {
         static if (slices[0].shape.length == 1)
-            seed = fun(seed, staticMap!(frontOf, slices));
+            seed = mixin("fun(seed, " ~ frontOf!(Slices.length) ~ ")");
         else
-            seed = .reduceImpl!fun(seed, staticMap!(frontOf, slices));
+            seed = mixin(".reduceImpl!fun(seed," ~ frontOf!(Slices.length) ~ ")");
         foreach(ref slice; slices)
             slice.popFront;
     }
@@ -267,9 +275,9 @@ void eachImpl(alias fun, Slices...)(Slices slices)
     do
     {
         static if (slices[0].shape.length == 1)
-            fun(staticMap!(frontOf, slices));
+            mixin("fun(" ~ frontOf!(Slices.length) ~ ");");
         else
-            .eachImpl!fun(staticMap!(frontOf, slices));
+            mixin(".eachImpl!fun(" ~ frontOf!(Slices.length) ~ ");");
         foreach(ref slice; slices)
             slice.popFront;
     }
@@ -385,7 +393,7 @@ size_t findImpl(alias fun, size_t N, Slices...)(ref size_t[N] backwardIndex, Sli
     {
         static if (slices[0].shape.length == 1)
         {
-            if (fun(staticMap!(frontOf, slices)))
+            if (mixin("fun(" ~ frontOf!(Slices.length) ~ ")"))
             {
                 backwardIndex[0] = slices[0].length;
                 return 1;
@@ -393,7 +401,7 @@ size_t findImpl(alias fun, size_t N, Slices...)(ref size_t[N] backwardIndex, Sli
         }
         else
         {
-            if (findImpl!fun(backwardIndex[1 .. $], staticMap!(frontOf, slices)))
+            if (mixin("findImpl!fun(backwardIndex[1 .. $], " ~ frontOf!(Slices.length) ~ ")"))
             {
                 backwardIndex[0] = slices[0].length;
                 return 1;
@@ -556,12 +564,12 @@ size_t anyImpl(alias fun, Slices...)(Slices slices)
     {
         static if (slices[0].shape.length == 1)
         {
-            if (fun(staticMap!(frontOf, slices)))
+            if (mixin("fun(" ~ frontOf!(Slices.length) ~ ")"))
                 return true;
         }
         else
         {
-            if (anyImpl!fun(staticMap!(frontOf, slices)))
+            if (mixin("anyImpl!fun(" ~ frontOf!(Slices.length) ~ ")"))
                 return true;
         }
         foreach(ref slice; slices)
@@ -676,12 +684,12 @@ size_t allImpl(alias fun, Slices...)(Slices slices)
     {
         static if (slices[0].shape.length == 1)
         {
-            if (!fun(staticMap!(frontOf, slices)))
+            if (!mixin("fun(" ~ frontOf!(Slices.length) ~ ")"))
                 return false;
         }
         else
         {
-            if (!allImpl!fun(staticMap!(frontOf, slices)))
+            if (!mixin("allImpl!fun(" ~ frontOf!(Slices.length) ~ ")"))
                 return false;
         }
         foreach(ref slice; slices)
@@ -1096,11 +1104,11 @@ size_t countImpl(alias fun, Slices...)(Slices slices)
     {
         static if (slices[0].shape.length == 1)
         {
-            if(fun(staticMap!(frontOf, slices)))
+            if(mixin("fun(" ~ frontOf!(Slices.length) ~ ")"))
                 ret++;
         }
         else
-            ret += .countImpl!fun(staticMap!(frontOf, slices));
+            ret += mixin(".countImpl!fun(" ~ frontOf!(Slices.length) ~ ")");
         foreach(ref slice; slices)
             slice.popFront;
     }
diff --git a/source/mir/ndslice/allocation.d b/source/mir/ndslice/allocation.d
index 278248bb..a5cedc8d 100644
--- a/source/mir/ndslice/allocation.d
+++ b/source/mir/ndslice/allocation.d
@@ -28,6 +28,7 @@ module mir.ndslice.allocation;
 import std.traits;
 import mir.ndslice.slice;
 import mir.ndslice.internal;
+import mir.ndslice.stack;
 
 @fastmath:
 
@@ -94,6 +95,14 @@ pure nothrow unittest
     assert(tensor[1, 1] == 5);
 }
 
+/// ditto
+auto slice(size_t dim, Slices...)(Stack!(dim, Slices) stack)
+{
+    auto ret = .slice!(Unqual!(stack.DeepElemType))(stack.shape);
+    ret[] = stack;
+    return ret;
+}
+
 pure nothrow unittest
 {
     import mir.ndslice.topology : iota;
diff --git a/source/mir/ndslice/package.d b/source/mir/ndslice/package.d
index fb2a0a30..0b03de7f 100644
--- a/source/mir/ndslice/package.d
+++ b/source/mir/ndslice/package.d
@@ -157,6 +157,17 @@ $(TR $(TDNW $(SUBMODULE algorithm)
     )
 )
 
+$(TR $(TDNW $(SUBMODULE stack)
+        $(BR) $(SMALL Concatenation and algorithms))
+     $(TD
+        $(SUBREF stack, isStack)
+        $(SUBREF stack, stack)
+        $(SUBREF stack, Stack)
+        $(SUBREF stack, stackDimension)
+        $(SUBREF stack, until)
+    )
+)
+
 $(TR $(TDNW $(SUBMODULE dynamic)
         $(BR) $(SMALL Dynamic dimension manipulators))
      $(TD
@@ -378,6 +389,7 @@ public import mir.ndslice.algorithm;
 public import mir.ndslice.allocation;
 public import mir.ndslice.dynamic;
 public import mir.ndslice.slice;
+public import mir.ndslice.stack;
 public import mir.ndslice.topology;
 
 
diff --git a/source/mir/ndslice/slice.d b/source/mir/ndslice/slice.d
index 4d6df111..22a81b30 100644
--- a/source/mir/ndslice/slice.d
+++ b/source/mir/ndslice/slice.d
@@ -33,6 +33,7 @@ import std.meta;
 
 import mir.internal.utility;
 import mir.ndslice.internal;
+import mir.ndslice.stack;
 import mir.primitives;
 
 @fastmath:
@@ -307,6 +308,8 @@ auto slicedField(Field)(Field field)
 Returns the element type of a $(LREF Slice).
 +/
 alias DeepElementType(S : Slice!(kind, packs, Iterator), SliceKind kind, size_t[] packs, Iterator) = S.DeepElemType;
+/// ditto
+alias DeepElementType(S : Stack!(dim, Slices), size_t dim, Slices...) = S.DeepElemType;
 
 ///
 unittest
@@ -1775,6 +1778,41 @@ struct Slice(SliceKind kind, size_t[] packs, Iterator)
                      [2, 2, 3, 3]]);
         }
 
+
+        private void opIndexOpAssignImplStack(string op, T)(T value)
+        {
+            auto sl = this;
+            static if (stackDimension!T)
+            {
+                if (!sl.empty) do
+                {
+                    mixin(`sl.front[] ` ~ op ~ `= value.front;`);
+                    value.popFront;
+                    sl.popFront;
+                }
+                while(!sl.empty);
+            }
+            else
+            {
+                foreach (ref slice; value._slices)
+                {
+                    mixin("sl[0 .. slice.length][] " ~ op ~ "= slice;");
+                    sl = sl[slice.length .. $];
+                }
+                assert(sl.empty);
+            }
+        }
+
+        ///
+        void opIndexAssign(T, Slices...)(T stack, Slices slices)
+            if (isFullPureSlice!Slices && isStack!T)
+        {
+            import mir.ndslice.topology : unpack;
+            auto sl = this[slices].unpack;
+            static assert(isSlice!(typeof(sl))[0] == stack.N);
+            sl.opIndexOpAssignImplStack!""(stack);
+        }
+
         /++
         Assignment of a value (e.g. a number) to a $(B fully defined slice).
 
@@ -1784,7 +1822,8 @@ struct Slice(SliceKind kind, size_t[] packs, Iterator)
         void opIndexAssign(T, Slices...)(T value, Slices slices)
             if (isFullPureSlice!Slices
                 && (!isDynamicArray!T || isDynamicArray!DeepElemType)
-                && !isSlice!T)
+                && !isSlice!T
+                && !isStack!T)
         {
             import mir.ndslice.topology : unpack;
             auto sl = this[slices].unpack;
@@ -2042,7 +2081,8 @@ struct Slice(SliceKind kind, size_t[] packs, Iterator)
         void opIndexOpAssign(string op, T, Slices...)(T value, Slices slices)
             if (isFullPureSlice!Slices
                 && (!isDynamicArray!T || isDynamicArray!DeepElemType)
-                && !isSlice!T)
+                && !isSlice!T
+                && !isStack!T)
         {
             import mir.ndslice.topology : unpack;
             auto sl = this[slices].unpack;
@@ -2067,6 +2107,16 @@ struct Slice(SliceKind kind, size_t[] packs, Iterator)
             assert(a[1] == [6, 6, 1]);
         }
 
+        ///
+        void opIndexOpAssign(string op,T, Slices...)(T stack, Slices slices)
+            if (isFullPureSlice!Slices && isStack!T)
+        {
+            import mir.ndslice.topology : unpack;
+            auto sl = this[slices].unpack;
+            static assert(isSlice!(typeof(sl))[0] == stack.N);
+            sl.opIndexOpAssignImplStack!op(stack);
+        }
+
         static if (doUnittest)
         /// Packed slices have the same behavior.
         pure nothrow unittest
diff --git a/source/mir/ndslice/stack.d b/source/mir/ndslice/stack.d
index 40e6058f..9300c433 100644
--- a/source/mir/ndslice/stack.d
+++ b/source/mir/ndslice/stack.d
@@ -2,7 +2,7 @@
 This is a submodule of $(MREF mir, ndslice).
 
 License:   $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).
-Copyright: Copyright © 2016-, Ilya Yaroshenko
+Copyright: Copyright © 2017-, Ilya Yaroshenko
 Authors:   Ilya Yaroshenko
 
 Macros:
@@ -12,49 +12,238 @@ T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
 module mir.ndslice.stack;
 
 import std.traits;
+import std.meta;
 
 import mir.ndslice.internal;
+import mir.ndslice.slice;
+import mir.internal.utility;
 
 @fastmath:
 
-struct Stack(Slices...)
-    if (Slices.lengths > 1)
+/++
+Creates a $(LREF Stack) view of multiple slices.
+
+Can be used in combination with itself, $(LREF until), $(SUBREF, allocation, slice),
+and $(SUBREF slice, Slice) assignment.
+until pred returns true.
+
+Returns: true if an element was 
+
+Params:
+    slices = tuple of slices and stacks. All slices and stacks must have the same dimension count.
+
+Returns: $(LREF Stack).
++/
+auto stack(size_t dim = 0, Slices...)(Slices slices)
 {
-    private Slices slices;
+    return Stack!(dim, Slices)(slices);
+}
 
-    private bool empty(size_t dimension)() @property
-        if (dimension)
-    {
+/// 1D
+unittest
+{
+    import mir.ndslice.allocation: slice;
+    import mir.ndslice.topology: iota;
 
-    }
+    size_t i;
+    auto a = 3.iota;
+    auto b = iota([6], a.length);
+    auto s = stack(a, b);
+    assert(s.length == a.length + b.length);
+    // iteration with until
+    s.until!((elem){ assert(elem == i++); return false; });
+    // allocation with slice
+    assert(s.slice == s.length.iota);
+    // assignment
+    auto d = slice!double(s.length);
+    d[] = s;
+    assert(d == s.length.iota);
+    d[] += s;
+    assert(d == iota([s.length], 0, 2));
+}
 
-    ///
-    void each(alias fun, size_t dimension = 0)()
+/// Multidimensional
+unittest
+{
+    import mir.ndslice.allocation: slice;
+    import mir.ndslice.topology: iota;
+
+    // 0, 1, 2
+    // 3, 4, 5
+    auto a = iota(2, 3);
+    // 0, 1
+    // 2, 3
+    auto b = iota(2, 2);
+    // 0, 1, 2, 3, 4
+    auto c = iota(1, 5);
+
+    // 0, 1, 2, | 0, 1
+    // 3, 4, 5, | 2, 3
+    // ---------------
+    // 0, 1, 2,   3, 4
+    auto s = stack(stack!1(a, b), c);
+    // allocation
+    auto d = s.slice;
+    assert(d == [
+        [0, 1, 2, 0, 1],
+        [3, 4, 5, 2, 3],
+        [0, 1, 2, 3, 4],
+        ]);
+}
+
+template frontOf(size_t N)
+{
+    static if (N == 0)
+        enum frontOf = "";
+    else
     {
-        static if (dimension)
-            foreach(i, ref slice; slices)
-                for (auto t = slice; !t.empty!dimension; t.popFront!dimension)
-                    t.front!dimension.each!fun;
-        else
-            foreach (i, ref slice; slices)
-                slice.each!fun;
+        enum i = N - 1;
+        enum frontOf = frontOf!i ~ "slices[" ~ i.stringof ~ "].front!d, ";
     }
 }
 
-struct Transposition(Slice, size_t dimension)
-    if (dimension)
+///
+enum bool isStack(T) = is(T : Stack!(dim, Slices), size_t dim, Slices...);
+///
+enum size_t stackDimension(T : Stack!(dim, Slices), size_t dim, Slices...) = dim; 
+
+///
+struct Stack(size_t dim, Slices...)
+    if (Slices.length > 1)
 {
-    private Slice slice;
+    /// Slices and sub-stacks
+    Slices _slices;
+
+    static if (isSlice!(Slices[0]))
+    /// Dimension count
+        enum N = isSlice!(Slices[0])[0];
+    else
+        enum N = Slices[0].N;
+
+    static assert(dim < N);
 
-    ///
-    void each(alias fun)()
+    package alias DeepElemType = CommonType!(staticMap!(DeepElementType, Slices));
+
+    /// Length primitive
+    size_t length(size_t d = 0)() const @property
     {
-        static if (is(Slice : Stack!(Slices), Slices...))
-            alias slices = slice.slices;
+        static if (d == dim)
+        {
+            size_t length;
+            foreach(i; Iota!(Slices.length))
+                length += _slices[i].length!d;
+            return length;
+        }
         else
-            alias slices = AliasSeq!slice;
-        foreach(i, ref slice; slices)
-            for (auto t = slice; !slice.empty!dimension; slice.popFront!dimension )
-            {}
+        {
+            return _slices[0].length!d;
+        }
+    }
+
+    /// Total elements count in the stack.
+    size_t elementsCount()() const @property
+    {
+        size_t count = 1;
+        foreach(i; Iota!N)
+            count *= length!i;
+        return count;
+    }
+
+    /// Shape of the stack.
+    size_t[N] shape()() const @property
+    {
+        typeof(return) ret = void;
+        foreach(i; Iota!N)
+            ret[i] = length!i;
+        return ret;
+    }
+
+    /// Multidimensional input range primitives
+    bool empty(size_t d = 0)() const @property
+        if (d != dim)
+    {
+        return _slices[0].empty!d;
+    }
+
+    /// ditto
+    void popFront(size_t d = 0)()
+        if (d != dim)
+    {
+        foreach_reverse (i, ref slice; _slices)
+            _slices[i].popFront!d;
+    }
+
+    /// ditto
+    auto front(size_t d = 0)()
+        if (d != dim)
+    {
+        enum elemDim = d < dim ? dim - 1 : dim;
+        alias slices = _slices;
+        return mixin(`stack!elemDim(` ~ frontOf!(Slices.length) ~ `)`);
+    }
+}
+
+/++
+Iterates elements in $(SUBREF slice, Slice) or $(LREF Stack)
+until pred returns true.
+
+Returns: false if pred returned false for all elements and true otherwise.
++/
+template until(alias pred)
+{
+    import mir.functional: naryFun;
+    static if (__traits(isSame, naryFun!pred, pred))
+    {
+        /++
+        Specialization for slices
+        Params:
+            sl = $(SUBREF slice, Slice)
+        +/
+        bool until(SliceKind kind, size_t[] packs, Iterator)(Slice!(kind, packs, Iterator) sl)
+        {
+            static if (packs[0] == 1)
+            {
+                pragma(inline, false);
+                alias f = pred;
+            }
+            else
+                alias f = .until!pred;
+            if (!sl.empty) do
+            {
+                if (f(sl.front))
+                    return true;
+                sl.popFront;
+            }
+            while(!sl.empty);
+            return false;
+        }
+
+        /++
+        Specialization for stacks
+        Params:
+            st = $(LREF Stack)
+        +/
+        bool until(size_t dim, Slices...)(Stack!(dim, Slices) st)
+        {
+            static if (dim == 0)
+            {
+               foreach (i, ref slice; st._slices)
+                    if (.until!pred(slice))
+                        return true;
+            }
+            else
+            {
+                if (!st.empty) do
+                {
+                    if (.until!pred(st.front))
+                        return true;
+                    st.popFront;
+                }
+                while(!st.empty);
+            }
+            return false;
+        }
     }
+    else
+        alias until = .until!(naryFun!pred);
 }