Golang matrix.FloatVector函数代码示例

func TestDaxpy(t *testing.T) {
	fmt.Printf("* L1 * test axpy: Y = alpha * X + Y\n")
	X := matrix.FloatVector([]float64{1, 1, 1})
	Y := matrix.FloatVector([]float64{0, 0, 0})
	fmt.Printf("before:\nX=\n%v\nY=\n%v\n", X, Y)
	Axpy(X, Y, matrix.FScalar(5.0))
	fmt.Printf("after:\nX=\n%v\nY=\n%v\n", X, Y)
}

func main() {

	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	gdata := [][]float64{
		[]float64{16., 7., 24., -8., 8., -1., 0., -1., 0., 0., 7.,
			-5., 1., -5., 1., -7., 1., -7., -4.},
		[]float64{-14., 2., 7., -13., -18., 3., 0., 0., -1., 0., 3.,
			13., -6., 13., 12., -10., -6., -10., -28.},
		[]float64{5., 0., -15., 12., -6., 17., 0., 0., 0., -1., 9.,
			6., -6., 6., -7., -7., -6., -7., -11.}}

	hdata := []float64{-3., 5., 12., -2., -14., -13., 10., 0., 0., 0., 68.,
		-30., -19., -30., 99., 23., -19., 23., 10.}

	c := matrix.FloatVector([]float64{-6., -4., -5.})
	G := matrix.FloatMatrixFromTable(gdata)
	h := matrix.FloatVector(hdata)

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{2})
	dims.Set("q", []int{4, 4})
	dims.Set("s", []int{3})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.ConeLp(c, G, h, nil, nil, dims, &solopts, nil, nil)
	if err == nil {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("Optimal\n")
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	} else {
		fmt.Printf("status: %s\n", err)
	}
}

func main() {
	flag.Parse()

	gdata0 := [][]float64{
		[]float64{12., 13., 12.},
		[]float64{6., -3., -12.},
		[]float64{-5., -5., 6.}}

	gdata1 := [][]float64{
		[]float64{3., 3., -1., 1.},
		[]float64{-6., -6., -9., 19.},
		[]float64{10., -2., -2., -3.}}

	c := matrix.FloatVector([]float64{-2.0, 1.0, 5.0})
	g0 := matrix.FloatMatrixFromTable(gdata0, matrix.ColumnOrder)
	g1 := matrix.FloatMatrixFromTable(gdata1, matrix.ColumnOrder)
	Ghq := sets.FloatSetNew("Gq", "hq")
	Ghq.Append("Gq", g0, g1)

	h0 := matrix.FloatVector([]float64{-12.0, -3.0, -2.0})
	h1 := matrix.FloatVector([]float64{27.0, 0.0, 3.0, -42.0})
	Ghq.Append("hq", h0, h1)

	var Gl, hl, A, b *matrix.FloatMatrix = nil, nil, nil, nil
	var solopts cvx.SolverOptions
	solopts.MaxIter = 30
	solopts.ShowProgress = true
	if maxIter > -1 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}

	sol, err := cvx.Socp(c, Gl, hl, A, b, Ghq, &solopts, nil, nil)
	fmt.Printf("status: %v\n", err)
	if sol != nil && sol.Status == cvx.Optimal {
		x := sol.Result.At("x")[0]
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		for i, m := range sol.Result.At("sq") {
			fmt.Printf("sq[%d]=\n%v\n", i, m.ToString("%.9f"))
		}
		for i, m := range sol.Result.At("zq") {
			fmt.Printf("zq[%d]=\n%v\n", i, m.ToString("%.9f"))
		}
		sq0 := sol.Result.At("sq")[0]
		sq1 := sol.Result.At("sq")[1]
		zq0 := sol.Result.At("zq")[0]
		zq1 := sol.Result.At("zq")[1]
		check(x, sq0, sq1, zq0, zq1)
	}
}

