func Mean(matrix [][]float64, axis int) (result []float64, ok bool) {
	var dim int
	sum, ok := Sum(matrix, axis)
	if !ok {
		return
	}
	result = make([]float64, len(sum))
	if axis == 0 {
		dim = len(matrix)
	} else {
		dim = len(matrix[0])
	}
	for ii, value := range sum {
		result[ii] = value / float64(dim)
	}
	return
}
Mean of matrix along with axis
axis=0 => row-wise
axis=1 => column-wise
package matrix
import (
	"math/cmplx"
	"sort"
)
func column(matrix [][]float64, columnNumber int) (result []float64) {
	result = make([]float64, len(matrix))
	for ii := 0; ii < len(matrix); ii++ {
		for jj := 0; jj < len(matrix[0]); jj++ {
			result[ii] = matrix[ii][columnNumber]
		}
	}
	return
}
func Add(firstMatrix [][]float64, secondMatrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(firstMatrix))
	for ii := 0; ii < len(firstMatrix); ii++ {
		result[ii] = make([]float64, len(firstMatrix[0]))
		for jj := 0; jj < len(firstMatrix[0]); jj++ {
			result[ii][jj] = firstMatrix[ii][jj] + secondMatrix[ii][jj]
		}
	}
	return
}
func Subtract(firstMatrix [][]float64, secondMatrix [][]float64) [][]float64 {
	secondMatrix = MultiplyByScalar(secondMatrix, -1.0)
	return Add(firstMatrix, secondMatrix)
}
func DotProduct(firstMatrix [][]float64, secondMatrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(firstMatrix))
	for ii := 0; ii < len(firstMatrix); ii++ {
		result[ii] = make([]float64, len(secondMatrix[0]))
		for jj := 0; jj < len(secondMatrix[0]); jj++ {
			for kk := 0; kk < len(secondMatrix); kk++ {
				result[ii][jj] += firstMatrix[ii][kk] * secondMatrix[kk][jj]
			}
		}
	}
	return
}
func Determinant(matrix [][]float64) (result float64) {
	matrixLength := len(matrix)
	sums := make([]float64, matrixLength*2)
	for ii := 0; ii < len(sums); ii++ {
		sums[ii] = 1
	}
	for ii := 0; ii < matrixLength; ii++ {
		for jj := 0; jj < matrixLength; jj++ {
			if ii-jj < 0 {
				sums[matrixLength+ii-jj] *= matrix[ii][jj]
			} else {
				sums[ii-jj] *= matrix[ii][jj]
			}
			if ii+jj >= matrixLength {
				sums[ii+jj] *= matrix[ii][jj]
			} else {
				sums[ii+jj+matrixLength] *= matrix[ii][jj]
			}
		}
	}
	dim := matrixLength * 2
	if matrixLength == 2 {
		dim = 2
		matrixLength = 1
	}
	for ii := 0; ii < dim; ii++ {
		if ii >= matrixLength {
			result -= sums[ii]
		} else {
			result += sums[ii]
		}
	}
	return
}
func MinorMatrix(matrix [][]float64) (result [][]float64) {
	var (
		matrixLength int
	)
	matrixLength = len(matrix)
	result = make([][]float64, matrixLength)
	for ii := 0; ii < matrixLength; ii++ {
		result[ii] = make([]float64, matrixLength)
		for jj := 0; jj < matrixLength; jj++ {
			auxM := [][]float64{}
			for iik := 0; iik < matrixLength; iik++ {
				if iik != ii {
					auxR := []float64{}
					for jjk := 0; jjk < matrixLength; jjk++ {
						if jjk != jj {
							auxR = append(auxR, matrix[iik][jjk])
						}
					}
					auxM = append(auxM, auxR)
				}
			}
			result[ii][jj] = Determinant(auxM)
		}
	}
	return
}
func CofactorMatrix(matrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(matrix))
	for ii := 0; ii < len(matrix); ii++ {
		result[ii] = make([]float64, len(matrix[0]))
		for jj := 0; jj < len(matrix[0]); jj++ {
			if (ii+jj)%2 == 0 {
				result[ii][jj] = matrix[ii][jj]
			} else {
				result[ii][jj] = -matrix[ii][jj]
			}
		}
	}
	return
}
func Inverse(matrix [][]float64) [][]float64 {
	determinant := Determinant(matrix)
	adj := Transpose(CofactorMatrix(MinorMatrix(matrix)))
	return MultiplyByScalar(adj, 1./determinant)
}
func Divide(firstMatrix [][]float64, secondMatrix [][]float64) [][]float64 {
	return DotProduct(firstMatrix, Inverse(secondMatrix))
}
func MultiplyByScalar(matrix [][]float64, scalar float64) [][]float64 {
	function := func(x float64) float64 {
		return x * scalar
	}
	return Apply(matrix, function)
}
func MultTranspose(firstMatrix [][]float64, secondMatrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(firstMatrix))
	for ii := 0; ii < len(firstMatrix); ii++ {
		result[ii] = make([]float64, len(secondMatrix))
		for jj := 0; jj < len(secondMatrix); jj++ {
			for kk := 0; kk < len(secondMatrix[0]); kk++ {
				result[ii][jj] += firstMatrix[ii][kk] * secondMatrix[jj][kk]
			}
		}
	}
	return
}
func Multiply(firstMatrix [][]float64, secondMatrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(firstMatrix))
	for ii := 0; ii < len(firstMatrix); ii++ {
		result[ii] = make([]float64, len(firstMatrix[0]))
		for jj := 0; jj < len(firstMatrix[0]); jj++ {
			result[ii][jj] = firstMatrix[ii][jj] * secondMatrix[ii][jj]
		}
	}
	return
}
func Transpose(matrix [][]float64) (result [][]float64) {
	result = make([][]float64, len(matrix[0]))
	for ii := 0; ii < len(matrix[0]); ii++ {
		result[ii] = make([]float64, len(matrix))
		for jj := 0; jj < len(matrix); jj++ {
			result[ii][jj] = matrix[jj][ii]
		}
	}
	return
}
func Sum(matrix [][]float64, axis int) (result []float64, ok bool) {
	rowSum := make([]float64, len(matrix))
	columnSum := make([]float64, len(matrix[0]))
	for ii := 0; ii < len(matrix); ii++ {
		for jj := 0; jj < len(matrix[0]); jj++ {
			rowSum[ii] += matrix[ii][jj]
			columnSum[jj] += matrix[ii][jj]
		}
	}
	if axis == 1 {
		result = rowSum
		ok = true
	} else {
		result = columnSum
		ok = true
	}
	return
}
func Max(matrix [][]float64, axis int) (result []float64, ok bool) {
	var (
		firstDimension  int
		secondDimension int
	)
	firstDimension = len(matrix)
	secondDimension = len(matrix[0])
	rowMax := make([]float64, firstDimension)
	columnMax := make([]float64, secondDimension)
	deepMatrix := deepCopyMatrix(matrix)
	for ii := 0; ii < firstDimension; ii++ {
		sort.Float64s(deepMatrix[ii])
		rowMax[ii] = deepMatrix[ii][secondDimension-1]
		if ii == firstDimension-1 {
			for jj := 0; jj < secondDimension; jj++ {
				deepCol := deepCopyArray(column(matrix, jj))
				sort.Float64s(deepCol)
				columnMax[jj] = deepCol[len(deepCol)-1]
			}
		}
	}
	if axis == 1 {
		result = rowMax
		ok = true
	} else {
		result = columnMax
		ok = true
	}
	return
}
func Min(matrix [][]float64, axis int) (result []float64, ok bool) {
	var (
		firstDimension  int
		secondDimension int
	)
	firstDimension = len(matrix)
	secondDimension = len(matrix[0])
	rowMin := make([]float64, firstDimension)
	columnMin := make([]float64, secondDimension)
	deepMatrix := deepCopyMatrix(matrix)
	for ii := 0; ii < firstDimension; ii++ {
		sort.Float64s(deepMatrix[ii])
		rowMin[ii] = deepMatrix[ii][0]
		if ii == firstDimension-1 {
			for jj := 0; jj < secondDimension; jj++ {
				deepCol := deepCopyArray(column(matrix, jj))
				sort.Float64s(deepCol)
				columnMin[jj] = deepCol[0]
			}
		}
	}
	if axis == 1 {
		result = rowMin
		ok = true
	} else {
		result = columnMin
		ok = true
	}
	return
}
func Median(matrix [][]float64, axis int) (result []float64, ok bool) {
	var (
		firstDimension       int
		secondDimension      int
		halfRow              int
		halfCol              int
		isFirstDivisibleBy2  bool
		isSecondDivisibleBy2 bool
	)
	firstDimension = len(matrix)
	secondDimension = len(matrix[0])
	if firstDimension%2 == 0 {
		halfRow = firstDimension / 2
		isFirstDivisibleBy2 = true
	} else {
		halfRow = (firstDimension - 1) / 2
	}
	if secondDimension%2 == 0 {
		halfCol = secondDimension / 2
		isSecondDivisibleBy2 = true
	} else {
		halfCol = (secondDimension - 1) / 2
	}
	rowMedian := make([]float64, firstDimension)
	columnMedian := make([]float64, secondDimension)
	deepMatrix := deepCopyMatrix(matrix)
	for ii := 0; ii < firstDimension; ii++ {
		sort.Float64s(deepMatrix[ii])
		if isSecondDivisibleBy2 && secondDimension > halfRow+1 {
			rowMedian[ii] = (deepMatrix[ii][halfRow] + deepMatrix[ii][halfRow+1]) / 2.
		} else if secondDimension == 2 {
			rowMedian[ii] = (deepMatrix[ii][0] + deepMatrix[ii][1]) / 2.
		} else {
			rowMedian[ii] = deepMatrix[ii][halfRow]
		}
		if ii == firstDimension-1 {
			for jj := 0; jj < secondDimension; jj++ {
				deepCol := deepCopyArray(column(matrix, jj))
				sort.Float64s(deepCol)
				if isFirstDivisibleBy2 && len(deepCol) > halfCol+1 {
					columnMedian[jj] = (deepCol[halfCol] + deepCol[halfCol+1]) / 2.
				} else if len(deepCol) == 2 {
					columnMedian[jj] = (deepCol[0] + deepCol[1]) / 2.
				} else {
					columnMedian[jj] = deepCol[halfCol]
				}
			}
		}
	}
	if axis == 1 {
		result = rowMedian
		ok = true
	} else {
		result = columnMedian
		ok = true
	}
	return
}