CVXOPT examples

593 days ago by mjfwest

import numpy from cvxopt.base import spmatrix from cvxopt.base import matrix as _matrix from cvxopt.blas import dot from cvxopt.solvers import qp, options RealNumber = float Integer=int 
       
n=4 n_S=numpy.matrix([[4e-2, 6e-3, -4e-3, 0.0],[6e-3, 1e-2, 0.0, 0.0],[-4e-3, 0.0, 2.5e-3, 0.0],[0.0, 0.0, 0.0, 0.0]]) 
       
S=_matrix(n_S) 
       
n_pbar=numpy.matrix([0.12,0.10,0.07,0.03]) n_pbar.shape=(4,1) pbar = _matrix(n_pbar) 
       
G=_matrix(0.0,(n,n)) G[::n+1]=-1.0 #print(G) 
       
h = _matrix(0.0, (n,1)) A = _matrix(1.0, (1,n)) b = _matrix(1.0) 
       
N = 100 mus = [ 10.0^(5.0*t/N-1.0) for t in xrange(N) ] options['show_progress'] = False xs = [qp(float(mu)*S, -pbar, G, h, A, b)['x'] for mu in mus ] 
       
returns = [ dot(pbar,x) for x in xs ] risks = [ sqrt(dot(x, S*x)) for x in xs ] 
       
P1=list_plot(zip(risks,returns), plotjoined=True) P1.axes_labels(['standard deviation','expected return']) P1.show() 
       
c1 = [ x[0] for x in xs ] c2 = [ x[0] + x[1] for x in xs ] c3 = [ x[0] + x[1] + x[2] for x in xs ] c4 = [ x[0] + x[1] + x[2] + x[3] for x in xs ] 
       
import pylab pylab.figure(1) pylab.fill(risks + [.20], c1 + [0.0], facecolor = '#F0F0F0') pylab.fill(risks[-1::-1] + risks, c2[-1::-1] + c1, facecolor = '#D0D0D0') pylab.fill(risks[-1::-1] + risks, c3[-1::-1] + c2, facecolor = '#F0F0F0') pylab.fill(risks[-1::-1] + risks, c4[-1::-1] + c3, facecolor = '#D0D0D0') pylab.axis([0.0, 0.2, 0.0, 1.0]) pylab.savefig('temp.png') 
       
from cvxopt import blas, lapack from cvxopt.solvers import cp def spdiag(x): return spmatrix(x,range(len(x)),range(len(x))) my_list=[1.0,1.0,1.0,1.0] print(len(my_list)) spdiag(my_list) 
       
4
<4x4 spmatrix, tc='d', nnz=4>
4
<4x4 spmatrix, tc='d', nnz=4>
V = _matrix([-2.1213, 2.1213, -2.2981, 1.9284, -2.4575, 1.7207, -2.5981, 1.5000, -2.7189, 1.2679, -2.8191, 1.0261, -2.8978, 0.7765, -2.9544, 0.5209, -2.9886, 0.2615, -3.0000, 0.0000, 1.5000, 0.0000, 1.4772, -0.2605, 1.4095, -0.5130, 1.2990, -0.7500, 1.1491, -0.9642, 0.9642, -1.1491, 0.7500, -1.2990, 0.5130, -1.4095, 0.2605, -1.4772, 0.0000, -1.5000 ], (2,20)) 
       
n = V.size[1] G = spmatrix(-1.0, range(n), range(n)) h = _matrix(0.0, (n,1)) A = _matrix(1.0, (1,n)) b = _matrix(1.0) 
       
# D-design # # minimize f(x) = -log det V*diag(x)*V' # subject to x >= 0 # sum(x) = 1 # # The gradient and Hessian of f are # # gradf = -diag(V' * X^-1 * V) # H = (V' * X^-1 * V)**2. # # where X = V * diag(x) * V'. def F(x=None, z=None): if x is None: return 0, _matrix(float(1.0), (n,1)) #print(x) #print(spdiag(x)) X = V * spdiag(x) * V.T L = +X try: lapack.potrf(L) except ArithmeticError: return None f = -2.0 * (log(L[0,0]) + log(L[1,1])) W = +V blas.trsm(L, W) gradf = _matrix(-1.0, (1,2)) * W**2 if z is None: return f, gradf H = _matrix(0.0,(n,n)) blas.syrk(W, H, trans='T') return f, gradf, z[0] * H**2 
       
xd = cp(F, G, h, A = A, b = b)['x'] 
       
pylab.figure(2, facecolor='w', figsize=(6,6)) pylab.plot(V[0,:], V[1,:],'ow', [0], [0], 'k+') I = [ k for k in xrange(n) if xd[k] > 1e-5 ] pylab.plot(V[0,I], V[1,I],'or') 
       
[<matplotlib.lines.Line2D object at 0x4c55450>,
<matplotlib.lines.Line2D object at 0x4c55490>]
[<matplotlib.lines.Line2D object at 0x4c55450>, <matplotlib.lines.Line2D object at 0x4c55490>]
from numpy import array as _array # Enclosing ellipse is {x | x' * (V*diag(xe)*V')^-1 * x = sqrt(2)} nopts = 1000 angles = _array(_matrix([ float(c) for c in [ a*2.0*pi/nopts for a in xrange(nopts)]] , (1,nopts) )) circle = _matrix(0.0, (2,nopts)) circle[0,:], circle[1,:] = cos(angles), sin(angles) 
       
W = V * spdiag(xd) * V.T lapack.potrf(W) ellipse = sqrt(2.0) * circle blas.trmm(W, ellipse) pylab.plot(ellipse[0,:].T, ellipse[1,:].T, 'k--') pylab.axis([-5, 5, -5, 5]) pylab.title('D-optimal design (fig. 7.9)') pylab.axis('off') pylab.savefig('temp2.png')