/**
Functions and types that manipulate multidimensional rectangular arrays.

Copyright: Denis Shelomovskij 2011-2014

License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).

Authors: Denis Shelomovskij

Macros:
	TR = <tr>$0</tr>
	TH = <th>$0</th>
	TD = <td>$0</td>
	TABLE = <table border=1 cellpadding=4 cellspacing=0>$0</table>
 */
module unstd.multidimarray;

import core.exception;
import std.exception;
import std.array;
import unstd.traits;
import unstd.generictuple;
import std.typecons;
import std.conv;
import std..string;
import std.range;
import std.algorithm;

import unstd.lifetime;


private struct R { size_t from, to; }

private template RCount(T...)
{
	static if(T.length)
		enum RCount = is(T[0] == R) + RCount!(T[1 .. $]);
	else
		enum RCount = 0u;
}

private enum isROrSize(T) = RCount!T || isImplicitlyConvertible!(T, const(size_t));

// Note: unittest can't be used as an example for `MultidimArray` member functions
// as it will lead to "forward reference to inferred return type of function call multidimArray".

/**
Implements multidimensional rectangular arrays.

Something like FORTRAN's one.
*/
struct MultidimArray(T, size_t n) if(n >= 1)
{
	/// Dimensions of this array.
	alias dimensions = n;

	private
	{
		size_t[n] _lengths = void;
		size_t[n] _strides = void;
		T[] _data; // _strides[0] * _lengths[0] <= _data.length
	}

	private this(in size_t[n] lengths, in size_t[n] strides, inout T[] data) inout @safe pure nothrow @nogc
	{
		_lengths = lengths;
		_strides = strides;
		_data = data;
	}

	private this(in size_t[n] lengths, inout T[] data) inout @safe pure nothrow
	in
	{
		if(data)
		{
			const count = lengths[].reduce!`a * b`(); // Ignore integral overflow here.
			assert(count == data.length,
				format("Data length for lengths %s must be %s, not %s.", lengths, count, data.length));
		}
	}
	body
	{
		if(n > 1 && lengths[].all!`a`()) // At least two dimensions and any elements.
		{
			size_t count = lengths[0];
			foreach(i; iotaTuple!(1, n)) count *= lengths[i];
			foreach(i; iotaTuple!(1, n)) count /= lengths[i];
			if(count != lengths[0])
				onOutOfMemoryError(); // Not enough address space to index all elements.
		}

		size_t[n] strides = void;
		strides[$ - 1] = 1;
		static if(n > 1)
		{
			strides[$ - 2] = lengths[$ - 1];
			foreach_reverse(i; iotaTuple!(n - 2))
				strides[i] = lengths[i + 1] * strides[i + 1];
		}
		this(lengths, strides, data ? data : new inout T[strides[0] * lengths[0]]);
	}

	this(in size_t[n] lengths...) @safe pure nothrow
	{
		this(lengths, null);
	}

	this(inout T[] data, in size_t[n] lengths...) inout @safe pure nothrow
	in { assert(data, "data must be non-null."); }
	body
	{
		this(lengths, data);
	}

	/// Returns the read only view at its lengths array.
	@property ref const(size_t)[n] lengths() const @safe pure nothrow @nogc
	{
		return _lengths;
	}

	/// Returns the elements count of the array.
	@property size_t elements() const @safe pure nothrow @nogc
	{
		size_t res = _lengths[0];
		foreach(i; iotaTuple!(1, n))
			res *= _lengths[i];
		return res;
	}

	/**
	Returns the maximum number of tail dimensions without pading. Note, that there can be no
	such dimensions.
	*/
	@property size_t packedDimensions() const @safe pure nothrow @nogc
	{
		static if(n > 1)
		{
			size_t packedStride = 1;
			foreach_reverse(i; iotaTuple!n)
			{
				static if(i < n - 1)
					packedStride *= _lengths[i + 1];
				assert(_strides[i] >= packedStride);
				if(_strides[i] > packedStride)
					return n - i - 1;
			}
		}
		return n;
	}

	/// Returns a forward range which has mutable elements and a length for iteration by an element.
	@property byElementForward() @safe pure nothrow @nogc
	{
		static struct Result
		{
		@safe pure nothrow @nogc:

			private
			{
				MultidimArray e; // entity
				size_t rest, shift;
				size_t[n] indices;
			}

			@property
			{
				auto save() inout
				{ return this; }

				size_t length() const
				{ return rest; }

				bool empty() const
				{ return !rest; }

				ref front() inout
				in { assert(!empty, "Trying to call front() on an empty MultidimArray.byElementForward"); }
				body
				{
					return e._data[shift];
				}
			}

			void popFront()
			in { assert(!empty, "Trying to call popFront() on an empty MultidimArray.byElementForward"); }
			body
			{
				--rest;
				foreach_reverse(i; iotaTuple!n)
				{
					shift += e._strides[i];
					if(++indices[i] < e._lengths[i])
					{
						break;
					}
					else
					{
						assert(i || !rest);
						shift -= indices[i] * e._strides[i];
						indices[i] = 0;
					}
				}
			}
		}

		return Result(this, elements, 0);
	}