func TestAcent(t *testing.T) {
	// matrix string in row order presentation
	Adata := [][]float64{
		[]float64{-7.44e-01, 1.11e-01, 1.29e+00, 2.62e+00, -1.82e+00},
		[]float64{4.59e-01, 7.06e-01, 3.16e-01, -1.06e-01, 7.80e-01},
		[]float64{-2.95e-02, -2.22e-01, -2.07e-01, -9.11e-01, -3.92e-01},
		[]float64{-7.75e-01, 1.03e-01, -1.22e+00, -5.74e-01, -3.32e-01},
		[]float64{-1.80e+00, 1.24e+00, -2.61e+00, -9.31e-01, -6.38e-01}}

	bdata := []float64{
		8.38e-01, 9.92e-01, 9.56e-01, 6.14e-01, 6.56e-01,
		3.57e-01, 6.36e-01, 5.08e-01, 8.81e-03, 7.08e-02}

	// these are solution obtained from running cvxopt acent.py with above data
	solData := []float64{-11.59728373909344512, -1.35196389161339936,
		7.21894899350256303, -3.29159917142051528, 4.90454147385329176}

	ntData := []float64{
		1.5163484265903457, 1.2433928210771914, 1.0562922103520955, 0.8816246051011607,
		0.7271128861543598, 0.42725003346248974, 0.0816777301914883, 0.0005458037072843131,
		1.6259980735305693e-10}

	b := matrix.FloatVector(bdata)
	Al := matrix.FloatMatrixFromTable(Adata, matrix.RowOrder)
	Au := matrix.Scale(Al, -1.0)
	A := matrix.FloatZeros(2*Al.Rows(), Al.Cols())
	A.SetSubMatrix(0, 0, Al)
	A.SetSubMatrix(Al.Rows(), 0, Au)

	x, nt, err := acent(A, b, 10)
	if err != nil {
		t.Logf("Acent error: %s", err)
		t.Fail()
	}
	solref := matrix.FloatVector(solData)
	ntref := matrix.FloatVector(ntData)
	soldf := matrix.Minus(x, solref)
	ntdf := matrix.Minus(matrix.FloatVector(nt), ntref)
	solNrm := blas.Nrm2Float(soldf)
	ntNrm := blas.Nrm2Float(ntdf)
	t.Logf("x  [diff=%.2e]:\n%v\n", solNrm, x)
	t.Logf("nt [diff=%.2e]:\n%v\n", ntNrm, nt)

	if solNrm > TOL {
		t.Log("solution deviates too much from expected\n")
		t.Fail()
	}
}

func _TestRankSmall(t *testing.T) {
	bM := 5
	bN := 5
	//bP := 5
	Adata := [][]float64{
		[]float64{1.0, 1.0, 1.0, 1.0, 1.0},
		[]float64{2.0, 2.0, 2.0, 2.0, 2.0},
		[]float64{3.0, 3.0, 3.0, 3.0, 3.0},
		[]float64{4.0, 4.0, 4.0, 4.0, 4.0},
		[]float64{5.0, 5.0, 5.0, 5.0, 5.0}}

	A := matrix.FloatMatrixFromTable(Adata, matrix.RowOrder)
	A0 := matrix.FloatMatrixFromTable(Adata, matrix.RowOrder)
	X := matrix.FloatVector([]float64{1.0, 2.0, 3.0, 4.0, 5.0})
	Y := matrix.FloatWithValue(bN, 1, 2.0)

	Ar := A.FloatArray()
	Xr := X.FloatArray()
	Yr := Y.FloatArray()

	blas.GerFloat(X, Y, A0, 1.0)

	DRankMV(Ar, Xr, Yr, 1.0, A.LeadingIndex(), 1, 1, 0, bN, 0, bM, 4, 4)
	ok := A0.AllClose(A)
	t.Logf("A0 == A1: %v\n", ok)
	if !ok {
		t.Logf("blas ger:\n%v\n", A0)
		t.Logf("A1: \n%v\n", A)
	}
}

// Dscal: X = alpha * X
func TestDscal(t *testing.T) {
	fmt.Printf("* L1 * test scal: X = alpha * X\n")
	alpha := matrix.FScalar(2.0)
	A := matrix.FloatVector([]float64{1.0, 1.0, 1.0, 1.0, 1.0, 1.0})
	Scal(A, alpha)
	fmt.Printf("Dscal 2.0 * A\n")
	fmt.Printf("%s\n", A)
	A = matrix.FloatVector([]float64{1.0, 1.0, 1.0, 1.0, 1.0, 1.0})
	Scal(A, alpha, &linalg.IOpt{"offset", 3})
	fmt.Printf("Dscal 2.0 * A[3:]\n")
	fmt.Printf("%s\n", A)
	A = matrix.FloatVector([]float64{1.0, 1.0, 1.0, 1.0, 1.0, 1.0})
	fmt.Printf("Dscal 2.0* A[::2]\n")
	Scal(A, alpha, &linalg.IOpt{"inc", 2})
	fmt.Printf("%s\n", A)
}

