Golang Matrix.Dims方法代码示例

//This is a temporal function. It returns the determinant of a 3x3 matrix. Panics if the matrix is not 3x3.
//It is also defined in the chem package which is not-so-clean.
func det(A mat64.Matrix) float64 {
	r, c := A.Dims()
	if r != 3 || c != 3 {
		panic(ErrDeterminant)
	}
	return (A.At(0, 0)*(A.At(1, 1)*A.At(2, 2)-A.At(2, 1)*A.At(1, 2)) - A.At(1, 0)*(A.At(0, 1)*A.At(2, 2)-A.At(2, 1)*A.At(0, 2)) + A.At(2, 0)*(A.At(0, 1)*A.At(1, 2)-A.At(1, 1)*A.At(0, 2)))
}

// MeanBatch predicts the mean at the set of locations specified by x. Stores in-place into yPred
// If yPred is nil new memory is allocated.
func (g *GP) MeanBatch(yPred []float64, x mat64.Matrix) []float64 {
	rx, cx := x.Dims()
	if cx != g.inputDim {
		panic(badInputLength)
	}
	if yPred == nil {
		yPred = make([]float64, rx)
	}
	ry := len(yPred)
	if rx != ry {
		panic(badOutputLength)
	}
	nSamples, _ := g.inputs.Dims()

	covariance := mat64.NewDense(nSamples, rx, nil)
	row := make([]float64, g.inputDim)
	for j := 0; j < rx; j++ {
		for k := 0; k < g.inputDim; k++ {
			row[k] = x.At(j, k)
		}
		for i := 0; i < nSamples; i++ {
			v := g.kernel.Distance(g.inputs.RawRowView(i), row)
			covariance.Set(i, j, v)
		}
	}
	yPredVec := mat64.NewVector(len(yPred), yPred)
	yPredVec.MulVec(covariance.T(), g.sigInvY)
	// Rescale the outputs
	for i, v := range yPred {
		yPred[i] = v*g.std + g.mean
	}
	return yPred
}

func columnMean(M mat.Matrix) mat.Matrix {
	r, c := M.Dims()

	SumMatrix := columnSum(M)

	switch t := SumMatrix.(type) {
	case *mat.Dense:
		M := mat.NewDense(1, c, nil)
		M.Scale(1/float64(r), SumMatrix)
		return M
	case mat.Mutable:
		_ = t
		V := SumMatrix.(mat.Mutable)
		_, cols := V.Dims()

		for i := 0; i < cols; i++ {
			V.Set(0, i, SumMatrix.At(0, i)/float64(r))
		}

		return V
	default:
		panic("M is of an unknown type")
	}

}

// findLinearlyIndependnt finds a set of linearly independent columns of A, and
// returns the column indexes of the linearly independent columns.
func findLinearlyIndependent(A mat64.Matrix) []int {
	m, n := A.Dims()
	idxs := make([]int, 0, m)
	columns := mat64.NewDense(m, m, nil)
	newCol := make([]float64, m)
	// Walk in reverse order because slack variables are typically the last columns
	// of A.
	for i := n - 1; i >= 0; i-- {
		if len(idxs) == m {
			break
		}
		mat64.Col(newCol, i, A)
		if len(idxs) == 0 {
			// A column is linearly independent from the null set.
			// This is what needs to be changed if zero columns are allowed, as
			// a column of all zeros is not linearly independent from itself.
			columns.SetCol(len(idxs), newCol)
			idxs = append(idxs, i)
			continue
		}
		if linearlyDependent(mat64.NewVector(m, newCol), columns.View(0, 0, m, len(idxs))) {
			continue
		}
		columns.SetCol(len(idxs), newCol)
		idxs = append(idxs, i)
	}
	return idxs
}

// Cov returns the covariance between a set of data points based on the current
// GP fit.
func (g *GP) Cov(m *mat64.SymDense, x mat64.Matrix) *mat64.SymDense {
	if m != nil {
		// TODO(btracey): Make this k**
		panic("resuing m not coded")
	}
	// The joint covariance matrix is
	// K(x_*, k_*) - k(x_*, x) k(x,x)^-1 k(x, x*)
	nSamp, nDim := x.Dims()
	if nDim != g.inputDim {
		panic(badInputLength)
	}

	// Compute K(x_*, x) K(x, x)^-1 K(x, x_*)
	kstar := g.formKStar(x)
	var tmp mat64.Dense
	tmp.SolveCholesky(g.cholK, kstar)
	var tmp2 mat64.Dense
	tmp2.Mul(kstar.T(), &tmp)

	// Compute k(x_*, x_*) and perform the subtraction.
	kstarstar := mat64.NewSymDense(nSamp, nil)
	for i := 0; i < nSamp; i++ {
		for j := i; j < nSamp; j++ {
			v := g.kernel.Distance(mat64.Row(nil, i, x), mat64.Row(nil, j, x))
			if i == j {
				v += g.noise
			}
			kstarstar.SetSym(i, j, v-tmp2.At(i, j))
		}
	}
	return kstarstar
}