	/// Returns a finite random-access range which has mutable elements and a length for iteration by an element.
	@property byElementRandomAccess() @safe pure nothrow @nogc
	{
		static struct Result
		{
		@safe pure nothrow /* dmd Issue 13118 @nogc*/:

			private
			{
				MultidimArray e; // entity
				size_t frontIndex, rest, frontShift, backShift;
				size_t[n] frontIndices, afterBackIndices;
			}

			@property
			{
				auto save() inout
				{ return this; }

				size_t length() const
				{ return rest; }

				bool empty() const
				{ return !rest; }

				ref front() inout
				in { assert(!empty, "Trying to call front() on an empty MultidimArray.byElementRandomAccess"); }
				body
				{
					return e._data[frontShift];
				}

				ref back() inout
				in { assert(!empty, "Trying to call back() on an empty MultidimArray.byElementRandomAccess"); }
				body
				{
					return e._data[backShift];
				}
			}

			ref opIndex(in size_t i) inout
			in
			{
				assert(!empty, format("Trying to call opIndex(%s) on an empty MultidimArray.byElementRandomAccess", i));
				assert(i < length, format("Index is out of bounds: trying to call opIndex(%s) on an MultidimArray.byElementRandomAccess with a length = %s", i, length));
			}
			body
			{
				size_t k = i + frontIndex;
				size_t shift = 0;
				foreach_reverse(j; iotaTuple!n)
				{
					// TODO: do % and / in one operation
					const currIndex = k % e._lengths[j];
					k /= e._lengths[j];
					shift += currIndex * e._strides[j];
				}

				return e._data[shift];
			}

			void popFront()
			in { assert(!empty, "Trying to call popFront() on an empty MultidimArray.byElementRandomAccess"); }
			body
			{
				
				++frontIndex;
				--rest;
				foreach_reverse(i; iotaTuple!n)
				{
					frontShift += e._strides[i];
					if(++frontIndices[i] < e._lengths[i])
						break;
					assert(i || !rest);
					frontShift -= frontIndices[i] * e._strides[i];
					frontIndices[i] = 0;
				}
			}

			void popBack()
			in { assert(!empty, "Trying to call popBack() on an empty MultidimArray.byElementRandomAccess"); }
			body
			{
				--rest;
				foreach_reverse(i; iotaTuple!n)
				{
					backShift -= e._strides[i];
					if(--afterBackIndices[i] > 0)
						break;
					assert(i || !rest);
					backShift += e._lengths[i] * e._strides[i];
					afterBackIndices[i] = e._lengths[i];
				}
			}
		}

		const els = elements;
		return Result(this, 0, els, 0, els - 1, 0, _lengths);
	}

	/// Returns a finite random-access range for iteration over the top dimension.
	/// It has mutable elements iff $(D dimensions) is 1.
	@property byTopDimension() @safe pure nothrow @nogc
	{
		static struct Result
		{
		@safe pure nothrow /* dmd Issue 13118 @nogc*/:

			private MultidimArray e; // entity

			@property
			{
				auto save() inout
				{ return this; }

				size_t length() const
				{ return e._lengths[0]; }

				bool empty() const
				{ return !e._lengths[0]; }

				auto ref front() inout
				in { assert(!empty, "Trying to call front() on an empty MultidimArray.byTopDimension"); }
				body
				{
					return opIndex(0);
				}

				auto ref back() inout
				in { assert(!empty, "Trying to call back() on an empty MultidimArray.byTopDimension"); }
				body
				{
					return opIndex(e._lengths[0] - 1);
				}
			}

			static if(n == 1)
			{
				ref opIndex(in size_t i) inout
				in
				{
					assert(!empty, format("Trying to call opIndex(%s) on an empty MultidimArray.byTopDimension", i));
					assert(i < length, format("Index is out of bounds: trying to call opIndex(%s) on an MultidimArray.byTopDimension with a length = %s", i, length));
				}
				body
				{
					return e._data[i];
				}
			}
			else
			{
				auto opIndex(in size_t i) inout @trusted
				in
				{
					assert(!empty, format("Trying to call opIndex(%s) on an empty MultidimArray.byTopDimension", i));
					assert(i < length, format("Index is out of bounds: trying to call opIndex(%s) on an MultidimArray.byTopDimension with a length = %s", i, length));
				}
				body
				{
					const shift = i * e._strides[0];
					return inout MultidimArray!(T, n - 1)(
						e._lengths[1 .. $],
						e._strides[1 .. $],
						e._data[shift .. shift + e._strides[0]]
					);
				}
			}

			void popFront()
			in { assert(!empty, "Trying to call popFront() on an empty MultidimArray.byTopDimension"); }
			body
			{
				e._data = e._data[e._strides[0] .. $];
				--e._lengths[0];
			}

			void popBack()
			in { assert(!empty, "Trying to call popBack() on an empty MultidimArray.byTopDimension"); }
			body
			{
				e._data = e._data[0 .. $ - e._strides[0]];
				--e._lengths[0];
			}
		}

		return Result(this);
	}