// v = X.T * Y
func TestDdot(t *testing.T) {
	fmt.Printf("* L1 * test dot: X.T*Y\n")
	A := matrix.FloatVector([]float64{1.0, 1.0, 1.0, 1.0, 1.0, 1.0})
	B := matrix.FloatVector([]float64{2.0, 2.0, 2.0, 2.0, 2.0, 2.0})
	v1 := Dot(A, B)
	v2 := Dot(A, B, &linalg.IOpt{"offset", 3})
	v3 := Dot(A, B, &linalg.IOpt{"inc", 2})
	fmt.Printf("Ddot: X.T * Y\n")
	fmt.Printf("%.3f\n", v1.Float())
	fmt.Printf("%.3f\n", v2.Float())
	fmt.Printf("%.3f\n", v3.Float())
	// Output:
	// 12.000
	// 6.000
	// 6.000
}

func (gp *gpConvexProg) F1(x *matrix.FloatMatrix) (f, Df *matrix.FloatMatrix, err error) {
	f = nil
	Df = nil
	err = nil
	f = matrix.FloatZeros(gp.mnl+1, 1)
	Df = matrix.FloatZeros(gp.mnl+1, gp.n)
	y := gp.g.Copy()
	blas.GemvFloat(gp.F, x, y, 1.0, 1.0)

	for i, s := range gp.ind {
		start := s[0]
		stop := s[1]
		// yi := exp(yi) = exp(Fi*x+gi)
		ymax := maxvec(y.FloatArray()[start:stop])
		// ynew = exp(y[start:stop] - ymax)
		ynew := matrix.Exp(matrix.FloatVector(y.FloatArray()[start:stop]).Add(-ymax))
		y.SetIndexesFromArray(ynew.FloatArray(), matrix.Indexes(start, stop)...)

		// fi = log sum yi = log sum exp(Fi*x+gi)
		ysum := blas.AsumFloat(y, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
		f.SetIndex(i, ymax+math.Log(ysum))

		blas.ScalFloat(y, 1.0/ysum, &la.IOpt{"n", stop - start}, &la.IOpt{"offset", start})
		blas.GemvFloat(gp.F, y, Df, 1.0, 0.0, la.OptTrans, &la.IOpt{"m", stop - start},
			&la.IOpt{"incy", gp.mnl + 1}, &la.IOpt{"offseta", start},
			&la.IOpt{"offsetx", start}, &la.IOpt{"offsety", i})
	}
	return
}

func main() {

	var A, b *matrix.FloatMatrix = nil, nil
	m, n := 20, 20

	blas.PanicOnError(true)
	matrix.PanicOnError(true)

	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}
	if len(AVal) > 0 {
		A, _ = matrix.FloatParse(AVal)
		if A == nil {
			fmt.Printf("could not parse:\n%s\n", AVal)
			return
		}
	} else {
		A = matrix.FloatNormal(m, n)
	}
	if len(bVal) > 0 {
		b, _ = matrix.FloatParse(bVal)
		if b == nil {
			fmt.Printf("could not parse:\n%s\n", bVal)
			return
		}
	} else {
		b = matrix.FloatNormal(m, 1)
	}

	sol, err := qcl1(A, b)
	if sol != nil {
		r := sol.Result.At("x")[0]
		x := matrix.FloatVector(r.FloatArray()[:A.Cols()])
		r = sol.Result.At("z")[0]
		z := matrix.FloatVector(r.FloatArray()[r.NumElements()-A.Rows():])
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, z)
	} else {
		fmt.Printf("status: %v\n", err)
	}
}

