func showLU(a *mat.DenseMatrix) {
	fmt.Printf("\na:\n%v\n", a)
	l, u, p := a.LU()
	fmt.Printf("l:\n%v\n", l)
	fmt.Printf("u:\n%v\n", u)
	fmt.Printf("p:\n%v\n", p)
}

func demo(m *mat.DenseMatrix) {
	fmt.Println("A:")
	fmt.Println(m)
	l, err := m.Cholesky()
	if err != nil {
		fmt.Println(err)
		return
	}
	fmt.Println("L:")
	fmt.Println(l)
}

func MatrixNormal(M, Omega, Sigma *mx.DenseMatrix) func() (X *mx.DenseMatrix) {
	checkMatrixNormal(M, Omega, Sigma)

	Mv := mx.Vectorize(M)
	Cov := mx.Kronecker(Omega, Sigma)
	normal := MVNormal(Mv, Cov)
	return func() (X *mx.DenseMatrix) {
		Xv := normal()
		X = mx.Unvectorize(Xv, M.Rows(), M.Cols())
		return
	}
}

func getNextLocations(piece byte, allMoves *matrix.DenseMatrix, ellipse *matrix.DenseMatrix, baseNode *node.Node) []*node.Node {
	result := make([]*node.Node, 0, 20)

	x := baseNode.GetX()
	y := baseNode.GetY()
	next := baseNode.GetStep() + 1

	// Start with the single move board
	var singleMove *matrix.DenseMatrix
	switch piece {
	case Pawn:
		singleMove = singlePawnMove
	case Rook:
		singleMove = singleRookMove
	case Knight:
		singleMove = singleKnightMove
	case Bishop:
		singleMove = singleBishopMove
	case Queen:
		singleMove = singleQueenMove
	case King:
		singleMove = singleKingMove
	case Puppy:
		singleMove = singlePuppyMove
	}
	singleMove = shiftMatrix(singleMove, x-8, y-8)
	singleMove = singleMove.GetMatrix(7, 0, 8, 8)

	// fmt.Println("(", x, ", ", y, ")")
	// fmt.Println("Single Move: \n", singleMove)

	// fmt.Println("\nAll Moves: \n", allMoves)
	// fmt.Println("\nEllipse: ", ellipse)

	for i := 0; i < 8; i++ {
		for j := 0; j < 8; j++ {
			// fmt.Println("(", i+1, ", ", 8-j, ") : single: ", singleMove.Get(j, i), "ellipse: ", ellipse.Get(j, i), "all: ", allMoves.Get(j, i))

			if singleMove.Get(j, i) != 0 && ellipse.Get(j, i) != 0 && allMoves.Get(j, i) == float64(next) {
				// fmt.Println("New Child Node: (", i+1, ", ", 8-j, ")")

				var newNode *node.Node
				newNode = new(node.Node)
				newNode.SetX(i + 1)
				newNode.SetY(8 - j)
				newNode.SetStep(next)

				baseNode.AddChild(newNode)
				result = append(result, newNode)
			}
		}
	}

	return result
}

func Wishart(n int, V *m.DenseMatrix) func() *m.DenseMatrix {
	p := V.Rows()
	zeros := m.Zeros(p, 1)
	rowGen := MVNormal(zeros, V)
	return func() *m.DenseMatrix {
		x := make([][]float64, n)
		for i := 0; i < n; i++ {
			x[i] = rowGen().Array()
		}
		X := m.MakeDenseMatrixStacked(x)
		S, _ := X.Transpose().TimesDense(X)
		return S
	}
}

//Computes the Euclidean distance between two vectors.
func (KNN *KNNClassifier) ComputeDistance(vector *mat.DenseMatrix, testrow *mat.DenseMatrix) float64 {
	var sum float64

	difference, err := testrow.MinusDense(vector)
	flat := difference.Array()

	if err != nil {
		fmt.Println(err)
	}

	for _, i := range flat {
		squared := math.Pow(i, 2)
		sum += squared
	}

	eucdistance := math.Sqrt(sum)
	return eucdistance
}

func addMovesToBoard(current *matrix.DenseMatrix, newMoves *matrix.DenseMatrix, steps int) *matrix.DenseMatrix {
	// fmt.Println("New Moves: \n", newMoves.String())

	result := current
	for i := 0; i < 15; i++ {
		for j := 0; j < 15; j++ {
			if newMoves.Get(i, j) != 0 && result.Get(i, j) == float64(0) {
				// fmt.Println("    Adding ", steps, " to (", i, ", ", j, ")")
				// if HoleBoard.Get(i, j) > 0 {
				// 	result.Set(i, j, 500)
				// } else {
				result.Set(i, j, float64(steps))
				// }

			}
		}
	}

	return result
}

