**
** Base matrix (n-dimensional array) class. The 'of' and 'shape' are constant for each Mat.
**
** TODO Support this technique for batched operations, if that turns out to matter for performance:
**
**   Mat a = Mat.batch {
**     b + c.trans * d
**   }
**
mixin Mat {

	**
	** A low-level tool when you want to specify the type and shape manually.
	**
	** Shape to Mat or 'Obj'?
	**
	static Mat createOf(Type of, Int[] shape, Obj?[] values) {
		// TODO Base on array type.
		ObjMat.dupOrGen(of, shape, values)
	}

	**
	** Clears the link to underlying data, and allows the data to go back into a pool
	** for reuse. Could make chunking through many large matrices faster.
	**
	** TODO Profile the effectiveness of pooling, if it's implemented.
	**
	** TODO Implement weak reference automated pooling on platforms that support it?
	**
	Void dispose() {
	}

	**
	** Mat instead of Int[].
	**
	Void eachCoords(|Int[]| handler) {
		Int[] coords := List(Int#, shape.size)
		shape.size.times {
			coords.add(0)
		}
		size.times {
			if (it > 0) {
				coords[-1]++
				// TODO Only go while == condition is true.
				(coords.size - 1 .. 0).each {
					if (coords[it] == shape[it]) {
						coords[it - 1]++
						coords[it] = 0
					}
				}
			}
			handler(coords)
		}
	}

	override Bool equals(Obj? that) {
		if (that isnot Mat) {
			return false
		}
		mat := that as Mat
		if (shape != mat.shape) {
			return false
		}
		if (of != mat.of) {
			// TODO Support of 'of' imperfect so far because of non-dynamic return types.
			// TODO Put this back in once we handle 'of' better. Maybe need explicit type literals.
			// return false
		}
		// Check the actual values.
		// TODO Break the value check out for optimizing in subclasses. This part is potentially slower than other parts of the test.
		result := true
		size.times |i| {
			if (get(i) != mat[i]) {
				result = false
				return
			}
		}
		// All clear.
		return result
	}

	static Mat eye(Int nrows, Int ncols := nrows) {
		// TODO FloatMat here instead of ObjMat.
		return gen(Float#, [nrows, ncols]) {
			it[0] == it[1] ? 1f : 0f
		}
	}

	**
	** TODO Change this to 'Obj' or 'Obj[]'? Should support n dims. Or lists of Mats, etc.
	**
	static Mat from(Obj[][] rows) {
		// TODO Base on the array type.
		return ObjMat(rows)
	}

	**
	** TODO Allow (or require) Mats for shape and for coords.
	**
	static Mat gen(Type of, Int[] shape, |Int[] coords -> Obj?| generator) {
		// TODO Delete this echo once typing is resolved!!!
		// TODO There was supposedly a fix for this. I need to look into it!!
		echo("Mat.gen: $generator")
		// TODO Check type and maybe even shape to determine class.
		return ObjMat.makeGen(of, shape, generator)
	}

	**
	** Pass in an 'Int[]' or 'Mat' for dimensional coords or an 'Int' for direct, undimensional array access.
	**
	abstract Obj? get(Obj coords)

	**
	** May avoid boxing or arrays for fast 'Float' access.
	** Will not convert. Only casts.
	**
	** Use c as -1 to imply r as an index into the full underlying data, ignoring dimensions.
	**
	virtual Float getf(Int r, Int c := Int.minVal) {
		get(c == Int.minVal ? r : [r, c])
	}

	**
	** May avoid boxing or arrays for fast 'Int' access.
	** Will not convert. Only casts.
	**
	** Use c as Int.minVal to imply r as an index into the full underlying data, ignoring dimensions.
	**
	virtual Int geti(Int r, Int c := Int.minVal) {
		get(c == Int.minVal ? r : [r, c])
	}

	**
	** TODO Standardize this! For now, not guaranteed to be the same in the future.
	**
	override Int hash() {
		toStr.hash
	}

	**
	** Convert the coords (Int[] or Mat or Int) to a raw index.
	**
	Int index(Obj coords) {
		index := 0
		if (coords is Int[]) {
			// TODO Method in Mat for sub2ind or whatever.
			// TODO Support negative indexes.
			coordsList := coords as Int[]
			coordsList.each |coord, dim| {
				if (dim > 0) {
					index *= shape[dim - 1]
				}
				index += coord
			}
		} else if (coords is Int) {
			index = coords
		} else {
			throw Err("coords not supported: $coords")
		}
		return index
	}

	**
	** Create a new matrix from this.
	** The matrix passed in is the coordinates.
	** TODO Change Int[] to Mat.
	**
	Mat map(Type mappedOf, |Obj?, Int[] -> Obj?| mapper) {
		// TODO Use the return type to match List#map, once we can control return types.
		// TODO Delete this echo once typing is resolved!!!
		// TODO There was supposedly a fix for this. I need to look into it!!
		echo("Mat.map: $mapper")
		i := 0
		return gen(mappedOf, shape) {
			mapper(this[i++], it)
		}
	}
	
	**
	** Create a new matrix from this.
	** The matrix passed in is the coordinates.
	** TODO Override in FloatMat for performance.
	** TODO Change Int[] to Mat.
	**
	virtual Mat mapf(|Float, Int[] -> Float| mapper) {
		// TODO Delete this echo once typing is resolved!!!
		// TODO There was supposedly a fix for this. I need to look into it!!
		echo("Mat.mapf: $mapper")
		return map(Float#, mapper)
	}

	**
	** Returns a new matrix, leaving this alone.
	**
	Mat max(Int dim) {
		return zeros(1)
	}

	**
	** Returns a new matrix, leaving this alone.
	**
	Mat mean(Int dim) {
		return zeros(1)
	}

	**
	** Returns a new matrix with the min values along dim.
	**
	Mat min(Int dim) {
		reduce(dim, null) |a, b| {a->compare(b) <= 0 ? a : b}
	}

	**
	** Returns a new matrix, leaving this alone.
	**
	Mat mult(Mat other) {
		// TODO Greater than 2 dims?
		if (ndims > 2) {
			throw Err("this.ndims > 2: $ndims")
		}
		if (other.ndims > 2) {
			throw Err("other.ndims > 2: ${other.ndims}")
		}
		if (ncols != other.nrows) {
			throw Err("mismatched inner sizes: $ncols vs. ${other.nrows}")
		}
		// TODO Deal with dims of size 0?
		result := gen(of, [nrows, other.ncols]) |coords| {
			sum := this[[coords[0], 0]]->mult(other[[0, coords[1]]])
			(1 ..< ncols).each |i| {
				sum = sum->plus(this[[coords[0], i]]->mult(other[[i, coords[1]]]))
			}
			return sum
		}
		return result
	}

	Mat multEach(Mat other) {
		// TODO Broadcasting like numpy?
		i := 0
		return map(of) {
			it->mult(other[i++])
		}
	}

	**
	** Returns a new matrix, leaving this alone.
	**
	Mat multFloat(Float float) {
		mapf{float * it}
	}

	Int ncols() {
		size(-1)
	}

	Int ndims() {
		// Or should this be how many dimensions have size > 1? No, unless we reverse order as in MATLAB.
		// Follow NumPy conventions here.
		// Is a 1x1 zero-dimensional? No. All unspecified sizes are implicitly 1.
		// Just trust the user to say how many dims they care about? Yes.
		// TODO Provide a method to collapse empty dimensions (to min of either 1 or 2)?
		shape.size
	}

	Int nrows() {
		size(-2)
	}

	**
	** Unlike List, actually require matching types on set! (Fan might start requiring soon too, though. It's TBD.)
	** If we can enforce type, and maybe size, that would allow consistent optimization for various types and sizes.
	**
	abstract Type of()

	Mat plus(Mat other) {
		i := 0
		return map(of) {
			it->plus(other[i++])
		}
	}

	**
	** Reduces the matrix along the given dimension.
	**
	** Reduces the Mat along one dimension perform item-wise operations.
	**
	** TODO Go back to supporting a reduced-matrix-wise form, too.
	**
	** TODO Should we support reducing multiple dimensions in one call?
	**
	** If init is null, then use the first item along dim.
	** In such cases, the reducer will only be called if size(dim) >= 2.
	**
	** The first 'Mat' arg in the callback is the current item at the given index.
	** The second is the next item to include in the reduction.
	** The third is the current coordinates in the original matrix.
	**
	** Returning a new Mat each time could be wasteful without a Mat pool.
	** Do we need modifiable Mats and just modify the current collapsed one?
	**
	Mat reduce(Int dim, Obj? init, |Obj?, Obj?, Int[] -> Obj?| reducer) {
		newShape := shape.map |v, i| {i == dim ? 1 : v}
		newSize := sizeFor(newShape)
		if (size(dim) > 0) {
			// TODO If init === null, then we need to fill with the first element.
			// TODO If init === null, then we need size(dim) > 1 to bother calling back.
			// TODO Use the reduceMat form for implementing this one???
			size(dim).times {
				// TODO Actually make any of this work.
				next := null
				result := reducer(init, next, Int[,])
			}
		}
		return init
	}

	**
	** Generates a selection for this matrix of size ndims.
	** For example, to create a selection of the first element of one dim, but all of the others, try:
	**
	**   sel{it == dim ? 0 : Step.all)
	**
	** The result can then be used by methods such as 'get' or 'viewSlice'.
	**
	Obj[] sel(|Int -> Obj| gen) {
		Step.mapTimes(Obj#, ndims, gen)
	}

	**
	** Returns the shape specified when creating the Mat as an immutable list.
	** TODO Maybe as an immutable Mat in the future.
	**
	** TODO Support specifying the size of the shape vector. Could prepend 1s or remove leading 1s in many cases.
	**
	abstract Int[] shape()

	**
	** If the dim is Int.minVal, then return the total number of elements.
	** Using Int.minVal allows for avoiding object creation (vs. using null).
	**
	virtual Int size(Int dim := Int.minVal) {
		dim == Int.minVal ? sizeFor(shape) : shape[dim]
	}

	static Int sizeFor(Int[] shape) {
		shape.reduce(1) |Int a, Int b -> Int| {return a * b}
	}

	**
	** TODO Better alignment of columns, ...
	**
	override Str toStr() {
		//
		buf := StrBuf()
		// TODO Normalize newlines?
		// TODO Inline for single row?
		// TODO Only head and tail and size for big matrices?
		buf.add("Mat <| -> $of")
		eachCoords {
			if (it[-1] == 0) {
				if (it.size > 2 && it[-2] == 0) {
					buf.add("\npage ${it[0..-3]}")
				}
				buf.add("\n\t")
			}
			buf.add("${this[it]} ")
		}
		buf.add("\n|>\n")
		return buf.toStr
	}

	**
	** Modifies this. No way!!!
	**
	Mat transpose() {
		return zeros(1)
	}

	**
	** Returns a new matrix, leaving this alone.
	**
	Mat transposed() {
		return zeros(1)
	}

	**
	** Returns a view to this Mat sliced along each dimension according to sel.
	**
	SliceMat viewSlice(Obj sel) {
		//
		SliceMat(this, sel)
	}

	static Mat zeros(Int nrows, Int ncols := nrows) {
		// TODO FloatMat here instead of ObjMat.
		return gen(Float#, [nrows, ncols]) {0f}
	}

}