# Minimum risk Markowitz portfolio selection example
#   minimize     f
#   subject to  [ r'  0   0 ] [ x ]   [ t-r'*w0 ] 
#               [ G   0  -I ]*[ f ] = [  -G*w0  ] 
#               [ e'  0   0 ] [ g ]   [    0    ]    
#               ||g||_2 <= f ,  and optionally w0 + x >= 0    
import sys, csv, mosek
def streamprinter(x): sys.stdout.write(x); sys.stdout.flush()
env = mosek.Env(); env.set_Stream(mosek.streamtype.log, streamprinter)
task = env.Task(0,0); task.set_Stream(mosek.streamtype.log, streamprinter)
SHORTSELLING = 1

def main ():
  f_itrT = open('../data/sp500-ret.csv','rb')  # expected returns
  rT     = csv.reader(f_itrT, quoting=csv.QUOTE_NONNUMERIC).next()
  f_itG  = open('../data/sp500-Gqr.csv','rb')  # covariance factor
  itG    = csv.reader(f_itG, quoting=csv.QUOTE_NONNUMERIC)
  G      = [[(i,v) for i,v in enumerate(row) if v!=0.0] for row in itG]
  n = len(rT);  p = 0;  w0 = [1.0/n]*n;  t  = 1.005;  inf = 0.0

  task.putobjsense(mosek.objsense.minimize)
  task.append(mosek.accmode.var, n+1+len(G)) # x, f and g
  task.append(mosek.accmode.con, 1+len(G)+1) # constraints
  task.putcj(n,1.0)

  bc = t - sum([rT[i]*w0[i] for i in range(n)])
  task.putavec(mosek.accmode.con, 0, range(n), rT) 
  task.putbound(mosek.accmode.con, 0, mosek.boundkey.fx, bc, bc)

  for row in G:
    Gidx = [v[0] for v in row];  Gval = [v[1] for v in row];
    bc = -sum([Gval[i]*w0[Gidx[i]] for i in range(len(Gval))])
    task.putavec(mosek.accmode.con, 1+p, Gidx+[n+1+p], Gval+[-1.0])
    task.putbound(mosek.accmode.con, 1+p, mosek.boundkey.fx, bc, bc)
    p += 1

  task.putavec(mosek.accmode.con, 1+p, range(n), [1.0]*n)
  task.putbound(mosek.accmode.con, 1+p, mosek.boundkey.fx, 0.0, 0.0)
  
  task.putboundlist(mosek.accmode.var, range(n+1+p), 
    [mosek.boundkey.fr]*(n+1+p), [-inf]*(n+1+p), [inf]*(n+1+p))
  if SHORTSELLING==0: task.putboundlist(mosek.accmode.var, range(n),
                        [mosek.boundkey.lo]*n, [-i for i in w0], [inf]*n)

  task.appendcone(mosek.conetype.quad, 0.0, range(n, n+1+p))

  task.optimize();  x = [0.0]*n;  f = [0.0]*1
  task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, 0, n, x)
  task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, n, n+1, f)
  [prosta,solsta] = task.getsolutionstatus(mosek.soltype.itr)
if __name__ == '__main__': main()