	/**
	Returns a forward range which has mutable elements for iteration
	using indices defined by $(D pred) starting from $(D a = 0) and
	incrementing it while indices are in valid range.

	Example:
	----
	auto matrix = multidimArray!char(30, 20);
	matrix[] = ' ';

	foreach(ref el; matrix.byFunction!`a, a`) // fills a diagonal
		el = 'X';

	foreach(ref el; matrix.byFunction!`a^^2 / 5, a`()) // fills a parabola points
		el = 'Y';

	import std.stdio;
	writeln(matrix);
	---
	*/
	@property auto byFunction(string pred)()
	{
		static struct Result
		{
			private
			{
				MultidimArray marr;
				size_t a;
				size_t[n] indices;
			}

			@property @safe pure
			{
				auto save() inout nothrow
				{ return this; }

				bool empty() const nothrow
				{ return !marr.elementExists(indices); }

				ref front() // FIXME `opIndex` isn't `inout nothrow`
				{ return marr[indices]; }
			}

			void popFront()
			{
				++a;
				indices = mixin('[' ~ pred ~ ']');
			}
		}

		auto res = Result(this, -1);
		res.popFront();
		return res;
	}

	/**
	Implements by-element iteration with inidces starting from the top dimension.

	Example:
	----
	auto matrix = multidimArray!int(2, 3, 4);
	foreach(z, y, x, ref el; matrices)
		el = z * 100 + y * 10 + x;
	----
	*/
	int opApply(int delegate(RepeatTuple!(n, size_t), ref T) dg)
	{
		if(!elements)
			return 0;

		RepeatTuple!(n, size_t) indices = 0;
		indices[$ - 1] = -1;

		for(;;)
		{
			foreach_reverse(const plane, ref index; indices)
			{
				if(++index < _lengths[plane])
					break;
				else if(plane)
					index = 0;
				else
					return 0;
			}

			if(const res = dg(indices, _data[getOffset(indices)]))
				return res;
		}
	}

	static if(isAssignable!T)
	{
		/**
		Implements elements initialisation with a $(D value), where $(D value) can be
		of type $(D T) or an input range which $(D front) can be assigned to an element.
		The range should contain exectly $(D elements) elements, otherwise an $(D Exception)
		will be thrown.

		Returns:
		If $(D value) is of type $(D T) or a forward range, returns $(D value).
		Otherwise ($(D value) is an input range but not a forward range) returns $(D void).

		Example:
		----
		auto a23 = multidimArray!int(2, 3);
		auto a46 = multidimArray!int(4, 6);
		auto a234 = multidimArray!int(2, 3, 4);

		a23[] = a234[] = 7;
		a23[] = take(a46[] = a234[] = iota(24), 6);
		----
		*/
		MultidimArray opSliceAssign(T value)
		{
			fill(byElementForward, value);
			return this;
		}

		/// ditto
		MultidimArray opSliceAssign(Range)(Range value)
		if(isInputRange!Range && isAssignable!(T, ElementType!Range))
		in
		{
			static if(hasLength!Range)
				assert(value.length == elements, format(
					"MultidimArray.opSliceAssign(Range): value length (%s) doesn't match array elements count (%s)",
					value.length, elements));
		}
		body
		{
			foreach(ref el; byElementForward)
			{
				assert(!value.empty, format(
					"MultidimArray.opSliceAssign(Range): value doesn't contain enough elements (< %s)",
					elements).assumeWontThrow());
				el = value.front;
				value.popFront();
			}

			assert(value.empty, format(
				"MultidimArray.opSliceAssign(Range): value contains too many elements (> %s)",
				elements).assumeWontThrow());

			return this;
		}

		/// ditto
		MultidimArray opSliceAssign(U)(MultidimArray!(U, n) value)
		if(isAssignable!(T, U))
		in
		{
			assert(value._lengths == _lengths, format(
				"MultidimArray.opSliceAssign(MultidimArray): value lengths %s aren't equal to this lengths %s",
				value._lengths, _lengths));
		}
		body
		{
			return opSliceAssign(value.byElementForward);
		}
	}

	private auto makeCopy(U)()
	{
		auto res = MultidimArray!(U, n)(lengths);
		auto r = byElementForward;
		foreach(ref el; res._data)
		{
			_assumeSafeCopyConstructFrom(el, r.front);
			r.popFront();
		}
		assert(r.empty);
		return res;
	}

	/**
	Support for $(D dup) and $(D idup) properties for MultidimArray.
	*/
	static if(__traits(compiles, { Unqual!T t = rvalueOf!T; }))
	@property dup()
	{
		return makeCopy!(Unqual!T)();
	}

	/// ditto
	static if(is(T == immutable) || __traits(compiles, { immutable t = rvalueOf!T; }))
	@property idup()
	{
		static if(is(T == immutable))
			return this;
		else
			return makeCopy!(immutable T)();
	}


	size_t opDollar(size_t dim)() const @safe pure nothrow @nogc
	{ return lengths[dim]; }

