changeset 65:3d53309236c4

Added cvxopt_lp.py. Some bugfixes in gap.py and l1eq_pd.py

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyCSalgos/BP/cvxopt_lp.py	Fri Mar 16 13:35:47 2012 +0200
@@ -0,0 +1,41 @@
+
+import numpy
+import cvxopt
+import cvxopt.solvers
+import cvxopt.msk
+import mosek
+
+def cvxopt_lp(y, A):
+
+    N = A.shape[1]
+    AA = numpy.hstack((A, -A))
+    c = numpy.ones((2*N, 1))
+    
+    G = numpy.vstack((numpy.zeros((2*N,2*N)), -numpy.eye(2*N)))
+    h = numpy.zeros((4*N,1))
+
+    # Convert numpy arrays to cvxopt matrices
+    cvx_c = cvxopt.matrix(c)
+    cvx_G = cvxopt.matrix(G)
+    cvx_h = cvxopt.matrix(h)
+    cvx_AA = cvxopt.matrix(AA)
+    cvx_y = cvxopt.matrix(y.reshape(y.size,1))
+    
+    # Options    
+    cvxopt.solvers.options['show_progress'] = False
+    #cvxopt.solvers.options['MOSEK'] = {mosek.iparam.log: 0}
+    
+    # Solve
+    #res = cvxopt.solvers.lp(cvx_c, cvx_G, cvx_h, A=cvx_AA, b=cvx_y, solver='mosek')
+    res = cvxopt.solvers.lp(cvx_c, cvx_G, cvx_h, A=cvx_AA, b=cvx_y)
+      
+    primal = numpy.squeeze(numpy.array(res['x']))
+    gamma = primal[:N] - primal[N:]
+    return gamma
+    
+    #lb = zeros(2*N,1);
+    #ub = Inf*ones(2*N,1);
+    ##[primal,obj,exitflag,output2,dual] = linprog(c,[],[],A,y,lb,ub,[],opt);
+    ##[primal,~,~,~,~] = linprog(c,[],[],A,aggy,lb,ub);
+    #[primal,obj,exitflag,output2,dual] = linprog(c,[],[],A,aggy,lb,ub);
+    #gamma = primal(1:N) - primal((N+1):(2*N));
\ No newline at end of file

--- a/pyCSalgos/BP/l1eq_pd.py	Tue Mar 13 16:18:11 2012 +0200
+++ b/pyCSalgos/BP/l1eq_pd.py	Fri Mar 16 13:35:47 2012 +0200
@@ -221,6 +221,11 @@
     # End of original Matab code
     #----------------------------    
     
+    # Nic: check if b is 0; if so, return 0
+    #    Otherwise it will break later
+    if numpy.linalg.norm(b) < 1e-16:
+        return numpy.zeros_like(x0)
+    
     #largescale = isa(A,'function_handle');
     if hasattr(A, '__call__'):
         largescale = True
@@ -257,7 +262,7 @@
         #x0 = A'*w;
         try:
             w = scipy.linalg.solve(numpy.dot(A,A.T), b, sym_pos=True)
-            hcond = 1.0/numpy.linalg.cond(np.dot(A,A.T))
+            hcond = 1.0/numpy.linalg.cond(numpy.dot(A,A.T))
         except scipy.linalg.LinAlgError:
             if verbose:
               print 'A*At is ill-conditioned: cannot find starting point'

--- a/pyCSalgos/GAP/GAP.py	Tue Mar 13 16:18:11 2012 +0200
+++ b/pyCSalgos/GAP/GAP.py	Fri Mar 16 13:35:47 2012 +0200
@@ -55,7 +55,8 @@
     #if strcmp(params.l2solver, 'pseudoinverse')
     if params['l2solver'] == 'pseudoinverse':
     #if strcmp(class(M), 'double') && strcmp(class(Omega), 'double')
-      if M.dtype == 'float64' and Omega.dtype == 'double':
+      #if M.dtype == 'float64' and Omega.dtype == 'double':
+      if 1:
         while True:
             alpha = math.sqrt(lagmult);
             xhat = numpy.linalg.lstsq(numpy.concatenate((M, alpha*Omega[Lambdahat,:])), numpy.concatenate((y, numpy.zeros(Lambdahat.size))))[0]