func calculate_even_fib_until(limit int64) (m []float64) {
	if limit <= 2 { // base case: This cannot be determined by the recurrence relation
		m = []float64{2, 0}
		return
	}
	done := false
	var a, b, c *matrix.DenseMatrix
	a = relation_mat // a holds older b before multiplication
	b = relation_mat // b holds result of matrix multiplication
	c = relation_mat // c holds matrix that is to be multiplied in current iteration
	dict[float64(1)] = relation_mat

	for !done {
		a = b
		b, _ = b.TimesDense(c)
		t, _ := b.TimesDense(init_mat) // calculate resultant vector for a certain power of a matrix
		if t.ColCopy(0)[1] < float64(limit) {
			if pow[1] == -1 { // pow[1] = -1 means we are exponentiating the matrix (first mode).
				pow[0] = 2 * pow[0]
				c = b
				dict[pow[0]] = b
			} else {
				pow[0] = math.Ceil((pow[1] + pow[0]) / float64(2))
				c = dict[math.Ceil((pow[1]-pow[0])/float64(2))]
			}
		} else {
			if pow[1] == -1 {
				pow[1] = pow[0]
				pow[0] = pow[0] / 2
				dict[pow[1]] = b
				if pow[1]-pow[0] < 2 {
					c = dict[float64(1)]
				} else {
					c = dict[math.Ceil((pow[1]-pow[0])/float64(2))]
				}
				b = a
			} else {
				pow[1] = math.Ceil((pow[1] + pow[0]) / float64(2))
				c = dict[math.Ceil((pow[1]-pow[0])/float64(2))]
				b = a
			}
		}
		t, _ = b.TimesDense(init_mat)
		if pow[1] != -1 && t.ColCopy(0)[1] < float64(limit) && t.ColCopy(0)[0] >= float64(limit) {
			done = true
			m = t.ColCopy(0)
		}
	}
	return
}

func Wishart_PDF(n int, V *m.DenseMatrix) func(W *m.DenseMatrix) float64 {
	p := V.Rows()
	Vdet := V.Det()
	Vinv, _ := V.Inverse()
	normalization := pow(2, -0.5*float64(n*p)) *
		pow(Vdet, -0.5*float64(n)) /
		Γ(0.5*float64(n))
	return func(W *m.DenseMatrix) float64 {
		VinvW, _ := Vinv.Times(W)
		return normalization * pow(W.Det(), 0.5*float64(n-p-1)) *
			exp(-0.5*VinvW.Trace())
	}
}

func Wishart_LnPDF(n int, V *m.DenseMatrix) func(W *m.DenseMatrix) float64 {

	p := V.Rows()
	Vdet := V.Det()
	Vinv, _ := V.Inverse()
	normalization := log(2)*(-0.5*float64(n*p)) +
		log(Vdet)*(-0.5*float64(n)) -
		LnΓ(0.5*float64(n))
	return func(W *m.DenseMatrix) float64 {
		VinvW, _ := Vinv.Times(W)
		return normalization +
			log(W.Det())*0.5*float64(n-p-1) -
			0.5*VinvW.Trace()
	}
}

func MatrixT(M, Omega, Sigma *mx.DenseMatrix, n int) func() (T *mx.DenseMatrix) {
	checkMatrixT(M, Omega, Sigma, n)

	fmt.Println("M:", M)
	fmt.Println("Sigma:", Sigma)
	fmt.Println("Omega:", Omega)

	p := M.Rows()
	m := M.Cols()

	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	Sdist := Wishart(n+p-1, OmegaInv)

	Xdist := MatrixNormal(mx.Zeros(p, m), mx.Eye(p), Sigma)

	return func() (T *mx.DenseMatrix) {
		S := Sdist()
		Sinv, err := S.Inverse()
		if err != nil {
			panic(err)
		}
		Sinvc, err := Sinv.Cholesky()
		if err != nil {
			panic(err)
		}
		X := Xdist()
		fmt.Println("Sinvc:", Sinvc)
		fmt.Println("X:", X)
		T, err = Sinvc.Transpose().TimesDense(X)
		if err != nil {
			panic(err)
		}
		err = T.AddDense(M)
		if err != nil {
			panic(err)
		}
		return
	}
}

func (this *KnownVarianceLRPosterior) Remove(x, y *mx.DenseMatrix) {
	xxt, _ := x.TimesDense(x.Transpose())
	this.XXt.Subtract(xxt)
	yxt, _ := y.TimesDense(x.Transpose())
	this.YXt.Subtract(yxt)
}