	R opSlice(size_t dim)(in size_t from, in size_t to) const @safe pure nothrow @nogc
	{ return R(from, to); }


	/**
	Inexing/slicing.

	A parameter can be:

$(TABLE
	$(TR $(TH type)           $(TH meaning)         $(TH effect on a resulting dimensions))
	$(TR $(TD $(D n))         $(TD a position)      $(TD -1))
	$(TR $(TD $(D m .. n))    $(TD a range)         $(TD 0))
)
	Examples:
	See $(D MultidimArray) examples.

	Bugs:
	A bit ugly syntax is used because dmd hasn't support for a better one yet (see  $(DBUGZILLA 6798)).
	*/
	ref opIndex(in size_t[n] indices...) inout @safe pure nothrow /* dmd Issue 13118 @nogc*/
	in { assert(elementExists(indices), format("MultidimArray.opIndex(size_t[n]): indices %s are out of bounds (lengths are %s)", indices, _lengths)); }
	body
	{
		return _data[getOffset(indices)];
	}

	/// ditto
	auto opIndex(A...)(A args) inout @trusted pure nothrow /* dmd Issue 13118 @nogc*/
	if(args.length == n && allTuple!(isROrSize, A) && RCount!A)
	in
	{
		string formatOutOfBounds(int i, string idx, string reason) @safe pure nothrow
		{
			return format("MultidimArray.opIndex: Index #%s = %s is out of bounds 0 .. %s (%s)",
				i + 1, idx, _lengths[i], reason).assumeWontThrow();
		}

		foreach(const i, const a; args)
		{
			alias U = Unqual!(A[i]);

			static if(is(U == R))
			{
				assert(a.from <= a.to, format("MultidimArray.opIndex: Index #%s = %s..%s is a range with from > to", i + 1, a.from, a.to));
				assert(a.from >= 0, formatOutOfBounds(i, format("%s..%s", a.from, a.to), "from < 0"));
				assert(a.to <= _lengths[i], formatOutOfBounds(i, format("%s..%s", a.from, a.to), "to > lengths[i]"));
			}
		}
	}
	body
	{
		size_t[n] firstIndices;

		MultidimArray!(T, RCount!A) res;
		static if(RCount!A == n)
			res._strides = _strides;

		foreach(const i, const a; args)
		{
			alias U = Unqual!(A[i]);

			static if(RCount!U)
			{
				enum j = RCount!(A[0 .. i]);
				static if(RCount!A != n)
					res._strides[j] = _strides[i];
			}

			static if(is(U == R))
			{
				firstIndices[i] = a.from;
				res._lengths[j] = a.to - a.from;
			}
			else
			{
				firstIndices[i] = a;
			}
		}

		res._data = cast(T[]) _data[getOffset(firstIndices) .. $]; // TODO $ -> actual bound
		return cast(inout) res;
	}

	/// ditto
	static if(isAssignable!T)
	auto opIndexAssign(U, A...)(U value, A args)
	if(args.length == n && allTuple!(isROrSize, A))
	in
	{
		static if(!RCount!A)
			assert(elementExists(args), format("MultidimArray.opIndexAssign: index out of bounds (lengths: %s, indices: %s)",
				_lengths, [args]));
	}
	body
	{
		static if(RCount!A)
			return this[args][] = value;
		else
			return _data[getOffset(args)] = value;
	}

	/**
	Creates a slice of this entire array with reordered indices. $(D newOrder[i] = n) means that
	$(D i)-th index of a resulting array will behave like $(D n)-th index of the original array.
	Every index sould be used once, otherwise an $(D Exception) will be thrown.

	Example:
	----
	auto matrix3x4 = multidimArray!int(3, 4);
	auto transposed = matrix3x4.reorderIndices(1, 0);
	assert(transposed.lengths == [4, 3]);
	assert(&matrix3x4[2, 3] == &transposed[3, 2]);
	----

	Example:
	----
	auto a = multidimArray!int(2, 3, 4);
	auto b = a.reorderIndices(2, 0, 1);
	assert(b.lengths == [4, 2, 3]);
	assert(&a[1, 2, 3] == &b[3, 1, 2]);
	----
	*/
	auto reorderIndices(in size_t[n] newOrder...) @safe pure nothrow /* dmd Issue 13118 @nogc*/
	in
	{
		bool[n] used = false;
		foreach(const newPlane, const oldPlane; newOrder)
		{
			assert(oldPlane >= 0 || oldPlane < n, format(
				"MultidimArray.reorderIndices: %s isn't a valid index number for a %s-dimnsional array",
				oldPlane, n));

			assert(!used[oldPlane], format(
				"MultidimArray.reorderIndices: Index number %s is used more than once in %s",
				oldPlane, newOrder));

			used[oldPlane] = true;
		}
	}
	body
	{
		MultidimArray res;
		foreach(const newPlane, const oldPlane; newOrder)
		{
			res._strides[newPlane] = _strides[oldPlane];
			res._lengths[newPlane] = _lengths[oldPlane];
		}
		res._data = _data;
		return res;
	}

