Golang DimensionSet.SumPacked方法代碼示例

/*
   Copy x to y using packed storage.

   The vector x is an element of S, with the 's' components stored in
   unpacked storage.  On return, x is copied to y with the 's' components
   stored in packed storage and the off-diagonal entries scaled by
   sqrt(2).
*/
func pack(x, y *matrix.FloatMatrix, dims *sets.DimensionSet, opts ...la_.Option) (err error) {
	/*DEBUGGED*/
	err = nil
	mnl := la_.GetIntOpt("mnl", 0, opts...)
	offsetx := la_.GetIntOpt("offsetx", 0, opts...)
	offsety := la_.GetIntOpt("offsety", 0, opts...)

	nlq := mnl + dims.At("l")[0] + dims.Sum("q")
	blas.Copy(x, y, &la_.IOpt{"n", nlq}, &la_.IOpt{"offsetx", offsetx},
		&la_.IOpt{"offsety", offsety})

	iu, ip := offsetx+nlq, offsety+nlq
	for _, n := range dims.At("s") {
		for k := 0; k < n; k++ {
			blas.Copy(x, y, &la_.IOpt{"n", n - k}, &la_.IOpt{"offsetx", iu + k*(n+1)},
				&la_.IOpt{"offsety", ip})
			y.SetIndex(ip, (y.GetIndex(ip) / math.Sqrt(2.0)))
			ip += n - k
		}
		iu += n * n
	}
	np := dims.SumPacked("s")
	blas.ScalFloat(y, math.Sqrt(2.0), &la_.IOpt{"n", np}, &la_.IOpt{"offset", offsety + nlq})
	return
}