func (this *KnownVarianceLRPosterior) Insert(x, y *mx.DenseMatrix) {
	xxt, _ := x.TimesDense(x.Transpose())
	this.XXt.Add(xxt)
	yxt, _ := y.TimesDense(x.Transpose())
	this.YXt.Add(yxt)
}

/*
 M is r x c, o x i
 Sigma is r x r, o x o
 Phi is c x c, i x i

 Sigma matches Y o x 1 output dimension
 Phi matches X i x 1 input dimension
*/
func NewKnownVarianceLRPosterior(M, Sigma, Phi *mx.DenseMatrix) (this *KnownVarianceLRPosterior) {
	if M.Rows() != Sigma.Rows() {
		panic("M.Rows != Sigma.Rows")
	}
	if M.Cols() != Phi.Cols() {
		panic("M.Cols != Phi.Cols")
	}
	if Sigma.Rows() != Sigma.Cols() {
		panic("Sigma is not square")
	}
	if Phi.Rows() != Phi.Cols() {
		panic("Phi is not square")
	}
	this = &KnownVarianceLRPosterior{
		M:     M,
		Sigma: Sigma,
		Phi:   Phi,
		XXt:   mx.Zeros(Phi.Cols(), Phi.Cols()),
		YXt:   mx.Zeros(Sigma.Cols(), Phi.Cols()),
	}

	return
}

/*
	If Y ~ N(AX, Sigma, I)
	and A ~ N(M, Sigma, Phi)
	this returns a sampler for P(A|X,Y,Sigma,M,Phi)
*/
func KnownVariancePosterior(Y, X, Sigma, M, Phi *mx.DenseMatrix) func() (A *mx.DenseMatrix) {
	o := Y.Rows()
	i := X.Rows()
	n := Y.Cols()
	if n != X.Cols() {
		panic("X and Y don't have the same number of columns")
	}
	if o != M.Rows() {
		panic("Y.Rows != M.Rows")
	}
	if i != M.Cols() {
		panic("Y.Rows != M.Cols")
	}
	if o != Sigma.Rows() {
		panic("Y.Rows != Sigma.Rows")
	}
	if Sigma.Cols() != Sigma.Rows() {
		panic("Sigma is not square")
	}
	if i != Phi.Rows() {
		panic("X.Rows != Phi.Rows")
	}
	if Phi.Cols() != Phi.Rows() {
		panic("Phi is not square")
	}

	Xt := X.Transpose()

	PhiInv, err := Phi.Inverse()
	if err != nil {
		panic(err)
	}

	XXt, err := X.TimesDense(Xt)
	if err != nil {
		panic(err)
	}

	XXtpPhiInv, err := XXt.PlusDense(PhiInv)
	if err != nil {
		panic(err)
	}

	Omega, err := XXtpPhiInv.Inverse()
	if err != nil {
		panic(err)
	}

	YXtpMPhiInv, err := Y.TimesDense(Xt)
	if err != nil {
		panic(err)
	}

	MPhiInv, err := M.TimesDense(PhiInv)
	if err != nil {
		panic(err)
	}

	err = YXtpMPhiInv.AddDense(MPhiInv)
	if err != nil {
		panic(err)
	}

	Mxy, err := YXtpMPhiInv.TimesDense(Omega)
	if err != nil {
		panic(err)
	}

	if false {
		fmt.Printf("Mxy:\n%v\n", Mxy)
		fmt.Printf("Sigma:\n%v\n", Sigma)
		fmt.Printf("Omega:\n%v\n", Omega)
	}

	return dst.MatrixNormal(Mxy, Sigma, Omega)
}

func MatrixTLnPDF(M, Omega, Sigma *mx.DenseMatrix, n int) func(T *mx.DenseMatrix) (ll float64) {
	checkMatrixT(M, Omega, Sigma, n)

	nf := float64(n)
	p := M.Rows()
	pf := float64(p)
	m := M.Cols()
	mf := float64(m)

	var norm float64 = 0

	norm += logΓpr(p, 0.5*(nf+mf+pf-1), 0.5*(nf+pf-1))
	norm += pow(π, -0.5*mf*pf)
	norm += pow(Omega.Det(), -0.5*mf)
	norm += pow(Sigma.Det(), -0.5*pf)

	SigmaInv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	return func(T *mx.DenseMatrix) (ll float64) {
		ll = norm

		diff, err := T.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := OmegaInv.Copy()
		inner, _ = inner.TimesDense(diff)
		inner, _ = inner.TimesDense(SigmaInv)
		inner, _ = inner.TimesDense(diff.Transpose())

		ll += log(inner.Det()) * -0.5 * (nf + mf + pf - 1)

		return
	}
}

