Data Structures exam

Readme.md

Orthogonal Range Query

Реализована следующая структура данных на языке D (dlang.org):

• Orthogonal Range Query: статическая задача в R^d c O(log^(d-1) n) на запрос

В качестве основного источника теории использовались эти слайды: http://www.cs.ucsb.edu/~suri/cs235/RangeSearching.pdf

Компиляция с помощью DMD (http://dlang.org/download.html#dmd , использовалась версия DMD64 D Compiler v2.072.2):

dmd -oforq -release -inline -O orq.d

Для достижения сложности запроса O(log^(d-1) n) вместо O(log^d n) реализован fractional cascading, который позволяет делать бинарный поиск на последнем уровне только один раз: там, где разделяются правые и левые пути.

Далее ORQ значит Orthogonal Range Query.

Пояснения по реализации

• Vec -- структура вектора со статическим размером. Позволяет получать компонент вектора по индексу.

• DumbQuery -- класс для проверки ответов. Умеет делать ORQ по заданному массиву методом линейного поиска.

• mergeArrays -- сливает два отсортированных массива. Почему-то реализация из стандартной библиотеки (std.algorithm : merge) работает в 1,5 раза медленней.

• computeCascade -- генерирует ссылки для последующего каскадирования вниз из массива и результата его слияния с другим массивом.

• RangeTree -- класс, реализующий ORQ. Далее всё внутри него:

• TreeNode -- интерфейс внутренней структуры поиска. Статически параметризуется "рангом" -- той координатой, по которой строится дерево на этом уровне.

• LeafNode -- специализация для листов дерева.

• InnerNode -- основная структура поиска. На последнем уровне просто содержит отсортированный массив, на остальных -- двоичное дерево поиска и указателем на дерево следующего уровня для всех точек, входящих в это поддерево. На предпоследнем уровне также хранит указатели для каскадного спуска и дополнительный аргумент в функциях, чтоб этот индекс прокидывать вниз.

orq.d

import std.stdio : writeln;
import std.array : appender, array, uninitializedArray;
import std.exception : enforce;
import std.algorithm : sort, isSorted, filter, until, equal;
import std.range : drop, take, put, assumeSorted, generate, isInputRange, isOutputRange, ElementType;
import std.functional : not;
import std.random : uniform;


struct Vec(size_t _dim = 2)
{
	static assert(_dim > 0, "Can't have less than 1 dimension");

	alias dim = _dim;
	double[dim] x;

	this(double[dim] args...) { x[] = args[]; }
	int opCmp(ref const Vec!dim v2) const
	{
		import std.math : cmp;
		foreach(d; 0 .. dim)
		{
			int c = cmp(x[d], v2[d]);
			if(c != 0) return c;
		}
		return 0;
	}

	alias x this;
}

alias Vec1 = Vec!1;
alias Vec2 = Vec!2;
alias Vec3 = Vec!3;


class DumbQuery(V)
{
	V[] storage;

	this(T)(T data) { prepareData(data); }

	void prepareData(Stuff)(Stuff data)
		if(isInputRange!Stuff && is(ElementType!Stuff == V))
	{
		storage = data.array;
	}

	auto rangeQuery(V min, V max)
	{
		foreach(d; 0 .. V.dim) enforce(min[d] < max[d], "Misformed query");
		return storage.filter!((element) {
			foreach(d; 0 .. V.dim) if(min[d] > element[d] || max[d] < element[d]) return false;
			return true;
		}).array;
	}
}

import std.functional : binaryFun;
E[] mergeArrays(alias less = "a < b", E)(E[] a, E[] b)
in
{
	assert(isSorted!(binaryFun!less)(a));
	assert(isSorted!(binaryFun!less)(b));
}
out(result)
{
	assert(isSorted!(binaryFun!less)(result));
}
body
{
	alias comp = binaryFun!less;
	size_t res_length = a.length + b.length;
	E[] res = uninitializedArray!(E[])(res_length);
	size_t index_a = 0, index_b = 0, index_res = 0;
	for(index_res = 0; index_res < res_length; index_res++)
	{
		bool a_empty = !(index_a < a.length), b_empty = !(index_b < b.length);
		if(a_empty) res[index_res] = b[index_b++];
		else if(b_empty) res[index_res] = a[index_a++];
		else
		{
			if(comp(a[index_a], b[index_b])) res[index_res] = a[index_a++];
			else res[index_res] = b[index_b++];
		}
	}
	return res;
}

size_t[] computeCascade(alias less = "a < b", E)(E[] from, E[] to)
in
{
	assert(isSorted!(binaryFun!less)(from));
	assert(isSorted!(binaryFun!less)(to));
	assert(from.length >= to.length);
}
out(result)
{
	assert(result.length == from.length + 1);
	foreach(i, e; from)
	{
		assert(result[i] == 0 || binaryFun!less(to[result[i] - 1], e));
	}
	assert(result[result.length - 1] == to.length);
}
body
{
	alias comp = binaryFun!less;
	size_t[] result = uninitializedArray!(size_t[])(from.length + 1);
	size_t to_index = 0;
	for(size_t i = 0; i < from.length; i++)
	{
		while(to_index < to.length && comp(to[to_index], from[i])) to_index++;
		result[i] = to_index;
	}
	result[from.length] = to.length;
	return result;
}


class RangeTree(V)
{
	import std.range : isInputRange, ElementType, isOutputRange;

	alias VElement = typeof(V.init[0]);
	enum MaxDimension = V.dim - 1;
	alias Output = typeof(appender!(V[])());

	private TreeNode!(0, Output) root = null;

	this(T)(T data) { prepareData(data); }



	void prepareData(Stuff)(Stuff data)
		if(isInputRange!Stuff && is(ElementType!Stuff == V))
	{
		V[] sorted = data.sort!(lessByDimension!0).array;
		root = makeTree!0(sorted);
	}

	auto rangeQuery(V min, V max)
	{
		foreach(d; 0 .. MaxDimension + 1) enforce(min[d] < max[d], "Misformed query");
		auto output = appender!(V[])();
		if(root !is null)
		{
			static if(MaxDimension == 0) root.rangeQuery(min, max, output, root.getCascadeIndex(min));
			else root.rangeQuery(min, max, output);
		}
		return output.data;
	}

	static bool lessByDimension(size_t compDim)(in V a, in V b)
	{
		static assert(compDim <= MaxDimension, "Unable to compare by this dimension");
		return a[compDim] < b[compDim];
	}

	static TreeNode!(rank, Output) makeTree(size_t rank)(ref V[] data)
	{
		if(data.length > 1) return new InnerNode!(rank, Output)(data);
		else return new LeafNode!(rank, Output)(data);
	}

	mixin template CommonTreeDeclarations(size_t rank, size_t maxRank, R)
	{
		static assert(rank <= MaxDimension);
		static assert(isOutputRange!(R, V));
		enum isLastRank = rank == MaxDimension;
		enum isCascadeRank = rank == MaxDimension - 1;
	}


	interface TreeNode(size_t rank, R)
	{
		mixin CommonTreeDeclarations!(rank, MaxDimension, R);

		static if(isLastRank)
		{
			void rangeQuery(V min, V max, R output, size_t cascade);
			size_t getCascadeIndex(V min);
		}
		else
		{
			void rangeQuery(V min, V max, R output);
			static if(isCascadeRank)
			{
				void leftPath(V min, V max, R output, size_t cascade);
				void rightPath(V min, V max, R output, size_t cascade);
				void rangeQueryInCannonicalSubtree(V min, V max, R output, size_t cascade);
			}
			else
			{
				void leftPath(V min, V max, R output);
				void rightPath(V min, V max, R output);
				void rangeQueryInCannonicalSubtree(V min, V max, R output);
			}
		}
	}

	static class InnerNode(size_t rank, R) : TreeNode!(rank, R)
	{
		mixin CommonTreeDeclarations!(rank, MaxDimension, R);

		VElement maxInLeftSubtree;
		TreeNode!(rank, R) left;
		TreeNode!(rank, R) right;

		static if(!isLastRank)
			TreeNode!(rank + 1, R) cannonicalSubtree;
		else
			V[] sortedPoints;

		static if(isCascadeRank)
		{
			size_t[] leftCascade;
			size_t[] rightCascade;
		}

		this(ref V[] data)
		in
		{
			assert(data.length > 1);
			assert(data.isSorted!(lessByDimension!rank));
		}
		out
		{
			static if(!isLastRank)
				assert(data.isSorted!(lessByDimension!(rank + 1)));
			else
				sortedPoints.isSorted!(lessByDimension!rank);
		}
		body
		{
			static if(!isLastRank)
			{
				size_t split = (data.length + 1) / 2;
				maxInLeftSubtree = data[split - 1][rank];
				V[] leftData = data[0 .. split], rightData = data[split .. $];
				left = makeTree!rank(leftData);
				right = makeTree!rank(rightData);
				data = mergeArrays!(lessByDimension!(rank + 1))(leftData, rightData);
				static if(isCascadeRank)
				{
					leftCascade = computeCascade!(lessByDimension!(rank + 1))(data, leftData);
					rightCascade = computeCascade!(lessByDimension!(rank + 1))(data, rightData);
					cannonicalSubtree = new InnerNode!(rank + 1, R)(data);
				}
				else
				{
					V[] copy = data.dup;
					cannonicalSubtree = new InnerNode!(rank + 1, R)(copy);
				}
			}
			else // isLastRank
			{
				sortedPoints = data;
			}
		}

		static if(isLastRank)
		{
			void rangeQuery(V min, V max, R output, size_t cascade)
			in
			{
				if(cascade > 0) assert(lessByDimension!rank(sortedPoints[cascade - 1], min[rank]));
				if(cascade < sortedPoints.lessByDimension) assert(not!(lessByDimension!rank(sortedPoints[cascade], min[rank])));
			}
			body
			{
				put(output, sortedPoints.drop(cascade).until!(not!(lessByDimension!rank))(max));
			}

			size_t getCascadeIndex(V min)
			{
				auto sorted = sortedPoints.assumeSorted!(lessByDimension!rank); // O(1)
				return sorted.lowerBound(min).length;
			}
		}
		else // !isLastRank
		{
			void rangeQuery(V min, V max, R output)
			{
				if(min[rank] > maxInLeftSubtree) right.rangeQuery(min, max, output);
				else if(max[rank] < maxInLeftSubtree) left.rangeQuery(min, max, output);
				else
				{
					static if(!isCascadeRank)
					{
						left.leftPath(min, max, output);
						right.rightPath(min, max, output);
					}
					else // isCascadeRank
					{
						size_t cascade = cannonicalSubtree.getCascadeIndex(min);
						left.leftPath(min, max, output, leftCascade[cascade]);
						right.rightPath(min, max, output, rightCascade[cascade]);
					}
				}
			}
		}

		static if(!isLastRank)
		{
			static if(isCascadeRank)
			{
				void leftPath(V min, V max, R output, size_t cascade)
				{
					if(min[rank] > maxInLeftSubtree) right.leftPath(min, max, output, rightCascade[cascade]);
					else
					{
						right.rangeQueryInCannonicalSubtree(min, max, output, rightCascade[cascade]);
						left.leftPath(min, max, output, leftCascade[cascade]);
					}
				}

				void rightPath(V min, V max, R output, size_t cascade)
				{
					if(max[rank] > maxInLeftSubtree)
					{
						left.rangeQueryInCannonicalSubtree(min, max, output, leftCascade[cascade]);
						right.rightPath(min, max, output, rightCascade[cascade]);
					}
					else left.rightPath(min, max, output, leftCascade[cascade]);
				}

				void rangeQueryInCannonicalSubtree(V min, V max, R output, size_t cascade)
				{
					return cannonicalSubtree.rangeQuery(min, max, output, cascade);
				}
			}
			else
			{
				void leftPath(V min, V max, R output)
				{
					if(min[rank] > maxInLeftSubtree) right.leftPath(min, max, output);
					else
					{
						right.rangeQueryInCannonicalSubtree(min, max, output);
						left.leftPath(min, max, output);
					}
				}

				void rightPath(V min, V max, R output)
				{
					if(max[rank] > maxInLeftSubtree)
					{
						left.rangeQueryInCannonicalSubtree(min, max, output);
						right.rightPath(min, max, output);
					}
					else left.rightPath(min, max, output);
				}

				void rangeQueryInCannonicalSubtree(V min, V max, R output)
				{
					return cannonicalSubtree.rangeQuery(min, max, output);
				}
			}
		}
	}

	static class LeafNode(size_t rank, R) : TreeNode!(rank, R)
	{
		mixin CommonTreeDeclarations!(rank, MaxDimension, R);

		private V data;

		this(V[] data)
		in
		{
			assert(data.length == 1);
		}
		body
		{
			this.data = data[0];
		}

		void rangeQuery(V min, V max, R output)
		{
			foreach(r; rank .. MaxDimension + 1) if(data[r] < min[r] || data[r] > max[r]) return;
			put(output, data);
		}

		static if(isLastRank)
		{
			void rangeQuery(V min, V max, R output, size_t cascade) { rangeQuery(min, max, output); }
			size_t getCascadeIndex(V min) { return 0; }
		}

		static if(!isLastRank)
		{
			static if(!isCascadeRank)
			{
				void leftPath(V min, V max, R output) { rangeQuery(min, max, output); }
				void rightPath(V min, V max, R output) { rangeQuery(min, max, output); }
				void rangeQueryInCannonicalSubtree(V min, V max, R output) { rangeQuery(min, max, output); }
			}
			else
			{
				void leftPath(V min, V max, R output, size_t cascade) { rangeQuery(min, max, output); }
				void rightPath(V min, V max, R output, size_t cascade) { rangeQuery(min, max, output); }
				void rangeQueryInCannonicalSubtree(V min, V max, R output, size_t cascade) { rangeQuery(min, max, output); }
			}
		}
	}

}

int main()
{
	import std.datetime;

	for(size_t size = 100_000; size <= 3_000_000; size += 100_000)
	{
		auto testData = generate!(() => Vec2(uniform(0.0, 10.0), uniform(0.0, 10.0))).take(size).array;
		auto q = new DumbQuery!Vec2(testData);
		RangeTree!Vec2 t;
		auto cons_time = benchmark!(() { t = new RangeTree!Vec2(testData); })(1);

		Vec2 a = Vec2(4.0, 4.0);
		Vec2 b = Vec2(4.1, 4.1);
		Vec2[] resultQ;
		Vec2[] resultT;
		auto time = comparingBenchmark!(() { resultQ = q.rangeQuery(a, b); }, () { resultT = t.rangeQuery(a, b); }, 1)();
		enforce(equal(sort(resultQ), sort(resultT)));
		writeln(size, " ", cons_time[0].length, " ", time.baseTime.length, " ", time.targetTime.length, " ", time.point, " ", resultQ.length);
	}

	writeln("\n3D\n");
	for(size_t size = 10_000; size <= 300_000; size += 10_000)
	{
		auto testData = generate!(() => Vec3(uniform(0.0, 10.0), uniform(0.0, 10.0), uniform(0.0, 10.0))).take(size).array;
		auto q = new DumbQuery!Vec3(testData);
		RangeTree!Vec3 t;
		auto cons_time = benchmark!(() { t = new RangeTree!Vec3(testData); })(1);

		Vec3 a = Vec3(4.0, 4.0, 4.0);
		Vec3 b = Vec3(4.4, 4.4, 4.4);
		Vec3[] resultQ;
		Vec3[] resultT;
		auto time = comparingBenchmark!(() { resultQ = q.rangeQuery(a, b); }, () { resultT = t.rangeQuery(a, b); }, 1)();
		enforce(equal(sort(resultQ), sort(resultT)));
		writeln(size, " ", cons_time[0].length, " ", time.baseTime.length, " ", time.targetTime.length, " ", time.point, " ", resultQ.length);
	}


	return 0;
}