	static if(n <= 3)
	{
		/**
		Conversion to string function for debugging purposes.

		Implemented for $(D dimensions <= 3).
		*/
		string toString()
		{ return toStringImpl(0); }

		private string toStringImpl(in size_t minElementLength)
		{
			size_t elementLength = minElementLength;

			static if(n == 1)
			{
				string res = "[";
				foreach(i, el; this)
					res ~= to!string(el).rightJustify(elementLength) ~ (i == _lengths[0] - 1 ? "" : ", ");
				res ~= ']';
			}
			else static if(n == 2)
			{
				foreach(el; byElementForward)
					elementLength = max(elementLength, to!string(el).length);

				string res;
				foreach(i; 0 .. _lengths[0])
					res ~= this[i, 0..$].toStringImpl(elementLength) ~ '\n';
			}
			else static if(n == 3)
			{
				foreach(el; byElementForward)
					elementLength = max(elementLength, to!string(el).length);

				string res = "[";
				foreach(i; 0 .. _lengths[0])
					res ~= '\n' ~ this[i, 0..$, 0..$].toStringImpl(elementLength);
				res ~= "]\n";
			}
			else
				static assert(0);

			return res;
		}
	}

private:

	size_t getOffset(in size_t[n] indices...) const @safe pure nothrow @nogc
	in
	{
		foreach(const plane, const index; indices)
			assert(index >= 0 && index <= _lengths[plane]);
	}
	body
	{
		size_t res = 0;
		foreach(const plane; iotaTuple!n)
			res += indices[plane] * _strides[plane];
		return res;
	}

	bool elementExists(in size_t[n] indices...) const @safe pure nothrow @nogc
	{
		foreach(const plane; iotaTuple!n)
			if(indices[plane] < 0 || indices[plane] >= _lengths[plane])
				return false;
		return true;
	}
}

///
unittest
{
	// Let's creates an GC allocated three-dimensional rectangular array from 2 matrices 3x4
	auto matrices = multidimArray!int(2, 3, 4); // matrices has a type MultidimArray!(int, 3)

	// Setting an element at intersection of the first column and
	// the third row of the secon matrix to seven:
	matrices[1, 0, 2] = 7;

	// Filling the whole first column of the secon matrix with sixes:
	matrices[1, 0, 0..$][] = 6;

	// Filling the whole array with fives:
	matrices[] = 5;

	// Iterating the array
	foreach(z, y, x, ref el; matrices) // using opApply
		el = cast(int) (z * 100 + y * 10 + x);

	int c = 0;
	foreach(ref el; matrices.byElementForward)
		el = c++;

	c = 0;
	foreach(i; 0 .. matrices.elements)
		matrices.byElementRandomAccess[i] = c++;

	c = 0;
	foreach(matrix; matrices.byTopDimension)       // for each of two matrices
		foreach(row; matrix.byTopDimension)        // for each row
			foreach(ref el; row.byTopDimension) // for each element
				el = c++;

	c = 0;
	foreach_reverse(ref el; matrices.byElementRandomAccess)
		el = c++;

	c = 0;
	foreach_reverse(i; 0 .. matrices.elements)
		matrices.byElementRandomAccess[i] = c++;

	// Inexing/slicing
	// * use <integer> to select a position
	// * use <integer> .. <integer> to select a range
	// E.g. use 0..$  to select the whole range
	matrices = matrices[0..$, 0..$, 0..$];  // the entire array, same as [0..2, 0..3, 0..4]
	auto array2d = matrices[0, 0..$, 0..$]; // the first matrix
	auto array1d = matrices[0, 1, 0..$];  // the second row of the first matrix
	array1d = matrices[0, 0..$, 1];       // the second column of the first matrix
	matrices[0, 1, 1] = 9;                // setting an element at a crossing of the row an the column

	// first two rows and three columns of the secon matrix
	array2d = matrices[1, 0 .. 2, 0 .. 3];
}

/**
Convenience function that returns an $(D MultidimArray!(T, n)) object.

Returns:
The first overload returns a $(D MultidimArray) with a newly allocated data
Others use an existing storage.

Params:
data = A memory storage for a resulting array of type $(D T[]).
array = An array to wrap. It can be a multidimensional static array or a
slice of it (has a dynamic top dimension).
lengths = Lengths of a resulting array.

Template_parameters:
$(D T) Element type of a resulting array. Should be explicitly defined only
for the first overload which has no memory storage.

$(D n) Dimensions of a resulting array. Can be explicitly defined to use only
first $(D n) of $(D array) dimensions.

$(D A) Type of a wrapping $(D array). It is inferred from the $(D array) argument
and should not be explicitly defined.

See_Also: MultidimArray

Throws: The first overload throws an $(D RangeError) in $(D debug) build if $(D data) length isn't equal to $(D lengths) prouct.
*/
// #1: allocate new
auto multidimArray(T, size_t n)(size_t[n] lengths...) @safe pure nothrow
if(n > 0)
{
	return MultidimArray!(T, n)(lengths);
}

