Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
C = cv.matrix(np.c_[np.ones(n), np.arange(1., n+1.)/n])
nC = C.size[1]
# Solve the problem:
# .5*(M*q + B*l + C*d - EDA_raw)^2 + alpha*sum(A,1)*p + .5*gamma*l'*l
# s.t. A*q >= 0
old_options = cv.solvers.options.copy()
cv.solvers.options.clear()
cv.solvers.options.update(options)
if solver == 'conelp':
# Use conelp
z = lambda m,n: cv.spmatrix([],[],[],(m,n))
G = cv.sparse([[-A,z(2,n),M,z(nB+2,n)],[z(n+2,nC),C,z(nB+2,nC)],
[z(n,1),-1,1,z(n+nB+2,1)],[z(2*n+2,1),-1,1,z(nB,1)],
[z(n+2,nB),B,z(2,nB),cv.spmatrix(1.0, range(nB), range(nB))]])
h = cv.matrix([z(n,1),.5,.5,EDA_raw,.5,.5,z(nB,1)])
c = cv.matrix([(cv.matrix(alpha, (1,n)) * A).T,z(nC,1),1,gamma,z(nB,1)])
res = cv.solvers.conelp(c, G, h, dims={'l':n,'q':[n+2,nB+2],'s':[]})
obj = res['primal objective']
else:
# Use qp
Mt, Ct, Bt = M.T, C.T, B.T
H = cv.sparse([[Mt*M, Ct*M, Bt*M], [Mt*C, Ct*C, Bt*C],
[Mt*B, Ct*B, Bt*B+gamma*cv.spmatrix(1.0, range(nB), range(nB))]])
f = cv.matrix([(cv.matrix(alpha, (1,n)) * A).T - Mt*EDA_raw, -(Ct*EDA_raw), -(Bt*EDA_raw)])
res = cv.solvers.qp(H, f, cv.spmatrix(-A.V, A.I, A.J, (n,len(f))),
cv.matrix(0., (n,1)), solver=solver)
obj = res['primal objective'] + .5 * (EDA_raw.T * EDA_raw)
cv.solvers.options.clear()
cv.solvers.options.update(old_options)
h = cvxopt.matrix([_cvx(n, 1), 0.5, 0.5, eda, 0.5, 0.5, _cvx(nB, 1)])
c = cvxopt.matrix([(cvxopt.matrix(alpha, (1, n)) * A).T, _cvx(nC, 1), 1, gamma, _cvx(nB, 1)])
res = cvxopt.solvers.conelp(c, G, h, dims={"l": n, "q": [n + 2, nB + 2], "s": []})
else:
# Use qp
Mt, Ct, Bt = M.T, C.T, B.T
H = cvxopt.sparse(
[
[Mt * M, Ct * M, Bt * M],
[Mt * C, Ct * C, Bt * C],
[Mt * B, Ct * B, Bt * B + gamma * cvxopt.spmatrix(1.0, range(nB), range(nB))],
]
)
f = cvxopt.matrix([(cvxopt.matrix(alpha, (1, n)) * A).T - Mt * eda, -(Ct * eda), -(Bt * eda)])
res = cvxopt.solvers.qp(
H, f, cvxopt.spmatrix(-A.V, A.I, A.J, (n, len(f))), cvxopt.matrix(0.0, (n, 1)), solver=solver
)
cvxopt.solvers.options.clear()
cvxopt.solvers.options.update(old_options)
tonic_splines = res["x"][-nB:]
drift = res["x"][n : n + nC]
tonic = B * tonic_splines + C * drift
q = res["x"][:n]
phasic = M * q
out = pd.DataFrame({"EDA_Tonic": np.array(tonic)[:, 0], "EDA_Phasic": np.array(phasic)[:, 0]})
return out
num_vars = len(self._var_list)
num_eq_cons = len(self._eq_cons_list)
num_leq_cons = len(self._leq_cons_list)
for var in self._var_list:
if var.lb is not None:
num_leq_cons += 1
if var.ub is not None:
num_leq_cons += 1
# Make the matrices (some need to be dense).
c = cvxopt.matrix([0.0] * num_vars)
h = cvxopt.matrix([0.0] * num_leq_cons)
g_mat = cvxopt.spmatrix([], [], [], (num_leq_cons, num_vars))
a_mat = None
b = None
if num_eq_cons > 0:
a_mat = cvxopt.spmatrix([], [], [], (num_eq_cons, num_vars))
b = cvxopt.matrix([0.0] * num_eq_cons)
# Objective coefficients: c
for var_label in self._obj_coeffs:
value = self._obj_coeffs[var_label]
vid = self._vars[var_label].vid
if self._objective == OBJ_MAX:
c[vid] = -value # negate the value because it's a max
else:
c[vid] = value # min objective matches cvxopt
# Inequality constraints: G, h
row = 0
for cons in self._leq_cons_list:
# If it's >= then need to negate all coeffs and the rhs
if cons.rhs is not None:
h[row] = cons.rhs if cons.ctype == CONS_TYPE_LEQ else -cons.rhs
for var_label in cons.coeffs:
def _cvx(m, n):
return cvxopt.spmatrix([], [], [], (m, n))
def identity(self, size):
return cvxopt.spmatrix(1, range(size), range(size))
L = np.zeros((n * (N + 1), 1))
for i in range(0, N):
ind1 = n + i * n + np.arange(n)
ind2x = i * n + np.arange(n)
ind2u = i * d + np.arange(d)
Gx[np.ix_(ind1, ind2x)] = -A[i]
Gu[np.ix_(ind1, ind2u)] = -B[i]
L[ind1, :] = C[i]
G = np.hstack((Gx, Gu))
if Controller.Solver == "CVX":
G_sparse = spmatrix(G[np.nonzero(G)], np.nonzero(G)[0].astype(int), np.nonzero(G)[1].astype(int), G.shape)
E_sparse = spmatrix(E[np.nonzero(E)], np.nonzero(E)[0].astype(int), np.nonzero(E)[1].astype(int), E.shape)
L_sparse = spmatrix(L[np.nonzero(L)], np.nonzero(L)[0].astype(int), np.nonzero(L)[1].astype(int), L.shape)
G_return = G_sparse
E_return = E_sparse
L_return = L_sparse
else:
G_return = G
E_return = E
L_return = L
return G_return, E_return, L_return
zeros_system_build_t = lambda shape: cvxopt.spmatrix( [], [], [], shape )
def find_nearest_valid_distribution(u_alpha, kernel, initial=None, reg=0):
""" (solution,distance_sqd)=find_nearest_valid_distribution(u_alpha,kernel):
Given a n-vector u_alpha summing to 1, with negative terms,
finds the distance (squared) to the nearest n-vector summing to 1,
with non-neg terms. Distance calculated using nxn matrix kernel.
Regularization parameter reg --
min_v (u_alpha - v)^\top K (u_alpha - v) + reg* v^\top v"""
P = matrix(2 * kernel)
n = kernel.shape[0]
q = matrix(np.dot(-2 * kernel, u_alpha))
A = matrix(np.ones((1, n)))
b = matrix(1.)
G = spmatrix(-1., range(n), range(n))
h = matrix(np.zeros(n))
dims = {'l': n, 'q': [], 's': []}
solvers.options['show_progress'] = False
solution = solvers.coneqp(
P,
q,
G,
h,
dims,
A,
b,
initvals=initial
)
distance_sqd = solution['primal objective'] + np.dot(u_alpha.T,
np.dot(kernel, u_alpha))[0, 0]
return (solution, distance_sqd)
I = - np.eye(2*N)
Zeros = np.zeros((2*N, F_hard.shape[1]))
Positivity = np.hstack((Zeros, I))
F = np.vstack((F_1, Positivity))
# np.savetxt('F.csv', F, delimiter=',', fmt='%f')
# pdb.set_trace()
# b vector
b = np.hstack((b_hard, np.zeros(2*N)))
if LMPC.Solver == "CVX":
F_sparse = spmatrix(F[np.nonzero(F)], np.nonzero(F)[0].astype(int), np.nonzero(F)[1].astype(int), F.shape)
F_return = F_sparse
else:
F_return = F
return F_return, b
# Original comments below
# Boyd & Vandenberghe, "Convex Optimization"
# Joelle Skaf - 08/23/05
#
# Solved a QCQP with 3 inequalities:
# minimize 1/2 x'*P0*x + q0'*r + r0
# s.t. 1/2 x'*Pi*x + qi'*r + ri <= 0 for i=1,2,3
# and verifies that strong duality holds.
# Input data
n = 6
eps = sys.float_info.epsilon
P0 = cvxopt.normal(n, n)
eye = cvxopt.spmatrix(1.0, range(n), range(n))
P0 = P0.T * P0 + eps * eye
print(P0)
P1 = cvxopt.normal(n, n)
P1 = P1.T*P1
P2 = cvxopt.normal(n, n)
P2 = P2.T*P2
P3 = cvxopt.normal(n, n)
P3 = P3.T*P3
q0 = cvxopt.normal(n, 1)
q1 = cvxopt.normal(n, 1)
q2 = cvxopt.normal(n, 1)
q3 = cvxopt.normal(n, 1)