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