How to use the cvxopt.base.matrix function in cvxopt
To help you get started, we’ve selected a few cvxopt examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
math1um / objects-invariants-properties / gt.py View on Github 
if n == 1:
return 1.0
#the definition of Xrow assumes that the vertices are integers from 0 to n-1, so we relabel the graph
gc.relabel()
d = m + n
c = -1 * cvxopt.base.matrix([0.0]*(n-1) + [2.0]*(d-n))
Xrow = [i*(1+n) for i in xrange(n-1)] + [b+a*n for (a, b) in gc.edge_iterator(labels=False)]
Xcol = range(n-1) + range(d-1)[n-1:]
X = cvxopt.base.spmatrix(1.0, Xrow, Xcol, (n*n, d-1))
for i in xrange(n-1):
X[n*n-1, i] = -1.0
sol = cvxopt.solvers.sdp(c, Gs=[-X], hs=[-cvxopt.base.matrix([0.0]*(n*n-1) + [-1.0], (n,n))])
v = 1.0 + cvxopt.base.matrix(-c, (1, d-1)) * sol['x']
# TODO: Rounding here is a total hack, sometimes it can come in slightly
# under the analytical answer, for example, 2.999998 instead of 3, which
# screws up the floor() call when checking difficult graphs.
return round(v[0], 3)def _update_w(self):
def updatesingleW(i):
# optimize alpha using qp solver from cvxopt
FA = base.matrix(np.float64(np.dot(-self.H, self.data[i,:].T)))
al = solvers.qp(HA, FA, INQa, INQb)
self.W[i,:] = np.array(al['x']).reshape((1,-1))
# float64 required for cvxopt
HA = base.matrix(np.float64(np.dot(self.H, self.H.T)))
INQa = base.matrix(-np.eye(self._num_bases))
INQb = base.matrix(0.0, (self._num_bases,1))
map(updatesingleW, xrange(self._data_dimension))
self.W = self.W/np.sum(self.W, axis=1)def update_w(self):
""" alternating least squares step, update W under the convexity
constraint """
def update_single_w(i):
""" compute single W[:,i] """
# optimize beta using qp solver from cvxopt
FB = base.matrix(np.float64(np.dot(-self.data.T, W_hat[:,i])))
be = solvers.qp(HB, FB, INQa, INQb, EQa, EQb)
self.beta[i,:] = np.array(be['x']).reshape((1, self._num_samples))
# float64 required for cvxopt
HB = base.matrix(np.float64(np.dot(self.data[:,:].T, self.data[:,:])))
EQb = base.matrix(1.0, (1, 1))
W_hat = np.dot(self.data, pinv(self.H))
INQa = base.matrix(-np.eye(self._num_samples))
INQb = base.matrix(0.0, (self._num_samples, 1))
EQa = base.matrix(1.0, (1, self._num_samples))
for i in xrange(self._num_bases):
update_single_w(i)
self.W = np.dot(self.beta, self.data.T).Tcdim_pckd = dims['l'] + sum(dims['q']) + sum([ int(k*(k+1)/2) for k in
dims['s'] ])
# A' = [Q1, Q2] * [R1; 0]
if type(A) is matrix:
QA = +A.T
else:
QA = matrix(A.T)
tauA = matrix(0.0, (p,1))
lapack.geqrf(QA, tauA)
Gs = matrix(0.0, (cdim, n))
tauG = matrix(0.0, (n-p,1))
u = matrix(0.0, (cdim_pckd, 1))
vv = matrix(0.0, (n,1))
w = matrix(0.0, (cdim_pckd, 1))
def factor(W):
# Gs = W^{-T}*G, in packed storage.
Gs[:,:] = G
scale(Gs, W, trans = 'T', inverse = 'I')
pack2(Gs, dims)
# Gs := [ Gs1, Gs2 ]
# = Gs * [ Q1, Q2 ]
lapack.ormqr(QA, tauA, Gs, side = 'R', m = cdim_pckd)
# QR factorization Gs2 := [ Q3, Q4 ] * [ R3; 0 ]
lapack.geqrf(Gs, tauG, n = n-p, m = cdim_pckd, offsetA =
Gs.size[0]*p)def _update_h(self):
""" alternating least squares step, update H enforcing a convexity
constraint.
"""
def update_single_h(i):
""" compute single H[:,i] """
# optimize alpha using qp solver from cvxopt
FA = base.matrix(np.float64(np.dot(-self.W.T, self.data[:,i])))
al = solvers.qp(HA, FA, INQa, INQb, EQa, EQb)
self.H[:,i] = np.array(al['x']).reshape((1, self._num_bases))
EQb = base.matrix(1.0, (1,1))
# float64 required for cvxopt
HA = base.matrix(np.float64(np.dot(self.W.T, self.W)))
INQa = base.matrix(-np.eye(self._num_bases))
INQb = base.matrix(0.0, (self._num_bases,1))
EQa = base.matrix(1.0, (1, self._num_bases))
for i in xrange(self._num_samples):
update_single_h(i)array(blx),
array(bux))
task.putobjsense(msk.objsense.minimize)
task.optimize()
task.solutionsummary (msk.streamtype.msg);
prosta, solsta = list(), list()
task.getsolutionstatus(msk.soltype.bas, prosta, solsta)
x, z = zeros(n, Float), zeros(m, Float)
task.getsolutionslice(msk.soltype.bas, msk.solitem.xx, 0, n, x)
task.getsolutionslice(msk.soltype.bas, msk.solitem.suc, 0, m, z)
x, z = matrix(x), matrix(z)
if p is not 0:
yu, yl = zeros(p, Float), zeros(p, Float)
task.getsolutionslice(msk.soltype.bas, msk.solitem.suc, m, m+p, yu)
task.getsolutionslice(msk.soltype.bas, msk.solitem.slc, m, m+p, yl)
y = matrix(yu) - matrix(yl)
else:
y = matrix(0.0, (0,1))
if (solsta[0] is msk.solsta.unknown):
return (solsta[0], None, None, None)
else:
return (solsta[0], x, z, y)ind += m**2
rx, ry = xnewcopy(x0), ynewcopy(b)
rznl, rzl = matrix(0.0, (mnl, 1)), matrix(0.0, (cdim, 1)),
dx, dy = xnewcopy(x), ynewcopy(y)
dz, ds = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1))
lmbda = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) +
sum(dims['s']), 1))
lmbdasq = matrix(0.0, (mnl + dims['l'] + sum(dims['q']) +
sum(dims['s']), 1))
sigs = matrix(0.0, (sum(dims['s']), 1))
sigz = matrix(0.0, (sum(dims['s']), 1))
dz2, ds2 = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1))
newx, newy = xnewcopy(x), ynewcopy(y)
newz, news = matrix(0.0, (mnl + cdim, 1)), matrix(0.0, (mnl + cdim, 1))
newrx = xnewcopy(x0)
newrznl = matrix(0.0, (mnl, 1))
rx0, ry0 = xnewcopy(x0), ynewcopy(b)
rznl0, rzl0 = matrix(0.0, (mnl, 1)), matrix(0.0, (cdim, 1)),
x0, dx0 = xnewcopy(x), xnewcopy(dx)
y0, dy0 = ynewcopy(y), ynewcopy(dy)
z0 = matrix(0.0, (mnl + cdim, 1))
dz0 = matrix(0.0, (mnl + cdim, 1))
dz20 = matrix(0.0, (mnl + cdim, 1))
s0 = matrix(0.0, (mnl + cdim, 1))
ds0 = matrix(0.0, (mnl + cdim, 1))
ds20 = matrix(0.0, (mnl + cdim, 1))# Solve for ux, uy, uz
f3(x, y, z)
# s := s - z
# = lambda o\ bs - uz.
blas.axpy(z, s, alpha=-1.0)
# f4(x, y, z, s) solves the same system as f4_no_ir, but applies
# iterative refinement.
if iters == 0:
if refinement:
wx, wy = matrix(q), matrix(b)
wz, ws = matrix(0.0, (cdim, 1)), matrix(0.0, (cdim, 1))
if refinement:
wx2, wy2 = matrix(q), matrix(b)
wz2, ws2 = matrix(0.0, (cdim, 1)), matrix(0.0, (cdim, 1))
def f4(x, y, z, s):
if refinement:
wx = matrix(np.copy(x))
wy = matrix(np.copy(y))
wz = matrix(np.copy(z))
ws = matrix(np.copy(s))
f4_no_ir(x, y, z, s)
for i in range(refinement):
wx2 = matrix(np.copy(wx))
wy2 = matrix(np.copy(wy))
wz2 = matrix(np.copy(wz))
ws2 = matrix(np.copy(ws))
res(x, y, z, s, wx2, wy2, wz2, ws2, W, lmbda)
f4_no_ir(wx2, wy2, wz2, ws2)if self.taps:
# Indices of branches with non-zero tap ratio
idxs = [branches.index(e) for e in branches if e.ratio != 0.0]
# Transformer off nominal turns ratio ( = 0 for lines ) (taps at
# "from" bus, impedance at 'to' bus, i.e. ratio = Vf / Vt)"
ratio = matrix([e.ratio for e in branches])
# Assign non-zero tap ratios
tap[idxs] = ratio[idxs]
# Phase shifters
# tap = tap .* exp(j*pi/180 * branch(:, SHIFT));
# Convert branch attribute in degrees to radians
if self.phase_shift:
phase_shift = matrix([e.phase_shift * pi / 180 for e in branches])
else:
phase_shift = matrix(0.0, (n_branches, 1))
tap = mul(tap, exp(j * phase_shift))
# Ytt = Ys + j*Bc/2;
# Yff = Ytt ./ (tap .* conj(tap));
# Yft = - Ys ./ conj(tap);
# Ytf = - Ys ./ tap;
Ytt = Ys + j*Bc/2
Yff = div(Ytt, (mul(tap, conj(tap))))
Yft = div(-Ys, conj(tap))
Ytf = div(-Ys, tap)
# Shunt admittance
# Ysh = (bus(:, GS) + j * bus(:, BS)) / baseMVA;
g_shunt = matrix([v.g_shunt for v in buses])
if self.bus_shunts:def f4(x, y, z, s):
if refinement:
wx = matrix(np.copy(x))
wy = matrix(np.copy(y))
wz = matrix(np.copy(z))
ws = matrix(np.copy(s))
f4_no_ir(x, y, z, s)
for i in range(refinement):
wx2 = matrix(np.copy(wx))
wy2 = matrix(np.copy(wy))
wz2 = matrix(np.copy(wz))
ws2 = matrix(np.copy(ws))
res(x, y, z, s, wx2, wy2, wz2, ws2, W, lmbda)
f4_no_ir(wx2, wy2, wz2, ws2)
y += wx2
y += wy2
blas.axpy(wz2, z)
blas.axpy(ws2, s)