/// ditto
// #2: use existing storage
auto multidimArray(size_t n, T)(T[] data, size_t[n] lengths...) @safe pure nothrow
if(n > 0)
{
	return MultidimArray!(T, n)(data, lengths);
}

/// ditto
// #3: use some dimensions of an existing static array
auto multidimArray(size_t n, A)(ref A array) pure nothrow
if(n > 0 && n <= staticArrayDims!A)
{
	alias ElementType = MultidimStaticArrayElementType!(A, n);
	return multidimArray!(n, ElementType)(cast(ElementType[]) array, multidimStaticArrayLengths!(A, n));
}

/// ditto
// #4: use all dimensions of an existing static array
auto multidimArray(A)(ref A array) pure nothrow
if(isStaticArray!A)
{
	return multidimArray!(staticArrayDims!A)(array);
}

/// ditto
// #5: use some dimensions of an existing dynamic array of static arrays
auto multidimArray(size_t n, A)(A array) @safe pure nothrow
if(isDynamicArray!A && n > 0 && n - 1 <= staticArrayDims!(ElementType!A))
{
	alias U = MultidimStaticArrayElementType!(ElementType!A, n - 1);
	return multidimArray!(n, U)(cast(U[]) array, array.length, multidimStaticArrayLengths!(ElementType!A, n - 1));
}

/// ditto
// #6: use all dimensions of an existing dynamic array of static arrays
auto multidimArray(A)(A array) @safe pure nothrow
if(isDynamicArray!A)
{
	return multidimArray!(1 + staticArrayDims!(ElementType!A))(array);
}

///
pure nothrow unittest
{
	// Let's create an GC allocated three-dimensional rectangular array from 2 matrices 3x4
	auto matrix1 = multidimArray!int(2, 3, 4);

	// Let's create the same array using an existing storage
	auto darr2 = new int[24]; // At least 24 elements are needed
	auto matrix2 = multidimArray(darr2, 2, 3, 4); // No need for explicit element type declaration

	// Let's create the same array using an existing static array as data storage
	int[4][3][2] sarr3; // or in a postfix form: int sarr[2][3][4];
	auto matrix3 = multidimArray(sarr3); // No need for any explicit template declarations

	// The head array can be dynamic
	int[4][3][] darr3 = sarr3[];
	auto matrix31 = multidimArray(darr3); // Works like previous one

	// Let's create an array of static arrays
	ubyte[4][4][3][2] sarr4; // a postfix form: ubyte[4] sarr[2][3][4];
	auto matrix4 = multidimArray!3(sarr4); // Only 3 major of 4 dimensions are indeces

	// The head array can also be dynamic
	auto matrix41 = multidimArray!3(sarr4[]); // Works like previous one
}

pure nothrow unittest // multidimArray
{
	void test234matrix(T)(ref T matrix)
	{
		static assert(isForwardRange!(typeof(matrix.byElementForward)));
		static assert(hasAssignableElements!(typeof(matrix.byElementForward)));
		static assert(hasLength!(typeof(matrix.byElementForward)));

		static assert(isBidirectionalRange!(typeof(matrix.byElementRandomAccess)));
		static assert(isRandomAccessRange!(typeof(matrix.byElementRandomAccess)));
		static assert(hasAssignableElements!(typeof(matrix.byElementRandomAccess)));
		static assert(hasLength!(typeof(matrix.byElementRandomAccess)));

		static assert(isBidirectionalRange!(typeof(matrix.byTopDimension)));
		static assert(isRandomAccessRange!(typeof(matrix.byTopDimension)));
		static assert(hasAssignableElements!(typeof(matrix.byTopDimension)) == (T.dimensions == 1));
		static assert(hasLength!(typeof(matrix.byTopDimension)));
		with(matrix)
		{
			assert(_data.length == 24);
			static assert(lengths.length == 3);
			assert(lengths == [2, 3, 4]);
		}
	}

	auto matrix1 = multidimArray!int(2, 3, 4); // #1
	test234matrix(matrix1);

	auto darr2 = new int[24];
	auto matrix2 = multidimArray(darr2, 2, 3, 4); // #2
	test234matrix(matrix2);

	int[4][3][2] sarr3;
	auto matrix3 = multidimArray(sarr3); // #4
	test234matrix(matrix3);

	int[4][3][] darr3 = sarr3[];
	auto matrix31 = multidimArray(darr3); // #6
	test234matrix(matrix31);

	ubyte[4][4][3][2] sarr4;
	auto matrix4 = multidimArray!3(sarr4); // #3
	test234matrix(matrix4);

	auto matrix41 = multidimArray!3(sarr4[]); // #5
	test234matrix(matrix41);
}