func GenerateMoveBoard(piece byte, x int, y int) *matrix.DenseMatrix {
	var singleMove *matrix.DenseMatrix
	var result *matrix.DenseMatrix

	// Start with the single move board
	switch piece {
	case Pawn:
		singleMove = singlePawnMove
	case Rook:
		singleMove = singleRookMove
	case Knight:
		singleMove = singleKnightMove
	case Bishop:
		singleMove = singleBishopMove
	case Queen:
		singleMove = singleQueenMove
	case King:
		singleMove = singleKingMove
	case Puppy:
		singleMove = singlePuppyMove
	}
	result = shiftMatrix(singleMove, x-8, y-8)

	for i := 0; i < 15; i++ {
		for j := 0; j < 15; j++ {
			if HoleBoard.Get(j, i) > float64(0) {
				result.Set(j, i, 500)
			}
		}
	}

	// fmt.Println(HoleBoard.String())

	// Get the secondary moves
	for n := 1; n < 20; n++ {
		// fmt.Println("Current State: \n", result.String())
		for i := 0; i < 15; i++ {
			for j := 0; j < 15; j++ {
				// Check if the current position needs to generate it's child moves
				if result.Get(i, j) == float64(n) {
					// Shift the single move matrix
					// fmt.Println("Generating moves from position (", j, ", ", i, ") as ", n)
					result = addMovesToBoard(result, shiftMatrix(singleMove, j-7, 15-(i+8)), n+1)
				}
			}
		}
	}

	// Remove the -1
	for i := 0; i < 15; i++ {
		for j := 0; j < 15; j++ {
			// Check if the current position needs to generate it's child moves
			if result.Get(i, j) == float64(-1) {
				result.Set(i, j, 0)
			}
		}
	}

	// fmt.Println(result.String())
	result = result.GetMatrix(7, 0, 8, 8)
	// for i := 0; i < 8; i++ {
	// 	for j := 0; j < 8; j++ {
	// 		if HoleBoard.Get(i, j) == 1 {
	// 			result.Set(j, i, 500)
	// 		}
	// 	}
	// }

	return result
}

func MatrixT_PDF(M, Omega, Sigma *mx.DenseMatrix, n int) func(T *mx.DenseMatrix) (l float64) {
	checkMatrixT(M, Omega, Sigma, n)

	nf := float64(n)
	p := M.Rows()
	pf := float64(p)
	m := M.Cols()
	mf := float64(m)

	var norm float64 = 1

	norm *= GammaPRatio(p, 0.5*(nf+mf+pf-1), 0.5*(nf+pf-1))
	norm *= math.Pow(math.Pi, -0.5*mf*pf)
	norm *= math.Pow(Omega.Det(), -0.5*mf)
	norm *= math.Pow(Sigma.Det(), -0.5*pf)

	SigmaInv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	return func(T *mx.DenseMatrix) (l float64) {
		l = norm

		diff, err := T.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := OmegaInv.Copy()
		inner, _ = inner.TimesDense(diff)
		inner, _ = inner.TimesDense(SigmaInv)
		inner, _ = inner.TimesDense(diff.Transpose())

		l *= math.Pow(inner.Det(), -0.5*(nf+mf+pf-1))

		return
	}
}

func checkMatrixNormal(M, Omega, Sigma *mx.DenseMatrix) {
	p := M.Rows()
	m := M.Cols()
	if Omega.Rows() != p {
		panic(fmt.Sprintf("Omega.Rows != M.Rows, %d != %d", Omega.Rows(), p))
	}
	if Omega.Cols() != Omega.Rows() {
		panic("Omega is not square")
	}
	if Sigma.Rows() != m {
		panic(fmt.Sprintf("Sigma.Cols != M.Cols, %d != %d", Sigma.Cols(), m))
	}
	if Sigma.Cols() != Sigma.Rows() {
		panic("Sigma is not square")
	}
}

func MatrixNormal_LnPDF(M, Omega, Sigma *mx.DenseMatrix) func(A *mx.DenseMatrix) float64 {
	checkMatrixNormal(M, Omega, Sigma)

	pf := float64(M.Rows())
	mf := float64(M.Cols())

	sinv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	oinv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	norm := (2 * math.Pi) * (-0.5 * mf * pf)
	norm += Omega.Det() * (-0.5 * mf)
	norm += Sigma.Det() * (-0.5 * pf)

	return func(X *mx.DenseMatrix) (lp float64) {
		lp = norm
		diff, err := X.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := oinv

		inner, err = inner.TimesDense(diff.Transpose())
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(sinv)
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(diff)
		if err != nil {
			panic(err)
		}

		innerTrace := inner.Trace()

		lp += -0.5 * innerTrace

		return
	}
}