Python cvxopt.matrix() Examples

def radius(K):
    """evaluate the radius of the MEB (Minimum Enclosing Ball) of examples in
    feature space.

    Parameters
    ----------
    K : (n,n) ndarray,
        the kernel that represents the data.

    Returns
    -------
    r : np.float64,
        the radius of the minimum enclosing ball of examples in feature space.
    """
    K = validation.check_K(K).numpy()
    n = K.shape[0]
    P = 2 * matrix(K)
    p = -matrix(K.diagonal())
    G = -spdiag([1.0] * n)
    h = matrix([0.0] * n)
    A = matrix([1.0] * n).T
    b = matrix([1.0])
    solvers.options['show_progress']=False
    sol = solvers.qp(P,p,G,h,A,b)
    return abs(sol['primal objective'])**.5 

def _update_grad(self,Kc, YY, beta, alpha, gamma):
        n = self.n_kernels
        
        eb = np.exp(beta)
        a,b = [], []
        gammaY = gamma['x'].T*YY
        for K in self.KL:   #optimized for generators
            K = matrix(K.numpy().astype(np.double))
            a.append( 1.-(alpha['x'].T*matrix(K)*alpha['x'])[0] )
            b.append( (gammaY*matrix(K)*gammaY.T)[0] )
        ebb, eba = np.dot(eb,b), np.dot(eb,a)
        den = np.dot(eb,b)**2
        num = [eb[r] * (a[r]*ebb - b[r]*eba) for r in range(n)]
        
        new_beta = np.array([beta[k] - self.learning_rate * (num[k]/den) for k in range(n)])
        return new_beta 

def _evalPolicyMatrix(self):
        # Evaluate the value function of the policy using linear equations.
        #
        # Arguments
        # ---------
        # Let S = number of states, A = number of actions
        # P(SxSxA) = transition matrix
        #      P could be an array with 3 dimensions or a cell array (1xA),
        #      each cell containing a matrix (SxS) possibly sparse
        # R(SxSxA) or (SxA) = reward matrix
        #      R could be an array with 3 dimensions (SxSxA) or
        #      a cell array (1xA), each cell containing a sparse matrix (SxS)
        #      or a 2D array(SxA) possibly sparse
        # discount = discount rate in ]0; 1[
        # policy(S) = a policy
        #
        # Evaluation
        # ----------
        # Vpolicy(S) = value function of the policy
        #
        Ppolicy, Rpolicy = self._computePpolicyPRpolicy()
        # V = PR + gPV  => (I-gP)V = PR  => V = inv(I-gP)* PR
        self.V = _np.linalg.solve(
            (_sp.eye(self.S, self.S) - self.discount * Ppolicy), Rpolicy) 