@safe pure nothrow unittest // MultidimArray properties: dimensions, lengths, elements, packedDimensions
{
	auto marr3 = multidimArray!int(3, 4, 5);

	with(marr3)
	{
		static assert(dimensions == 3);
		assert(lengths == [3, 4, 5]);
		assert(elements == 60);
		assert(packedDimensions == 3);
	}

	auto marr3s = marr3[0..$, 0..$, 0..$];
	with(marr3)
	{
		static assert(dimensions == 3);
		assert(lengths == [3, 4, 5]);
		assert(elements == 60);
		assert(packedDimensions == 3);
	}

	marr3s = marr3[0..2, 0..$, 0..$];
	with(marr3s)
	{
		static assert(dimensions == 3);
		assert(lengths == [2, 4, 5]);
		assert(elements == 40);
		assert(packedDimensions == 3);
	}

	marr3s = marr3[0..1, 0..$, 0..$];
	with(marr3s)
	{
		static assert(dimensions == 3);
		assert(lengths == [1, 4, 5]);
		assert(elements == 20);
		assert(packedDimensions == 3);
	}

	foreach(i; 0 .. 4)
	{
		marr3s = marr3[0..$, 0..i, 0..$];
		with(marr3s)
		{
			static assert(dimensions == 3);
			assert(lengths == [3, i, 5]);
			assert(elements == 15 * i);
			assert(packedDimensions == 2);
		}
	}

	foreach(i; 0 .. 5)
	{
		marr3s = marr3[0..$, 0..$, 0..i];
		with(marr3s)
		{
			static assert(dimensions == 3);
			assert(lengths == [3, 4, i]);
			assert(elements == 12 * i);
			assert(packedDimensions == 1);
		}
	}

	auto marr2 = marr3[1, 0..$, 0..$];
	with(marr2)
	{
		static assert(dimensions == 2);
		assert(lengths == [4, 5]);
		assert(elements == 20);
		assert(packedDimensions == 2);
	}

	marr2 = marr3[0..$, 1, 0..$];
	with(marr2)
	{
		static assert(dimensions == 2);
		assert(lengths == [3, 5]);
		assert(elements == 15);
		assert(packedDimensions == 1);
	}

	marr2 = marr3[0..$, 0..$, 1];
	with(marr2)
	{
		static assert(dimensions == 2);
		assert(lengths == [3, 4]);
		assert(elements == 12);
		assert(packedDimensions == 0);
	}
}

unittest // MultidimArray iterations: byElementForward, byElementRandomAccess, byTopDimension, opApply
{
	void test(T)()
	{
		T sarr;
		auto darr = cast(int[])sarr;
		auto matrix = multidimArray(sarr);

		darr[] = -1;
		int c;
		foreach(ref el; matrix.byElementForward)
		{
			el = c;
			assert(darr[c] == c && el == c && &el == &darr[c]);
			++c;
		}

		darr[] = -1;
		c = 0;
		foreach(ref el; matrix.byElementRandomAccess)
		{
			el = c;
			assert(darr[c] == c && el == c && &el == &darr[c]);
			++c;
		}

		darr[] = -1;
		c = cast(int) matrix.elements;
		foreach_reverse(ref el; matrix.byElementRandomAccess)
		{
			--c;
			el = c;
			assert(darr[c] == c && el == c && &el == &darr[c]);
		}

		darr[] = -1;
		foreach(i; 0 .. matrix.elements)
		{
			auto ptr = &matrix.byElementRandomAccess[i];
			*ptr = cast(int) i;
			assert(darr[i] == i && ptr == &darr[i]);
		}

		darr[] = -1;
		{
			auto r = matrix.byElementRandomAccess;

			foreach(i; 0 .. matrix.elements)
			{
				assert(&r[matrix.elements - i - 1] == &darr[$ - 1]);
				auto ptr = &r[0];
				// &r.front() instead of &r.front because of bad property syntax
				assert(ptr == &darr[i] && ptr == &r.front());
				r.popFront();
			}
		}

		darr[] = -1;
		c = 0;
		static if(matrix.dimensions == 3)
			foreach(x, y, z, ref el; matrix)
			{
				el = c;
				assert(c == x * matrix._strides[0] + y * matrix._strides[1] + z * matrix._strides[2]);
				assert(darr[c] == c && el == c && &el == &darr[c]);
				++c;
			}
		else static if(matrix.dimensions == 2)
			foreach(x, y, ref el; matrix)
			{
				el = c;
				assert(c == x * matrix._strides[0] + y * matrix._strides[1]);
				assert(darr[c] == c && el == c && &el == &darr[c]);
				++c;
			}

		darr[] = -1;
		c = 0;
		static if(matrix.dimensions == 3)
			foreach(subMatrix; matrix.byTopDimension)
				foreach(row; subMatrix.byTopDimension)
					foreach(ref el; row.byTopDimension)
					{
						el = c;
						assert(darr[c] == c && el == c && &el == &darr[c]);
						++c;
					}
		else static if(matrix.dimensions == 2)
			foreach(row; matrix.byTopDimension)
				foreach(ref el; row.byTopDimension)
				{
					el = c;
					assert(darr[c] == c && el == c && &el == &darr[c]);
					++c;
				}
		else static if(matrix.dimensions == 1)
			foreach(ref el; matrix.byTopDimension)
			{
				el = c;
				assert(darr[c] == c && el == c && &el == &darr[c]);
				++c;
			}

		/*darr[] = -1;
		c = 0;
		foreach(size_t[matrix.dimensions] indices, ref el; matrix)
		{
			el = c;
			assert(darr[c] == c && el == c && &el == &darr[c]);
			++c;
		}*/
	}
	test!(int[0])();
	test!(int[2])();
	test!(int[2][0])();
	test!(int[0][2])();
	test!(int[1][1])();
	test!(int[7][3])();
	test!(int[4][3][2])();
	test!(int[1][2][3][4])();
	test!(int[7][3][7][3])();
}

