mir.algorithm.iteration source code

1 // Written in the D programming language.
2 /**
3 This module contains generic _iteration algorithms.
4 $(SCRIPT inhibitQuickIndex = 1;)
5 
6 $(BOOKTABLE $(H2 Function),
7 $(TR $(TH Function Name) $(TH Description))
8 $(T2 all, Checks if all elements satisfy to a predicate.)
9 $(T2 any, Checks if at least one element satisfy to a predicate.)
10 $(T2 cmp, Compares two slices.)
11 $(T2 count, Counts elements in a slices according to a predicate.)
12 $(T2 each, Iterates elements.)
13 $(T2 eachLower, Iterates lower triangle of matrix.)
14 $(T2 eachOnBorder, Iterates elementes on tensors borders and corners.)
15 $(T2 eachUploPair, Iterates upper and lower pairs of elements in square matrix.)
16 $(T2 eachUpper, Iterates upper triangle of matrix.)
17 $(T2 equal, Compares two slices for equality.)
18 $(T2 filter, Filters elements in a range or an ndslice.)
19 $(T2 find, Finds backward index.)
20 $(T2 findIndex, Finds index.)
21 $(T2 fold, Accumulates all elements (different parameter order than `reduce`).)
22 $(T2 isSymmetric, Checks if the matrix is symmetric.)
23 $(T2 maxIndex, Finds index of the maximum.)
24 $(T2 maxPos, Finds backward index of the maximum.)
25 $(T2 minIndex, Finds index of the minimum.)
26 $(T2 minmaxIndex, Finds indices of the minimum and the maximum.)
27 $(T2 minmaxPos, Finds backward indices of the minimum and the maximum.)
28 $(T2 minPos, Finds backward index of the minimum.)
29 $(T2 nBitsToCount, Сount bits until set bit count is reached.)
30 $(T2 reduce, Accumulates all elements.)
31 $(T2 Chequer, Chequer color selector to work with $(LREF each) .)
32 $(T2 uniq, Iterates over the unique elements in a range or an ndslice, which is assumed sorted.)
33 )
34 
35 Transform function is represented by $(NDSLICEREF topology, map).
36 
37 All operators are suitable to change slices using `ref` argument qualification in a function declaration.
38 Note, that string lambdas in Mir are `auto ref` functions.
39 
40 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
41 Copyright: 2020 Ilia Ki, Kaleidic Associates Advisory Limited, Symmetry Investments
42 Authors: Ilia Ki, John Michael Hall, Andrei Alexandrescu (original Phobos code)
43 
44 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
45 Copyright: 2020 Ilia Ki, Kaleidic Associates Advisory Limited, Symmetry Investments
46 
47 Authors: , Ilia Ki (Mir & BetterC rework).
48 Source: $(PHOBOSSRC std/algorithm/_iteration.d)
49 Macros:
50     NDSLICEREF = $(REF_ALTTEXT $(TT $2), $2, mir, ndslice, $1)$(NBSP)
51     T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
52  */
53 module mir.algorithm.iteration;
54 
55 import mir.functional: naryFun;
56 import mir.internal.utility;
57 import mir.math.common: fmamath;
58 import mir.ndslice.field: BitField;
59 import mir.ndslice.internal;
60 import mir.ndslice.iterator: FieldIterator, RetroIterator;
61 import mir.ndslice.slice;
62 import mir.primitives;
63 import mir.qualifier;
64 import std.meta;
65 import std.traits;
66 
67 /++
68 Chequer color selector to work with $(LREF each)
69 +/
70 enum Chequer : bool
71 {
72     /// Main diagonal color
73     black,
74     /// First sub-diagonal color
75     red,
76 }
77 
78 ///
79 version(mir_test) unittest
80 {
81     import mir.ndslice.allocation: slice;
82     auto s = [5, 4].slice!int;
83 
84     Chequer.black.each!"a = 1"(s);
85     assert(s == [
86         [1, 0, 1, 0],
87         [0, 1, 0, 1],
88         [1, 0, 1, 0],
89         [0, 1, 0, 1],
90         [1, 0, 1, 0],
91     ]);
92 
93     Chequer.red.each!((ref b) => b = 2)(s);
94     assert(s == [
95         [1, 2, 1, 2],
96         [2, 1, 2, 1],
97         [1, 2, 1, 2],
98         [2, 1, 2, 1],
99         [1, 2, 1, 2],
100     ]);
101 
102 }
103 
104 @fmamath:
105 
106 /+
107 Bitslice representation for accelerated bitwise algorithm.
108 1-dimensional contiguousitslice can be split into three chunks: head bits, body chunks, and tail bits.
109 
110 Bitslice can have head bits because it has slicing and the zero bit may not be aligned to the zero of a body chunk.
111 +/
112 private struct BitSliceAccelerator(Field, I = typeof(Field.init[size_t.init]))
113     if (__traits(isUnsigned, I))
114 {
115     import mir.bitop;
116     import mir.qualifier: lightConst;
117     import mir.ndslice.traits: isIterator;
118     import mir.ndslice.iterator: FieldIterator;
119     import mir.ndslice.field: BitField;
120 
121     ///
122     alias U = typeof(I + 1u);
123     /// body bits chunks
124     static if (isIterator!Field)
125         Slice!Field bodyChunks;
126     else
127         Slice!(FieldIterator!Field) bodyChunks;
128     /// head length
129     int headLength;
130     /// tail length
131     int tailLength;
132 
133 @fmamath:
134 
135     this(Slice!(FieldIterator!(BitField!(Field, I))) slice)
136     {
137         enum mask = bitShiftMask!I;
138         enum shift = bitElemShift!I;
139         size_t length = slice.length;
140         size_t index = slice._iterator._index;
141         if (auto hlen = index & mask)
142         {
143             auto l = I.sizeof * 8 - hlen;
144             if (l > length)
145             {
146                 // central problem
147                 headLength = -cast(int) length;
148                 tailLength = cast(int) hlen;
149                 goto F;
150             }
151             else
152             {
153                 headLength = cast(uint) l;
154                 length -= l;
155                 index += l;
156             }
157         }
158         tailLength = cast(int) (length & mask);
159     F:
160         length >>= shift;
161         index >>= shift;
162         bodyChunks._lengths[0] = length;
163         static if (isIterator!Field)
164         {
165             bodyChunks._iterator = slice._iterator._field._field;
166             bodyChunks._iterator += index;
167         }
168         else
169         {
170             bodyChunks._iterator._index = index;
171             bodyChunks._iterator._field = slice._iterator._field._field;
172         }
173     }
174 
175 scope const:
176 
177     bool isCentralProblem()
178     {
179         return headLength < 0;
180     }
181 
182     U centralBits()
183     {
184         assert(isCentralProblem);
185         return *bodyChunks._iterator.lightConst >>> tailLength;
186     }
187 
188     uint centralLength()
189     {
190         assert(isCentralProblem);
191         return -headLength;
192     }
193 
194     /// head bits (last `headLength` bits are valid).
195     U headBits()
196     {
197         assert(!isCentralProblem);
198         if (headLength == 0)
199             return U.init;
200         static if (isIterator!Field)
201             return bodyChunks._iterator.lightConst[-1];
202         else
203             return bodyChunks._iterator._field.lightConst[bodyChunks._iterator._index - 1];
204     }
205 
206     /// tail bits (first `tailLength` bits are valid).
207     U tailBits()
208     {
209         assert(!isCentralProblem);
210         if (tailLength == 0)
211             return U.init;
212         static if (isIterator!Field)
213             return bodyChunks._iterator.lightConst[bodyChunks.length];
214         else
215             return bodyChunks._iterator._field.lightConst[bodyChunks._iterator._index + bodyChunks.length];
216     }
217 
218     U negCentralMask()
219     {
220         return U.max << centralLength;
221     }
222 
223     U negHeadMask()
224     {
225         return U.max << headLength;
226     }
227 
228     U negTailMask()
229     {
230         return U.max << tailLength;
231     }
232 
233     U negCentralMaskS()
234     {
235         return U.max >> centralLength;
236     }
237 
238     U negHeadMaskS()
239     {
240         return U.max >> headLength;
241     }
242 
243     U negTailMaskS()
244     {
245         return U.max >> tailLength;
246     }
247 
248     U centralBitsWithRemainingZeros()
249     {
250         return centralBits & ~negCentralMask;
251     }
252 
253     U centralBitsWithRemainingZerosS()
254     {
255         return centralBits << (U.sizeof * 8 - centralLength);
256     }
257 
258     U headBitsWithRemainingZeros()
259     {
260         return headBits >>> (I.sizeof * 8 - headLength);
261     }
262 
263     U headBitsWithRemainingZerosS()
264     {
265         static if (U.sizeof > I.sizeof)
266             return (headBits << (U.sizeof - I.sizeof) * 8) & ~negTailMaskS;
267         else
268             return headBits & ~negTailMaskS;
269     }
270 
271     U tailBitsWithRemainingZeros()
272     {
273         return tailBits & ~negTailMask;
274     }
275 
276     U tailBitsWithRemainingZerosS()
277     {
278         return tailBits << (U.sizeof * 8 - tailLength);
279     }
280 
281     U centralBitsWithRemainingOnes()
282     {
283         return centralBits | negCentralMask;
284     }
285 
286     U centralBitsWithRemainingOnesS()
287     {
288         return centralBitsWithRemainingZerosS | negCentralMaskS;
289     }
290 
291     U headBitsWithRemainingOnes()
292     {
293         return headBitsWithRemainingZeros | negHeadMask;
294     }
295 
296     U headBitsWithRemainingOnesS()
297     {
298         return headBitsWithRemainingZerosS | negHeadMaskS;
299     }
300 
301     U tailBitsWithRemainingOnes()
302     {
303         return tailBits | negTailMask;
304     }
305 
306     U tailBitsWithRemainingOnesS()
307     {
308         return tailBitsWithRemainingZerosS | negTailMaskS;
309     }
310 
311     size_t ctpop()
312     {
313         import mir.bitop: ctpop;
314         if (isCentralProblem)
315             return centralBitsWithRemainingZeros.ctpop;
316         size_t ret;
317         if (headLength)
318             ret = cast(size_t) headBitsWithRemainingZeros.ctpop;
319         if (bodyChunks.length)
320         {
321             auto bc = bodyChunks.lightConst;
322             do
323             {
324                 ret += cast(size_t) bc.front.ctpop;
325                 bc.popFront;
326             }
327             while(bc.length);
328         }
329         if (tailBits)
330             ret += cast(size_t) tailBitsWithRemainingZeros.ctpop;
331         return ret;
332     }
333 
334     bool any()
335     {
336         if (isCentralProblem)
337             return centralBitsWithRemainingZeros != 0;
338         if (headBitsWithRemainingZeros != 0)
339             return true;
340         if (bodyChunks.length)
341         {
342             auto bc = bodyChunks.lightConst;
343             do
344             {
345                 if (bc.front != 0)
346                     return true;
347                 bc.popFront;
348             }
349             while(bc.length);
350         }
351         if (tailBitsWithRemainingZeros != 0)
352             return true;
353         return false;
354     }
355 
356     bool all()
357     {
358         if (isCentralProblem)
359             return centralBitsWithRemainingOnes != U.max;
360         size_t ret;
361         if (headBitsWithRemainingOnes != U.max)
362             return false;
363         if (bodyChunks.length)
364         {
365             auto bc = bodyChunks.lightConst;
366             do
367             {
368                 if (bc.front != I.max)
369                     return false;
370                 bc.popFront;
371             }
372             while(bc.length);
373         }
374         if (tailBitsWithRemainingOnes != U.max)
375             return false;
376         return true;
377     }
378 
379     size_t cttz()
380     {
381         U v;
382         size_t ret;
383         if (isCentralProblem)
384         {
385             v = centralBitsWithRemainingOnes;
386             if (v)
387                 goto R;
388             ret = centralLength;
389             goto L;
390         }
391         v = headBitsWithRemainingOnes;
392         if (v)
393             goto R;
394         ret = headLength;
395         if (bodyChunks.length)
396         {
397             auto bc = bodyChunks.lightConst;
398             do
399             {
400                 v = bc.front;
401                 if (v)
402                     goto R;
403                 ret += I.sizeof * 8;
404                 bc.popFront;
405             }
406             while(bc.length);
407         }
408         v = tailBitsWithRemainingOnes;
409         if (v)
410             goto R;
411         ret += tailLength;
412         goto L;
413     R:
414         ret += v.cttz;
415     L:
416         return ret;
417     }
418 
419     size_t ctlz()
420     {
421         U v;
422         size_t ret;
423         if (isCentralProblem)
424         {
425             v = centralBitsWithRemainingOnes;
426             if (v)
427                 goto R;
428             ret = centralLength;
429             goto L;
430         }
431         v = tailBitsWithRemainingOnesS;
432         if (v)
433             goto R;
434         ret = tailLength;
435         if (bodyChunks.length)
436         {
437             auto bc = bodyChunks.lightConst;
438             do
439             {
440                 v = bc.back;
441                 if (v)
442                     goto R;
443                 ret += I.sizeof * 8;
444                 bc.popBack;
445             }
446             while(bc.length);
447         }
448         v = headBitsWithRemainingOnesS;
449         if (v)
450             goto R;
451         ret += headLength;
452         goto L;
453     R:
454         ret += v.ctlz;
455     L:
456         return ret;
457     }
458 
459     sizediff_t nBitsToCount(size_t count)
460     {
461         pragma(inline, false);
462         size_t ret;
463         if (count == 0)
464             return count;
465         U v, cnt;
466         if (isCentralProblem)
467         {
468             v = centralBitsWithRemainingZeros;
469             goto E;
470         }
471         v = headBitsWithRemainingZeros;
472         cnt = v.ctpop;
473         if (cnt >= count)
474             goto R;
475         ret += headLength;
476         count -= cast(size_t) cnt;
477         if (bodyChunks.length)
478         {
479             auto bc = bodyChunks.lightConst;
480             do
481             {
482                 v = bc.front;
483                 cnt = v.ctpop;
484                 if (cnt >= count)
485                     goto R;
486                 ret += I.sizeof * 8;
487                 count -= cast(size_t) cnt;
488                 bc.popFront;
489             }
490             while(bc.length);
491         }
492         v = tailBitsWithRemainingZeros;
493     E:
494         cnt = v.ctpop;
495         if (cnt >= count)
496             goto R;
497         return -1;
498     R:
499         return ret + v.nTrailingBitsToCount(count);
500     }
501 
502     sizediff_t retroNBitsToCount(size_t count)
503     {
504         if (count == 0)
505             return count;
506         size_t ret;
507         U v, cnt;
508         if (isCentralProblem)
509         {
510             v = centralBitsWithRemainingZerosS;
511             goto E;
512         }
513         v = tailBitsWithRemainingZerosS;
514         cnt = v.ctpop;
515         if (cnt >= count)
516             goto R;
517         ret += tailLength;
518         count -= cast(size_t) cnt;
519         if (bodyChunks.length)
520         {
521             auto bc = bodyChunks.lightConst;
522             do
523             {
524                 v = bc.back;
525                 cnt = v.ctpop;
526                 if (cnt >= count)
527                     goto R;
528                 ret += I.sizeof * 8;
529                 count -= cast(size_t) cnt;
530                 bc.popBack;
531             }
532             while(bc.length);
533         }
534         v = headBitsWithRemainingZerosS;
535     E:
536         cnt = v.ctpop;
537         if (cnt >= count)
538             goto R;
539         return -1;
540     R:
541         return ret + v.nLeadingBitsToCount(count);
542     }
543 }
544 
545 /++
546 Сount bits until set bit count is reached. Works with ndslices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
547 Returns: bit count if set bit count is reached or `-1` otherwise.
548 +/
549 sizediff_t nBitsToCount(Field, I)(Slice!(FieldIterator!(BitField!(Field, I))) bitSlice, size_t count)
550 {
551     return BitSliceAccelerator!(Field, I)(bitSlice).nBitsToCount(count);
552 }
553 
554 ///ditto
555 sizediff_t nBitsToCount(Field, I)(Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))) bitSlice, size_t count)
556 {
557     import mir.ndslice.topology: retro;
558     return BitSliceAccelerator!(Field, I)(bitSlice.retro).retroNBitsToCount(count);
559 }
560 
561 /++
562 Сount bits starting from the end until set bit count is reached. Works with ndslices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
563 Returns: bit count if set bit count is reached or `-1` otherwise.
564 +/
565 sizediff_t retroNBitsToCount(Field, I)(Slice!(FieldIterator!(BitField!(Field, I))) bitSlice, size_t count)
566 {
567     return BitSliceAccelerator!(Field, I)(bitSlice).retroNBitsToCount(count);
568 }
569 
570 ///ditto
571 sizediff_t retroNBitsToCount(Field, I)(Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))) bitSlice, size_t count)
572 {
573     import mir.ndslice.topology: retro;
574     return BitSliceAccelerator!(Field, I)(bitSlice.retro).nBitsToCount(count);
575 }
576 
577 ///
578 version(mir_test)
579 pure unittest
580 {
581     import mir.ndslice.allocation: bitSlice;
582     import mir.ndslice.topology: retro;
583     auto s = bitSlice(1000);
584     s[50] = true;
585     s[100] = true;
586     s[200] = true;
587     s[300] = true;
588     s[400] = true;
589     assert(s.nBitsToCount(4) == 301);
590     assert(s.retro.nBitsToCount(4) == 900);
591 }
592 
593 private void checkShapesMatch(
594     string fun = __FUNCTION__,
595     string pfun = __PRETTY_FUNCTION__,
596     Slices...)
597     (Slices slices)
598     if (Slices.length > 1)
599 {
600     enum msgShape = "all slices must have the same shape"  ~ tailErrorMessage!(fun, pfun);
601     enum N = slices[0].shape.length;
602     foreach (i, Slice; Slices)
603     {
604         static if (i == 0)
605             continue;
606         else
607         static if (slices[i].shape.length == N)
608             assert(slices[i].shape == slices[0].shape, msgShape);
609         else
610         {
611             import mir.ndslice.fuse: fuseShape;
612             static assert(slices[i].fuseShape.length >= N);
613             assert(cast(size_t[N])slices[i].fuseShape[0 .. N] == slices[0].shape, msgShape);
614         }
615     }
616 }
617 
618 
619 package(mir) template allFlattened(args...)
620 {
621     static if (args.length)
622     {
623         alias arg = args[0];
624         @fmamath @property allFlattenedMod()()
625         {
626             import mir.ndslice.topology: flattened;
627             return flattened(arg);
628         }
629         alias allFlattened = AliasSeq!(allFlattenedMod, allFlattened!(args[1..$]));
630     }
631     else
632         alias allFlattened = AliasSeq!();
633 }
634 
635 private template areAllContiguousSlices(Slices...)
636 {
637     import mir.ndslice.traits: isContiguousSlice;
638      static if (allSatisfy!(isContiguousSlice, Slices))
639         enum areAllContiguousSlices = Slices[0].N > 1 && areAllContiguousSlicesImpl!(Slices[0].N, Slices[1 .. $]);
640      else
641         enum areAllContiguousSlices = false;
642 }
643 
644 private template areAllContiguousSlicesImpl(size_t N, Slices...)
645 {
646      static if (Slices.length == 0)
647         enum areAllContiguousSlicesImpl = true;
648      else
649         enum areAllContiguousSlicesImpl = Slices[0].N == N && areAllContiguousSlicesImpl!(N, Slices[1 .. $]);
650 }
651 
652 version(LDC) {}
653 else version(GNU) {}
654 else version (Windows) {}
655 else version (X86_64)
656 {
657     //Compiling with DMD for x86-64 for Linux & OS X with optimizations enabled,
658     //"Tensor mutation on-the-fly" unittest was failing. Disabling inlining
659     //caused it to succeed.
660     //TODO: Rework so this is unnecessary!
661     version = Mir_disable_inlining_in_reduce;
662 }
663 
664 version(Mir_disable_inlining_in_reduce)
665 {
666     private enum Mir_disable_inlining_in_reduce = true;
667 
668     private template _naryAliases(size_t n)
669     {
670         static if (n == 0)
671             enum _naryAliases = "";
672         else
673         {
674             enum i = n - 1;
675             enum _naryAliases = _naryAliases!i ~ "alias " ~ cast(char)('a' + i) ~ " = args[" ~ i.stringof ~ "];\n";
676         }
677     }
678 
679     private template nonInlinedNaryFun(alias fun)
680     {
681         import mir.math.common : fmamath;
682         static if (is(typeof(fun) : string))
683         {
684             /// Specialization for string lambdas
685             @fmamath auto ref nonInlinedNaryFun(Args...)(auto ref Args args)
686                 if (args.length <= 26)
687             {
688                 pragma(inline,false);
689                 mixin(_naryAliases!(Args.length));
690                 return mixin(fun);
691             }
692         }
693         else static if (is(typeof(fun.opCall) == function))
694         {
695             @fmamath auto ref nonInlinedNaryFun(Args...)(auto ref Args args)
696                 if (is(typeof(fun.opCall(args))))
697             {
698                 pragma(inline,false);
699                 return fun.opCall(args);
700             }
701         }
702         else
703         {
704             @fmamath auto ref nonInlinedNaryFun(Args...)(auto ref Args args)
705                 if (is(typeof(fun(args))))
706             {
707                 pragma(inline,false);
708                 return fun(args);
709             }
710         }
711     }
712 }
713 else
714 {
715     private enum Mir_disable_inlining_in_reduce = false;
716 }
717 
718 S reduceImpl(alias fun, S, Slices...)(S seed, Slices slices)
719 {
720     do
721     {
722         static if (DimensionCount!(Slices[0]) == 1)
723             mixin(`seed = fun(seed,` ~ frontOf!(Slices.length) ~ `);`);
724         else
725             mixin(`seed = .reduceImpl!fun(seed,` ~ frontOf!(Slices.length) ~ `);`);
726         foreach_reverse(ref slice; slices)
727             slice.popFront;
728     }
729     while(!slices[0].empty);
730     return seed;
731 }
732 
733 /++
734 Implements the homonym function (also known as `accumulate`,
735 `compress`, `inject`, or `fold`) present in various programming
736 languages of functional flavor. The call `reduce!(fun)(seed, slice1, ..., sliceN)`
737 first assigns `seed` to an internal variable `result`,
738 also called the accumulator. Then, for each set of element `x1, ..., xN` in
739 `slice1, ..., sliceN`, `result = fun(result, x1, ..., xN)` gets evaluated. Finally,
740 `result` is returned.
741 
742 `reduce` allows to iterate multiple slices in the lockstep.
743 
744 Note:
745     $(NDSLICEREF topology, pack) can be used to specify dimensions.
746 Params:
747     fun = A function.
748 See_Also:
749     $(HTTP llvm.org/docs/LangRef.html#fast-math-flags, LLVM IR: Fast Math Flags)
750 
751     $(HTTP en.wikipedia.org/wiki/Fold_(higher-order_function), Fold (higher-order function))
752 +/
753 template reduce(alias fun)
754 {
755     import mir.functional: naryFun;
756     static if (__traits(isSame, naryFun!fun, fun)
757         && !Mir_disable_inlining_in_reduce)
758     /++
759     Params:
760         seed = An initial accumulation value.
761         slices = One or more slices, range, and arrays.
762     Returns:
763         the accumulated `result`
764     +/
765     @fmamath auto reduce(S, Slices...)(S seed, Slices slices)
766         if (Slices.length)
767     {
768         static if (Slices.length > 1)
769             slices.checkShapesMatch;
770         static if (areAllContiguousSlices!Slices)
771         {
772             import mir.ndslice.topology: flattened;
773             return .reduce!fun(seed, allFlattened!(allLightScope!slices));
774         }
775         else
776         {
777             if (slices[0].anyEmpty)
778                 return cast(Unqual!S) seed;
779             static if (is(S : Unqual!S))
780                 alias UT = Unqual!S;
781             else
782                 alias UT = S;
783             return reduceImpl!(fun, UT, Slices)(seed, allLightScope!slices);
784         }
785     }
786     else version(Mir_disable_inlining_in_reduce)
787     //As above, but with inlining disabled.
788     @fmamath auto reduce(S, Slices...)(S seed, Slices slices)
789         if (Slices.length)
790     {
791         static if (Slices.length > 1)
792             slices.checkShapesMatch;
793         static if (areAllContiguousSlices!Slices)
794         {
795             import mir.ndslice.topology: flattened;
796             return .reduce!fun(seed, allFlattened!(allLightScope!slices));
797         }
798         else
799         {
800             if (slices[0].anyEmpty)
801                 return cast(Unqual!S) seed;
802             static if (is(S : Unqual!S))
803                 alias UT = Unqual!S;
804             else
805                 alias UT = S;
806             return reduceImpl!(nonInlinedNaryFun!fun, UT, Slices)(seed, allLightScope!slices);
807         }
808     }
809     else
810         alias reduce = .reduce!(naryFun!fun);
811 }
812 
813 /// Ranges and arrays
814 version(mir_test)
815 unittest
816 {
817     auto ar = [1, 2, 3];
818     auto s = 0.reduce!"a + b"(ar);
819     assert (s == 6);
820 }
821 
822 /// Single slice
823 version(mir_test)
824 unittest
825 {
826     import mir.ndslice.topology : iota;
827 
828     //| 0 1 2 | => 3  |
829     //| 3 4 5 | => 12 | => 15
830     auto sl = iota(2, 3);
831 
832     // sum of all element in the slice
833     auto res = size_t(0).reduce!"a + b"(sl);
834 
835     assert(res == 15);
836 }
837 
838 /// Multiple slices, dot product
839 version(mir_test)
840 unittest
841 {
842     import mir.ndslice.allocation : slice;
843     import mir.ndslice.topology : as, iota;
844 
845     //| 0 1 2 |
846     //| 3 4 5 |
847     auto a = iota([2, 3], 0).as!double.slice;
848     //| 1 2 3 |
849     //| 4 5 6 |
850     auto b = iota([2, 3], 1).as!double.slice;
851 
852     alias dot = reduce!"a + b * c";
853     auto res = dot(0.0, a, b);
854 
855     // check the result:
856     import mir.ndslice.topology : flattened;
857     import std.numeric : dotProduct;
858     assert(res == dotProduct(a.flattened, b.flattened));
859 }
860 
861 /// Zipped slices, dot product
862 pure
863 version(mir_test) unittest
864 {
865     import std.typecons : Yes;
866     import std.numeric : dotProduct;
867     import mir.ndslice.allocation : slice;
868     import mir.ndslice.topology : as, iota, zip, universal;
869     import mir.math.common : fmamath;
870 
871     static @fmamath T fmuladd(T, Z)(const T a, Z z)
872     {
873         return a + z.a * z.b;
874     }
875 
876     // 0 1 2
877     // 3 4 5
878     auto sl1 = iota(2, 3).as!double.slice.universal;
879     // 1 2 3
880     // 4 5 6
881     auto sl2 = iota([2, 3], 1).as!double.slice;
882 
883     // slices must have the same strides for `zip!true`.
884     assert(sl1.strides == sl2.strides);
885 
886     auto z = zip!true(sl1, sl2);
887 
888     auto dot = reduce!fmuladd(0.0, z);
889 
890     assert(dot == dotProduct(iota(6), iota([6], 1)));
891 }
892 
893 /// Tensor mutation on-the-fly
894 version(mir_test)
895 unittest
896 {
897     import mir.ndslice.allocation : slice;
898     import mir.ndslice.topology : as, iota;
899     import mir.math.common : fmamath;
900 
901     static @fmamath T fun(T)(const T a, ref T b)
902     {
903         return a + b++;
904     }
905 
906     //| 0 1 2 |
907     //| 3 4 5 |
908     auto sl = iota(2, 3).as!double.slice;
909 
910     auto res = reduce!fun(double(0), sl);
911 
912     assert(res == 15);
913 
914     //| 1 2 3 |
915     //| 4 5 6 |
916     assert(sl == iota([2, 3], 1));
917 }
918 
919 /++
920 Packed slices.
921 
922 Computes minimum value of maximum values for each row.
923 +/
924 version(mir_test)
925 unittest
926 {
927     import mir.math.common;
928     import mir.ndslice.allocation : slice;
929     import mir.ndslice.dynamic : transposed;
930     import mir.ndslice.topology : as, iota, pack, map, universal;
931 
932     alias maxVal = (a) => reduce!fmax(-double.infinity, a);
933     alias minVal = (a) => reduce!fmin(double.infinity, a);
934     alias minimaxVal = (a) => minVal(a.pack!1.map!maxVal);
935 
936     auto sl = iota(2, 3).as!double.slice;
937 
938     // Vectorized computation: row stride equals 1.
939     //| 0 1 2 | => | 2 |
940     //| 3 4 5 | => | 5 | => 2
941     auto res = minimaxVal(sl);
942     assert(res == 2);
943 
944     // Common computation: row stride does not equal 1.
945     //| 0 1 2 |    | 0 3 | => | 3 |
946     //| 3 4 5 | => | 1 4 | => | 4 |
947     //             | 2 5 | => | 5 | => 3
948     auto resT = minimaxVal(sl.universal.transposed);
949     assert(resT == 3);
950 }
951 
952 /// Dlang Range API support.
953 version(mir_test)
954 unittest
955 {
956     import mir.algorithm.iteration: each;
957     import std.range: phobos_iota = iota;
958 
959     int s;
960     // 0 1 2 3
961     4.phobos_iota.each!(i => s += i);
962     assert(s == 6);
963 }
964 
965 @safe pure nothrow @nogc
966 version(mir_test) unittest
967 {
968     import mir.ndslice.topology : iota;
969     auto a = reduce!"a + b"(size_t(7), iota([0, 1], 1));
970     assert(a == 7);
971 }
972 
973 void eachImpl(alias fun, Slices...)(Slices slices)
974 {
975     foreach(ref slice; slices)
976         assert(!slice.empty);
977     do
978     {
979         static if (DimensionCount!(Slices[0]) == 1)
980             mixin(`fun(`~ frontOf!(Slices.length) ~ `);`);
981         else
982             .eachImpl!fun(frontOf2!slices);
983         foreach_reverse(i; Iota!(Slices.length))
984             slices[i].popFront;
985     }
986     while(!slices[0].empty);
987 }
988 
989 void chequerEachImpl(alias fun, Slices...)(Chequer color, Slices slices)
990 {
991     foreach(ref slice; slices)
992         assert(!slice.empty);
993     static if (DimensionCount!(Slices[0]) == 1)
994     {
995         if (color)
996         {
997             foreach_reverse(i; Iota!(Slices.length))
998                 slices[i].popFront;
999             if (slices[0].empty)
1000                 return;
1001         }
1002         eachImpl!fun(strideOf!slices);
1003     }
1004     else
1005     {
1006         do
1007         {
1008             mixin(`.chequerEachImpl!fun(color,` ~ frontOf!(Slices.length) ~ `);`);
1009             color = cast(Chequer)!color;
1010             foreach_reverse(i; Iota!(Slices.length))
1011                 slices[i].popFront;
1012         }
1013         while(!slices[0].empty);
1014     }
1015 }
1016 
1017 /++
1018 The call `each!(fun)(slice1, ..., sliceN)`
1019 evaluates `fun` for each set of elements `x1, ..., xN` in
1020 the borders of `slice1, ..., sliceN` respectively.
1021 
1022 `each` allows to iterate multiple slices in the lockstep.
1023 
1024 Params:
1025     fun = A function.
1026 Note:
1027     $(NDSLICEREF dynamic, transposed) and
1028     $(NDSLICEREF topology, pack) can be used to specify dimensions.
1029 +/
1030 template eachOnBorder(alias fun)
1031 {
1032     import mir.functional: naryFun;
1033     static if (__traits(isSame, naryFun!fun, fun))
1034     /++
1035     Params:
1036         slices = One or more slices.
1037     +/
1038     @fmamath void eachOnBorder(Slices...)(Slices slices)
1039         if (allSatisfy!(isSlice, Slices))
1040     {
1041         import mir.ndslice.traits: isContiguousSlice;
1042         static if (Slices.length > 1)
1043             slices.checkShapesMatch;
1044         if (!slices[0].anyEmpty)
1045         {
1046             alias N = DimensionCount!(Slices[0]);
1047             static if (N == 1)
1048             {
1049                 mixin(`fun(`~ frontOf!(Slices.length) ~ `);`);
1050                 if (slices[0].length > 1)
1051                     fun(backOf!slices);
1052             }
1053             else
1054             static if (anySatisfy!(isContiguousSlice, Slices))
1055             {
1056                 import mir.ndslice.topology: canonical;
1057                 template f(size_t i)
1058                 {
1059                     static if (isContiguousSlice!(Slices[i]))
1060                         auto f () { return canonical(slices[i]); }
1061                     else
1062                         alias f = slices[i];
1063                 }
1064                 eachOnBorder(staticMap!(f, Iota!(Slices.length)));
1065             }
1066             else
1067             {
1068                 foreach (dimension; Iota!N)
1069                 {
1070                     eachImpl!fun(frontOfD!(dimension, slices));
1071                     foreach_reverse(ref slice; slices)
1072                         slice.popFront!dimension;
1073                     if (slices[0].empty!dimension)
1074                         return;
1075                     eachImpl!fun(backOfD!(dimension, slices));
1076                     foreach_reverse(ref slice; slices)
1077                         slice.popBack!dimension;
1078                     if (slices[0].empty!dimension)
1079                         return;
1080                 }
1081             }
1082         }
1083     }
1084     else
1085         alias eachOnBorder = .eachOnBorder!(naryFun!fun);
1086 }
1087 
1088 ///
1089 @safe pure nothrow 
1090 version(mir_test) unittest
1091 {
1092     import mir.ndslice.allocation : slice;
1093     import mir.ndslice.topology : repeat, iota;
1094 
1095     auto sl = [3, 4].iota.slice;
1096     auto zeros = repeat(0, [3, 4]);
1097 
1098     sl.eachOnBorder!"a = b"(zeros);
1099 
1100     assert(sl == 
1101         [[0, 0, 0 ,0],
1102          [0, 5, 6, 0],
1103          [0, 0, 0 ,0]]);
1104     
1105     sl.eachOnBorder!"a = 1";
1106     sl[0].eachOnBorder!"a = 2";
1107 
1108     assert(sl == 
1109         [[2, 1, 1, 2],
1110          [1, 5, 6, 1],
1111          [1, 1, 1 ,1]]);
1112 }
1113 
1114 /++
1115 The call `each!(fun)(slice1, ..., sliceN)`
1116 evaluates `fun` for each set of elements `x1, ..., xN` in
1117 `slice1, ..., sliceN` respectively.
1118 
1119 `each` allows to iterate multiple slices in the lockstep.
1120 Params:
1121     fun = A function.
1122 Note:
1123     $(NDSLICEREF dynamic, transposed) and
1124     $(NDSLICEREF topology, pack) can be used to specify dimensions.
1125 See_Also:
1126     This is functionally similar to $(LREF reduce) but has not seed.
1127 +/
1128 template each(alias fun)
1129 {
1130     import mir.functional: naryFun;
1131     static if (__traits(isSame, naryFun!fun, fun))
1132     {
1133         /++
1134         Params:
1135             slices = One or more slices, ranges, and arrays.
1136         +/
1137         @fmamath auto each(Slices...)(Slices slices)
1138             if (Slices.length && !is(Slices[0] : Chequer))
1139         {
1140             static if (Slices.length > 1)
1141                 slices.checkShapesMatch;
1142             static if (areAllContiguousSlices!Slices)
1143             {
1144                 import mir.ndslice.topology: flattened;
1145                 .each!fun(allFlattened!(allLightScope!slices));
1146             }
1147             else
1148             {
1149                 if (slices[0].anyEmpty)
1150                     return;
1151                 eachImpl!fun(allLightScope!slices);
1152             }
1153         }
1154 
1155         /++
1156         Iterates elements of selected $(LREF Chequer) color.
1157         Params:
1158             color = $(LREF Chequer).
1159             slices = One or more slices.
1160         +/
1161         @fmamath auto each(Slices...)(Chequer color, Slices slices)
1162             if (Slices.length && allSatisfy!(isSlice, Slices))
1163         {
1164             static if (Slices.length > 1)
1165                 slices.checkShapesMatch;
1166             if (slices[0].anyEmpty)
1167                 return;
1168             chequerEachImpl!fun(color, allLightScope!slices);
1169         }
1170     }
1171     else
1172         alias each = .each!(naryFun!fun);
1173 }
1174 
1175 /// Ranges and arrays
1176 version(mir_test)
1177 unittest
1178 {
1179     auto ar = [1, 2, 3];
1180     ar.each!"a *= 2";
1181     assert (ar == [2, 4, 6]);
1182 }
1183 
1184 /// Single slice, multiply-add
1185 version(mir_test)
1186 unittest
1187 {
1188     import mir.ndslice.allocation : slice;
1189     import mir.ndslice.topology : as, iota;
1190 
1191     //| 0 1 2 |
1192     //| 3 4 5 |
1193     auto sl = iota(2, 3).as!double.slice;
1194 
1195     sl.each!((ref a) { a = a * 10 + 5; });
1196 
1197     assert(sl ==
1198         [[ 5, 15, 25],
1199          [35, 45, 55]]);
1200 }
1201 
1202 /// Swap two slices
1203 version(mir_test)
1204 unittest
1205 {
1206     import mir.utility : swap;
1207     import mir.ndslice.allocation : slice;
1208     import mir.ndslice.topology : as, iota;
1209 
1210     //| 0 1 2 |
1211     //| 3 4 5 |
1212     auto a = iota([2, 3], 0).as!double.slice;
1213     //| 10 11 12 |
1214     //| 13 14 15 |
1215     auto b = iota([2, 3], 10).as!double.slice;
1216 
1217     each!swap(a, b);
1218 
1219     assert(a == iota([2, 3], 10));
1220     assert(b == iota([2, 3], 0));
1221 }
1222 
1223 /// Swap two zipped slices
1224 version(mir_test)
1225 unittest
1226 {
1227     import mir.utility : swap;
1228     import mir.ndslice.allocation : slice;
1229     import mir.ndslice.topology : as, zip, iota;
1230 
1231     //| 0 1 2 |
1232     //| 3 4 5 |
1233     auto a = iota([2, 3], 0).as!double.slice;
1234     //| 10 11 12 |
1235     //| 13 14 15 |
1236     auto b = iota([2, 3], 10).as!double.slice;
1237 
1238     auto z = zip(a, b);
1239 
1240     z.each!(z => swap(z.a, z.b));
1241 
1242     assert(a == iota([2, 3], 10));
1243     assert(b == iota([2, 3], 0));
1244 }
1245 
1246 @safe pure nothrow
1247 version(mir_test) unittest
1248 {
1249     import mir.ndslice.topology : iota;
1250     size_t i;
1251     iota(0, 2).each!((a){i++;});
1252     assert(i == 0);
1253 }
1254 
1255 /++
1256 The call `eachUploPair!(fun)(matrix)`
1257 evaluates `fun` for each pair (`matrix[j, i]`, `matrix[i, j]`),
1258 for i <= j (default) or i < j (if includeDiagonal is false).
1259 
1260 Params:
1261     fun = A function.
1262     includeDiagonal = true if applying function to diagonal,
1263                       false (default) otherwise.
1264 +/
1265 template eachUploPair(alias fun, bool includeDiagonal = false)
1266 {
1267     import mir.functional: naryFun;
1268     static if (__traits(isSame, naryFun!fun, fun))
1269     {
1270         /++
1271         Params:
1272             matrix = Square matrix.
1273         +/
1274         auto eachUploPair(Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) matrix)
1275         in
1276         {
1277             assert(matrix.length!0 == matrix.length!1, "matrix must be square.");
1278         }
1279         do
1280         {
1281             static if (kind == Contiguous)
1282             {
1283                 import mir.ndslice.topology: canonical;
1284                 .eachUploPair!(fun, includeDiagonal)(matrix.canonical);
1285             }
1286             else
1287             {
1288                 static if (includeDiagonal == true)
1289                 {
1290                     if (matrix.length) do
1291                     {
1292                         eachImpl!fun(matrix.lightScope.front!0, matrix.lightScope.front!1);
1293                         matrix.popFront!1;
1294                         matrix.popFront!0;
1295                         // hint for optimizer
1296                         matrix._lengths[1] = matrix._lengths[0];
1297                     }
1298                     while (matrix.length);
1299                 }
1300                 else
1301                 {
1302                     if (matrix.length) for(;;)
1303                     {
1304                         assert(!matrix.empty!0);
1305                         assert(!matrix.empty!1);
1306                         auto l = matrix.lightScope.front!1;
1307                         auto u = matrix.lightScope.front!0;
1308                         matrix.popFront!1;
1309                         matrix.popFront!0;
1310                         l.popFront;
1311                         u.popFront;
1312                         // hint for optimizer
1313                         matrix._lengths[1] = matrix._lengths[0] = l._lengths[0] = u._lengths[0];
1314                         if (u.length == 0)
1315                             break;
1316                         eachImpl!fun(u, l);
1317                     }
1318                 }
1319             }
1320         }
1321     }
1322     else
1323     {
1324         alias eachUploPair = .eachUploPair!(naryFun!fun, includeDiagonal);
1325     }
1326 }
1327 
1328 /// Transpose matrix in place.
1329 version(mir_test)
1330 unittest
1331 {
1332     import mir.ndslice.allocation: slice;
1333     import mir.ndslice.topology: iota, universal;
1334     import mir.ndslice.dynamic: transposed;
1335     import mir.utility: swap;
1336 
1337     auto m = iota(4, 4).slice;
1338 
1339     m.eachUploPair!swap;
1340 
1341     assert(m == iota(4, 4).universal.transposed);
1342 }
1343 
1344 /// Reflect Upper matrix part to lower part.
1345 version(mir_test)
1346 unittest
1347 {
1348     import mir.ndslice.allocation: slice;
1349     import mir.ndslice.topology: iota, universal;
1350     import mir.ndslice.dynamic: transposed;
1351     import mir.utility: swap;
1352 
1353     // 0 1 2
1354     // 3 4 5
1355     // 6 7 8
1356     auto m = iota(3, 3).slice;
1357 
1358     m.eachUploPair!((u, ref l) { l = u; });
1359 
1360     assert(m == [
1361         [0, 1, 2],
1362         [1, 4, 5],
1363         [2, 5, 8]]);
1364 }
1365 
1366 /// Fill lower triangle and diagonal with zeroes.
1367 version(mir_test)
1368 unittest
1369 {
1370     import mir.ndslice.allocation: slice;
1371     import mir.ndslice.topology: iota;
1372 
1373     // 1 2 3
1374     // 4 5 6
1375     // 7 8 9
1376     auto m = iota([3, 3], 1).slice;
1377 
1378     m.eachUploPair!((u, ref l) { l = 0; }, true);
1379 
1380     assert(m == [
1381         [0, 2, 3],
1382         [0, 0, 6],
1383         [0, 0, 0]]);
1384 }
1385 
1386 version(mir_test)
1387 unittest
1388 {
1389     import mir.ndslice.allocation: slice;
1390     import mir.ndslice.topology: iota;
1391 
1392     // 0 1 2
1393     // 3 4 5
1394     // 6 7 8
1395     auto m = iota(3, 3).slice;
1396     m.eachUploPair!((u, ref l) { l = l + 1; }, true);
1397     assert(m == [
1398         [1, 1, 2],
1399         [4, 5, 5],
1400         [7, 8, 9]]);
1401 }
1402 
1403 version(mir_test)
1404 unittest
1405 {
1406     import mir.ndslice.allocation: slice;
1407     import mir.ndslice.topology: iota;
1408 
1409     // 0 1 2
1410     // 3 4 5
1411     // 6 7 8
1412     auto m = iota(3, 3).slice;
1413     m.eachUploPair!((u, ref l) { l = l + 1; }, false);
1414 
1415     assert(m == [
1416         [0, 1, 2],
1417         [4, 4, 5],
1418         [7, 8, 8]]);
1419 }
1420 
1421 /++
1422 Checks if the matrix is symmetric.
1423 +/
1424 template isSymmetric(alias fun = "a == b")
1425 {
1426     import mir.functional: naryFun;
1427     static if (__traits(isSame, naryFun!fun, fun))
1428     /++
1429     Params:
1430         matrix = 2D ndslice.
1431     +/
1432     bool isSymmetric(Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) matrix)
1433     {
1434         static if (kind == Contiguous)
1435         {
1436             import mir.ndslice.topology: canonical;
1437             return .isSymmetric!fun(matrix.canonical);
1438         }
1439         else
1440         {
1441             if (matrix.length!0 != matrix.length!1)
1442                 return false;
1443             if (matrix.length) do
1444             {
1445                 if (!allImpl!fun(matrix.lightScope.front!0, matrix.lightScope.front!1))
1446                 {
1447                     return false;
1448                 }
1449                 matrix.popFront!1;
1450                 matrix.popFront!0;
1451                 matrix._lengths[1] = matrix._lengths[0];
1452             }
1453             while (matrix.length);
1454             return true;
1455         }
1456     }
1457     else
1458         alias isSymmetric = .isSymmetric!(naryFun!fun);
1459 }
1460 
1461 ///
1462 version(mir_test)
1463 unittest
1464 {
1465     import mir.ndslice.slice: sliced;
1466     import mir.ndslice.topology: iota;
1467     assert(iota(2, 2).isSymmetric == false);
1468 
1469     assert(
1470         [1, 2,
1471          2, 3].sliced(2, 2).isSymmetric == true);
1472 }
1473 
1474 bool minPosImpl(alias fun, Iterator, size_t N, SliceKind kind)(scope ref size_t[N] backwardIndex, scope ref Iterator iterator, Slice!(Iterator, N, kind) slice)
1475 {
1476     bool found;
1477     do
1478     {
1479         static if (slice.shape.length == 1)
1480         {
1481             if (fun(*slice._iterator, *iterator))
1482             {
1483                 backwardIndex[0] = slice.length;
1484                 iterator = slice._iterator;
1485                 found = true;
1486             }
1487         }
1488         else
1489         {
1490             if (minPosImpl!(fun, LightScopeOf!Iterator, N - 1, kind)(backwardIndex[1 .. $], iterator, lightScope(slice).front))
1491             {
1492                 backwardIndex[0] = slice.length;
1493                 found = true;
1494             }
1495         }
1496         slice.popFront;
1497     }
1498     while(!slice.empty);
1499     return found;
1500 }
1501 
1502 bool[2] minmaxPosImpl(alias fun, Iterator, size_t N, SliceKind kind)(scope ref size_t[2][N] backwardIndex, scope ref Iterator[2] iterator, Slice!(Iterator, N, kind) slice)
1503 {
1504     bool[2] found;
1505     do
1506     {
1507         static if (slice.shape.length == 1)
1508         {
1509             if (fun(*slice._iterator, *iterator[0]))
1510             {
1511                 backwardIndex[0][0] = slice.length;
1512                 iterator[0] = slice._iterator;
1513                 found[0] = true;
1514             }
1515             else
1516             if (fun(*iterator[1], *slice._iterator))
1517             {
1518                 backwardIndex[0][1] = slice.length;
1519                 iterator[1] = slice._iterator;
1520                 found[1] = true;
1521             }
1522         }
1523         else
1524         {
1525             auto r = minmaxPosImpl!(fun, LightScopeOf!Iterator, N - 1, kind)(backwardIndex[1 .. $], iterator, lightScope(slice).front);
1526             if (r[0])
1527             {
1528                 backwardIndex[0][0] = slice.length;
1529             }
1530             if (r[1])
1531             {
1532                 backwardIndex[0][1] = slice.length;
1533             }
1534         }
1535         slice.popFront;
1536     }
1537     while(!slice.empty);
1538     return found;
1539 }
1540 
1541 /++
1542 Finds a positions (ndslices) such that
1543 `position[0].first` is minimal and `position[1].first` is maximal elements in the slice.
1544 
1545 Position is sub-ndslice of the same dimension in the right-$(RPAREN)down-$(RPAREN)etc$(LPAREN)$(LPAREN) corner.
1546 
1547 Params:
1548     pred = A predicate.
1549 
1550 See_also:
1551     $(LREF minmaxIndex),
1552     $(LREF minPos),
1553     $(LREF maxPos),
1554     $(NDSLICEREF slice, Slice.backward).
1555 +/
1556 template minmaxPos(alias pred = "a < b")
1557 {
1558     import mir.functional: naryFun;
1559     static if (__traits(isSame, naryFun!pred, pred))
1560     /++
1561     Params:
1562         slice = ndslice.
1563     Returns:
1564         2 subslices with minimal and maximal `first` elements.
1565     +/
1566     @fmamath Slice!(Iterator, N, kind == Contiguous && N > 1 ? Canonical : kind)[2]
1567         minmaxPos(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
1568     {
1569         typeof(return) pret;
1570         if (!slice.anyEmpty)
1571         {
1572             size_t[2][N] ret;
1573             auto scopeSlice = lightScope(slice);
1574             auto it = scopeSlice._iterator;
1575             LightScopeOf!Iterator[2] iterator = [it, it];
1576             minmaxPosImpl!(pred, LightScopeOf!Iterator, N, kind)(ret, iterator, scopeSlice);
1577             foreach (i; Iota!N)
1578             {
1579                 pret[0]._lengths[i] = ret[i][0];
1580                 pret[1]._lengths[i] = ret[i][1];
1581             }
1582             pret[0]._iterator = slice._iterator + (iterator[0] - scopeSlice._iterator);
1583             pret[1]._iterator = slice._iterator + (iterator[1] - scopeSlice._iterator);
1584         }
1585         auto strides = slice.strides;
1586         foreach(i; Iota!(0, pret[0].S))
1587         {
1588             pret[0]._strides[i] = strides[i];
1589             pret[1]._strides[i] = strides[i];
1590         }
1591         return pret;
1592     }
1593     else
1594         alias minmaxPos = .minmaxPos!(naryFun!pred);
1595 }
1596 
1597 ///
1598 version(mir_test)
1599 unittest
1600 {
1601     import mir.ndslice.slice: sliced;
1602     auto s = [
1603         2, 6, 4, -3,
1604         0, -4, -3, 3,
1605         -3, -2, 7, 2,
1606         ].sliced(3, 4);
1607 
1608     auto pos = s.minmaxPos;
1609 
1610     assert(pos[0] == s[$ - 2 .. $, $ - 3 .. $]);
1611     assert(pos[1] == s[$ - 1 .. $, $ - 2 .. $]);
1612 
1613     assert(pos[0].first == -4);
1614     assert(s.backward(pos[0].shape) == -4);
1615     assert(pos[1].first ==  7);
1616     assert(s.backward(pos[1].shape) ==  7);
1617 }
1618 
1619 /++
1620 Finds a backward indices such that
1621 `slice[indices[0]]` is minimal and `slice[indices[1]]` is maximal elements in the slice.
1622 
1623 Params:
1624     pred = A predicate.
1625 
1626 See_also:
1627     $(LREF minmaxIndex),
1628     $(LREF minPos),
1629     $(LREF maxPos),
1630     $(REF Slice.backward, mir,ndslice,slice).
1631 +/
1632 template minmaxIndex(alias pred = "a < b")
1633 {
1634     import mir.functional: naryFun;
1635     static if (__traits(isSame, naryFun!pred, pred))
1636     /++
1637     Params:
1638         slice = ndslice.
1639     Returns:
1640         Subslice with minimal (maximal) `first` element.
1641     +/
1642     @fmamath size_t[N][2] minmaxIndex(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
1643     {
1644         typeof(return) pret = size_t.max;
1645         if (!slice.anyEmpty)
1646         {
1647             auto shape = slice.shape;
1648             size_t[2][N] ret;
1649             foreach (i; Iota!N)
1650             {
1651                 ret[i][1] = ret[i][0] = shape[i];
1652             }
1653             auto scopeSlice = lightScope(slice);
1654             auto it = scopeSlice._iterator;
1655             LightScopeOf!Iterator[2] iterator = [it, it];
1656             minmaxPosImpl!(pred, LightScopeOf!Iterator, N, kind)(ret, iterator, scopeSlice);
1657             foreach (i; Iota!N)
1658             {
1659                 pret[0][i] = slice._lengths[i] - ret[i][0];
1660                 pret[1][i] = slice._lengths[i] - ret[i][1];
1661             }
1662         }
1663         return pret;
1664     }
1665     else
1666         alias minmaxIndex = .minmaxIndex!(naryFun!pred);
1667 }
1668 
1669 ///
1670 version(mir_test)
1671 unittest
1672 {
1673     import mir.ndslice.slice: sliced;
1674     auto s = [
1675         2, 6, 4, -3,
1676         0, -4, -3, 3,
1677         -3, -2, 7, 8,
1678         ].sliced(3, 4);
1679 
1680     auto indices = s.minmaxIndex;
1681 
1682     assert(indices == [[1, 1], [2, 3]]);
1683     assert(s[indices[0]] == -4);
1684     assert(s[indices[1]] ==  8);
1685 }
1686 
1687 /++
1688 Finds a backward index such that
1689 `slice.backward(index)` is minimal(maximal).
1690 
1691 Params:
1692     pred = A predicate.
1693 
1694 See_also:
1695     $(LREF minIndex),
1696     $(LREF maxPos),
1697     $(LREF maxIndex),
1698     $(REF Slice.backward, mir,ndslice,slice).
1699 +/
1700 template minPos(alias pred = "a < b")
1701 {
1702     import mir.functional: naryFun;
1703     static if (__traits(isSame, naryFun!pred, pred))
1704     /++
1705     Params:
1706         slice = ndslice.
1707     Returns:
1708         Multidimensional backward index such that element is minimal(maximal).
1709         Backward index equals zeros, if slice is empty.
1710     +/
1711     @fmamath Slice!(Iterator, N, kind == Contiguous && N > 1 ? Canonical : kind)
1712         minPos(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
1713     {
1714         typeof(return) ret;
1715         auto iterator = slice.lightScope._iterator;
1716         if (!slice.anyEmpty)
1717         {
1718             minPosImpl!(pred, LightScopeOf!Iterator, N, kind)(ret._lengths, iterator, lightScope(slice));
1719             ret._iterator = slice._iterator + (iterator - slice.lightScope._iterator);
1720         }
1721         auto strides = slice.strides;
1722         foreach(i; Iota!(0, ret.S))
1723         {
1724             ret._strides[i] = strides[i];
1725         }
1726         return ret;
1727     }
1728     else
1729         alias minPos = .minPos!(naryFun!pred);
1730 }
1731 
1732 /// ditto
1733 template maxPos(alias pred = "a < b")
1734 {
1735     import mir.functional: naryFun, reverseArgs;
1736     alias maxPos = minPos!(reverseArgs!(naryFun!pred));
1737 }
1738 
1739 ///
1740 version(mir_test)
1741 unittest
1742 {
1743     import mir.ndslice.slice: sliced;
1744     auto s = [
1745         2, 6, 4, -3,
1746         0, -4, -3, 3,
1747         -3, -2, 7, 2,
1748         ].sliced(3, 4);
1749 
1750     auto pos = s.minPos;
1751 
1752     assert(pos == s[$ - 2 .. $, $ - 3 .. $]);
1753     assert(pos.first == -4);
1754     assert(s.backward(pos.shape) == -4);
1755 
1756     pos = s.maxPos;
1757 
1758     assert(pos == s[$ - 1 .. $, $ - 2 .. $]);
1759     assert(pos.first == 7);
1760     assert(s.backward(pos.shape) == 7);
1761 }
1762 
1763 /++
1764 Finds an index such that
1765 `slice[index]` is minimal(maximal).
1766 
1767 Params:
1768     pred = A predicate.
1769 
1770 See_also:
1771     $(LREF minIndex),
1772     $(LREF maxPos),
1773     $(LREF maxIndex).
1774 +/
1775 template minIndex(alias pred = "a < b")
1776 {
1777     import mir.functional: naryFun;
1778     static if (__traits(isSame, naryFun!pred, pred))
1779     /++
1780     Params:
1781         slice = ndslice.
1782     Returns:
1783         Multidimensional index such that element is minimal(maximal).
1784         Index elements equal to `size_t.max`, if slice is empty.
1785     +/
1786     @fmamath size_t[N] minIndex(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
1787     {
1788         size_t[N] ret = size_t.max;
1789         if (!slice.anyEmpty)
1790         {
1791             ret = slice.shape;
1792             auto scopeSlice = lightScope(slice);
1793             auto iterator = scopeSlice._iterator;
1794             minPosImpl!(pred, LightScopeOf!Iterator, N, kind)(ret, iterator, scopeSlice);
1795             foreach (i; Iota!N)
1796                 ret[i] = slice._lengths[i] - ret[i];
1797         }
1798         return ret;
1799     }
1800     else
1801         alias minIndex = .minIndex!(naryFun!pred);
1802 }
1803 
1804 /// ditto
1805 template maxIndex(alias pred = "a < b")
1806 {
1807     import mir.functional: naryFun, reverseArgs;
1808     alias maxIndex = minIndex!(reverseArgs!(naryFun!pred));
1809 }
1810 
1811 ///
1812 version(mir_test)
1813 unittest
1814 {
1815     import mir.ndslice.slice: sliced;
1816     auto s = [
1817         2, 6, 4, -3,
1818         0, -4, -3, 3,
1819         -3, -2, 7, 8,
1820         ].sliced(3, 4);
1821 
1822     auto index = s.minIndex;
1823 
1824     assert(index == [1, 1]);
1825     assert(s[index] == -4);
1826 
1827     index = s.maxIndex;
1828 
1829     assert(index == [2, 3]);
1830     assert(s[index] == 8);
1831 }
1832 
1833 ///
1834 version(mir_test)
1835 unittest
1836 {
1837     import mir.ndslice.slice: sliced;
1838     auto s = [
1839         -8, 6, 4, -3,
1840         0, -4, -3, 3,
1841         -3, -2, 7, 8,
1842         ].sliced(3, 4);
1843 
1844     auto index = s.minIndex;
1845 
1846     assert(index == [0, 0]);
1847     assert(s[index] == -8);
1848 }
1849 
1850 version(mir_test)
1851 unittest
1852 {
1853     import mir.ndslice.slice: sliced;
1854     auto s = [
1855             0, 1, 2, 3,
1856             4, 5, 6, 7,
1857             8, 9, 10, 11
1858             ].sliced(3, 4);
1859 
1860     auto index = s.minIndex;
1861     assert(index == [0, 0]);
1862     assert(s[index] == 0);
1863 
1864     index = s.maxIndex;
1865     assert(index == [2, 3]);
1866     assert(s[index] == 11);
1867 }
1868 
1869 bool findImpl(alias fun, size_t N, Slices...)(scope ref size_t[N] backwardIndex, Slices slices)
1870     if (Slices.length)
1871 {
1872     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(FieldIterator!(BitField!(Field, I))), Field, I))
1873     {
1874         auto cnt = BitSliceAccelerator!(Field, I)(slices[0]).cttz;
1875         if (cnt = -1)
1876             return false;
1877         backwardIndex[0] = slices[0].length - cnt;
1878     }
1879     else
1880     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))), Field, I))
1881     {
1882         import mir.ndslice.topology: retro;
1883         auto cnt = BitSliceAccelerator!(Field, I)(slices[0].retro).ctlz;
1884         if (cnt = -1)
1885             return false;
1886         backwardIndex[0] = slices[0].length - cnt;
1887     }
1888     else
1889     {
1890         do
1891         {
1892             static if (DimensionCount!(Slices[0]) == 1)
1893             {
1894                 if (mixin(`fun(`~ frontOf!(Slices.length) ~ `)`))
1895                 {
1896                     backwardIndex[0] = slices[0].length;
1897                     return true;
1898                 }
1899             }
1900             else
1901             {
1902                 if (mixin(`findImpl!fun(backwardIndex[1 .. $],` ~ frontOf!(Slices.length) ~ `)`))
1903                 {
1904                     backwardIndex[0] = slices[0].length;
1905                     return true;
1906                 }
1907             }
1908             foreach_reverse(ref slice; slices)
1909                 slice.popFront;
1910         }
1911         while(!slices[0].empty);
1912         return false;
1913     }
1914 }
1915 
1916 /++
1917 Finds an index such that
1918 `pred(slices[0][index], ..., slices[$-1][index])` is `true`.
1919 
1920 Params:
1921     pred = A predicate.
1922 
1923 See_also:
1924     $(LREF find),
1925     $(LREF any).
1926 Optimization:
1927     `findIndex!"a"` has accelerated specialization for slices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
1928 +/
1929 template findIndex(alias pred)
1930 {
1931     import mir.functional: naryFun;
1932     static if (__traits(isSame, naryFun!pred, pred))
1933     /++
1934     Params:
1935         slices = One or more slices.
1936     Returns:
1937         Multidimensional index such that the predicate is true.
1938         Index equals `size_t.max`, if the predicate evaluates `false` for all indices.
1939     Constraints:
1940         All slices must have the same shape.
1941     +/
1942     @fmamath Select!(DimensionCount!(Slices[0]) > 1, size_t[DimensionCount!(Slices[0])], size_t) findIndex(Slices...)(Slices slices)
1943         if (Slices.length)
1944     {
1945         static if (Slices.length > 1)
1946             slices.checkShapesMatch;
1947         size_t[DimensionCount!(Slices[0])] ret = -1;
1948         auto lengths = slices[0].shape;
1949         if (!slices[0].anyEmpty && findImpl!pred(ret, allLightScope!slices))
1950             foreach (i; Iota!(DimensionCount!(Slices[0])))
1951                 ret[i] = lengths[i] - ret[i];
1952         static if (DimensionCount!(Slices[0]) > 1)
1953             return ret;
1954         else
1955             return ret[0];
1956     }
1957     else
1958         alias findIndex = .findIndex!(naryFun!pred);
1959 }
1960 
1961 /// Ranges and arrays
1962 version(mir_test)
1963 unittest
1964 {
1965     import std.range : iota;
1966     // 0 1 2 3 4 5
1967     auto sl = iota(5);
1968     size_t index = sl.findIndex!"a == 3";
1969 
1970     assert(index == 3);
1971     assert(sl[index] == 3);
1972 
1973     assert(sl.findIndex!(a => a == 8) == size_t.max);
1974 }
1975 
1976 ///
1977 @safe pure nothrow @nogc
1978 version(mir_test) unittest
1979 {
1980     import mir.ndslice.topology : iota;
1981     // 0 1 2
1982     // 3 4 5
1983     auto sl = iota(2, 3);
1984     size_t[2] index = sl.findIndex!(a => a == 3);
1985 
1986     assert(sl[index] == 3);
1987 
1988     index = sl.findIndex!"a == 6";
1989     assert(index[0] == size_t.max);
1990     assert(index[1] == size_t.max);
1991 }
1992 
1993 /++
1994 Finds a backward index such that
1995 `pred(slices[0].backward(index), ..., slices[$-1].backward(index))` is `true`.
1996 
1997 Params:
1998     pred = A predicate.
1999 
2000 Optimization:
2001     To check if any element was found
2002     use the last dimension (row index).
2003     This will slightly optimize the code.
2004 --------
2005 if (backwardIndex)
2006 {
2007     auto elem1 = slice1.backward(backwardIndex);
2008     //...
2009     auto elemK = sliceK.backward(backwardIndex);
2010 }
2011 else
2012 {
2013     // not found
2014 }
2015 --------
2016 
2017 See_also:
2018     $(LREF findIndex),
2019     $(LREF any),
2020     $(REF Slice.backward, mir,ndslice,slice).
2021 
2022 Optimization:
2023     `find!"a"` has accelerated specialization for slices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
2024 +/
2025 template find(alias pred)
2026 {
2027     import mir.functional: naryFun;
2028     static if (__traits(isSame, naryFun!pred, pred))
2029     /++
2030     Params:
2031         slices = One or more slices.
2032     Returns:
2033         Multidimensional backward index such that the predicate is true.
2034         Backward index equals zeros, if the predicate evaluates `false` for all indices.
2035     Constraints:
2036         All slices must have the same shape.
2037     +/
2038     @fmamath Select!(DimensionCount!(Slices[0]) > 1, size_t[DimensionCount!(Slices[0])], size_t) find(Slices...)(Slices slices)
2039         if (Slices.length && allSatisfy!(hasShape, Slices))
2040     {
2041         static if (Slices.length > 1)
2042             slices.checkShapesMatch;
2043         size_t[DimensionCount!(Slices[0])] ret;
2044         if (!slices[0].anyEmpty)
2045             findImpl!pred(ret, allLightScope!slices);
2046         static if (DimensionCount!(Slices[0]) > 1)
2047             return ret;
2048         else
2049             return ret[0];
2050     }
2051     else
2052         alias find = .find!(naryFun!pred);
2053 }
2054 
2055 /// Ranges and arrays
2056 version(mir_test)
2057 unittest
2058 {
2059     import std.range : iota;
2060 
2061     auto sl = iota(10);
2062     size_t index = sl.find!"a == 3";
2063 
2064     assert(sl[$ - index] == 3);
2065 }
2066 
2067 ///
2068 @safe pure nothrow @nogc
2069 version(mir_test) unittest
2070 {
2071     import mir.ndslice.topology : iota;
2072     // 0 1 2
2073     // 3 4 5
2074     auto sl = iota(2, 3);
2075     size_t[2] bi = sl.find!"a == 3";
2076     assert(sl.backward(bi) == 3);
2077     assert(sl[$ - bi[0], $ - bi[1]] == 3);
2078 
2079     bi = sl.find!"a == 6";
2080     assert(bi[0] == 0);
2081     assert(bi[1] == 0);
2082 }
2083 
2084 /// Multiple slices
2085 @safe pure nothrow @nogc
2086 version(mir_test) unittest
2087 {
2088     import mir.ndslice.topology : iota;
2089 
2090     // 0 1 2
2091     // 3 4 5
2092     auto a = iota(2, 3);
2093     // 10 11 12
2094     // 13 14 15
2095     auto b = iota([2, 3], 10);
2096 
2097     size_t[2] bi = find!((a, b) => a * b == 39)(a, b);
2098     assert(a.backward(bi) == 3);
2099     assert(b.backward(bi) == 13);
2100 }
2101 
2102 /// Zipped slices
2103 @safe pure nothrow
2104 version(mir_test) unittest
2105 {
2106     import mir.ndslice.topology : iota, zip;
2107 
2108     // 0 1 2
2109     // 3 4 5
2110     auto a = iota(2, 3);
2111     // 10 11 12
2112     // 13 14 15
2113     auto b = iota([2, 3], 10);
2114 
2115     size_t[2] bi = zip!true(a, b).find!"a.a * a.b == 39";
2116 
2117     assert(a.backward(bi) == 3);
2118     assert(b.backward(bi) == 13);
2119 }
2120 
2121 /// Mutation on-the-fly
2122 pure nothrow
2123 version(mir_test) unittest
2124 {
2125     import mir.ndslice.allocation : slice;
2126     import mir.ndslice.topology : as, iota;
2127 
2128     // 0 1 2
2129     // 3 4 5
2130     auto sl = iota(2, 3).as!double.slice;
2131 
2132     static bool pred(T)(ref T a)
2133     {
2134         if (a == 5)
2135             return true;
2136         a = 8;
2137         return false;
2138     }
2139 
2140     size_t[2] bi = sl.find!pred;
2141 
2142     assert(bi == [1, 1]);
2143     assert(sl.backward(bi) == 5);
2144 
2145     import mir.test;
2146     // sl was changed
2147     sl.should == [[8, 8, 8],
2148                   [8, 8, 5]];
2149 }
2150 
2151 @safe pure nothrow
2152 version(mir_test) unittest
2153 {
2154     import mir.ndslice.topology : iota;
2155     size_t i;
2156     size_t[2] bi = iota(2, 0).find!((elem){i++; return true;});
2157     assert(i == 0);
2158     assert(bi == [0, 0]);
2159 }
2160 
2161 size_t anyImpl(alias fun, Slices...)(Slices slices)
2162     if (Slices.length)
2163 {
2164     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(FieldIterator!(BitField!(Field, I))), Field, I))
2165     {
2166         return BitSliceAccelerator!(Field, I)(slices[0]).any;
2167     }
2168     else
2169     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))), Field, I))
2170     {
2171         // pragma(msg, S);
2172         import mir.ndslice.topology: retro;
2173         return .anyImpl!fun(lightScope(slices[0]).retro);
2174     }
2175     else
2176     {
2177         do
2178         {
2179             static if (DimensionCount!(Slices[0]) == 1)
2180             {
2181                 if (mixin(`fun(`~ frontOf!(Slices.length) ~ `)`))
2182                     return true;
2183             }
2184             else
2185             {
2186                 if (anyImpl!fun(frontOf2!slices))
2187                     return true;
2188             }
2189             foreach_reverse(ref slice; slices)
2190                 slice.popFront;
2191         }
2192         while(!slices[0].empty);
2193         return false;
2194     }
2195 }
2196 
2197 /++
2198 Like $(LREF find), but only returns whether or not the search was successful.
2199 
2200 Params:
2201     pred = The predicate.
2202 Optimization:
2203     `any!"a"` has accelerated specialization for slices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
2204 +/
2205 template any(alias pred = "a")
2206 {
2207     import mir.functional: naryFun;
2208     static if (__traits(isSame, naryFun!pred, pred))
2209     /++
2210     Params:
2211         slices = One or more slices, ranges, and arrays.
2212     Returns:
2213         `true` if the search was successful and `false` otherwise.
2214     Constraints:
2215         All slices must have the same shape.
2216     +/
2217     @fmamath bool any(Slices...)(Slices slices)
2218         if ((Slices.length == 1 || !__traits(isSame, pred, "a")) && Slices.length)
2219     {
2220         static if (Slices.length > 1)
2221             slices.checkShapesMatch;
2222         static if (areAllContiguousSlices!Slices)
2223         {
2224             import mir.ndslice.topology: flattened;
2225             return .any!pred(allFlattened!(allLightScope!slices));
2226         }
2227         else
2228         {
2229             return !slices[0].anyEmpty && anyImpl!pred(allLightScope!slices);
2230         }
2231     }
2232     else
2233         alias any = .any!(naryFun!pred);
2234 }
2235 
2236 /// Ranges and arrays
2237 @safe pure nothrow @nogc
2238 version(mir_test) unittest
2239 {
2240     import std.range : iota;
2241     // 0 1 2 3 4 5
2242     auto r = iota(6);
2243 
2244     assert(r.any!"a == 3");
2245     assert(!r.any!"a == 6");
2246 }
2247 
2248 ///
2249 @safe pure nothrow @nogc
2250 version(mir_test) unittest
2251 {
2252     import mir.ndslice.topology : iota;
2253     // 0 1 2
2254     // 3 4 5
2255     auto sl = iota(2, 3);
2256 
2257     assert(sl.any!"a == 3");
2258     assert(!sl.any!"a == 6");
2259 }
2260 
2261 /// Multiple slices
2262 @safe pure nothrow @nogc
2263 version(mir_test) unittest
2264 {
2265     import mir.ndslice.topology : iota;
2266 
2267     // 0 1 2
2268     // 3 4 5
2269     auto a = iota(2, 3);
2270     // 10 11 12
2271     // 13 14 15
2272     auto b = iota([2, 3], 10);
2273 
2274     assert(any!((a, b) => a * b == 39)(a, b));
2275 }
2276 
2277 /// Zipped slices
2278 @safe pure nothrow
2279 version(mir_test) unittest
2280 {
2281     import mir.ndslice.topology : iota, zip;
2282 
2283     // 0 1 2
2284     // 3 4 5
2285     auto a = iota(2, 3);
2286     // 10 11 12
2287     // 13 14 15
2288     auto b = iota([2, 3], 10);
2289 
2290     // slices must have the same strides
2291 
2292     assert(zip!true(a, b).any!"a.a * a.b == 39");
2293 }
2294 
2295 /// Mutation on-the-fly
2296 pure nothrow
2297 version(mir_test) unittest
2298 {
2299     import mir.ndslice.allocation : slice;
2300     import mir.ndslice.topology : as, iota;
2301 
2302     // 0 1 2
2303     // 3 4 5
2304     auto sl = iota(2, 3).as!double.slice;
2305 
2306     static bool pred(T)(ref T a)
2307     {
2308         if (a == 5)
2309             return true;
2310         a = 8;
2311         return false;
2312     }
2313 
2314     assert(sl.any!pred);
2315 
2316     // sl was changed
2317     assert(sl == [[8, 8, 8],
2318                   [8, 8, 5]]);
2319 }
2320 
2321 size_t allImpl(alias fun, Slices...)(Slices slices)
2322     if (Slices.length)
2323 {
2324     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(FieldIterator!(BitField!(Field, I))), Field, I))
2325     {
2326         return BitSliceAccelerator!(LightScopeOf!Field, I)(lightScope(slices[0])).all;
2327     }
2328     else
2329     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))), Field, I))
2330     {
2331         // pragma(msg, S);
2332         import mir.ndslice.topology: retro;
2333         return .allImpl!fun(lightScope(slices[0]).retro);
2334     }
2335     else
2336     {
2337         do
2338         {
2339             static if (DimensionCount!(Slices[0]) == 1)
2340             {
2341                 if (!mixin(`fun(`~ frontOf!(Slices.length) ~ `)`))
2342                     return false;
2343             }
2344             else
2345             {
2346                 if (!allImpl!fun(frontOf2!slices))
2347                     return false;
2348             }
2349             foreach_reverse(ref slice; slices)
2350                 slice.popFront;
2351         }
2352         while(!slices[0].empty);
2353         return true;
2354     }
2355 }
2356 
2357 /++
2358 Checks if all of the elements verify `pred`.
2359 
2360 Params:
2361     pred = The predicate.
2362 Optimization:
2363     `all!"a"` has accelerated specialization for slices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
2364 +/
2365 template all(alias pred = "a")
2366 {
2367     import mir.functional: naryFun;
2368     static if (__traits(isSame, naryFun!pred, pred))
2369     /++
2370     Params:
2371         slices = One or more slices.
2372     Returns:
2373         `true` all of the elements verify `pred` and `false` otherwise.
2374     Constraints:
2375         All slices must have the same shape.
2376     +/
2377     @fmamath bool all(Slices...)(Slices slices)
2378         if ((Slices.length == 1 || !__traits(isSame, pred, "a")) && Slices.length)
2379     {
2380         static if (Slices.length > 1)
2381             slices.checkShapesMatch;
2382         static if (areAllContiguousSlices!Slices)
2383         {
2384             import mir.ndslice.topology: flattened;
2385             return .all!pred(allFlattened!(allLightScope!slices));
2386         }
2387         else
2388         {
2389             return slices[0].anyEmpty || allImpl!pred(allLightScope!slices);
2390         }
2391     }
2392     else
2393         alias all = .all!(naryFun!pred);
2394 }
2395 
2396 /// Ranges and arrays
2397 @safe pure nothrow @nogc
2398 version(mir_test) unittest
2399 {
2400     import std.range : iota;
2401     // 0 1 2 3 4 5
2402     auto r = iota(6);
2403 
2404     assert(r.all!"a < 6");
2405     assert(!r.all!"a < 5");
2406 }
2407 
2408 ///
2409 @safe pure nothrow
2410 version(mir_test) unittest
2411 {
2412     import mir.ndslice.topology : iota;
2413 
2414     // 0 1 2
2415     // 3 4 5
2416     auto sl = iota(2, 3);
2417 
2418     assert(sl.all!"a < 6");
2419     assert(!sl.all!"a < 5");
2420 }
2421 
2422 /// Multiple slices
2423 @safe pure nothrow
2424 version(mir_test) unittest
2425 {
2426     import mir.ndslice.topology : iota;
2427 
2428     // 0 1 2
2429     // 3 4 5
2430     auto sl = iota(2, 3);
2431 
2432     assert(all!"a - b == 0"(sl, sl));
2433 }
2434 
2435 /// Zipped slices
2436 @safe pure nothrow
2437 version(mir_test) unittest
2438 {
2439     import mir.ndslice.topology : iota, zip;
2440 
2441     // 0 1 2
2442     // 3 4 5
2443     auto sl = iota(2, 3);
2444 
2445 
2446     assert(zip!true(sl, sl).all!"a.a - a.b == 0");
2447 }
2448 
2449 /// Mutation on-the-fly
2450 pure nothrow
2451 version(mir_test) unittest
2452 {
2453     import mir.ndslice.allocation : slice;
2454     import mir.ndslice.topology : as, iota;
2455 
2456     // 0 1 2
2457     // 3 4 5
2458     auto sl = iota(2, 3).as!double.slice;
2459 
2460     static bool pred(T)(ref T a)
2461     {
2462         if (a < 4)
2463         {
2464             a = 8;
2465             return true;
2466         }
2467         return false;
2468     }
2469 
2470     assert(!sl.all!pred);
2471 
2472     // sl was changed
2473     assert(sl == [[8, 8, 8],
2474                   [8, 4, 5]]);
2475 }
2476 
2477 @safe pure nothrow
2478 version(mir_test) unittest
2479 {
2480     import mir.ndslice.topology : iota;
2481     size_t i;
2482     assert(iota(2, 0).all!((elem){i++; return true;}));
2483     assert(i == 0);
2484 }
2485 
2486 /++
2487 Counts elements in slices according to the `fun`.
2488 Params:
2489     fun = A predicate.
2490 
2491 Optimization:
2492     `count!"a"` has accelerated specialization for slices created with $(REF bitwise, mir,ndslice,topology), $(REF bitSlice, mir,ndslice,allocation).
2493 +/
2494 template count(alias fun)
2495 {
2496     import mir.functional: naryFun;
2497     static if (__traits(isSame, naryFun!fun, fun))
2498     /++
2499     Params:
2500         slices = One or more slices, ranges, and arrays.
2501 
2502     Returns: The number elements according to the `fun`.
2503 
2504     Constraints:
2505         All slices must have the same shape.
2506     +/
2507     @fmamath size_t count(Slices...)(Slices slices)
2508         if (Slices.length)
2509     {
2510         static if (Slices.length > 1)
2511             slices.checkShapesMatch;
2512         static if (__traits(isSame, fun, naryFun!"true"))
2513         {
2514             return slices[0].elementCount;
2515         }
2516         else
2517         static if (areAllContiguousSlices!Slices)
2518         {
2519             import mir.ndslice.topology: flattened;
2520             return .count!fun(allFlattened!(allLightScope!slices));
2521         }
2522         else
2523         {
2524             if (slices[0].anyEmpty)
2525                 return 0;
2526             return countImpl!(fun)(allLightScope!slices);
2527         }
2528     }
2529     else
2530         alias count = .count!(naryFun!fun);
2531 
2532 }
2533 
2534 /// Ranges and arrays
2535 @safe pure nothrow @nogc
2536 version(mir_test) unittest
2537 {
2538     import std.range : iota;
2539     // 0 1 2 3 4 5
2540     auto r = iota(6);
2541 
2542     assert(r.count!"true" == 6);
2543     assert(r.count!"a" == 5);
2544     assert(r.count!"a % 2" == 3);
2545 }
2546 
2547 /// Single slice
2548 version(mir_test)
2549 unittest
2550 {
2551     import mir.ndslice.topology : iota;
2552 
2553     //| 0 1 2 |
2554     //| 3 4 5 |
2555     auto sl = iota(2, 3);
2556 
2557     assert(sl.count!"true" == 6);
2558     assert(sl.count!"a" == 5);
2559     assert(sl.count!"a % 2" == 3);
2560 }
2561 
2562 /// Accelerated set bit count
2563 version(mir_test)
2564 unittest
2565 {
2566     import mir.ndslice.topology: retro, iota, bitwise;
2567     import mir.ndslice.allocation: slice;
2568 
2569     //| 0 1 2 |
2570     //| 3 4 5 |
2571     auto sl = iota!size_t(2, 3).bitwise;
2572 
2573     assert(sl.count!"true" == 6 * size_t.sizeof * 8);
2574 
2575     assert(sl.slice.count!"a" == 7);
2576 
2577     // accelerated
2578     assert(sl.count!"a" == 7);
2579     assert(sl.retro.count!"a" == 7);
2580 
2581     auto sl2 = iota!ubyte([6], 128).bitwise;
2582     // accelerated
2583     assert(sl2.count!"a" == 13);
2584     assert(sl2[4 .. $].count!"a" == 13);
2585     assert(sl2[4 .. $ - 1].count!"a" == 12);
2586     assert(sl2[4 .. $ - 1].count!"a" == 12);
2587     assert(sl2[41 .. $ - 1].count!"a" == 1);
2588 }
2589 
2590 version(mir_test)
2591 unittest
2592 {
2593     import mir.ndslice.allocation: slice;
2594     import mir.ndslice.topology: bitwise, assumeFieldsHaveZeroShift;
2595     auto sl = slice!uint([6]).bitwise;
2596     auto slb = slice!ubyte([6]).bitwise;
2597     slb[4] = true;
2598     auto d = slb[4];
2599     auto c = assumeFieldsHaveZeroShift(slb & ~slb);
2600     // pragma(msg, typeof(c));
2601     assert(!sl.any);
2602     assert((~sl).all);
2603     // pragma(msg, typeof(~slb));
2604     // pragma(msg, typeof(~slb));
2605     // assert(sl.findIndex);
2606 }
2607 
2608 /++
2609 Compares two or more slices for equality, as defined by predicate `pred`.
2610 
2611 See_also: $(NDSLICEREF slice, Slice.opEquals)
2612 
2613 Params:
2614     pred = The predicate.
2615 +/
2616 template equal(alias pred = "a == b")
2617 {
2618     import mir.functional: naryFun;
2619     static if (__traits(isSame, naryFun!pred, pred))
2620     {
2621         /++
2622         Params:
2623             slices = Two or more ndslices, ranges, and arrays.
2624 
2625         Returns:
2626             `true` any of the elements verify `pred` and `false` otherwise.
2627         +/
2628         bool equal(Slices...)(Slices slices) @safe
2629             if (Slices.length >= 2)
2630         {
2631             import mir.internal.utility;
2632             static if (allSatisfy!(hasShape, Slices))
2633             {
2634                 auto shape0 = slices[0].shape;
2635                 enum N = DimensionCount!(Slices[0]);
2636                 foreach (ref slice; slices[1 .. $])
2637                 {
2638                     if (slice.shape != shape0)
2639                         goto False;
2640                 }
2641                 return all!pred(allLightScope!slices);
2642             }
2643             else
2644             {
2645                 for(;;)
2646                 {
2647                     auto empty = slices[0].empty;
2648                     foreach (ref slice; slices[1 .. $])
2649                     {
2650                         if (slice.empty != empty)
2651                             goto False;
2652                     }
2653                     if (empty)
2654                         return true;
2655                     if (!mixin(`pred(`~ frontOf!(Slices.length) ~ `)`))
2656                         goto False;
2657                     foreach (ref slice; slices)
2658                         slice.popFront;
2659                 }
2660             }
2661             False: return false;
2662         }
2663     }
2664     else
2665         alias equal = .equal!(naryFun!pred);
2666 }
2667 
2668 /// Ranges and arrays
2669 @safe pure nothrow
2670 version(mir_test) unittest
2671 {
2672     import std.range : iota;
2673     auto r = iota(6);
2674     assert(r.equal([0, 1, 2, 3, 4, 5]));
2675 }
2676 
2677 ///
2678 @safe pure nothrow @nogc
2679 version(mir_test) unittest
2680 {
2681     import mir.ndslice.allocation : slice;
2682     import mir.ndslice.topology : iota;
2683 
2684     // 0 1 2
2685     // 3 4 5
2686     auto sl1 = iota(2, 3);
2687     // 1 2 3
2688     // 4 5 6
2689     auto sl2 = iota([2, 3], 1);
2690 
2691     assert(equal(sl1, sl1));
2692     assert(sl1 == sl1); //can also use opEquals for two Slices
2693     assert(equal!"2 * a == b + c"(sl1, sl1, sl1));
2694 
2695     assert(equal!"a < b"(sl1, sl2));
2696 
2697     assert(!equal(sl1[0 .. $ - 1], sl1));
2698     assert(!equal(sl1[0 .. $, 0 .. $ - 1], sl1));
2699 }
2700 
2701 @safe pure nothrow @nogc
2702 version(mir_test) unittest
2703 {
2704     import mir.math.common: approxEqual;
2705     import mir.ndslice.allocation: rcslice;
2706     import mir.ndslice.topology: as, iota;
2707 
2708     auto x = 5.iota.as!double.rcslice;
2709     auto y = x.rcslice;
2710 
2711     assert(equal(x, y));
2712     assert(equal!approxEqual(x, y));
2713 }
2714 
2715 ptrdiff_t cmpImpl(alias pred, A, B)
2716     (scope A sl1, scope B sl2)
2717     if (DimensionCount!A == DimensionCount!B)
2718 {
2719     for (;;)
2720     {
2721         static if (DimensionCount!A == 1)
2722         {
2723             import mir.functional : naryFun;
2724             if (naryFun!pred(sl1.front, sl2.front))
2725                 return -1;
2726             if (naryFun!pred(sl2.front, sl1.front))
2727                 return 1;
2728         }
2729         else
2730         {
2731             if (auto res = .cmpImpl!pred(sl1.front, sl2.front))
2732                 return res;
2733         }
2734         sl1.popFront;
2735         if (sl1.empty)
2736             return -cast(ptrdiff_t)(sl2.length > 1);
2737         sl2.popFront;
2738         if (sl2.empty)
2739             return 1;
2740     }
2741 }
2742 
2743 /++
2744 Performs three-way recursive lexicographical comparison on two slices according to predicate `pred`.
2745 Iterating `sl1` and `sl2` in lockstep, `cmp` compares each `N-1` dimensional element `e1` of `sl1`
2746 with the corresponding element `e2` in `sl2` recursively.
2747 If one of the slices has been finished,`cmp` returns a negative value if `sl1` has fewer elements than `sl2`,
2748 a positive value if `sl1` has more elements than `sl2`,
2749 and `0` if the ranges have the same number of elements.
2750 
2751 Params:
2752     pred = The predicate.
2753 +/
2754 template cmp(alias pred = "a < b")
2755 {
2756     import mir.functional: naryFun;
2757     static if (__traits(isSame, naryFun!pred, pred))
2758     /++
2759     Params:
2760         sl1 = First slice, range, or array.
2761         sl2 = Second slice, range, or array.
2762 
2763     Returns:
2764         `0` if both ranges compare equal.
2765         Negative value if the first differing element of `sl1` is less than the corresponding
2766         element of `sl2` according to `pred`.
2767         Positive value if the first differing element of `sl2` is less than the corresponding
2768         element of `sl1` according to `pred`.
2769     +/
2770     auto cmp(A, B)
2771         (scope A sl1, scope B sl2)
2772         if (DimensionCount!A == DimensionCount!B)
2773     {
2774         auto b = sl2.anyEmpty;
2775         if (sl1.anyEmpty)
2776         {
2777             if (!b)
2778                 return -1;
2779             auto sh1 = sl1.shape;
2780             auto sh2 = sl2.shape;
2781             foreach (i; Iota!(DimensionCount!A))
2782                 if (sh1[i] != sh2[i])
2783                     return sh1[i] > sh2[i] ? 1 : -1;
2784             return 0;
2785         }
2786         if (b)
2787             return 1;
2788         return cmpImpl!pred(lightScope(sl1), lightScope(sl2));
2789     }
2790     else
2791         alias cmp = .cmp!(naryFun!pred);
2792 }
2793 
2794 /// Ranges and arrays
2795 @safe pure nothrow
2796 version(mir_test) unittest
2797 {
2798     import std.range : iota;
2799 
2800     // 0 1 2 3 4 5
2801     auto r1 = iota(0, 6);
2802     // 1 2 3 4 5 6
2803     auto r2 = iota(1, 7);
2804 
2805     assert(cmp(r1, r1) == 0);
2806     assert(cmp(r1, r2) < 0);
2807     assert(cmp!"a >= b"(r1, r2) > 0);
2808 }
2809 
2810 ///
2811 @safe pure nothrow @nogc
2812 version(mir_test) unittest
2813 {
2814     import mir.ndslice.topology : iota;
2815 
2816     // 0 1 2
2817     // 3 4 5
2818     auto sl1 = iota(2, 3);
2819     // 1 2 3
2820     // 4 5 6
2821     auto sl2 = iota([2, 3], 1);
2822 
2823     assert(cmp(sl1, sl1) == 0);
2824     assert(cmp(sl1, sl2) < 0);
2825     assert(cmp!"a >= b"(sl1, sl2) > 0);
2826 }
2827 
2828 @safe pure nothrow @nogc
2829 version(mir_test) unittest
2830 {
2831     import mir.ndslice.topology : iota;
2832 
2833     auto sl1 = iota(2, 3);
2834     auto sl2 = iota([2, 3], 1);
2835 
2836     assert(cmp(sl1[0 .. $ - 1], sl1) < 0);
2837     assert(cmp(sl1, sl1[0 .. $, 0 .. $ - 1]) > 0);
2838 
2839     assert(cmp(sl1[0 .. $ - 2], sl1) < 0);
2840     assert(cmp(sl1, sl1[0 .. $, 0 .. $ - 3]) > 0);
2841     assert(cmp(sl1[0 .. $, 0 .. $ - 3], sl1[0 .. $, 0 .. $ - 3]) == 0);
2842     assert(cmp(sl1[0 .. $, 0 .. $ - 3], sl1[0 .. $ - 1, 0 .. $ - 3]) > 0);
2843     assert(cmp(sl1[0 .. $ - 1, 0 .. $ - 3], sl1[0 .. $, 0 .. $ - 3]) < 0);
2844 }
2845 
2846 size_t countImpl(alias fun, Slices...)(Slices slices)
2847 {
2848     size_t ret;
2849     alias S = Slices[0];
2850     import mir.functional: naryFun;
2851     import mir.ndslice.iterator: FieldIterator, RetroIterator;
2852     import mir.ndslice.field: BitField;
2853     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(FieldIterator!(BitField!(Field, I))), Field, I))
2854     {
2855         ret = BitSliceAccelerator!(Field, I)(slices[0]).ctpop;
2856     }
2857     else
2858     static if (__traits(isSame, fun, naryFun!"a") && is(S : Slice!(RetroIterator!(FieldIterator!(BitField!(Field, I)))), Field, I))
2859     {
2860         // pragma(msg, S);
2861         import mir.ndslice.topology: retro;
2862         ret = .countImpl!fun(lightScope(slices[0]).retro);
2863     }
2864     else
2865     do
2866     {
2867         static if (DimensionCount!(Slices[0]) == 1)
2868         {
2869             if(mixin(`fun(`~ frontOf!(Slices.length) ~ `)`))
2870                 ret++;
2871         }
2872         else
2873             ret += .countImpl!fun(frontOf2!slices);
2874         foreach_reverse(ref slice; slices)
2875             slice.popFront;
2876     }
2877     while(!slices[0].empty);
2878     return ret;
2879 }
2880 
2881 /++
2882 Returns: max length across all dimensions.
2883 +/
2884 size_t maxLength(S)(auto ref scope S s)
2885  if (hasShape!S)
2886 {
2887     auto shape = s.shape;
2888     size_t length = 0;
2889     foreach(i; Iota!(shape.length))
2890         if (shape[i] > length)
2891             length = shape[i];
2892     return length;
2893 }
2894 
2895 /++
2896 The call `eachLower!(fun)(slice1, ..., sliceN)` evaluates `fun` on the lower
2897 triangle in `slice1, ..., sliceN` respectively.
2898 
2899 `eachLower` allows iterating multiple slices in the lockstep.
2900 
2901 Params:
2902     fun = A function
2903 See_Also:
2904     This is functionally similar to $(LREF each).
2905 +/
2906 template eachLower(alias fun)
2907 {
2908     import mir.functional : naryFun;
2909 
2910     static if (__traits(isSame, naryFun!fun, fun))
2911     {
2912         /++
2913         Params:
2914             inputs = One or more two-dimensional slices and an optional
2915                      integer, `k`.
2916 
2917             The value `k` determines which diagonals will have the function
2918             applied:
2919             For k = 0, the function is also applied to the main diagonal
2920             For k = 1 (default), only the non-main diagonals below the main
2921             diagonal will have the function applied.
2922             For k > 1, fewer diagonals below the main diagonal will have the
2923             function applied.
2924             For k < 0, more diagonals above the main diagonal will have the
2925             function applied.
2926         +/
2927         void eachLower(Inputs...)(scope Inputs inputs)
2928             if (((Inputs.length > 1) &&
2929                  (isIntegral!(Inputs[$ - 1]))) ||
2930                 (Inputs.length))
2931         {
2932             import mir.ndslice.traits : isMatrix;
2933 
2934             size_t val;
2935 
2936             static if ((Inputs.length > 1) && (isIntegral!(Inputs[$ - 1])))
2937             {
2938                 immutable(sizediff_t) k = inputs[$ - 1];
2939                 alias Slices = Inputs[0..($ - 1)];
2940                 alias slices = inputs[0..($ - 1)];
2941             }
2942             else
2943             {
2944                 enum sizediff_t k = 1;
2945                 alias Slices = Inputs;
2946                 alias slices = inputs;
2947             }
2948 
2949             static assert (allSatisfy!(isMatrix, Slices),
2950                 "eachLower: Every slice input must be a two-dimensional slice");
2951             static if (Slices.length > 1)
2952                 slices.checkShapesMatch;
2953             if (slices[0].anyEmpty)
2954                 return;
2955 
2956             foreach(ref slice; slices)
2957                 assert(!slice.empty);
2958 
2959             immutable(size_t) m = slices[0].length!0;
2960             immutable(size_t) n = slices[0].length!1;
2961 
2962             if ((n + k) < m)
2963             {
2964                 val = m - (n + k);
2965                 .eachImpl!fun(selectBackOf!(val, slices));
2966             }
2967 
2968             size_t i;
2969 
2970             if (k > 0)
2971             {
2972                 foreach(ref slice; slices)
2973                     slice.popFrontExactly!0(k);
2974                 i = k;
2975             }
2976 
2977             do
2978             {
2979                 val = i - k + 1;
2980                 .eachImpl!fun(frontSelectFrontOf!(val, slices));
2981 
2982                 foreach(ref slice; slices)
2983                         slice.popFront!0;
2984                 i++;
2985             } while ((i < (n + k)) && (i < m));
2986         }
2987     }
2988     else
2989     {
2990         alias eachLower = .eachLower!(naryFun!fun);
2991     }
2992 }
2993 
2994 pure nothrow
2995 version(mir_test) unittest
2996 {
2997     import mir.ndslice.allocation: slice;
2998     import mir.ndslice.topology: iota, canonical, universal;
2999     alias AliasSeq(T...) = T;
3000 
3001     pure nothrow
3002     void test(alias func)()
3003     {
3004         //| 1 2 3 |
3005         //| 4 5 6 |
3006         //| 7 8 9 |
3007         auto m = func(iota([3, 3], 1).slice);
3008         m.eachLower!"a = 0"(0);
3009         assert(m == [
3010             [0, 2, 3],
3011             [0, 0, 6],
3012             [0, 0, 0]]);
3013     }
3014 
3015     @safe pure nothrow @nogc
3016     T identity(T)(T x)
3017     {
3018         return x;
3019     }
3020 
3021     alias kinds = AliasSeq!(identity, canonical, universal);
3022     test!(kinds[0]);
3023     test!(kinds[1]);
3024     test!(kinds[2]);
3025 }
3026 
3027 ///
3028 pure nothrow
3029 version(mir_test) unittest
3030 {
3031     import mir.ndslice.allocation: slice;
3032     import mir.ndslice.topology: iota;
3033 
3034     //| 1 2 3 |
3035     //| 4 5 6 |
3036     //| 7 8 9 |
3037     auto m = iota([3, 3], 1).slice;
3038     m.eachLower!"a = 0";
3039     assert(m == [
3040         [1, 2, 3],
3041         [0, 5, 6],
3042         [0, 0, 9]]);
3043 }
3044 
3045 pure nothrow
3046 version(mir_test) unittest
3047 {
3048     import mir.ndslice.allocation: slice;
3049     import mir.ndslice.topology: iota;
3050 
3051     //| 1 2 3 |
3052     //| 4 5 6 |
3053     //| 7 8 9 |
3054     auto m = iota([3, 3], 1).slice;
3055     m.eachLower!"a = 0"(-1);
3056     assert(m == [
3057         [0, 0, 3],
3058         [0, 0, 0],
3059         [0, 0, 0]]);
3060 }
3061 
3062 pure nothrow
3063 version(mir_test) unittest
3064 {
3065     import mir.ndslice.allocation: slice;
3066     import mir.ndslice.topology: iota;
3067 
3068     //| 1 2 3 |
3069     //| 4 5 6 |
3070     //| 7 8 9 |
3071     auto m = iota([3, 3], 1).slice;
3072     m.eachLower!"a = 0"(2);
3073     assert(m == [
3074         [1, 2, 3],
3075         [4, 5, 6],
3076         [0, 8, 9]]);
3077 }
3078 
3079 pure nothrow
3080 version(mir_test) unittest
3081 {
3082     import mir.ndslice.allocation: slice;
3083     import mir.ndslice.topology: iota;
3084 
3085     //| 1 2 3 |
3086     //| 4 5 6 |
3087     //| 7 8 9 |
3088     auto m = iota([3, 3], 1).slice;
3089     m.eachLower!"a = 0"(-2);
3090     assert(m == [
3091         [0, 0, 0],
3092         [0, 0, 0],
3093         [0, 0, 0]]);
3094 }
3095 
3096 pure nothrow
3097 version(mir_test) unittest
3098 {
3099     import mir.ndslice.allocation: slice;
3100     import mir.ndslice.topology: iota;
3101 
3102     //| 1  2  3  4 |
3103     //| 5  6  7  8 |
3104     //| 9 10 11 12 |
3105     auto m = iota([3, 4], 1).slice;
3106     m.eachLower!"a = 0"(0);
3107     assert(m == [
3108         [0, 2, 3, 4],
3109         [0, 0, 7, 8],
3110         [0, 0, 0, 12]]);
3111 }
3112 
3113 pure nothrow
3114 version(mir_test) unittest
3115 {
3116     import mir.ndslice.allocation: slice;
3117     import mir.ndslice.topology: iota;
3118 
3119     //| 1  2  3  4 |
3120     //| 5  6  7  8 |
3121     //| 9 10 11 12 |
3122     auto m = iota([3, 4], 1).slice;
3123     m.eachLower!"a = 0";
3124     assert(m == [
3125         [1, 2, 3, 4],
3126         [0, 6, 7, 8],
3127         [0, 0, 11, 12]]);
3128 }
3129 
3130 pure nothrow
3131 version(mir_test) unittest
3132 {
3133     import mir.ndslice.allocation: slice;
3134     import mir.ndslice.topology: iota;
3135 
3136     //| 1  2  3  4 |
3137     //| 5  6  7  8 |
3138     //| 9 10 11 12 |
3139     auto m = iota([3, 4], 1).slice;
3140     m.eachLower!"a = 0"(-1);
3141     assert(m == [
3142         [0, 0, 3, 4],
3143         [0, 0, 0, 8],
3144         [0, 0, 0, 0]]);
3145 }
3146 
3147 pure nothrow
3148 version(mir_test) unittest
3149 {
3150     import mir.ndslice.allocation: slice;
3151     import mir.ndslice.topology: iota;
3152 
3153     //| 1  2  3  4 |
3154     //| 5  6  7  8 |
3155     //| 9 10 11 12 |
3156     auto m = iota([3, 4], 1).slice;
3157     m.eachLower!"a = 0"(2);
3158     assert(m == [
3159         [1, 2, 3, 4],
3160         [5, 6, 7, 8],
3161         [0, 10, 11, 12]]);
3162 }
3163 
3164 pure nothrow
3165 version(mir_test) unittest
3166 {
3167     import mir.ndslice.allocation: slice;
3168     import mir.ndslice.topology: iota;
3169 
3170     //| 1  2  3  4 |
3171     //| 5  6  7  8 |
3172     //| 9 10 11 12 |
3173     auto m = iota([3, 4], 1).slice;
3174     m.eachLower!"a = 0"(-2);
3175     assert(m == [
3176         [0, 0, 0, 4],
3177         [0, 0, 0, 0],
3178         [0, 0, 0, 0]]);
3179 }
3180 
3181 pure nothrow
3182 version(mir_test) unittest
3183 {
3184     import mir.ndslice.allocation: slice;
3185     import mir.ndslice.topology: iota;
3186 
3187     //|  1  2  3 |
3188     //|  4  5  6 |
3189     //|  7  8  9 |
3190     //| 10 11 12 |
3191     auto m = iota([4, 3], 1).slice;
3192     m.eachLower!"a = 0"(0);
3193     assert(m == [
3194         [0, 2, 3],
3195         [0, 0, 6],
3196         [0, 0, 0],
3197         [0, 0, 0]]);
3198 }
3199 
3200 pure nothrow
3201 version(mir_test) unittest
3202 {
3203     import mir.ndslice.allocation: slice;
3204     import mir.ndslice.topology: iota;
3205 
3206     //|  1  2  3 |
3207     //|  4  5  6 |
3208     //|  7  8  9 |
3209     //| 10 11 12 |
3210     auto m = iota([4, 3], 1).slice;
3211     m.eachLower!"a = 0";
3212     assert(m == [
3213         [1, 2, 3],
3214         [0, 5, 6],
3215         [0, 0, 9],
3216         [0, 0, 0]]);
3217 }
3218 
3219 pure nothrow
3220 version(mir_test) unittest
3221 {
3222     import mir.ndslice.allocation: slice;
3223     import mir.ndslice.topology: iota;
3224 
3225     //|  1  2  3 |
3226     //|  4  5  6 |
3227     //|  7  8  9 |
3228     //| 10 11 12 |
3229     auto m = iota([4, 3], 1).slice;
3230     m.eachLower!"a = 0"(-1);
3231     assert(m == [
3232         [0, 0, 3],
3233         [0, 0, 0],
3234         [0, 0, 0],
3235         [0, 0, 0]]);
3236 }
3237 
3238 pure nothrow
3239 version(mir_test) unittest
3240 {
3241     import mir.ndslice.allocation: slice;
3242     import mir.ndslice.topology: iota;
3243 
3244     //|  1  2  3 |
3245     //|  4  5  6 |
3246     //|  7  8  9 |
3247     //| 10 11 12 |
3248     auto m = iota([4, 3], 1).slice;
3249     m.eachLower!"a = 0"(2);
3250     assert(m == [
3251         [1, 2, 3],
3252         [4, 5, 6],
3253         [0, 8, 9],
3254         [0, 0, 12]]);
3255 }
3256 
3257 pure nothrow
3258 version(mir_test) unittest
3259 {
3260     import mir.ndslice.allocation: slice;
3261     import mir.ndslice.topology: iota;
3262 
3263     //|  1  2  3 |
3264     //|  4  5  6 |
3265     //|  7  8  9 |
3266     //| 10 11 12 |
3267     auto m = iota([4, 3], 1).slice;
3268     m.eachLower!"a = 0"(-2);
3269     assert(m == [
3270         [0, 0, 0],
3271         [0, 0, 0],
3272         [0, 0, 0],
3273         [0, 0, 0]]);
3274 }
3275 
3276 /// Swap two slices
3277 pure nothrow
3278 version(mir_test) unittest
3279 {
3280     import mir.utility : swap;
3281     import mir.ndslice.allocation : slice;
3282     import mir.ndslice.topology : as, iota;
3283 
3284     //| 0 1 2 |
3285     //| 3 4 5 |
3286     //| 6 7 8 |
3287     auto a = iota([3, 3]).as!double.slice;
3288     //| 10 11 12 |
3289     //| 13 14 15 |
3290     //| 16 17 18 |
3291     auto b = iota([3, 3], 10).as!double.slice;
3292 
3293     eachLower!swap(a, b);
3294 
3295     assert(a == [
3296         [ 0,  1, 2],
3297         [13,  4, 5],
3298         [16, 17, 8]]);
3299     assert(b == [
3300         [10, 11, 12],
3301         [ 3, 14, 15],
3302         [ 6,  7, 18]]);
3303 }
3304 
3305 /// Swap two zipped slices
3306 pure nothrow
3307 version(mir_test) unittest
3308 {
3309     import mir.utility : swap;
3310     import mir.ndslice.allocation : slice;
3311     import mir.ndslice.topology : as, zip, iota;
3312 
3313     //| 0 1 2 |
3314     //| 3 4 5 |
3315     //| 6 7 8 |
3316     auto a = iota([3, 3]).as!double.slice;
3317     //| 10 11 12 |
3318     //| 13 14 15 |
3319     //| 16 17 18 |
3320     auto b = iota([3, 3], 10).as!double.slice;
3321 
3322     auto z = zip(a, b);
3323 
3324     z.eachLower!(z => swap(z.a, z.b));
3325 
3326     assert(a == [
3327         [ 0,  1, 2],
3328         [13,  4, 5],
3329         [16, 17, 8]]);
3330     assert(b == [
3331         [10, 11, 12],
3332         [ 3, 14, 15],
3333         [ 6,  7, 18]]);
3334 }
3335 
3336 /++
3337 The call `eachUpper!(fun)(slice1, ..., sliceN)` evaluates `fun` on the upper
3338 triangle in `slice1, ..., sliceN`, respectively.
3339 
3340 `eachUpper` allows iterating multiple slices in the lockstep.
3341 
3342 Params:
3343     fun = A function
3344 See_Also:
3345     This is functionally similar to $(LREF each).
3346 +/
3347 template eachUpper(alias fun)
3348 {
3349     import mir.functional: naryFun;
3350 
3351     static if (__traits(isSame, naryFun!fun, fun))
3352     {
3353         /++
3354         Params:
3355             inputs = One or more two-dimensional slices and an optional
3356                      integer, `k`.
3357 
3358             The value `k` determines which diagonals will have the function
3359             applied:
3360             For k = 0, the function is also applied to the main diagonal
3361             For k = 1 (default), only the non-main diagonals above the main
3362             diagonal will have the function applied.
3363             For k > 1, fewer diagonals below the main diagonal will have the
3364             function applied.
3365             For k < 0, more diagonals above the main diagonal will have the
3366             function applied.
3367         +/
3368         void eachUpper(Inputs...)(scope Inputs inputs)
3369             if (((Inputs.length > 1) &&
3370                  (isIntegral!(Inputs[$ - 1]))) ||
3371                 (Inputs.length))
3372         {
3373             import mir.ndslice.traits : isMatrix;
3374 
3375             size_t val;
3376 
3377             static if ((Inputs.length > 1) && (isIntegral!(Inputs[$ - 1])))
3378             {
3379                 immutable(sizediff_t) k = inputs[$ - 1];
3380                 alias Slices = Inputs[0..($ - 1)];
3381                 alias slices = inputs[0..($ - 1)];
3382             }
3383             else
3384             {
3385                 enum sizediff_t k = 1;
3386                 alias Slices = Inputs;
3387                 alias slices = inputs;
3388             }
3389 
3390             static assert (allSatisfy!(isMatrix, Slices),
3391                 "eachUpper: Every slice input must be a two-dimensional slice");
3392             static if (Slices.length > 1)
3393                 slices.checkShapesMatch;
3394             if (slices[0].anyEmpty)
3395                 return;
3396 
3397             foreach(ref slice; slices)
3398                 assert(!slice.empty);
3399 
3400             immutable(size_t) m = slices[0].length!0;
3401             immutable(size_t) n = slices[0].length!1;
3402 
3403             size_t i;
3404 
3405             if (k < 0)
3406             {
3407                 val = -k;
3408                 .eachImpl!fun(selectFrontOf!(val, slices));
3409 
3410                 foreach(ref slice; slices)
3411                     slice.popFrontExactly!0(-k);
3412                 i = -k;
3413             }
3414 
3415             do
3416             {
3417                 val = (n - k) - i;
3418                 .eachImpl!fun(frontSelectBackOf!(val, slices));
3419 
3420                 foreach(ref slice; slices)
3421                     slice.popFront;
3422                 i++;
3423             } while ((i < (n - k)) && (i < m));
3424         }
3425     }
3426     else
3427     {
3428         alias eachUpper = .eachUpper!(naryFun!fun);
3429     }
3430 }
3431 
3432 pure nothrow
3433 version(mir_test) unittest
3434 {
3435     import mir.ndslice.allocation: slice;
3436     import mir.ndslice.topology: iota, canonical, universal;
3437 
3438     pure nothrow
3439     void test(alias func)()
3440     {
3441         //| 1 2 3 |
3442         //| 4 5 6 |
3443         //| 7 8 9 |
3444         auto m = func(iota([3, 3], 1).slice);
3445         m.eachUpper!"a = 0"(0);
3446         assert(m == [
3447             [0, 0, 0],
3448             [4, 0, 0],
3449             [7, 8, 0]]);
3450     }
3451 
3452     @safe pure nothrow @nogc
3453     T identity(T)(T x)
3454     {
3455         return x;
3456     }
3457 
3458     alias kinds = AliasSeq!(identity, canonical, universal);
3459     test!(kinds[0]);
3460     test!(kinds[1]);
3461     test!(kinds[2]);
3462 }
3463 
3464 ///
3465 pure nothrow
3466 version(mir_test) unittest
3467 {
3468     import mir.ndslice.allocation: slice;
3469     import mir.ndslice.topology: iota;
3470 
3471     //| 1 2 3 |
3472     //| 4 5 6 |
3473     //| 7 8 9 |
3474     auto m = iota([3, 3], 1).slice;
3475     m.eachUpper!"a = 0";
3476     assert(m == [
3477         [1, 0, 0],
3478         [4, 5, 0],
3479         [7, 8, 9]]);
3480 }
3481 
3482 pure nothrow
3483 version(mir_test) unittest
3484 {
3485     import mir.ndslice.allocation: slice;
3486     import mir.ndslice.topology: iota;
3487 
3488     //| 1 2 3 |
3489     //| 4 5 6 |
3490     //| 7 8 9 |
3491     auto m = iota([3, 3], 1).slice;
3492     m.eachUpper!"a = 0"(-1);
3493     assert(m == [
3494         [0, 0, 0],
3495         [0, 0, 0],
3496         [7, 0, 0]]);
3497 }
3498 
3499 pure nothrow
3500 version(mir_test) unittest
3501 {
3502     import mir.ndslice.allocation: slice;
3503     import mir.ndslice.topology: iota;
3504 
3505     //| 1 2 3 |
3506     //| 4 5 6 |
3507     //| 7 8 9 |
3508     auto m = iota([3, 3], 1).slice;
3509     m.eachUpper!"a = 0"(2);
3510     assert(m == [
3511         [1, 2, 0],
3512         [4, 5, 6],
3513         [7, 8, 9]]);
3514 }
3515 
3516 pure nothrow
3517 version(mir_test) unittest
3518 {
3519     import mir.ndslice.allocation: slice;
3520     import mir.ndslice.topology: iota;
3521 
3522     //| 1 2 3 |
3523     //| 4 5 6 |
3524     //| 7 8 9 |
3525     auto m = iota([3, 3], 1).slice;
3526     m.eachUpper!"a = 0"(-2);
3527     assert(m == [
3528         [0, 0, 0],
3529         [0, 0, 0],
3530         [0, 0, 0]]);
3531 }
3532 
3533 pure nothrow
3534 version(mir_test) unittest
3535 {
3536     import mir.ndslice.allocation: slice;
3537     import mir.ndslice.topology: iota;
3538 
3539     //| 1  2  3  4 |
3540     //| 5  6  7  8 |
3541     //| 9 10 11 12 |
3542     auto m = iota([3, 4], 1).slice;
3543     m.eachUpper!"a = 0"(0);
3544     assert(m == [
3545         [0, 0, 0, 0],
3546         [5, 0, 0, 0],
3547         [9, 10, 0, 0]]);
3548 }
3549 
3550 pure nothrow
3551 version(mir_test) unittest
3552 {
3553     import mir.ndslice.allocation: slice;
3554     import mir.ndslice.topology: iota;
3555 
3556     //| 1  2  3  4 |
3557     //| 5  6  7  8 |
3558     //| 9 10 11 12 |
3559     auto m = iota([3, 4], 1).slice;
3560     m.eachUpper!"a = 0";
3561     assert(m == [
3562         [1, 0, 0, 0],
3563         [5, 6, 0, 0],
3564         [9, 10, 11, 0]]);
3565 }
3566 
3567 pure nothrow
3568 version(mir_test) unittest
3569 {
3570     import mir.ndslice.allocation: slice;
3571     import mir.ndslice.topology: iota;
3572 
3573     //| 1  2  3  4 |
3574     //| 5  6  7  8 |
3575     //| 9 10 11 12 |
3576     auto m = iota([3, 4], 1).slice;
3577     m.eachUpper!"a = 0"(-1);
3578     assert(m == [
3579         [0, 0, 0, 0],
3580         [0, 0, 0, 0],
3581         [9, 0, 0, 0]]);
3582 }
3583 
3584 pure nothrow
3585 version(mir_test) unittest
3586 {
3587     import mir.ndslice.allocation: slice;
3588     import mir.ndslice.topology: iota;
3589 
3590     //| 1  2  3  4 |
3591     //| 5  6  7  8 |
3592     //| 9 10 11 12 |
3593     auto m = iota([3, 4], 1).slice;
3594     m.eachUpper!"a = 0"(2);
3595     assert(m == [
3596         [1, 2, 0, 0],
3597         [5, 6, 7, 0],
3598         [9, 10, 11, 12]]);
3599 }
3600 
3601 pure nothrow
3602 version(mir_test) unittest
3603 {
3604     import mir.ndslice.allocation: slice;
3605     import mir.ndslice.topology: iota;
3606 
3607     //| 1  2  3  4 |
3608     //| 5  6  7  8 |
3609     //| 9 10 11 12 |
3610     auto m = iota([3, 4], 1).slice;
3611     m.eachUpper!"a = 0"(-2);
3612     assert(m == [
3613         [0, 0, 0, 0],
3614         [0, 0, 0, 0],
3615         [0, 0, 0, 0]]);
3616 }
3617 
3618 pure nothrow
3619 version(mir_test) unittest
3620 {
3621     import mir.ndslice.allocation: slice;
3622     import mir.ndslice.topology: iota;
3623 
3624     //|  1  2  3 |
3625     //|  4  5  6 |
3626     //|  7  8  9 |
3627     //| 10 11 12 |
3628     auto m = iota([4, 3], 1).slice;
3629     m.eachUpper!"a = 0"(0);
3630     assert(m == [
3631         [0, 0, 0],
3632         [4, 0, 0],
3633         [7, 8, 0],
3634         [10, 11, 12]]);
3635 }
3636 
3637 pure nothrow
3638 version(mir_test) unittest
3639 {
3640     import mir.ndslice.allocation: slice;
3641     import mir.ndslice.topology: iota;
3642 
3643     //|  1  2  3 |
3644     //|  4  5  6 |
3645     //|  7  8  9 |
3646     //| 10 11 12 |
3647     auto m = iota([4, 3], 1).slice;
3648     m.eachUpper!"a = 0";
3649     assert(m == [
3650         [1, 0, 0],
3651         [4, 5, 0],
3652         [7, 8, 9],
3653         [10, 11, 12]]);
3654 }
3655 
3656 pure nothrow
3657 version(mir_test) unittest
3658 {
3659     import mir.ndslice.allocation: slice;
3660     import mir.ndslice.topology: iota;
3661 
3662     //|  1  2  3 |
3663     //|  4  5  6 |
3664     //|  7  8  9 |
3665     //| 10 11 12 |
3666     auto m = iota([4, 3], 1).slice;
3667     m.eachUpper!"a = 0"(-1);
3668     assert(m == [
3669         [0, 0, 0],
3670         [0, 0, 0],
3671         [7, 0, 0],
3672         [10, 11, 0]]);
3673 }
3674 
3675 pure nothrow
3676 version(mir_test) unittest
3677 {
3678     import mir.ndslice.allocation: slice;
3679     import mir.ndslice.topology: iota;
3680 
3681     //|  1  2  3 |
3682     //|  4  5  6 |
3683     //|  7  8  9 |
3684     //| 10 11 12 |
3685     auto m = iota([4, 3], 1).slice;
3686     m.eachUpper!"a = 0"(2);
3687     assert(m == [
3688         [1, 2, 0],
3689         [4, 5, 6],
3690         [7, 8, 9],
3691         [10, 11, 12]]);
3692 }
3693 
3694 pure nothrow
3695 version(mir_test) unittest
3696 {
3697     import mir.ndslice.allocation: slice;
3698     import mir.ndslice.topology: iota;
3699 
3700     //|  1  2  3 |
3701     //|  4  5  6 |
3702     //|  7  8  9 |
3703     //| 10 11 12 |
3704     auto m = iota([4, 3], 1).slice;
3705     m.eachUpper!"a = 0"(-2);
3706     assert(m == [
3707         [0, 0, 0],
3708         [0, 0, 0],
3709         [0, 0, 0],
3710         [10, 0, 0]]);
3711 }
3712 
3713 /// Swap two slices
3714 pure nothrow
3715 version(mir_test) unittest
3716 {
3717     import mir.utility : swap;
3718     import mir.ndslice.allocation : slice;
3719     import mir.ndslice.topology : as, iota;
3720 
3721     //| 0 1 2 |
3722     //| 3 4 5 |
3723     //| 6 7 8 |
3724     auto a = iota([3, 3]).as!double.slice;
3725     //| 10 11 12 |
3726     //| 13 14 15 |
3727     //| 16 17 18 |
3728     auto b = iota([3, 3], 10).as!double.slice;
3729 
3730     eachUpper!swap(a, b);
3731 
3732     assert(a == [
3733         [0, 11, 12],
3734         [3,  4, 15],
3735         [6,  7,  8]]);
3736     assert(b == [
3737         [10,  1,  2],
3738         [13, 14,  5],
3739         [16, 17, 18]]);
3740 }
3741 
3742 /// Swap two zipped slices
3743 pure nothrow
3744 version(mir_test) unittest
3745 {
3746     import mir.utility : swap;
3747     import mir.ndslice.allocation : slice;
3748     import mir.ndslice.topology : as, zip, iota;
3749 
3750     //| 0 1 2 |
3751     //| 3 4 5 |
3752     //| 6 7 8 |
3753     auto a = iota([3, 3]).as!double.slice;
3754     //| 10 11 12 |
3755     //| 13 14 15 |
3756     //| 16 17 18 |
3757     auto b = iota([3, 3], 10).as!double.slice;
3758 
3759     auto z = zip(a, b);
3760 
3761     z.eachUpper!(z => swap(z.a, z.b));
3762 
3763     assert(a == [
3764         [0, 11, 12],
3765         [3,  4, 15],
3766         [6,  7,  8]]);
3767     assert(b == [
3768         [10,  1,  2],
3769         [13, 14,  5],
3770         [16, 17, 18]]);
3771 }
3772 
3773 // uniq
3774 /**
3775 Lazily iterates unique consecutive elements of the given range (functionality
3776 akin to the $(HTTP wikipedia.org/wiki/_Uniq, _uniq) system
3777 utility). Equivalence of elements is assessed by using the predicate
3778 $(D pred), by default $(D "a == b"). The predicate is passed to
3779 $(REF nary, mir,functional), and can either accept a string, or any callable
3780 that can be executed via $(D pred(element, element)). If the given range is
3781 bidirectional, $(D uniq) also yields a
3782 `std,range,primitives`.
3783 Params:
3784     pred = Predicate for determining equivalence between range elements.
3785 */
3786 template uniq(alias pred = "a == b")
3787 {
3788     static if (__traits(isSame, naryFun!pred, pred))
3789     {
3790         /++
3791         Params:
3792             r = An input range of elements to filter.
3793         Returns:
3794             An input range of
3795             consecutively unique elements in the original range. If `r` is also a
3796             forward range or bidirectional range, the returned range will be likewise.
3797         +/
3798         Uniq!(naryFun!pred, Range) uniq(Range)(Range r)
3799          if (isInputRange!Range && !isSlice!Range)
3800         {
3801             import core.lifetime: move;
3802             return typeof(return)(r.move);
3803         }
3804 
3805         /// ditto
3806         auto uniq(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
3807         {
3808             import mir.ndslice.topology: flattened; 
3809             import core.lifetime: move;
3810             auto r = slice.move.flattened;
3811             return Uniq!(pred, typeof(r))(move(r));
3812         }
3813     }
3814     else
3815         alias uniq = .uniq!(naryFun!pred);
3816 }
3817 
3818 ///
3819 @safe version(mir_test) unittest
3820 {
3821     int[] arr = [ 1, 2, 2, 2, 2, 3, 4, 4, 4, 5 ];
3822     assert(equal(uniq(arr), [ 1, 2, 3, 4, 5 ]));
3823 
3824     import std.algorithm.mutation : copy;
3825     // Filter duplicates in-place using copy
3826     arr.length -= arr.uniq.copy(arr).length;
3827     assert(arr == [ 1, 2, 3, 4, 5 ]);
3828 
3829     // Note that uniqueness is only determined consecutively; duplicated
3830     // elements separated by an intervening different element will not be
3831     // eliminated:
3832     assert(equal(uniq([ 1, 1, 2, 1, 1, 3, 1]), [1, 2, 1, 3, 1]));
3833 }
3834 
3835 /// N-dimensional case
3836 version(mir_test)
3837 @safe pure unittest
3838 {
3839     import mir.ndslice.fuse;
3840     import mir.ndslice.topology: byDim, map, iota;
3841 
3842     auto matrix = [ [1, 2, 2], [2, 2, 3], [4, 4, 4] ].fuse;
3843 
3844     assert(matrix.uniq.equal([ 1, 2, 3, 4 ]));
3845 
3846     // unique elements for each row
3847     assert(matrix.byDim!0.map!uniq.equal!equal([ [1, 2], [2, 3], [4] ]));
3848 }
3849 
3850 /++
3851 Authros: $(HTTP erdani.com, Andrei Alexandrescu) (original Phobos code), Ilia Ki (betterC rework)
3852 +/
3853 struct Uniq(alias pred, Range)
3854 {
3855     Range _input;
3856 
3857     void popFront() scope
3858     {
3859         assert(!empty, "Attempting to popFront an empty uniq.");
3860         auto last = _input.front;
3861         do
3862         {
3863             _input.popFront();
3864         }
3865         while (!_input.empty && pred(last, _input.front));
3866     }
3867 
3868     auto ref front() @property
3869     {
3870         assert(!empty, "Attempting to fetch the front of an empty uniq.");
3871         return _input.front;
3872     }
3873 
3874     import std.range.primitives: isBidirectionalRange;
3875 
3876     static if (isBidirectionalRange!Range)
3877     {
3878         void popBack() scope
3879         {
3880             assert(!empty, "Attempting to popBack an empty uniq.");
3881             auto last = _input.back;
3882             do
3883             {
3884                 _input.popBack();
3885             }
3886             while (!_input.empty && pred(last, _input.back));
3887         }
3888 
3889         auto ref back() return scope @property
3890         {
3891             assert(!empty, "Attempting to fetch the back of an empty uniq.");
3892             return _input.back;
3893         }
3894     }
3895 
3896     static if (isInfinite!Range)
3897     {
3898         enum bool empty = false;  // Propagate infiniteness.
3899     }
3900     else
3901     {
3902         @property bool empty() const { return _input.empty; }
3903     }
3904 
3905     ref opIndex()() scope return
3906     {
3907         return this;
3908     }
3909 
3910     auto opIndex()() const return scope
3911     {
3912         return Filter!(typeof(_input[]))(_input[]);
3913     }
3914 
3915     import std.range.primitives: isForwardRange;
3916 
3917     static if (isForwardRange!Range)
3918     {
3919         @property typeof(this) save() return scope
3920         {
3921             return typeof(this)(_input.save);
3922         }
3923     }
3924 }
3925 
3926 version(none)
3927 @safe version(mir_test) unittest
3928 {
3929     import std.internal.test.dummyrange;
3930     import std.range;
3931 
3932     int[] arr = [ 1, 2, 2, 2, 2, 3, 4, 4, 4, 5 ];
3933     auto r = uniq(arr);
3934     static assert(isForwardRange!(typeof(r)));
3935 
3936     assert(equal(r, [ 1, 2, 3, 4, 5 ][]));
3937     assert(equal(retro(r), retro([ 1, 2, 3, 4, 5 ][])));
3938 
3939     foreach (DummyType; AllDummyRanges)
3940     {
3941         DummyType d;
3942         auto u = uniq(d);
3943         assert(equal(u, [1,2,3,4,5,6,7,8,9,10]));
3944 
3945         static assert(d.rt == RangeType.Input || isForwardRange!(typeof(u)));
3946 
3947         static if (d.rt >= RangeType.Bidirectional)
3948         {
3949             assert(equal(retro(u), [10,9,8,7,6,5,4,3,2,1]));
3950         }
3951     }
3952 }
3953 
3954 @safe version(mir_test) unittest // https://issues.dlang.org/show_bug.cgi?id=17264
3955 {
3956     const(int)[] var = [0, 1, 1, 2];
3957     assert(var.uniq.equal([0, 1, 2]));
3958 }
3959 
3960 @safe version(mir_test) unittest {
3961     import mir.ndslice.allocation;
3962     import mir.math.common: approxEqual;
3963     auto x = rcslice!double(2);
3964     auto y = rcslice!double(2);
3965     x[] = [2, 3];
3966     y[] = [2, 3];
3967     assert(equal!approxEqual(x,y));
3968 }
3969 
3970 /++
3971 Implements the higher order filter function. The predicate is passed to
3972 `mir.functional.naryFun`, and can either accept a string, or any callable
3973 that can be executed via `pred(element)`.
3974 Params:
3975     pred = Function to apply to each element of range
3976 Returns:
3977     `filter!(pred)(range)` returns a new range containing only elements `x` in `range` for
3978     which `pred(x)` returns `true`.
3979 See_Also:
3980     $(HTTP en.wikipedia.org/wiki/Filter_(higher-order_function), Filter (higher-order function))
3981 Note:
3982     $(RED User and library code MUST call `empty` method ahead each call of pair or one of `front` and `popFront` methods.)
3983 +/
3984 template filter(alias pred = "a")
3985 {
3986     static if (__traits(isSame, naryFun!pred, pred))
3987     {
3988         /++
3989         Params:
3990             r = An input range of elements to filter.
3991         Returns:
3992             A new range containing only elements `x` in `range` for which `predicate(x)` returns `true`.
3993         +/
3994         Filter!(naryFun!pred, Range) filter(Range)(Range r)
3995             if (isInputRange!Range && !isSlice!Range)
3996         {
3997             import core.lifetime: move;
3998             return typeof(return)(r.move);
3999         }
4000 
4001         /// ditto
4002         auto filter(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
4003         {
4004             import mir.ndslice.topology: flattened; 
4005             import core.lifetime: move;
4006             auto r = slice.move.flattened;
4007             return Filter!(pred, typeof(r))(move(r));
4008         }
4009     }
4010     else
4011         alias filter = .filter!(naryFun!pred);
4012 }
4013 
4014 /// ditto
4015 struct Filter(alias pred, Range)
4016 {
4017     Range _input;
4018     version(assert) bool _freshEmpty;
4019 
4020     void popFront() scope
4021     {
4022         assert(!_input.empty, "Attempting to popFront an empty Filter.");
4023         version(assert) assert(_freshEmpty, "Attempting to pop the front of a Filter without calling '.empty' method ahead.");
4024         version(assert) _freshEmpty = false;
4025         _input.popFront;
4026     }
4027 
4028     auto ref front() @safe @property return scope
4029     {
4030         assert(!_input.empty, "Attempting to fetch the front of an empty Filter.");
4031         version(assert) assert(_freshEmpty, "Attempting to fetch the front of a Filter without calling '.empty' method ahead.");
4032         return _input.front;
4033     }
4034 
4035     bool empty() @safe scope @property
4036     {
4037         version(assert) _freshEmpty = true;
4038         for (;;)
4039         {
4040             if (auto r = _input.empty)
4041                 return true;
4042             if (pred(_input.front))
4043                 return false;
4044             _input.popFront;
4045         }
4046     }
4047 
4048     import std.range.primitives: isForwardRange;
4049     static if (isForwardRange!Range)
4050     {
4051         @property typeof(this) save() return scope
4052         {
4053             return typeof(this)(_input.save);
4054         }
4055     }
4056 
4057     ref opIndex()() scope return
4058     {
4059         return this;
4060     }
4061 
4062     auto opIndex()() const return scope
4063     {
4064         return Filter!(pred, typeof(_input[]))(_input[]);
4065     }
4066 }
4067 
4068 ///
4069 version(mir_test)
4070 @safe pure nothrow unittest
4071 {
4072     int[] arr = [ 0, 1, 2, 3, 4, 5 ];
4073 
4074     // Filter below 3
4075     auto small = filter!(a => a < 3)(arr);
4076     assert(equal(small, [ 0, 1, 2 ]));
4077 
4078     // Filter again, but with Uniform Function Call Syntax (UFCS)
4079     auto sum = arr.filter!(a => a < 3);
4080     assert(equal(sum, [ 0, 1, 2 ]));
4081 
4082     // Filter with the default predicate
4083     auto nonZeros = arr.filter;
4084     assert(equal(nonZeros, [ 1, 2, 3, 4, 5 ]));
4085 
4086     // In combination with concatenation() to span multiple ranges
4087     import mir.ndslice.concatenation;
4088 
4089     int[] a = [ 3, -2, 400 ];
4090     int[] b = [ 100, -101, 102 ];
4091     auto r = concatenation(a, b).filter!(a => a > 0);
4092     assert(equal(r, [ 3, 400, 100, 102 ]));
4093 
4094     // Mixing convertible types is fair game, too
4095     double[] c = [ 2.5, 3.0 ];
4096     auto r1 = concatenation(c, a, b).filter!(a => cast(int) a != a);
4097     assert(equal(r1, [ 2.5 ]));
4098 }
4099 
4100 /// N-dimensional filtering
4101 version(mir_test)
4102 @safe pure unittest
4103 {
4104     import mir.ndslice.fuse;
4105     import mir.ndslice.topology: byDim, map;
4106 
4107     auto matrix =
4108         [[   3,   -2, 400 ],
4109          [ 100, -101, 102 ]].fuse;
4110 
4111     alias filterPositive = filter!"a > 0";
4112 
4113     // filter all elements in the matrix
4114     auto r = filterPositive(matrix);
4115     assert(equal(r, [ 3, 400, 100, 102 ]));
4116 
4117     // filter all elements for each row
4118     auto rr = matrix.byDim!0.map!filterPositive;
4119     assert(equal!equal(rr, [ [3, 400], [100, 102] ]));
4120 
4121     // filter all elements for each column
4122     auto rc = matrix.byDim!1.map!filterPositive;
4123     assert(equal!equal(rc, [ [3, 100], [], [400, 102] ]));
4124 }
4125 
4126 /// N-dimensional filtering based on value in specific row/column
4127 version(mir_test)
4128 @safe pure
4129 unittest
4130 {
4131     import mir.ndslice.fuse;
4132     import mir.ndslice.topology: byDim;
4133 
4134     auto matrix =
4135     [[   3,   2, 400 ],
4136      [ 100, -101, 102 ]].fuse;
4137 
4138     // filter row based on value in index 1
4139     auto r1 = matrix.byDim!0.filter!("a[1] > 0");
4140     assert(equal!equal(r1, [ [3, 2, 400] ]));
4141 
4142     // filter column based on value in index 1
4143     auto r2 = matrix.byDim!1.filter!("a[1] > 0");
4144     assert(equal!equal(r2, [ [3, 100], [400, 102] ]));
4145 }
4146 
4147 /// Filter out NaNs
4148 version(mir_test)
4149 @safe pure nothrow
4150 unittest {
4151     import mir.algorithm.iteration: equal, filter;
4152     import mir.ndslice.slice: sliced;
4153     import std.math.traits: isNaN;
4154    
4155     static immutable result1 = [1.0, 2];
4156    
4157     double x;
4158     auto y = [1.0, 2, x].sliced;
4159     auto z = y.filter!(a => !isNaN(a));
4160     assert(z.equal(result1));
4161 }
4162 
4163 /// Filter out NaNs by row and by column
4164 version(mir_test)
4165 @safe pure
4166 unittest {
4167     import mir.algorithm.iteration: equal, filter;
4168     import mir.ndslice.fuse: fuse;
4169     import mir.ndslice.topology: byDim, map;
4170     import std.math.traits: isNaN;
4171    
4172     static immutable result1 = [[1.0, 2], [3.0, 4, 5]];
4173     static immutable result2 = [[1.0, 3], [2.0, 4], [5.0]];
4174    
4175     double x;
4176     auto y = [[1.0, 2, x], [3.0, 4, 5]].fuse;
4177 
4178     // by row
4179     auto z1 = y.byDim!0.map!(filter!(a => !isNaN(a)));
4180     assert(z1.equal!equal(result1));
4181     // by column
4182     auto z2 = y.byDim!1.map!(filter!(a => !isNaN(a)));
4183     assert(z2.equal!equal(result2));
4184 }
4185 
4186 /// Filter entire rows/columns that have NaNs
4187 version(mir_test)
4188 @safe pure
4189 unittest {
4190     import mir.algorithm.iteration: equal, filter;
4191     import mir.ndslice.fuse: fuse;
4192     import mir.ndslice.topology: byDim, map;
4193     import std.math.traits: isNaN;
4194    
4195     static immutable result1 = [[3.0, 4, 5]];
4196     static immutable result2 = [[1.0, 3], [2.0, 4]];
4197    
4198     double x;
4199     auto y = [[1.0, 2, x], [3.0, 4, 5]].fuse;
4200 
4201     // by row
4202     auto z1 = y.byDim!0.filter!(a => a.all!(b => !isNaN(b)));
4203     assert(z1.equal!equal(result1));
4204     // by column
4205     auto z2 = y.byDim!1.filter!(a => a.all!(b => !isNaN(b)));
4206     assert(z2.equal!equal(result2));
4207 }
4208 
4209 // Ensure filter works with each
4210 version(mir_test)
4211 @safe pure nothrow
4212 unittest {
4213     import mir.ndslice.slice: sliced;
4214 
4215     int[] result = [3, 4];
4216 
4217     int[] x = [1, 2, 3];
4218     auto y = x.filter!"a >= 2";
4219     y.each!"a++";
4220 
4221     assert(y.equal(result));
4222 }
4223 
4224 /++
4225 Implements the higher order filter and map function. The predicate and map functions are passed to
4226 `mir.functional.naryFun`, and can either accept a string, or any callable
4227 that can be executed via `pred(element)` and `map(element)`.
4228 Params:
4229     pred = Filter function to apply to each element of range (optional)
4230     map = Map function to apply to each  element of range
4231 Returns:
4232     `rcfilter!(pred)(range)` returns a new RCArray containing only elements `map(x)` in `range` for
4233     which `pred(x)` returns `true`.
4234 See_Also:
4235     $(HTTP en.wikipedia.org/wiki/Filter_(higher-order_function), Filter (higher-order function))
4236 +/
4237 template rcfilter(alias pred = "a", alias map = "a")
4238 {
4239     static if (__traits(isSame, naryFun!pred, pred) && __traits(isSame, naryFun!map, map))
4240     {
4241         /++ 
4242         Params:
4243             r = An input range of elements to filter.
4244         Returns:
4245             A new range containing only elements `x` in `range` for which `predicate(x)` returns `true`.
4246         +/
4247         auto rcfilter(Range)(Range r)
4248             if (isIterable!Range && (!isSlice!Range || DimensionCount!Range == 1))
4249         {
4250             import core.lifetime: forward;
4251             import mir.appender: scopedBuffer;
4252             import mir.primitives: isInputRange;
4253             import mir.rc.array: RCArray;
4254 
4255             alias T = typeof(map(r.front));
4256             auto buffer = scopedBuffer!T;
4257             foreach (ref e; r)
4258             {
4259                 if (pred(e))
4260                 {
4261                     static if (__traits(isSame, naryFun!"a", map))
4262                         buffer.put(forward!e);
4263                     else
4264                         buffer.put(map(forward!e));
4265                 }
4266             }
4267             return () @trusted
4268             {
4269                 auto ret = RCArray!T(buffer.length);
4270                 buffer.moveDataAndEmplaceTo(ret[]);
4271                 return ret;
4272             } ();
4273         }
4274 
4275         /// ditto
4276         auto rcfilter(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice)
4277             if (N > 1)
4278         {
4279             import mir.ndslice.topology: flattened; 
4280             import core.lifetime: move;
4281             return rcfilter(slice.move.flattened);
4282         }
4283     }
4284     else
4285         alias rcfilter = .rcfilter!(naryFun!pred, naryFun!map);
4286 }
4287 
4288 ///
4289 version(mir_test)
4290 @safe pure nothrow @nogc unittest
4291 {
4292     import mir.ndslice.topology: iota;
4293 
4294     auto val = 3;
4295     auto factor = 5;
4296     // Filter iota 2x3 matrix below 3
4297     assert(iota(2, 3).rcfilter!(a => a < val).moveToSlice.equal(val.iota));
4298     // Filter and map below 3
4299     assert(6.iota.rcfilter!(a => a < val, a => a * factor).moveToSlice.equal(val.iota * factor));
4300 }
4301 
4302 version(mir_test)
4303 @safe pure nothrow @nogc unittest
4304 {
4305     import mir.ndslice.topology: iota, as;
4306 
4307     auto val = 3;
4308     auto factor = 5;
4309     // Filter iota 2x3 matrix below 3
4310     assert(iota(2, 3).as!(const int).rcfilter!(a => a < val).moveToSlice.equal(val.iota));
4311     // Filter and map below 3
4312     assert(6.iota.as!(immutable int).rcfilter!(a => a < val, a => a * factor).moveToSlice.equal(val.iota * factor));
4313 }
4314 
4315 /++
4316 Implements the homonym function (also known as `accumulate`, $(D
4317 compress), `inject`, or `foldl`) present in various programming
4318 languages of functional flavor. The call `fold!(fun)(slice, seed)`
4319 first assigns `seed` to an internal variable `result`,
4320 also called the accumulator. Then, for each element `x` in $(D
4321 slice), `result = fun(result, x)` gets evaluated. Finally, $(D
4322 result) is returned.
4323 
4324 Params:
4325     fun = the predicate function to apply to the elements
4326 
4327 See_Also:
4328     $(HTTP en.wikipedia.org/wiki/Fold_(higher-order_function), Fold (higher-order function))
4329     $(LREF sum) is similar to `fold!((a, b) => a + b)` that offers
4330     precise summing of floating point numbers.
4331     This is functionally equivalent to $(LREF reduce) with the argument order
4332     reversed.
4333 +/
4334 template fold(alias fun)
4335 {
4336     /++
4337     Params:
4338         slice = A slice, range, and array.
4339         seed = An initial accumulation value.
4340     Returns:
4341         the accumulated result
4342     +/
4343     @fmamath auto fold(Slice, S)(scope Slice slice, S seed)
4344     {
4345         import core.lifetime: move;
4346         return reduce!fun(seed, slice.move);
4347     }
4348 }
4349 
4350 ///
4351 version(mir_test)
4352 @safe pure nothrow
4353 unittest
4354 {
4355     import mir.ndslice.slice: sliced;
4356     import mir.ndslice.topology: map;
4357 
4358     auto arr = [1, 2, 3, 4, 5].sliced;
4359 
4360     // Sum all elements
4361     assert(arr.fold!((a, b) => a + b)(0) == 15);
4362     assert(arr.fold!((a, b) => a + b)(6) == 21);
4363 
4364     // Can be used in a UFCS chain
4365     assert(arr.map!(a => a + 1).fold!((a, b) => a + b)(0) == 20);
4366 
4367     // Return the last element of any range
4368     assert(arr.fold!((a, b) => b)(0) == 5);
4369 }
4370 
4371 /// Works for matrices
4372 version(mir_test)
4373 @safe pure
4374 unittest
4375 {
4376     import mir.ndslice.fuse: fuse;
4377 
4378     auto arr = [
4379         [1, 2, 3], 
4380         [4, 5, 6]
4381     ].fuse;
4382 
4383     assert(arr.fold!((a, b) => a + b)(0) == 21);
4384 }
4385 
4386 version(mir_test)
4387 @safe pure nothrow
4388 unittest
4389 {
4390     import mir.ndslice.topology: map;
4391 
4392     int[] arr = [1, 2, 3, 4, 5];
4393 
4394     // Sum all elements
4395     assert(arr.fold!((a, b) => a + b)(0) == 15);
4396     assert(arr.fold!((a, b) => a + b)(6) == 21);
4397 
4398     // Can be used in a UFCS chain
4399     assert(arr.map!(a => a + 1).fold!((a, b) => a + b)(0) == 20);
4400 
4401     // Return the last element of any range
4402     assert(arr.fold!((a, b) => b)(0) == 5);
4403 }
4404 
4405 version(mir_test)
4406 @safe pure nothrow 
4407 unittest
4408 {
4409     int[] arr = [1];
4410     static assert(!is(typeof(arr.fold!()(0))));
4411     static assert(!is(typeof(arr.fold!(a => a)(0))));
4412     static assert(is(typeof(arr.fold!((a, b) => a)(0))));
4413     assert(arr.length == 1);
4414 }
4415 
4416 version(mir_test)
4417 unittest
4418 {
4419     import mir.rc.array: RCArray;
4420     import mir.algorithm.iteration: minmaxPos, minPos, maxPos, minmaxIndex, minIndex, maxIndex;
4421 
4422     static immutable a = [0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11];
4423 
4424     auto x = RCArray!double(12);
4425     foreach(i, ref e; x)
4426         e = a[i];
4427     auto y = x.asSlice;
4428     auto z0 = y.minmaxPos;
4429     auto z1 = y.minPos;
4430     auto z2 = y.maxPos;
4431     auto z3 = y.minmaxIndex;
4432     auto z4 = y.minIndex;
4433     auto z5 = y.maxIndex;
4434 }
The OpenD Programming Language