func TestDgemv(t *testing.T) {
	fmt.Printf("* L2 * test gemv: Y = alpha * A * X + beta * Y\n")
	A := matrix.FloatNew(3, 2, []float64{1, 1, 1, 2, 2, 2})
	X := matrix.FloatVector([]float64{1, 1})
	Y := matrix.FloatVector([]float64{0, 0, 0})
	alpha := matrix.FScalar(1.0)
	beta := matrix.FScalar(0.0)
	fmt.Printf("before: alpha=1.0, beta=0.0\nA=\n%v\nX=\n%v\nY=\n%v\n", A, X, Y)
	err := Gemv(A, X, Y, alpha, beta)
	fmt.Printf("after:\nA=\n%v\nX=\n%v\nY=\n%v\n", A, X, Y)
	fmt.Printf("* L2 * test gemv: X = alpha * A.T * Y + beta * X\n")
	err = Gemv(A, Y, X, alpha, beta, linalg.OptTrans)
	if err != nil {
		fmt.Printf("error: %s\n", err)
	}
	fmt.Printf("after:\nA=\n%v\nX=\n%v\nY=\n%v\n", A, X, Y)
}

func main() {
	flag.Parse()
	if len(spPath) > 0 {
		checkpnt.Reset(spPath)
		checkpnt.Activate()
		checkpnt.Verbose(spVerbose)
		checkpnt.Format("%.17f")
	}

	adata := [][]float64{
		[]float64{0.3, -0.4, -0.2, -0.4, 1.3},
		[]float64{0.6, 1.2, -1.7, 0.3, -0.3},
		[]float64{-0.3, 0.0, 0.6, -1.2, -2.0}}

	A := matrix.FloatMatrixFromTable(adata, matrix.ColumnOrder)
	b := matrix.FloatVector([]float64{1.5, 0.0, -1.2, -0.7, 0.0})

	_, n := A.Size()
	N := n + 1 + n

	h := matrix.FloatZeros(N, 1)
	h.SetIndex(n, 1.0)

	I0 := matrix.FloatDiagonal(n, -1.0)
	I1 := matrix.FloatIdentity(n)
	G, _ := matrix.FloatMatrixStacked(matrix.StackDown, I0, matrix.FloatZeros(1, n), I1)

	At := A.Transpose()
	P := At.Times(A)
	q := At.Times(b).Scale(-1.0)

	dims := sets.NewDimensionSet("l", "q", "s")
	dims.Set("l", []int{n})
	dims.Set("q", []int{n + 1})

	var solopts cvx.SolverOptions
	solopts.MaxIter = 20
	solopts.ShowProgress = true
	if maxIter > 0 {
		solopts.MaxIter = maxIter
	}
	if len(solver) > 0 {
		solopts.KKTSolverName = solver
	}
	sol, err := cvx.ConeQp(P, q, G, h, nil, nil, dims, &solopts, nil)
	if err == nil {
		x := sol.Result.At("x")[0]
		s := sol.Result.At("s")[0]
		z := sol.Result.At("z")[0]
		fmt.Printf("Optimal\n")
		fmt.Printf("x=\n%v\n", x.ToString("%.9f"))
		fmt.Printf("s=\n%v\n", s.ToString("%.9f"))
		fmt.Printf("z=\n%v\n", z.ToString("%.9f"))
		check(x, s, z)
	}

}

func (p *floorPlan) F2(x, z *matrix.FloatMatrix) (f, Df, H *matrix.FloatMatrix, err error) {
	f, Df, err = p.F1(x)
	x17 := matrix.FloatVector(x.FloatArray()[17:])
	tmp := matrix.Div(p.Amin, matrix.Pow(x17, 3.0))
	tmp = matrix.Mul(z, tmp).Scale(2.0)
	diag := matrix.FloatDiagonal(5, tmp.FloatArray()...)
	H = matrix.FloatZeros(22, 22)
	H.SetSubMatrix(17, 17, diag)
	return
}

