1 /**
2     Defines a Region and common operation with these. A Region is a set of
3     tagged intervals where differently tagged intervals are distinct.
4 
5     Copyright: © 2018 Arne Ludwig <arne.ludwig@posteo.de>
6     License: Subject to the terms of the MIT license, as written in the
7              included LICENSE file.
8     Authors: Arne Ludwig <arne.ludwig@posteo.de>
9 */
10 module dentist.util.region;
11 
12 import std.algorithm : all, cmp, filter, map, max, min, sort, sum;
13 import std.array : appender, array, join;
14 import std.exception : assertThrown;
15 import std.format : format;
16 import std.functional : unaryFun;
17 import std.range : assumeSorted, dropExactly, ElementType, isInputRange, only, retro;
18 import std.traits : isNumeric, Unqual;
19 
20 /// Returns the type of the property `tag` of `T`.
21 template TagType(T)
22 {
23     private T instance;
24 
25     alias TagType = typeof(instance.tag);
26 }
27 
28 /// Checks if T has a property `tag` implicitly convertible to `Tag` – if given.
29 template isTaggable(T)
30 {
31     enum isTaggable = is(TagType!T);
32 }
33 
34 unittest
35 {
36     struct Taggable
37     {
38         int tag;
39     }
40 
41     static assert(isTaggable!Taggable);
42     static assert(!isTaggable!int);
43 }
44 
45 unittest
46 {
47     struct Taggable
48     {
49         int tag;
50     }
51 
52     static assert(isTaggable!Taggable);
53     static assert(!isTaggable!int);
54 }
55 
56 /// Thrown if two operands require the same tag but different were provided.
57 static class MismatchingTagsException(Tag) : Exception
58 {
59     const(Tag[2]) tags;
60 
61     this(in Tag tagA, in Tag tagB)
62     {
63         this.tags = [tagA, tagB];
64 
65         super(format!"mismatching interval tags: %s"(this.tags));
66     }
67 }
68 
69 /*
70     Throws an exception if tags do not match.
71 
72     Throws: MismatchingTagsException if tags do not match.
73 */
74 void enforceMatchingTags(Taggable)(in Taggable taggableA, in Taggable taggableB) pure
75         if (isTaggable!Taggable)
76 {
77     alias Tag = TagType!Taggable;
78 
79     if (taggableA.tag != taggableB.tag)
80     {
81         throw new MismatchingTagsException!Tag(taggableA.tag, taggableB.tag);
82     }
83 }
84 
85 /// Thrown if regions is unexpectedly empty.
86 static class EmptyRegionException : Exception
87 {
88     this()
89     {
90         super("empty region");
91     }
92 }
93 
94 /*
95     Throws an exception if region is empty.
96 
97     Throws: EmptyRegionException if region is empty.
98 */
99 void enforceNonEmpty(R)(in R region) if (is(R : Region!Args, Args...))
100 {
101     if (region.empty)
102     {
103         throw new EmptyRegionException();
104     }
105 }
106 
107 /**
108     A Region is a set of tagged intervals where differently tagged intervals are distinct.
109 */
110 struct Region(Number, Tag, string tagAlias = null, Tag emptyTag = Tag.init)
111 {
112     static assert(isNumeric!Number, "interval limits must be numeric");
113 
114     static if (__traits(compiles, Number.infinity))
115     {
116         static enum numberSup = Number.infinity;
117     }
118     else
119     {
120         static enum numberSup = Number.max;
121     }
122 
123     /**
124         This represents a single tagged point.
125     */
126     static struct TaggedPoint
127     {
128         Tag tag = emptyTag;
129         Number value;
130         static if (!(tagAlias is null) && tagAlias != "tag")
131         {
132             mixin("alias " ~ tagAlias ~ " = tag;");
133         }
134     }
135 
136     /**
137         This is a right-open interval `[begin, end)` tagged with `tag`.
138         If `tagAlias` is given then the tag may be access as a property of
139         that name.
140     */
141     static struct TaggedInterval
142     {
143         Tag tag = emptyTag;
144         Number begin;
145         Number end;
146         static if (!(tagAlias is null) && tagAlias != "tag")
147         {
148             mixin("alias " ~ tagAlias ~ " = tag;");
149         }
150 
151         invariant
152         {
153             assert(begin <= end, "begin must be less than or equal to end");
154         }
155 
156         /// Returns the size of this interval.
157         @property Number size() pure const nothrow
158         {
159             return this.end - this.begin;
160         }
161 
162         ///
163         unittest
164         {
165             assert(Region!(int, int).TaggedInterval().size == 0);
166             assert(TaggedInterval(1, 10, 20).size == 10);
167             assert(TaggedInterval(2, 20, 40).size == 20);
168             assert(TaggedInterval(3, 30, 60).size == 30);
169         }
170 
171         /// Returns true iff the interval is empty. An interval is empty iff
172         /// `begin == end`.
173         @property bool empty() pure const nothrow
174         {
175             return begin >= end;
176         }
177 
178         ///
179         unittest
180         {
181             assert(TaggedInterval().empty);
182             assert(!TaggedInterval(1, 10, 20).empty);
183             assert(!TaggedInterval(2, 20, 40).empty);
184             assert(TaggedInterval(3, 60, 60).empty);
185         }
186 
187         /**
188             Returns the convex hull of the intervals.
189 
190             Throws: MismatchingTagsException if `tag`s differ.
191         */
192         static TaggedInterval convexHull(in TaggedInterval[] intervals...) pure
193         {
194             if (intervals.length == 0)
195             {
196                 return TaggedInterval();
197             }
198 
199             TaggedInterval convexHullInterval = intervals[0];
200 
201             foreach (interval; intervals[1 .. $])
202             {
203                 enforceMatchingTags(convexHullInterval, interval);
204 
205                 if (convexHullInterval.empty)
206                 {
207                     convexHullInterval = interval;
208                 }
209                 else if (interval.empty)
210                 {
211                     continue;
212                 }
213                 else
214                 {
215                     convexHullInterval.begin = min(convexHullInterval.begin, interval.begin);
216                     convexHullInterval.end = max(convexHullInterval.end, interval.end);
217                 }
218             }
219 
220             return convexHullInterval.empty
221                 ? TaggedInterval()
222                 : convexHullInterval;
223         }
224 
225         ///
226         unittest
227         {
228             alias R = Region!(int, int);
229             alias TI = R.TaggedInterval;
230 
231             assert(TI.convexHull(TI(0, 10, 20), TI(0, 0, 5)) == TI(0, 0, 20));
232             assert(TI.convexHull(TI(0, 10, 20), TI(0, 5, 15)) == TI(0, 5, 20));
233             assert(TI.convexHull(TI(0, 10, 20), TI(0, 12, 18)) == TI(0, 10, 20));
234             assert(TI.convexHull(TI(0, 10, 20), TI(0, 10, 20)) == TI(0, 10, 20));
235             assert(TI.convexHull(TI(0, 10, 20), TI(0, 15, 25)) == TI(0, 10, 25));
236             assert(TI.convexHull(TI(0, 10, 20), TI(0, 25, 30)) == TI(0, 10, 30));
237             assertThrown!(MismatchingTagsException!int)(TI.convexHull(TI(0, 10, 20), TI(1, 25, 30)));
238         }
239 
240         /// Returns the intersection of both intervals; empty if `tag`s differ.
241         TaggedInterval opBinary(string op)(in TaggedInterval other) const pure nothrow
242                 if (op == "&")
243         {
244             if (this.tag != other.tag)
245             {
246                 return TaggedInterval();
247             }
248 
249             auto newBegin = max(this.begin, other.begin);
250             auto newEnd = min(this.end, other.end);
251 
252             if (newBegin > newEnd)
253             {
254                 return TaggedInterval();
255             }
256 
257             return TaggedInterval(
258                 tag,
259                 newBegin,
260                 newEnd,
261             );
262         }
263 
264         /// ditto
265         TaggedInterval opOpAssign(string op)(in TaggedInterval other) if (op == "&")
266         {
267             {
268                 auto tmp = this & other;
269 
270                 this.tag = tmp.tag;
271                 this.begin = tmp.begin;
272                 this.end = tmp.end;
273 
274                 return this;
275             }
276         }
277         ///
278         unittest
279         {
280             alias R = Region!(int, int);
281             alias TI = R.TaggedInterval;
282 
283             assert((TI(0, 10, 20) & TI(0, 0, 5)).empty);
284             assert((TI(0, 10, 20) & TI(0, 5, 15)) == TI(0, 10, 15));
285             assert((TI(0, 10, 20) & TI(0, 12, 18)) == TI(0, 12, 18));
286             assert((TI(0, 10, 20) & TI(0, 10, 20)) == TI(0, 10, 20));
287             assert((TI(0, 10, 20) & TI(0, 15, 25)) == TI(0, 15, 20));
288             assert((TI(0, 10, 20) & TI(0, 25, 30)).empty);
289             assert((TI(0, 10, 20) & TI(1, 25, 30)).empty);
290         }
291 
292         /// Returns the difference of both intervals.
293         Region opBinary(string op)(in TaggedInterval other) const if (op == "-")
294         {
295             auto intersection = this & other;
296 
297             if (intersection.empty)
298             {
299                 TaggedInterval thisCopy = this;
300 
301                 return Region(this.empty ? [] : [thisCopy]);
302             }
303 
304             return Region(only(
305                 TaggedInterval(
306                     tag,
307                     this.begin,
308                     intersection.begin,
309                 ),
310                 TaggedInterval(
311                     tag,
312                     intersection.end,
313                     this.end,
314                 ),
315             ).filter!"!a.empty".array);
316         }
317 
318         ///
319         unittest
320         {
321             alias R = Region!(int, int);
322             alias TI = R.TaggedInterval;
323 
324             assert(TI(0, 10, 20) - TI(0, 0, 5) == R([TI(0, 10, 20)]));
325             assert(TI(0, 10, 20) - TI(0, 5, 15) == R([TI(0, 15, 20)]));
326             assert(TI(0, 10, 20) - TI(0, 12, 18) == R([TI(0, 10, 12), TI(0, 18, 20)]));
327             assert(TI(0, 10, 20) - TI(0, 10, 20) == R([]));
328             assert(TI(0, 10, 20) - TI(0, 15, 25) == R([TI(0, 10, 15)]));
329             assert(TI(0, 10, 20) - TI(0, 25, 30) == R([TI(0, 10, 20)]));
330             assert(TI(0, 10, 20) - TI(1, 25, 30) == R([TI(0, 10, 20)]));
331         }
332 
333         int opCmp(in TaggedInterval other) const pure nothrow
334         {
335             return cmp(
336                 only(this.tag, this.begin, this.end),
337                 only(other.tag, other.begin, other.end),
338             );
339         }
340 
341         ///
342         unittest
343         {
344             alias R = Region!(int, int);
345             alias TI = R.TaggedInterval;
346 
347             assert(TI(0, 10, 20) > TI(0, 0, 5));
348             assert(TI(0, 10, 20) > TI(0, 5, 15));
349             assert(TI(0, 10, 20) < TI(0, 12, 18));
350             assert(TI(0, 10, 20) < TI(0, 15, 25));
351             assert(TI(0, 10, 20) == TI(0, 10, 20));
352             assert(TI(0, 10, 20) < TI(0, 25, 30));
353             assert(TI(0, 10, 20) < TI(1, 25, 30));
354         }
355 
356         /// Returns true iff the tagged intervals intersect.
357         bool intersects(in TaggedInterval other) const pure nothrow
358         {
359             return !(this & other).empty;
360         }
361 
362         ///
363         unittest
364         {
365             alias R = Region!(int, int);
366             alias TI = R.TaggedInterval;
367 
368             assert(!TI(0, 10, 20).intersects(TI(0, 0, 5)));
369             assert(TI(0, 10, 20).intersects(TI(0, 5, 15)));
370             assert(TI(0, 10, 20).intersects(TI(0, 12, 18)));
371             assert(TI(0, 10, 20).intersects(TI(0, 15, 25)));
372             assert(TI(0, 10, 20).intersects(TI(0, 10, 20)));
373             assert(!TI(0, 10, 20).intersects(TI(0, 25, 30)));
374             assert(!TI(0, 10, 20).intersects(TI(1, 25, 30)));
375         }
376 
377         /// Returns true iff `this` is a subset of `other`, ie. fully included _in_.
378         bool opBinary(string op)(in TaggedInterval other) const pure nothrow if (op == "in")
379         {
380             return (this & other) == this;
381         }
382 
383         ///
384         unittest
385         {
386             alias R = Region!(int, int);
387             alias TI = R.TaggedInterval;
388 
389             assert(TI(0, 0, 5) !in TI(0, 10, 20));
390             assert(TI(0, 5, 15) !in TI(0, 10, 20));
391             assert(TI(0, 12, 18) in TI(0, 10, 20));
392             assert(TI(0, 15, 25) !in TI(0, 10, 20));
393             assert(TI(0, 10, 20) in TI(0, 10, 20));
394             assert(TI(0, 25, 30) !in TI(0, 10, 20));
395             assert(TI(1, 12, 18) !in TI(0, 10, 20));
396         }
397 
398         /// Returns true iff the tagged intervals do not intersect and `this < other`.
399         bool isStrictlyBefore(in TaggedInterval other) const pure nothrow
400         {
401             return !this.intersects(other) && this < other;
402         }
403 
404         ///
405         unittest
406         {
407             alias R = Region!(int, int);
408             alias TI = R.TaggedInterval;
409 
410             assert(!TI(0, 10, 20).isStrictlyBefore(TI(0, 0, 5)));
411             assert(!TI(0, 10, 20).isStrictlyBefore(TI(0, 5, 15)));
412             assert(!TI(0, 10, 20).isStrictlyBefore(TI(0, 12, 18)));
413             assert(!TI(0, 10, 20).isStrictlyBefore(TI(0, 15, 25)));
414             assert(!TI(0, 10, 20).isStrictlyBefore(TI(0, 10, 20)));
415             assert(TI(0, 10, 20).isStrictlyBefore(TI(0, 25, 30)));
416             assert(TI(0, 10, 20).isStrictlyBefore(TI(1, 25, 30)));
417         }
418 
419         /// Returns true iff the tagged intervals do not intersect and `this > other`.
420         bool isStrictlyAfter(in TaggedInterval other) const pure nothrow
421         {
422             return !this.intersects(other) && this > other;
423         }
424 
425         ///
426         unittest
427         {
428             alias R = Region!(int, int);
429             alias TI = R.TaggedInterval;
430 
431             assert(TI(0, 10, 20).isStrictlyAfter(TI(0, 0, 5)));
432             assert(!TI(0, 10, 20).isStrictlyAfter(TI(0, 5, 15)));
433             assert(!TI(0, 10, 20).isStrictlyAfter(TI(0, 12, 18)));
434             assert(!TI(0, 10, 20).isStrictlyAfter(TI(0, 15, 25)));
435             assert(!TI(0, 10, 20).isStrictlyAfter(TI(0, 10, 20)));
436             assert(!TI(0, 10, 20).isStrictlyAfter(TI(0, 25, 30)));
437             assert(!TI(0, 10, 20).isStrictlyAfter(TI(1, 25, 30)));
438         }
439     }
440 
441     ///
442     unittest
443     {
444         static enum emptyTag = 42;
445         alias R = Region!(int, int, "bucketId", emptyTag);
446         alias TI = R.TaggedInterval;
447 
448         TI emptyInterval;
449 
450         // Default constructor produces empty interval.
451         assert((emptyInterval).empty);
452         assert(emptyInterval.tag == emptyTag);
453 
454         auto ti1 = TI(1, 0, 10);
455 
456         // The tag can be aliased:
457         assert(ti1.tag == ti1.bucketId);
458 
459         auto ti2 = TI(1, 5, 15);
460 
461         // Tagged intervals with the same tag behave like regular intervals:
462         assert((ti1 & ti2) == TI(1, 5, 10));
463 
464         auto ti3 = TI(2, 0, 10);
465 
466         // Tagged intervals with different tags are distinct:
467         assert((ti1 & ti3).empty);
468     }
469 
470     TaggedInterval[] _intervals;
471 
472     this(TaggedInterval[] intervals)
473     {
474         this._intervals = intervals;
475         this.normalize();
476     }
477 
478     ///
479     unittest
480     {
481         alias R = Region!(int, int);
482         alias TI = R.TaggedInterval;
483 
484         auto region = R([TI(0, 15, 20), TI(0, 5, 10), TI(0, 0, 10)]);
485 
486         // Intervals get implicitly normalized.
487         assert(region.intervals == [TI(0, 0, 10), TI(0, 15, 20)]);
488     }
489 
490     this(TaggedInterval interval)
491     {
492         this._intervals = [interval];
493     }
494 
495     ///
496     unittest
497     {
498         alias R = Region!(int, int);
499         alias TI = R.TaggedInterval;
500 
501         R region = TI(0, 15, 20);
502 
503         assert(region.intervals == [TI(0, 15, 20)]);
504     }
505 
506     this(Tag tag, Number begin, Number end)
507     {
508         this._intervals = [TaggedInterval(tag, begin, end)];
509     }
510 
511     ///
512     unittest
513     {
514         alias R = Region!(int, int);
515         alias TI = R.TaggedInterval;
516 
517         auto region = R(0, 15, 20);
518 
519         assert(region.intervals == [TI(0, 15, 20)]);
520     }
521 
522     unittest
523     {
524         alias R = Region!(int, int);
525         alias TI = R.TaggedInterval;
526 
527         auto tis = [TI(0, 15, 20)];
528         auto region = R(tis);
529         auto regionDup = region;
530 
531         assert(region == regionDup);
532 
533         // Changing the original does not affect duplicate
534         region |= R(TI(0, 25, 30));
535 
536         assert(region != regionDup);
537     }
538 
539     /// Return a list of the tagged intervals in this region.
540     @property const(TaggedInterval)[] intervals() const pure nothrow
541     {
542         return _intervals;
543     }
544 
545     ///
546     unittest
547     {
548         alias R = Region!(int, int);
549         alias TI = R.TaggedInterval;
550 
551         R emptyRegion;
552         auto region1 = R(0, 5, 10);
553         auto region2 = R([TI(0, 5, 10), TI(0, 15, 20)]);
554 
555         assert(emptyRegion.intervals == []);
556         assert(region1.intervals == [TI(0, 5, 10)]);
557         assert(region2.intervals == [TI(0, 5, 10), TI(0, 15, 20)]);
558     }
559 
560     /// Returns the size of this region.
561     Number size() pure const nothrow
562     {
563         return _intervals.map!"a.size".sum;
564     }
565 
566     ///
567     unittest
568     {
569         alias R = Region!(int, int);
570         alias TI = R.TaggedInterval;
571 
572         R emptyRegion;
573         auto region1 = R(0, 0, 10);
574         auto region2 = R([TI(0, 0, 10), TI(0, 20, 30)]);
575         auto region3 = R([TI(0, 0, 20), TI(0, 10, 30)]);
576 
577         assert(emptyRegion.size == 0);
578         assert(region1.size == 10);
579         assert(region2.size == 20);
580         assert(region3.size == 30);
581     }
582 
583     /// Returns true iff the region is empty.
584     bool empty() pure const nothrow
585     {
586         return _intervals.length == 0;
587     }
588 
589     ///
590     unittest
591     {
592         alias R = Region!(int, int);
593         alias TI = R.TaggedInterval;
594 
595         R emptyRegion1;
596         auto emptyRegion2 = R([TI(0, 0, 0), TI(0, 10, 10)]);
597         auto emptyRegion3 = R([TI(0, 0, 0), TI(1, 0, 0)]);
598         auto region1 = R(0, 0, 10);
599 
600         assert(emptyRegion1.empty);
601         assert(emptyRegion2.empty);
602         assert(emptyRegion3.empty);
603         assert(!region1.empty);
604     }
605 
606     protected void normalize()
607     {
608         if (_intervals.length == 0)
609         {
610             return;
611         }
612 
613         _intervals.sort();
614         TaggedInterval accInterval = _intervals[0];
615         size_t insertIdx = 0;
616 
617         foreach (i, intervalB; _intervals[1 .. $])
618         {
619             if (intervalB.empty)
620             {
621                 continue;
622             }
623             else if (accInterval.intersects(intervalB))
624             {
625                 // If two intervals intersect their union is the same as the convex hull of both.
626                 accInterval = TaggedInterval.convexHull(accInterval, intervalB);
627             }
628             else
629             {
630                 if (!accInterval.empty)
631                 {
632                     _intervals[insertIdx++] = accInterval;
633                 }
634                 accInterval = intervalB;
635             }
636         }
637 
638         if (!accInterval.empty)
639         {
640             _intervals[insertIdx++] = accInterval;
641         }
642         _intervals.length = insertIdx;
643     }
644 
645     unittest
646     {
647         alias R = Region!(int, int);
648         alias TI = R.TaggedInterval;
649 
650         auto region = R([
651             TI(3, 8349, 8600),
652             TI(3, 8349, 8349),
653             TI(3, 8349, 8850),
654             TI(3, 8349, 9100),
655             TI(3, 8349, 8349),
656             TI(3, 8349, 9350),
657             TI(3, 8349, 9600),
658             TI(3, 8349, 8349),
659             TI(3, 8349, 9900),
660             TI(3, 8349, 10150),
661             TI(3, 8349, 8349),
662             TI(3, 8349, 10400),
663             TI(3, 8349, 10650),
664             TI(3, 8349, 10800),
665             TI(3, 8499, 10800),
666             TI(3, 8749, 10800),
667             TI(3, 8749, 8749),
668             TI(3, 8999, 10800),
669             TI(3, 9249, 10800),
670             TI(3, 9549, 10800),
671             TI(3, 9799, 10800),
672             TI(3, 10049, 10800),
673             TI(3, 10299, 10800),
674             TI(3, 10549, 10800),
675         ]);
676         auto normalizedRegion = R([
677             TI(3, 8349, 10800),
678         ]);
679 
680         assert(region == normalizedRegion);
681     }
682 
683     /// Computes the union of all tagged intervals.
684     Region opBinary(string op)(in Region other) const if (op == "|")
685     {
686         return Region(this._intervals.dup ~ other._intervals.dup);
687     }
688 
689     /// ditto
690     Region opBinary(string op)(in TaggedInterval other) const if (op == "|")
691     {
692         return Region(this._intervals.dup ~ [cast(TaggedInterval) other]);
693     }
694 
695     ///
696     unittest
697     {
698         alias R = Region!(int, int);
699         alias TI = R.TaggedInterval;
700 
701         assert((R(0, 10, 20) | R(0, 0, 5)) == R([TI(0, 10, 20), TI(0, 0, 5)]));
702         assert((R(0, 10, 20) | R(0, 5, 15)) == R(0, 5, 20));
703         assert((R(0, 10, 20) | R(0, 12, 18)) == R(0, 10, 20));
704         assert((R(0, 10, 20) | R(0, 10, 20)) == R(0, 10, 20));
705         assert((R(0, 10, 20) | R(0, 15, 25)) == R(0, 10, 25));
706         assert((R(0, 10, 20) | R(0, 25, 30)) == R([TI(0, 10, 20), TI(0, 25, 30)]));
707         assert((R(0, 10, 20) | R(1, 25, 30)) == R([TI(0, 10, 20), TI(1, 25, 30)]));
708     }
709 
710     /// Computes the intersection of the two regions.
711     Region opBinary(string op)(in Region other) const if (op == "&")
712     {
713         if (this.empty || other.empty)
714         {
715             return Region();
716         }
717 
718         auto intersectionAcc = appender!(TaggedInterval[]);
719         intersectionAcc.reserve(this._intervals.length + other._intervals.length);
720 
721         size_t lhsIdx = 0;
722         size_t lhsLength = this._intervals.length;
723         TaggedInterval lhsInterval;
724         size_t rhsIdx = 0;
725         size_t rhsLength = other._intervals.length;
726         TaggedInterval rhsInterval;
727 
728         while (lhsIdx < lhsLength && rhsIdx < rhsLength)
729         {
730             lhsInterval = this._intervals[lhsIdx];
731             rhsInterval = other._intervals[rhsIdx];
732 
733             if (lhsInterval.isStrictlyBefore(rhsInterval))
734             {
735                 ++lhsIdx;
736             }
737             else if (rhsInterval.isStrictlyBefore(lhsInterval))
738             {
739                 ++rhsIdx;
740             }
741             else
742             {
743                 assert(lhsInterval.intersects(rhsInterval));
744 
745                 intersectionAcc ~= lhsInterval & rhsInterval;
746 
747                 if (lhsIdx + 1 < lhsLength && rhsIdx + 1 < rhsLength)
748                 {
749                     if (this._intervals[lhsIdx + 1].begin < other._intervals[rhsIdx + 1].begin)
750                     {
751                         ++lhsIdx;
752                     }
753                     else
754                     {
755                         ++rhsIdx;
756                     }
757                 }
758                 else if (lhsIdx + 1 < lhsLength)
759                 {
760                     ++lhsIdx;
761                 }
762                 else
763                 {
764                     ++rhsIdx;
765                 }
766             }
767         }
768 
769         return Region(intersectionAcc.data);
770     }
771 
772     /// ditto
773     Region opBinary(string op)(in TaggedInterval other) const if (op == "&")
774     {
775         return this & Region([other]);
776     }
777 
778     ///
779     unittest
780     {
781         alias R = Region!(int, int);
782         alias TI = R.TaggedInterval;
783 
784         assert((R(0, 10, 20) & R(0, 0, 5)) == R([]));
785         assert((R(0, 10, 20) & R(0, 5, 15)) == R(0, 10, 15));
786         assert((R(0, 10, 20) & R(0, 12, 18)) == R(0, 12, 18));
787         assert((R(0, 10, 20) & R(0, 10, 20)) == R(0, 10, 20));
788         assert((R(0, 10, 20) & R(0, 15, 25)) == R(0, 15, 20));
789         assert((R(0, 10, 20) & R(0, 25, 30)) == R([]));
790         assert((R(0, 10, 20) & R(1, 25, 30)) == R([]));
791         // R1:       [-------)   [-------)   [-------)
792         // R2:             [-------)   [-------)   [-------)
793         // R1 & R2:        [-)   [-)   [-)   [-)   [-)
794         assert((R([
795             TI(0, 0, 30),
796             TI(0, 40, 70),
797             TI(0, 80, 110),
798         ]) & R([
799             TI(0, 20, 50),
800             TI(0, 60, 90),
801             TI(0, 100, 130),
802         ])) == R([
803             TI(0, 20, 30),
804             TI(0, 40, 50),
805             TI(0, 60, 70),
806             TI(0, 80, 90),
807             TI(0, 100, 110),
808         ]));
809     }
810 
811     Region opBinary(string op)(in TaggedInterval interval) const if (op == "-")
812     {
813         if (interval.empty)
814         {
815             return Region(_intervals.dup);
816         }
817 
818         auto differenceAcc = appender!(TaggedInterval[]);
819         differenceAcc.reserve(_intervals.length + 1);
820 
821         auto trisection = _intervals
822             .assumeSorted!"a.isStrictlyBefore(b)"
823             .trisect(interval);
824         auto copyHeadEnd = trisection[0].length;
825         auto copyTailBegin = trisection[0].length + trisection[1].length;
826 
827         differenceAcc ~= _intervals[0 .. copyHeadEnd];
828         foreach (lhsInterval; trisection[1])
829         {
830             assert(lhsInterval.intersects(interval));
831             auto tmpDifference = lhsInterval - interval;
832 
833             differenceAcc ~= tmpDifference._intervals;
834         }
835         differenceAcc ~= _intervals[copyTailBegin .. $];
836 
837         return Region(differenceAcc.data);
838     }
839 
840     Region opBinary(string op)(in Region other) const if (op == "-")
841     {
842         if (other.empty)
843         {
844             return Region(this._intervals.dup);
845         }
846 
847         auto otherDifferenceCandidates = getDifferenceCandidates(other).assumeSorted!"a.isStrictlyBefore(b)";
848         auto differenceAcc = appender!(TaggedInterval[]);
849         differenceAcc.reserve(this._intervals.length + otherDifferenceCandidates.length);
850 
851         foreach (lhsInterval; this._intervals)
852         {
853             auto intersectingRhs = otherDifferenceCandidates.equalRange(lhsInterval);
854             auto tmpDifference = Region(lhsInterval);
855 
856             foreach (rhsInterval; intersectingRhs)
857             {
858                 if (tmpDifference.empty
859                         || tmpDifference._intervals[$ - 1].isStrictlyBefore(rhsInterval))
860                 {
861                     // Remaining rhsItervals will not intersect anymore.
862                     break;
863                 }
864 
865                 tmpDifference -= rhsInterval;
866             }
867 
868             if (!tmpDifference.empty)
869             {
870                 differenceAcc ~= tmpDifference._intervals;
871             }
872         }
873 
874         return Region(differenceAcc.data);
875     }
876 
877     ///
878     unittest
879     {
880         alias R = Region!(int, int);
881         alias TI = R.TaggedInterval;
882 
883         assert((R(0, 10, 20) - R(0, 0, 5)) == R(0, 10, 20));
884         assert((R(0, 10, 20) - R(0, 5, 15)) == R(0, 15, 20));
885         assert((R(0, 10, 20) - R(0, 12, 18)) == R([TI(0, 10, 12), TI(0, 18, 20)]));
886         assert((R(0, 10, 20) - R(0, 10, 20)).empty);
887         assert((R(0, 10, 20) - R(0, 15, 25)) == R(0, 10, 15));
888         assert((R(0, 10, 20) - R(0, 25, 30)) == R(0, 10, 20));
889         assert((R(0, 10, 20) - R(1, 25, 30)) == R(0, 10, 20));
890     }
891 
892     private auto getDifferenceCandidates(in Region other) const pure nothrow
893     {
894         auto otherIntervals = other._intervals.assumeSorted!"a.isStrictlyBefore(b)";
895         auto sliceBegin = otherIntervals.lowerBound(this._intervals[0]).length;
896         auto sliceEnd = other._intervals.length - otherIntervals.upperBound(this._intervals[$ - 1]).length;
897 
898         if (sliceBegin < sliceEnd)
899         {
900             return other._intervals[sliceBegin .. sliceEnd];
901         }
902         else
903         {
904             return cast(typeof(other._intervals)) [];
905         }
906     }
907 
908     Region opOpAssign(string op, T)(in T other)
909             if (is(T : Region) || is(T : TaggedInterval))
910     {
911         static if (op == "|" || op == "&" || op == "-")
912         {
913             mixin("auto tmp = this " ~ op ~ " other;");
914 
915             this._intervals = tmp._intervals;
916 
917             return this;
918         }
919         else
920         {
921             static assert(0, "unsupported operator: " ~ op);
922         }
923     }
924 
925     unittest
926     {
927         alias R = Region!(int, int);
928         alias TI = R.TaggedInterval;
929 
930         R accRegion;
931         auto inputRegion1 = R([
932             TI(3, 8349, 8600),
933             TI(3, 8349, 8850),
934             TI(3, 8349, 9100),
935             TI(3, 8349, 9350),
936             TI(3, 8349, 9600),
937             TI(3, 8349, 9900),
938             TI(3, 8349, 10150),
939             TI(3, 8349, 10400),
940             TI(3, 8349, 10650),
941             TI(3, 8349, 10800),
942             TI(3, 8499, 10800),
943             TI(3, 8749, 10800),
944             TI(3, 8999, 10800),
945             TI(3, 9249, 10800),
946             TI(3, 9549, 10800),
947             TI(3, 9799, 10800),
948             TI(3, 10049, 10800),
949             TI(3, 10299, 10800),
950             TI(3, 10549, 10800),
951         ]);
952         auto inputRegion2 = R([
953             TI(3, 2297, 11371),
954         ]);
955         auto expectedResult = R([
956             TI(3, 2297, 11371),
957         ]);
958 
959         accRegion |= inputRegion1;
960 
961         assert(accRegion == inputRegion1);
962 
963         accRegion |= inputRegion2;
964 
965         assert(accRegion == expectedResult);
966         assert((inputRegion1 | inputRegion2) == expectedResult);
967     }
968 
969     /// Returns true iff point is in this region.
970     bool opBinaryRight(string op)(in TaggedPoint point) const pure nothrow
971             if (op == "in")
972     {
973         if (this.empty)
974         {
975             return false;
976         }
977 
978         auto needle = TaggedInterval(point.tag, point.value, numberSup);
979         auto sortedIntervals = intervals.assumeSorted;
980         auto candidateIntervals = sortedIntervals.lowerBound(needle).retro;
981 
982         if (candidateIntervals.empty)
983         {
984             return false;
985         }
986 
987         auto candidateInterval = candidateIntervals.front;
988 
989         return candidateInterval.begin <= point.value && point.value < candidateInterval.end;
990     }
991 
992     /// ditto
993     alias includes = opBinaryRight!"in";
994 
995     ///
996     unittest
997     {
998         alias R = Region!(int, int);
999         alias TI = R.TaggedInterval;
1000         alias TP = R.TaggedPoint;
1001 
1002         R emptyRegion;
1003         auto region = R([TI(0, 0, 10), TI(1, 0, 10)]);
1004 
1005         assert(TP(0, 0) !in emptyRegion);
1006         assert(TP(0, 5) !in emptyRegion);
1007         assert(TP(0, 10) !in emptyRegion);
1008         assert(TP(0, 20) !in emptyRegion);
1009         assert(TP(0, 0) in region);
1010         assert(TP(0, 5) in region);
1011         assert(TP(0, 10) !in region);
1012         assert(TP(0, 20) !in region);
1013         assert(TP(1, 0) in region);
1014         assert(TP(1, 5) in region);
1015         assert(TP(1, 10) !in region);
1016         assert(TP(1, 20) !in region);
1017     }
1018 }
1019 
1020 unittest
1021 {
1022     Region!(int, int) r;
1023 }
1024 
1025 /**
1026     Returns true iff `thing` is empty
1027 
1028     See_Also: Region.empty, Region.TaggedInterval.empty
1029 */
1030 bool empty(T)(in T thing) pure nothrow
1031         if (is(T : Region!Args, Args...) || is(T : Region!Args.TaggedInterval, Args...))
1032 {
1033     return thing.empty;
1034 }
1035 
1036 ///
1037 unittest
1038 {
1039     alias R = Region!(int, int);
1040     alias TI = R.TaggedInterval;
1041 
1042     R emptyRegion;
1043     TI emptyTI;
1044     auto region = R(0, 0, 10);
1045     auto ti = TI(0, 0, 10);
1046 
1047     assert(empty(emptyRegion));
1048     assert(empty(emptyTI));
1049     assert(!empty(region));
1050     assert(!empty(ti));
1051 }
1052 
1053 /**
1054     Returns the union of all elements.
1055 
1056     See_Also: Region.opBinary!"|", Region.TaggedInterval.opBinary!"|"
1057 */
1058 auto union_(Range)(Range regions)
1059         if (isInputRange!Range && is(ElementType!Range : Region!Args, Args...))
1060 {
1061     alias Region = Unqual!(ElementType!Range);
1062 
1063     return Region(regions
1064         .map!"a._intervals.dup"
1065         .join);
1066 }
1067 
1068 ///
1069 unittest
1070 {
1071     alias R = Region!(int, int);
1072     alias TI = R.TaggedInterval;
1073 
1074     R emptyRegion;
1075     TI emptyTI;
1076     auto region = R(0, 0, 10);
1077     auto ti = TI(0, 0, 10);
1078 
1079     assert(empty(emptyRegion));
1080     assert(empty(emptyTI));
1081     assert(!empty(region));
1082     assert(!empty(ti));
1083 }
1084 
1085 /**
1086     Returns the minimum/supremum point of the intervals. Both values are
1087     undefined for empty regions.
1088 
1089     Throws: MismatchingTagsException if `tag`s differ.
1090     Throws: EmptyRegionException if region is empty.
1091 */
1092 auto min(R)(R region) if (is(R : Region!Args, Args...))
1093 {
1094     alias TI = R.TaggedInterval;
1095 
1096     enforceNonEmpty(region);
1097     auto convexHull = TI.convexHull(region.intervals[0], region.intervals[$ - 1]);
1098 
1099     return convexHull.begin;
1100 }
1101 
1102 /// ditto
1103 auto sup(R)(R region) if (is(R : Region!Args, Args...))
1104 {
1105     alias TI = R.TaggedInterval;
1106 
1107     enforceNonEmpty(region);
1108     auto convexHull = TI.convexHull(region.intervals[0], region.intervals[$ - 1]);
1109 
1110     return convexHull.end;
1111 }
1112 
1113 ///
1114 unittest
1115 {
1116     alias R = Region!(int, int);
1117     alias TI = R.TaggedInterval;
1118 
1119     R emptyRegion;
1120     auto region1 = R([TI(0, 0, 10), TI(0, 20, 30)]);
1121     auto region2 = R([TI(0, 0, 10), TI(1, 0, 10)]);
1122 
1123     assert(min(region1) == 0);
1124     assert(sup(region1) == 30);
1125     assertThrown!EmptyRegionException(min(emptyRegion));
1126     assertThrown!EmptyRegionException(sup(emptyRegion));
1127     assertThrown!(MismatchingTagsException!int)(min(region2));
1128     assertThrown!(MismatchingTagsException!int)(sup(region2));
1129 }
1130 
1131 struct Tiling(R, N, T)
1132 {
1133     R region;
1134     N totalOverlap;
1135     T[] elements;
1136 }
1137 
1138 auto findTilings(alias toInterval, T, N)(T[] elements, in N maxLocalOverlap, in N maxGlobalOverlap = N.max)
1139 {
1140     alias interval = unaryFun!toInterval;
1141     alias Region = typeof(interval(elements[0]) - interval(elements[0]));
1142     alias Num = typeof(interval(elements[0]).size);
1143     alias overlapEachOther = (lhs, rhs) => interval(lhs).intersects(interval(rhs));
1144 
1145     auto tilingsAcc = appender!(Tiling!(Region, Num, T)[]);
1146 
1147     void findAllowedTilings(T[] candidates, Region region, T[] included, in Num globalOverlap)
1148     {
1149         if (globalOverlap > maxGlobalOverlap)
1150             return;
1151 
1152         tilingsAcc ~= Tiling!(Region, Num, T)(
1153             region,
1154             globalOverlap,
1155             included,
1156         );
1157 
1158         size_t numExpansions;
1159         foreach (i, candidate; candidates)
1160         {
1161             auto candidateInterval = interval(candidate);
1162             auto newRegion = region | candidateInterval;
1163             auto localOverlap = region.size + candidateInterval.size - newRegion.size;
1164 
1165             if (localOverlap <= maxLocalOverlap)
1166             {
1167                 findAllowedTilings(
1168                     candidates.dropExactly(i + 1),
1169                     newRegion,
1170                     included ~ [candidate],
1171                     globalOverlap + localOverlap,
1172                 );
1173                 ++numExpansions;
1174             }
1175         }
1176     }
1177 
1178     findAllowedTilings(
1179         elements,
1180         Region(),
1181         [],
1182         Num.init,
1183     );
1184 
1185     assert(tilingsAcc.data.length > 0);
1186 
1187     return tilingsAcc.data;
1188 }
1189 
1190 unittest
1191 {
1192     alias R = Region!(int, int);
1193     alias TI = R.TaggedInterval;
1194     auto elements = [
1195         [1, 5],
1196         [3, 9],
1197         [7, 10],
1198         [1, 10],
1199     ];
1200     auto tilings = findTilings!(
1201         interval => TI(0, interval[0], interval[1])
1202     )(
1203         elements,
1204         2,
1205         4,
1206     );
1207 
1208     import std.stdio;
1209 
1210     assert(tilings == [
1211         Tiling!(R, int, int[])(
1212             R([]),
1213             0,
1214             [],
1215         ),
1216         Tiling!(R, int, int[])(
1217             R([TI(0, 1, 5)]),
1218             0,
1219             [[1, 5]],
1220         ),
1221         Tiling!(R, int, int[])(
1222             R([TI(0, 1, 9)]),
1223             2,
1224             [[1, 5], [3, 9]],
1225         ),
1226         Tiling!(R, int, int[])(
1227             R([TI(0, 1, 10)]),
1228             4,
1229             [[1, 5], [3, 9], [7, 10]],
1230         ),
1231         Tiling!(R, int, int[])(
1232             R([TI(0, 1, 5), TI(0, 7, 10)]),
1233             0,
1234             [[1, 5], [7, 10]],
1235         ),
1236         Tiling!(R, int, int[])(
1237             R([TI(0, 3, 9)]),
1238             0,
1239             [[3, 9]],
1240         ),
1241         Tiling!(R, int, int[])(
1242             R([TI(0, 3, 10)]),
1243             2,
1244             [[3, 9], [7, 10]],
1245         ),
1246         Tiling!(R, int, int[])(
1247             R([TI(0, 7, 10)]),
1248             0,
1249             [[7, 10]],
1250         ),
1251         Tiling!(R, int, int[])(
1252             R([TI(0, 1, 10)]),
1253             0,
1254             [[1, 10]],
1255         ),
1256 ]);
1257 }