func benchmarkCovarianceMatrixInPlace(b *testing.B, m mat64.Matrix) {
	_, c := m.Dims()
	res := mat64.NewDense(c, c, nil)
	b.ResetTimer()
	for i := 0; i < b.N; i++ {
		CovarianceMatrix(res, m, nil)
	}
}

func MatrixToImage(src mat64.Matrix) image.Image {
	width, height := src.Dims()
	img := image.NewRGBA(image.Rect(0, 0, width, height))
	for x := 0; x < width; x++ {
		for y := 0; y < height; y++ {
			img.Set(x, y, Float64ToColor(src.At(x, y)))
		}
	}
	return img
}

// linearlyDependent returns whether the vector is linearly dependent
// with the columns of A. It assumes that A is a full-rank matrix.
func linearlyDependent(vec *mat64.Vector, A mat64.Matrix) bool {
	// Add vec to the columns of A, and see if the condition number is reasonable.
	m, n := A.Dims()
	aNew := mat64.NewDense(m, n+1, nil)
	aNew.Copy(A)
	col := mat64.Col(nil, 0, vec)
	aNew.SetCol(n, col)
	cond := mat64.Cond(aNew, 1)
	return cond > 1e12
}

func rowSum(M mat.Matrix) mat.Matrix {
	rows, _ := M.Dims()

	floatRes := make([]float64, rows)
	for i := 0; i < rows; i++ {
		floatRes[i] = mat.Sum(getRowVector(i, M))
	}

	return mat.NewDense(rows, 1, floatRes)
}

func columnSum(M mat.Matrix) mat.Matrix {
	_, cols := M.Dims()

	floatRes := make([]float64, cols)
	for i := 0; i < cols; i++ {
		floatRes[i] = mat.Sum(getColumnVector(i, M))
	}

	return mat.NewDense(1, cols, floatRes)
}

// extractColumns creates a new matrix out of the columns of A specified by cols.
// TODO(btracey): Allow this to take a receiver.
func extractColumns(A mat64.Matrix, cols []int) *mat64.Dense {
	r, _ := A.Dims()
	sub := mat64.NewDense(r, len(cols), nil)
	col := make([]float64, r)
	for j, idx := range cols {
		mat64.Col(col, idx, A)
		sub.SetCol(j, col)
	}
	return sub
}

func verifyInputs(initialBasic []int, c []float64, A mat64.Matrix, b []float64) error {
	m, n := A.Dims()
	if len(c) != n {
		panic("lp: c vector incorrect length")
	}
	if len(b) != m {
		panic("lp: b vector incorrect length")
	}
	if len(c) != n {
		panic("lp: c vector incorrect length")
	}
	if len(initialBasic) != 0 && len(initialBasic) != m {
		panic("lp: initialBasic incorrect length")
	}

	// Do some sanity checks so that ab does not become singular during the
	// simplex solution. If the ZeroRow checks are removed then the code for
	// finding a set of linearly indepent columns must be improved.

	// Check that if a row of A only has zero elements that corresponding
	// element in b is zero, otherwise the problem is infeasible.
	// Otherwise return ErrZeroRow.
	for i := 0; i < m; i++ {
		isZero := true
		for j := 0; j < n; j++ {
			if A.At(i, j) != 0 {
				isZero = false
				break
			}
		}
		if isZero && b[i] != 0 {
			// Infeasible
			return ErrInfeasible
		} else if isZero {
			return ErrZeroRow
		}
	}
	// Check that if a column only has zero elements that the respective C vector
	// is positive (otherwise unbounded). Otherwise return ErrZeroColumn.
	for j := 0; j < n; j++ {
		isZero := true
		for i := 0; i < m; i++ {
			if A.At(i, j) != 0 {
				isZero = false
				break
			}
		}
		if isZero && c[j] < 0 {
			return ErrUnbounded
		} else if isZero {
			return ErrZeroColumn
		}
	}
	return nil
}

// AssignCentroid assigns all of the examples in X to one of the groups
// in Mu
// X -> (m*n), Mu -> (K*n)
// returns (m*1)
func AssignCentroid(X, Mu mat.Matrix) *mat.Vector {
	m, _ := X.Dims()
	idx := mat.NewVector(m, nil)

	for i := 0; i < m; i++ {
		x := getRowVector(i, X)
		idx.SetVec(i, float64(NearestCentroid(x, Mu)))
	}

	return idx
}

func getRowVector(index int, M mat.Matrix) *mat.Vector {
	_, cols := M.Dims()
	var rowData []float64

	if cols == 0 {
		rowData = []float64{}
	} else {
		rowData = mat.Row(nil, index, M)
	}
	return mat.NewVector(cols, rowData)
}

func J(idx *mat.Vector, X, Mu mat.Matrix) float64 {
	Mux := ConstructXCentroidMatrix(idx, Mu)
	xRows, xCols := X.Dims()

	Diff := mat.NewDense(xRows, xCols, nil)
	Diff.Sub(X, Mux)
	Diff.MulElem(Diff, Diff)
	Diff = rowSum(Diff).(*mat.Dense)

	return columnSum(Diff).At(0, 0) / float64(xRows)
}