version(unittest) // MultidimArray unittest helper functions
{
	bool isSame(MArr1, MArr2)(MArr1 marr1, MArr2 marr2) pure nothrow
	in { assert(marr1.lengths == marr2.lengths); }
	body
	{
		return &marr1[0, 0, 0] == &marr2[0, 0, 0];
	}

	bool isCopy(MArr1, MArr2)(MArr1 marr1, MArr2 marr2)
	{
		return !isSame(marr1, marr2) && equal(marr1.byElementForward, marr2.byElementForward);
	}

	bool equalRange(MArr, Range)(MArr marr, Range r)
	{
		return equal(marr.byElementForward, r);
	}
}

pure nothrow unittest // MultidimArray copying: opSliceAssign, dup, idup
{
	alias repeat = std.range.repeat; // std.string.repeat will be removed in February 2012

	auto a23 = multidimArray!int(2, 3);
	auto a46 = multidimArray!int(4, 6);
	auto a234 = multidimArray!int(2, 3, 4);
	assert(equalRange(a23, repeat(0, 6)));
	assert(equalRange(a46, repeat(0, 24)));
	assert(equalRange(a234, repeat(0, 24)));

	a23[] = 7;
	assert(equalRange(a23, repeat(7, 6)));
	assert(equalRange(a23[] = iota(6), iota(6)));

	a234[] = iota(24);
	assert(equalRange(a234, iota(24)));
	assert((a23[] = a234[0..$, 0..$, 1]) is a23);
	assert(equalRange(a23, [1, 5, 9, 13, 17, 21]));

	auto b234 = a234.dup;
	assert(isCopy(a234, b234));

	b234[] = -1;
	assert(equalRange(a234, iota(24)));

	b234 = a234;
	assert(isSame(a234, b234));

	b234 = multidimArray!int(2, 3, 4);
	b234[] = a234;
	assert(isCopy(a234, b234));

	auto ia234 = a234.idup;
	static assert(is(typeof(ia234) == MultidimArray!(immutable int, 3u)));
	static assert(!__traits(compiles, (ia234[] = 7)));
	static assert(!__traits(compiles, (ia234[] = new int[24])));
	assert(isCopy(a234, ia234));
	assert(isSame(ia234, ia234.idup));
	assert(isCopy(ia234, ia234.dup));

	const(int)[4][3][2] carr;
	auto ca234 = multidimArray(carr);
	static assert(is(typeof(ca234) == MultidimArray!(const int, 3u)));
	static assert(!__traits(compiles, (ca234[] = 7)));
	static assert(!__traits(compiles, (ca234[] = new int[24])));
	assert(isCopy(ca234, ca234.idup));
	assert(isCopy(ca234, ca234.dup));


	b234[] = ia234;
	assert(isCopy(a234, b234));


	auto a123 = multidimArray!int(1, 2, 3);
	assert(equalRange(a123, repeat(0, 6)));

	a123[] = 7;
	assert(equalRange(a123, repeat(7, 6)));

	a123[] = [8, 8, 8,    8, 8, 8];
	assert(equalRange(a123, repeat(8, 6)));

	a123[] = repeat(9, 6);
	assert(equalRange(a123, repeat(9, 6)));

	auto matrix2 = multidimArray!(int[])(1, 2, 3);
	matrix2[] = [9];
}

@safe pure nothrow unittest // MultidimArray duplication: dup, idup
{
	static struct S { int i; }
	auto s1 = multidimArray!S(1);
	++s1[0].i;
	assert(s1.dup[0].i == 1);
	assert(s1.idup[0].i == 1);
}

@safe pure nothrow unittest // MultidimArray reordering: reorderIndices
{
	auto matrix3x4 = multidimArray!int(3, 4);
	auto transposed = matrix3x4.reorderIndices(1, 0);
	assert(transposed.lengths == [4, 3]);
	assert(&matrix3x4[2, 3] == &transposed[3, 2]);

	auto a = multidimArray!int(2, 3, 4);
	auto b = a.reorderIndices(2, 0, 1);
	assert(b.lengths == [4, 2, 3]);
	assert(&a[1, 2, 3] == &b[3, 1, 2]);
}

@safe pure nothrow unittest
{
	auto matrix = multidimArray!int(3, 4);

	foreach(ref el; matrix.byFunction!`a + 1, a + 1`)
		el = 1;
	foreach(ref el; matrix.byFunction!`a, (a + 1)^^2 - 1`())
		el = 2;

	assert(equalRange(matrix, [
		2, 0, 0, 0,
		0, 1, 0, 2,
		0, 0, 1, 0]));
}

// TODO unittests for: opIndexAssign, opIndex