// Solves a pair of primal and dual cone programs  using custom KKT solver.
//
func ConeLpCustomKKT(c, G, h, A, b *matrix.FloatMatrix, dims *sets.DimensionSet,
	kktsolver KKTConeSolver, solopts *SolverOptions, primalstart,
	dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	if c == nil || c.Cols() > 1 {
		err = errors.New("'c' must be matrix with 1 column")
		return
	}
	if h == nil {
		h = matrix.FloatZeros(0, 1)
	}
	if h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if dims == nil {
		dims = sets.NewDimensionSet("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}
	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")
	//cdim_diag := dims.Sum("l", "q", "s")

	if G == nil {
		G = matrix.FloatZeros(0, c.Rows())
	}
	if !G.SizeMatch(cdim, c.Rows()) {
		estr := fmt.Sprintf("'G' must be of size (%d,%d)", cdim, c.Rows())
		err = errors.New(estr)
		return
	}

	// Check A and set defaults if it is nil
	if A == nil {
		// zeros rows reduces Gemv to vector products
		A = matrix.FloatZeros(0, c.Rows())
	}
	if A.Cols() != c.Rows() {
		estr := fmt.Sprintf("'A' must have %d columns", c.Rows())
		err = errors.New(estr)
		return
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}
	if b.Rows() != A.Rows() {
		estr := fmt.Sprintf("'b' must have length %d", A.Rows())
		err = errors.New(estr)
		return
	}

	if b.Rows() > c.Rows() || b.Rows()+cdim_pckd < c.Rows() {
		err = errors.New("Rank(A) < p or Rank([G; A]) < n")
		return
	}

	mA := &matrixVarA{A}
	mG := &matrixVarG{G, dims}
	mc := &matrixVar{c}
	mb := &matrixVar{b}

	return conelp_problem(mc, mG, h, mA, mb, dims, kktsolver, solopts, primalstart, dualstart)
}

// Solves a pair of primal and dual cone programs using custom KKT solver and constraint
// interfaces MatrixG and MatrixA
//
func ConeLpCustomMatrix(c *matrix.FloatMatrix, G MatrixG, h *matrix.FloatMatrix,
	A MatrixA, b *matrix.FloatMatrix, dims *sets.DimensionSet, kktsolver KKTConeSolver,
	solopts *SolverOptions, primalstart, dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	err = nil

	if c == nil || c.Cols() > 1 {
		err = errors.New("'c' must be matrix with 1 column")
		return
	}
	if h == nil || h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if err = checkConeLpDimensions(dims); err != nil {
		return
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")
	//cdim_diag := dims.Sum("l", "q", "s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	// Data for kth 'q' constraint are found in rows indq[k]:indq[k+1] of G.
	indq := make([]int, 0)
	indq = append(indq, dims.At("l")[0])
	for _, k := range dims.At("q") {
		indq = append(indq, indq[len(indq)-1]+k)
	}

	// Data for kth 's' constraint are found in rows inds[k]:inds[k+1] of G.
	inds := make([]int, 0)
	inds = append(inds, indq[len(indq)-1])
	for _, k := range dims.At("s") {
		inds = append(inds, inds[len(inds)-1]+k*k)
	}

	// Check b and set defaults if it is nil
	if b == nil {
		b = matrix.FloatZeros(0, 1)
	}
	if b.Cols() != 1 {
		estr := fmt.Sprintf("'b' must be a matrix with 1 column")
		err = errors.New(estr)
		return
	}
	if b.Rows() > c.Rows() || b.Rows()+cdim_pckd < c.Rows() {
		err = errors.New("Rank(A) < p or Rank([G; A]) < n")
		return
	}

	if kktsolver == nil {
		err = errors.New("nil kktsolver not allowed.")
		return
	}

	var mA MatrixVarA
	var mG MatrixVarG
	if G == nil {
		mG = &matrixVarG{matrix.FloatZeros(0, c.Rows()), dims}
	} else {
		mG = &matrixIfG{G}
	}
	if A == nil {
		mA = &matrixVarA{matrix.FloatZeros(0, c.Rows())}
	} else {
		mA = &matrixIfA{A}
	}
	var mc = &matrixVar{c}
	var mb = &matrixVar{b}

	return conelp_problem(mc, mG, h, mA, mb, dims, kktsolver, solopts, primalstart, dualstart)
}

// Solves a pair of primal and dual cone programs
//
//        minimize    c'*x
//        subject to  G*x + s = h
//                    A*x = b
//                    s >= 0
//
//        maximize    -h'*z - b'*y
//        subject to  G'*z + A'*y + c = 0
//                    z >= 0.
//
// The inequalities are with respect to a cone C defined as the Cartesian
// product of N + M + 1 cones:
//
//        C = C_0 x C_1 x .... x C_N x C_{N+1} x ... x C_{N+M}.
//
// The first cone C_0 is the nonnegative orthant of dimension ml.
// The next N cones are second order cones of dimension r[0], ..., r[N-1].
// The second order cone of dimension m is defined as
//
//        { (u0, u1) in R x R^{m-1} | u0 >= ||u1||_2 }.
//
// The next M cones are positive semidefinite cones of order t[0], ..., t[M-1] >= 0.
//
// The structure of C is specified by DimensionSet dims which holds following sets
//
//   dims.At("l")  l, the dimension of the nonnegative orthant (array of length 1)
//   dims.At("q")  r[0], ... r[N-1], list with the dimesions of the second-order cones
//   dims.At("s")  t[0], ... t[M-1], array with the dimensions of the positive
//                 semidefinite cones
//
// The default value for dims is l: []int{G.Rows()}, q: []int{}, s: []int{}.
//
// Arguments primalstart, dualstart are optional starting points for primal and
// dual problems. If non-nil then primalstart is a FloatMatrixSet having two entries.
//
//  primalstart.At("x")[0]  starting point for x
//  primalstart.At("s")[0]  starting point for s
//  dualstart.At("y")[0]    starting point for y
//  dualstart.At("z")[0]    starting point for z
//
// On exit Solution contains the result and information about the accurancy of the
// solution. if SolutionStatus is Optimal then Solution.Result contains solutions
// for the problems.
//
//   Result.At("x")[0]  solution for x
//   Result.At("y")[0]  solution for y
//   Result.At("s")[0]  solution for s
//   Result.At("z")[0]  solution for z
//
func ConeLp(c, G, h, A, b *matrix.FloatMatrix, dims *sets.DimensionSet, solopts *SolverOptions,
	primalstart, dualstart *sets.FloatMatrixSet) (sol *Solution, err error) {

	if c == nil || c.Cols() > 1 {
		err = errors.New("'c' must be matrix with 1 column")
		return
	}
	if c.Rows() < 1 {
		err = errors.New("No variables, 'c' must have at least one row")
		return

	}
	if h == nil || h.Cols() > 1 {
		err = errors.New("'h' must be matrix with 1 column")
		return
	}

	if dims == nil {
		dims = sets.NewDimensionSet("l", "q", "s")
		dims.Set("l", []int{h.Rows()})
	}

	cdim := dims.Sum("l", "q") + dims.SumSquared("s")
	cdim_pckd := dims.Sum("l", "q") + dims.SumPacked("s")

	if h.Rows() != cdim {
		err = errors.New(fmt.Sprintf("'h' must be float matrix of size (%d,1)", cdim))
		return
	}

	if G == nil {
		G = matrix.FloatZeros(0, c.Rows())
	}
	if !G.SizeMatch(cdim, c.Rows()) {
		estr := fmt.Sprintf("'G' must be of size (%d,%d)", cdim, c.Rows())
		err = errors.New(estr)
		return
	}

	// Check A and set defaults if it is nil
	if A == nil {
		// zeros rows reduces Gemv to vector products
		A = matrix.FloatZeros(0, c.Rows())
	}
	if A.Cols() != c.Rows() {
		estr := fmt.Sprintf("'A' must have %d columns", c.Rows())
		err = errors.New(estr)
		return
	}

//.........這裏部分代碼省略.........