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