def lk(E=1.):
        """element stiffness matrix"""
        nu = 0.3
        k = np.array([0.5 - nu / 6., 0.125 + nu / 8., -0.25 - nu / 12.,
            -0.125 + 0.375 * nu, -0.25 + nu / 12., -0.125 - nu / 8., nu / 6.,
            0.125 - 0.375 * nu])
        KE = E / (1 - nu**2) * np.array([
            [k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
            [k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
            [k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
            [k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
            [k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
            [k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
            [k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
            [k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]]])
        return KE 

def __init__(self, transitions, reward, discount, skip_check=False):
        # Initialise a linear programming MDP.
        # import some functions from cvxopt and set them as object methods
        try:
            from cvxopt import matrix, solvers
            self._linprog = solvers.lp
            self._cvxmat = matrix
        except ImportError:
            raise ImportError("The python module cvxopt is required to use "
                              "linear programming functionality.")
        # initialise the MDP. epsilon and max_iter are not needed
        MDP.__init__(self, transitions, reward, discount, None, None,
                     skip_check=skip_check)
        # Set the cvxopt solver to be quiet by default, but ...
        # this doesn't do what I want it to do c.f. issue #3
        if not self.verbose:
            solvers.options['show_progress'] = False 

def fit(self, Xs, Xt):
        '''
        Fit source and target using KMM (compute the coefficients)
        :param Xs: ns * dim
        :param Xt: nt * dim
        :return: Coefficients (Pt / Ps) value vector (Beta in the paper)
        '''
        ns = Xs.shape[0]
        nt = Xt.shape[0]
        if self.eps == None:
            self.eps = self.B / np.sqrt(ns)
        K = kernel(self.kernel_type, Xs, None, self.gamma)
        kappa = np.sum(kernel(self.kernel_type, Xs, Xt, self.gamma) * float(ns) / float(nt), axis=1)

        K = matrix(K)
        kappa = matrix(kappa)
        G = matrix(np.r_[np.ones((1, ns)), -np.ones((1, ns)), np.eye(ns), -np.eye(ns)])
        h = matrix(np.r_[ns * (1 + self.eps), ns * (self.eps - 1), self.B * np.ones((ns,)), np.zeros((ns,))])

        sol = solvers.qp(K, -kappa, G, h)
        beta = np.array(sol['x'])
        return beta 

def _computeReward(self, reward, transition):
        # Compute the reward for the system in one state chosing an action.
        # Arguments
        # Let S = number of states, A = number of actions
        # P could be an array with 3 dimensions or  a cell array (1xA),
        # each cell containing a matrix (SxS) possibly sparse
        # R could be an array with 3 dimensions (SxSxA) or  a cell array
        # (1xA), each cell containing a sparse matrix (SxS) or a 2D
        # array(SxA) possibly sparse
        try:
            if reward.ndim == 1:
                return self._computeVectorReward(reward)
            elif reward.ndim == 2:
                return self._computeArrayReward(reward)
            else:
                r = tuple(map(self._computeMatrixReward, reward, transition))
                return r
        except (AttributeError, ValueError):
            if len(reward) == self.A:
                r = tuple(map(self._computeMatrixReward, reward, transition))
                return r
            else:
                return self._computeVectorReward(reward) 

def lk(E=1.):
        """element stiffness matrix"""
        nu = 0.3
        k = np.array([0.5 - nu / 6., 0.125 + nu / 8., -0.25 - nu / 12.,
            -0.125 + 0.375 * nu, -0.25 + nu / 12., -0.125 - nu / 8., nu / 6.,
            0.125 - 0.375 * nu])
        KE = E / (1 - nu**2) * np.array([
            [k[0], k[1], k[2], k[3], k[4], k[5], k[6], k[7]],
            [k[1], k[0], k[7], k[6], k[5], k[4], k[3], k[2]],
            [k[2], k[7], k[0], k[5], k[6], k[3], k[4], k[1]],
            [k[3], k[6], k[5], k[0], k[7], k[2], k[1], k[4]],
            [k[4], k[5], k[6], k[7], k[0], k[1], k[2], k[3]],
            [k[5], k[4], k[3], k[2], k[1], k[0], k[7], k[6]],
            [k[6], k[3], k[4], k[1], k[2], k[7], k[0], k[5]],
            [k[7], k[2], k[1], k[4], k[3], k[6], k[5], k[0]]])
        return KE 

def __init__(self, x, y, model, opts, sample_negative=False):
        """

        :param x: np.ndarray
            The instance matrix with ALL instances
        :param y: np.array
        :param model: Aad
        :param opts: AadOpts
        :param sample_negative: bool
        """
        self.x = x
        self.y = y
        self.model = model
        self.opts = opts
        self.sample_negative = sample_negative

        self.meta = None 

def test_ilp():
    """
    Problem defined in:
        https://en.wikipedia.org/wiki/Integer_programming#Example
    Solution from:
        https://stackoverflow.com/questions/33785396/python-the-integer-linear-programming-ilp-function-in-cvxopt-is-not-generati
    """
    import cvxopt
    from cvxopt import glpk

    glpk.options['msg_lev'] = 'GLP_MSG_OFF'
    c = cvxopt.matrix([0, -1], tc='d')
    G = cvxopt.matrix([[-1, 1], [3, 2], [2, 3], [-1, 0], [0, -1]], tc='d')
    h = cvxopt.matrix([1, 12, 12, 0, 0], tc='d')
    (status, x) = cvxopt.glpk.ilp(c, G.T, h, I=set([0, 1]))
    print (status)
    print ("%s, %s" % (str(x[0]), str(x[1])))
    print (sum(c.T * x)) 

def cvxopt_matrix(M):
    """
    Convert matrix to CVXOPT format.

    Parameters
    ----------
    M : numpy.array
        Matrix in NumPy format.

    Returns
    -------
    N : cvxopt.matrix
        Matrix in CVXOPT format.
    """
    if type(M) is ndarray:
        return matrix(M)
    elif type(M) is spmatrix or type(M) is matrix:
        return M
    coo = M.tocoo()
    return spmatrix(
        coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=M.shape) 

def min_singular(forces, torques, normals, soft_fingers=False, params=None):
        """ Min singular value of grasp matrix - measure of wrench that grasp is "weakest" at resisting.

        Parameters
        ----------
        forces : 3xN :obj:`numpy.ndarray`
            set of forces on object in object basis
        torques : 3xN :obj:`numpy.ndarray`
            set of torques on object in object basis
        normals : 3xN :obj:`numpy.ndarray`
            surface normals at the contact points
        soft_fingers : bool
            whether or not to use the soft finger contact model
        params : :obj:`GraspQualityConfig`
            set of parameters for grasp matrix and contact model

        Returns
        -------
        float : value of smallest singular value
        """
        G = PointGraspMetrics3D.grasp_matrix(forces, torques, normals, soft_fingers)
        _, S, _ = np.linalg.svd(G)
        min_sig = S[5]
        return min_sig 

def _QP(self, x, y):
        # In QP formulation (dual): m variables, 2m+1 constraints (1 equation, 2m inequations)
        m = len(y)
        print x.shape, y.shape
        P = self._kernel(x) * np.outer(y, y)
        P, q = matrix(P, tc='d'), matrix(-np.ones((m, 1)), tc='d')
        G = matrix(np.r_[-np.eye(m), np.eye(m)], tc='d')
        h = matrix(np.r_[np.zeros((m,1)), np.zeros((m,1)) + self.C], tc='d')
        A, b = matrix(y.reshape((1,-1)), tc='d'), matrix([0.0])
        # print "P, q:"
        # print P, q
        # print "G, h"
        # print G, h
        # print "A, b"
        # print A, b
        solution = solvers.qp(P, q, G, h, A, b)
        if solution['status'] == 'unknown':
            print 'Not PSD!'
            exit(2)
        else:
            self.alphas = np.array(solution['x']).squeeze() 

def calculateA(input_,label,C,K):          #using cvxopt bag figure out the convex quadratic programming
    num = input_.shape[0]
    P = cvxopt.matrix(numpy.outer(label,label)*K)
    q = cvxopt.matrix(numpy.ones(num)*-1)
    A = cvxopt.matrix(label,(1,num))
    b = cvxopt.matrix(0.0)
    
    if C is None:
        G = cvxopt.matrix(numpy.diag(numpy.ones(num)*-1))
        h = cvxopt.matrix(numpy.zeros(num))
    else:
        temp1 = numpy.diag(numpy.ones(num)*-1)
        temp2 = bp.identity(num)
        G = cvxopt.matrix(numpy.vstack(temp1,temp2))
        temp1 = numpy.zeros(num)
        temp2 = numpy.ones(num)*self.C
        h = cvxopt.matrix(numpy.hstack(temp1,temp2))        #P\q\A\b\G\h are parameters of cvxopt.solvers.qp() function

    solution = cvxopt.solvers.qp(P,q,G,h,A,b)               #figure out the model
    
    a = numpy.ravel(solution['x'])      #transfer the 'a' into a vector
    return a 

def projection_in_norm(self, x, M):
        """
        Projection of x to simplex induced by matrix M. Uses quadratic programming.
        """
        m = M.shape[0]

        # Constrains matrices
        P = opt.matrix(2 * M)
        q = opt.matrix(-2 * M * x)
        G = opt.matrix(-np.eye(m))
        h = opt.matrix(np.zeros((m, 1)))
        A = opt.matrix(np.ones((1, m)))
        b = opt.matrix(1.)

        # Solve using quadratic programming
        sol = opt.solvers.qp(P, q, G, h, A, b)
        return np.squeeze(sol['x']) 

def rebalance(self, obs):
        """
        Performs portfolio rebalance within environment
        :param obs: pandas DataFrame: Environment observation
        :return: numpy array: Portfolio vector
        """
        n_pairs = obs.columns.levels[0].shape[0]
        if self.step:
            prev_posit = self.get_portfolio_vector(obs, index=self.reb)
            price_relative = self.predict(obs)
            return self.update(prev_posit, price_relative)
        else:
            action = np.ones(n_pairs)
            action[-1] = 0
            self.sigma = np.matrix(np.eye(n_pairs) / n_pairs ** 2)
            return array_normalize(action) 

def decide_by_history(self, x, last_b):
        '''
        :param x: input matrix
        :param last_b: last portfolio
        '''
        x = self.get_last_rpv(x)
        if self.A is None:
            self.init_portfolio(x)

        # calculate gradient
        grad = np.mat(x / np.dot(last_b, x)).T
        # update A
        self.A += grad * grad.T
        # update b
        self.b += (1 + 1./self.beta) * grad

        # projection of p induced by norm A
        pp = self.projection_in_norm(self.delta * self.A.I * self.b, self.A)
        return pp * (1 - self.eta) + np.ones(len(x)) / float(len(x)) * self.eta 

def test1():
    import cvxopt
    from cvxopt import matrix

    P = matrix(np.diag([1.0, 0.0]))
    q = matrix(np.array([3.0, 4.0]).T)
    G = matrix(np.array([[-1.0, 0.0], [0, -1.0], [-1.0, -3.0], [2.0, 5.0], [3.0, 4.0]]))
    h = matrix(np.array([0.0, 0.0, -15.0, 100.0, 80.0]).T)

    sol = cvxopt.solvers.qp(P, q, G, h)

    #  print(sol)
    print(sol["x"])
    #  print(sol["primal objective"])

    assert sol["x"][0] - 0.0, "Error1"
    assert sol["x"][1] - 5.0, "Error2"

    P = np.diag([1.0, 0.0])
    q = np.matrix([3.0, 4.0]).T
    G = np.matrix([[-1.0, 0.0], [0, -1.0], [-1.0, -3.0], [2.0, 5.0], [3.0, 4.0]])
    h = np.matrix([0.0, 0.0, -15.0, 100.0, 80.0]).T

    sol2 = ecosqp(P, q, G, h)

    for i in range(len(sol["x"])):
        assert (sol["x"][i] - sol2["x"][i]) <= 0.001, "Error1" 

def plot_forest_contours_2D(x, y, xx, yy, budget, forest, pdfpath_contours, dash_xy, dash_wh):
    # Original detector contours
    tm = Timer()
    tm.start()
    baseline_scores = 0.5 - forest.decision_function(x)
    queried = np.argsort(-baseline_scores)
    # logger.debug("baseline scores:%s\n%s" % (str(baseline_scores.shape), str(list(baseline_scores))))

    # n_found = np.cumsum(y[queried[np.arange(budget)]])
    # print n_found

    Z_if = 0.5 - forest.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z_if = Z_if.reshape(xx.shape)

    dp = DataPlotter(pdfpath=pdfpath_contours, rows=1, cols=1)
    pl = dp.get_next_plot()
    pl.contourf(xx, yy, Z_if, 20, cmap=plt.cm.get_cmap('jet'))

    dp.plot_points(x, pl, labels=y, lbl_color_map={0: "grey", 1: "red"})

    dp.plot_points(matrix(x[queried[np.arange(budget)], :], nrow=budget),
                   pl, labels=y[queried[np.arange(budget)]], defaultcol="red",
                   lbl_color_map={0: "green", 1: "red"}, edgecolor="black",
                   marker=matplotlib.markers.MarkerStyle('o', fillstyle=None), s=35)

    # plot the sidebar
    anom_idxs = np.where(y[queried] == 1)[0]
    dash = 1 - (anom_idxs * 1.0 / x.shape[0])
    plot_sidebar(dash, dash_xy, dash_wh, pl)

    dp.close()

    logger.debug(tm.message("plotted contours")) 

def test1():
    import cvxopt
    from cvxopt import matrix

    P = matrix(np.diag([1.0, 0.0]))
    q = matrix(np.array([3.0, 4.0]).T)
    G = matrix(np.array([[-1.0, 0.0], [0, -1.0], [-1.0, -3.0], [2.0, 5.0], [3.0, 4.0]]))
    h = matrix(np.array([0.0, 0.0, -15.0, 100.0, 80.0]).T)

    sol = cvxopt.solvers.qp(P, q, G, h)

    #  print(sol)
    print(sol["x"])
    #  print(sol["primal objective"])

    assert sol["x"][0] - 0.0, "Error1"
    assert sol["x"][1] - 5.0, "Error2"

    P = np.diag([1.0, 0.0])
    q = np.matrix([3.0, 4.0]).T
    G = np.matrix([[-1.0, 0.0], [0, -1.0], [-1.0, -3.0], [2.0, 5.0], [3.0, 4.0]])
    h = np.matrix([0.0, 0.0, -15.0, 100.0, 80.0]).T

    sol2 = ecosqp(P, q, G, h)

    for i in range(len(sol["x"])):
        assert (sol["x"][i] - sol2["x"][i]) <= 0.001, "Error1" 

def force_closure_qp(forces, torques, normals, soft_fingers=False,
                         wrench_norm_thresh=1e-3, wrench_regularizer=1e-10,
                         params=None):
        """ Checks force closure by solving a quadratic program (whether or not zero is in the convex hull)

        Parameters
        ----------
        forces : 3xN :obj:`numpy.ndarray`
            set of forces on object in object basis
        torques : 3xN :obj:`numpy.ndarray`
            set of torques on object in object basis
        normals : 3xN :obj:`numpy.ndarray`
            surface normals at the contact points
        soft_fingers : bool
            whether or not to use the soft finger contact model
        wrench_norm_thresh : float
            threshold to use to determine equivalence of target wrenches
        wrench_regularizer : float
            small float to make quadratic program positive semidefinite
        params : :obj:`GraspQualityConfig`
            set of parameters for grasp matrix and contact model

        Returns
        -------
        int : 1 if in force closure, 0 otherwise
        """
        if params is not None:
            if 'wrench_norm_thresh' in list(params.keys()):
                wrench_norm_thresh = params.wrench_norm_thresh
            if 'wrench_regularizer' in list(params.keys()):
                wrench_regularizer = params.wrench_regularizer

        G = PointGraspMetrics3D.grasp_matrix(forces, torques, normals, soft_fingers, params=params)
        min_norm, _ = PointGraspMetrics3D.min_norm_vector_in_facet(G, wrench_regularizer=wrench_regularizer)
        return 1 * (min_norm < wrench_norm_thresh)  # if greater than wrench_norm_thresh, 0 is outside of hull 

def fit(self, x, y):
        # x, y: np.ndarray
        # x.shape: (m, n), where m = # samples, n = # features.
        # y.shape: (m,), m labels which range from {-1.0, +1.0}.
        assert type(x) == np.ndarray
        # print type(x), type(y)
        print x.shape, y.shape
        # In the design matrix x: m examples, n features
        # In QP formulation (dual): m variables, 2m+1 constraints (1 equation, 2m inequations)
        # print 'x = ', x
        # print 'y = ', y
        if self.kernel == 'rbf' and self.gamma is None:
            self.gamma = 1.0 / x.shape[1]
            print 'gamma = ', self.gamma
        if self.alg == 'dual':
            self._QP(x, y)
        else:
            assert self.alg == 'SMO'
            K = self._kernel(x)
            self._SMO5(K, y)

        logging.info('self.alphas = ' + str(self.alphas))
        sv_indices = filter(lambda i:self.alphas[i] > self.eps, xrange(len(y)))
        self.sv_x, self.sv_y, self.alphas = x[sv_indices], y[sv_indices], self.alphas[sv_indices]
        self.n_sv = len(sv_indices)
        logging.info('sv_indices:' + str(sv_indices))
        print self.n_sv, 'SVs!'
        logging.info(str(np.c_[self.sv_x, self.sv_y]))
        if self.kernel == 'linear':
            self.w = np.dot(self.alphas * self.sv_y, self.sv_x)
        if self.alg == 'dual':
            # calculate b: average over all support vectors
            sv_boundary = self.alphas < self.C - self.eps
            self.b = np.mean(self.sv_y[sv_boundary] - np.dot(self.alphas * self.sv_y,
                                                             self._kernel(self.sv_x, self.sv_x[sv_boundary]))) 

def fit(self,input_,label):         #fit data
        samples,features = input_.shape
        
        K = makeKernelMatrix(input_,self.kernel)          #calculate the kernel matrix
        a = calculateA(input_,label,self.C,K)             #calculate the process parameter 'a'
        supportVector = a > 1e-5                          #if a < 1e-5,then think a is 1
        indexis = numpy.arange(len(a))[supportVector]     #the support vectors' indexis
        self.a = a[supportVector]                         #the support vectors' a
        self.supportVector = input_[supportVector]        #the support vectors
        self.supportVectorLabel = label[supportVector]    #the support vectorss' label
        
        print(len(self.a),' support vectors out of ',samples,' points')
        
        self.b = calculateB(self.a,self.supportVectorLabel,supportVector,K,indexis)   #calculate the model's risdual
        self.w = calculateWeight(self.kernel,features,self.a,self.supportVector,self.supportVectorLabel)    #calculate the model's weights 

def run(self):
        # Run the policy iteration algorithm.
        self._startRun()

        while True:
            self.iter += 1
            # these _evalPolicy* functions will update the classes value
            # attribute
            if self.eval_type == "matrix":
                self._evalPolicyMatrix()
            elif self.eval_type == "iterative":
                self._evalPolicyIterative()
            # This should update the classes policy attribute but leave the
            # value alone
            policy_next, null = self._bellmanOperator()
            del null
            # calculate in how many places does the old policy disagree with
            # the new policy
            n_different = (policy_next != self.policy).sum()
            # if verbose then continue printing a table
            if self.verbose:
                _printVerbosity(self.iter, n_different)
            # Once the policy is unchanging of the maximum number of
            # of iterations has been reached then stop
            if n_different == 0:
                if self.verbose:
                    print(_MSG_STOP_UNCHANGING_POLICY)
                break
            elif self.iter == self.max_iter:
                if self.verbose:
                    print(_MSG_STOP_MAX_ITER)
                break
            else:
                self.policy = policy_next

        self._endRun() 

def makeKernelMatrix(input_,kernel,par = 3):      #set up kernel matrix according to the kernel
    num = input_.shape[0]
    K = numpy.zeros((num,num))
    for i in range(num):
        for j in range(num):
            if kernel == 'linearKernel':
                K[i,j] = linearKernel(input_[i],input_[j])
            else:
                K[i,j] = gaussianKernel(input_[i],input_[j],par)
    return K 

def SOCP1Test(solver_to_test):
        #first test (A optimality)
        primal=prob_multiresponse_A.copy()
        try:
                primal.solve(solver=solver_to_test,timelimit=1,maxit=50)
        except Exception as ex:
                traceback.print_exc(file=sys.stdout)
                return (False,repr(ex))
        
        primf = primal.check_current_value_feasibility()
        if not(primf[0]):
                return (False,'not primal feasible|{0:1.0e}'.format(primf[1]))
        obj=1.0759874194855403
        
        pgap = abs(primal.obj_value()-obj)/abs(obj)
        if pgap>0.1:
                return (False,'failed')        
        elif pgap>1e-5:
                return (False,'not primal optimal|{0:1.0e}'.format(pgap))
        
        Zvar=prob_dual_A.get_variable('Z')
        muvar=prob_dual_A.get_variable('mu')
        
        try:
                Z=[cvx.matrix(cs.dual[1:],Zvar[i].size) for i,cs in enumerate(primal.constraints)]
                mu=[cs.dual[0] for cs in primal.constraints]
        except TypeError:
                return (False,'no dual computed')
        
        muvar.value=mu
        for i,zi in enumerate(Z):
                Zvar[i].value=zi
                
        dualf=prob_dual_A.check_current_value_feasibility(tol=1e-5)
        if not(dualf[0]):
                return (False,'not dual feasible|{0:1.0e}'.format(dualf[1]))
        dgap = abs(prob_dual_A.obj_value()-obj)/abs(obj)
        if dgap >1e-5:
                return (False,'not dual optimal|{0:1.0e}'.format(dgap))
        return (True,primal.status) 

def projection_in_norm(self, x, M):
        """ Projection of x to simplex indiced by matrix M. Uses quadratic programming.
        """
        m = M.shape[0]

        P = matrix(2*M)
        q = matrix(-2 * M * x)
        G = matrix(-np.eye(m))
        h = matrix(np.zeros((m,1)))
        A = matrix(np.ones((1,m)))
        b = matrix(1.)

        sol = solvers.qp(P, q, G, h, A, b)
        return np.squeeze(sol['x']) 

def slack_var(self):
        return cvx.matrix([1 - norm(self.exp, 'inf').value,
                           self.radius - norm(self.exp, 1).value]) 

def fit_nnl2reg(self):
        try:
            from cvxopt import matrix, solvers
        except ImportError:
            raise ImportError("To infer titer models, you need a working installation of cvxopt")
        n_params = self.design_matrix.shape[1]
        P = matrix(np.dot(self.design_matrix.T, self.design_matrix) + self.lam_drop*np.eye(n_params))
        q = matrix( -np.dot( self.titer_dist, self.design_matrix))
        h = matrix(np.zeros(n_params)) # Gw <=h
        G = matrix(-np.eye(n_params))
        W = solvers.qp(P,q,G,h)
        return np.array([x for x in W['x']]) 

def wrench_volume(forces, torques, normals, soft_fingers=False, params=None):
        """ Volume of grasp matrix singular values - score of all wrenches that the grasp can resist.

        Parameters
        ----------
        forces : 3xN :obj:`numpy.ndarray`
            set of forces on object in object basis
        torques : 3xN :obj:`numpy.ndarray`
            set of torques on object in object basis
        normals : 3xN :obj:`numpy.ndarray`
            surface normals at the contact points
        soft_fingers : bool
            whether or not to use the soft finger contact model
        params : :obj:`GraspQualityConfig`
            set of parameters for grasp matrix and contact model

        Returns
        -------
        float : value of wrench volume
        """
        k = 1
        if params is not None and 'k' in list(params.keys()):
            k = params.k

        G = PointGraspMetrics3D.grasp_matrix(forces, torques, normals, soft_fingers)
        _, S, _ = np.linalg.svd(G)
        sig = S
        return k * np.sqrt(np.prod(sig))