func main() {
	flag.Parse()

	x := floorplan(matrix.FloatWithValue(5, 1, 100.0))
	if x != nil {
		W := x.GetIndex(0)
		H := x.GetIndex(1)
		xs := matrix.FloatVector(x.FloatArray()[2:7])
		ys := matrix.FloatVector(x.FloatArray()[7:12])
		ws := matrix.FloatVector(x.FloatArray()[12:17])
		hs := matrix.FloatVector(x.FloatArray()[17:])
		fmt.Printf("W = %.5f, H = %.5f\n", W, H)
		fmt.Printf("x = \n%v\n", xs.ToString("%.5f"))
		fmt.Printf("y = \n%v\n", ys.ToString("%.5f"))
		fmt.Printf("w = \n%v\n", ws.ToString("%.5f"))
		fmt.Printf("h = \n%v\n", hs.ToString("%.5f"))
		check(x)
	}
}

func sinv(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, mnl int) (err error) {
	/*DEBUGGED*/

	err = nil

	// For the nonlinear and 'l' blocks:
	//
	//     yk o\ xk = yk .\ xk.

	ind := mnl + dims.At("l")[0]
	blas.Tbsv(y, x, &la_.IOpt{"n", ind}, &la_.IOpt{"k", 0}, &la_.IOpt{"ldA", 1})

	// For the 'q' blocks:
	//
	//                        [ l0   -l1'              ]
	//     yk o\ xk = 1/a^2 * [                        ] * xk
	//                        [ -l1  (a*I + l1*l1')/l0 ]
	//
	// where yk = (l0, l1) and a = l0^2 - l1'*l1.

	for _, m := range dims.At("q") {
		aa := blas.Nrm2Float(y, &la_.IOpt{"n", m - 1}, &la_.IOpt{"offset", ind + 1})
		ee := y.GetIndex(ind)
		aa = (ee + aa) * (ee - aa)
		cc := x.GetIndex(ind)
		dd := blas.DotFloat(x, y, &la_.IOpt{"n", m - 1}, &la_.IOpt{"offsetx", ind + 1},
			&la_.IOpt{"offsety", ind + 1})
		x.SetIndex(ind, cc*ee-dd)
		blas.ScalFloat(x, aa/ee, &la_.IOpt{"n", m - 1}, &la_.IOpt{"offset", ind + 1})
		blas.AxpyFloat(y, x, dd/ee-cc, &la_.IOpt{"n", m - 1},
			&la_.IOpt{"offsetx", ind + 1}, &la_.IOpt{"offsety", ind + 1})
		blas.ScalFloat(x, 1.0/aa, &la_.IOpt{"n", m}, &la_.IOpt{"offset", ind})
		ind += m
	}

	// For the 's' blocks:
	//
	//     yk o\ xk =  xk ./ gamma
	//
	// where gammaij = .5 * (yk_i + yk_j).

	ind2 := ind
	for _, m := range dims.At("s") {
		for j := 0; j < m; j++ {
			u := matrix.FloatVector(y.FloatArray()[ind2+j : ind2+m])
			u.Add(y.GetIndex(ind2 + j))
			u.Scale(0.5)
			blas.Tbsv(u, x, &la_.IOpt{"n", m - j}, &la_.IOpt{"k", 0}, &la_.IOpt{"lda", 1},
				&la_.IOpt{"offsetx", ind + j*(m+1)})
		}
		ind += m * m
		ind2 += m
	}
	return
}

func (p *floorPlan) F1(x *matrix.FloatMatrix) (f, Df *matrix.FloatMatrix, err error) {
	err = nil
	mn := x.Min(-1, -2, -3, -4, -5)
	if mn <= 0.0 {
		f, Df = nil, nil
		return
	}
	zeros := matrix.FloatZeros(5, 12)
	dk1 := matrix.FloatDiagonal(5, -1.0)
	dk2 := matrix.FloatZeros(5, 5)
	x17 := matrix.FloatVector(x.FloatArray()[17:])
	// -( Amin ./ (x17 .* x17) )
	diag := matrix.Div(p.Amin, matrix.Mul(x17, x17)).Scale(-1.0)
	dk2.SetIndexesFromArray(diag.FloatArray(), matrix.MakeDiagonalSet(5)...)
	Df, _ = matrix.FloatMatrixStacked(matrix.StackRight, zeros, dk1, dk2)

	x12 := matrix.FloatVector(x.FloatArray()[12:17])
	// f = -x[12:17] + div(Amin, x[17:]) == div(Amin, x[17:]) - x[12:17]
	f = matrix.Minus(matrix.Div(p.Amin, x17), x12)
	return
}