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