How to use the cvxopt.spmatrix function in cvxopt

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)