import numpy as np; np.set_printoptions(precision=3, suppress=True)
import pandas as pd; from cvxopt import matrix; import cvxpy as cp;

num_assets = 8
x_orig = matrix(np.array([[ 0.   ], [ 0.   ], [ 0.069], [ 0.   ], [ 0.148], [ 0.086], [ 0.111], [ 0.585]]))
cov_mat = matrix(np.array([
 [ 0.154,  0.121,  0.025,  0.021,  0.028,  0.002,  0.054,  0.072],
 [ 0.121,  0.188,  0.013,  0.032,  0.068,  0.035,  0.053,  0.064],
 [ 0.025,  0.013,  0.067,  0.047, -0.023, -0.003,  0.007,  0.015],
 [ 0.021,  0.032,  0.047,  0.059,  0.012,  0.012,  0.006,  0.012],
 [ 0.028,  0.068, -0.023,  0.012,  0.153,  0.029,  0.043,  0.021],
 [ 0.002,  0.035, -0.003,  0.012,  0.029,  0.055, -0.035, -0.011],
 [ 0.054,  0.053,  0.007,  0.006,  0.043, -0.035,  0.103,  0.04 ],
 [ 0.072,  0.064,  0.015,  0.012,  0.021, -0.011,  0.04,   0.072]]))
e = matrix(np.ones((num_assets,1))); zeros = matrix(np.zeros((num_assets,1)))

def inverse_opt(xbar, cov, u, mu_bar, sig=0.2):
    alpha = cp.Variable(); beta = cp.Variable(num_assets); gamma = cp.Variable()
    mu = cp.Variable(num_assets);
    objective = cp.Minimize(cp.norm1(mu - mu_bar));
    constraints = [alpha >= 0, beta >= zeros, mu - 2*alpha*cov*xbar - beta - gamma*e <= zeros]
    constraints += [(u.T - xbar.T)*beta == 0, 2*sig**2*alpha + u.T*beta + gamma == mu.T*xbar]
    result = cp.Problem(objective, constraints).solve(solver=cp.ECOS)
    if cp.get_status(result) is cp.SOLVED:
        return mu.value
    else:
        print result
        return matrix([np.nan]*num_assets)

print inverse_opt(x_orig, cov_mat, e, zeros)