import theano
from theano import tensor, function
import numpy as np
import cvxopt



def l1regls(M_val, w_val):
    """
    Solves:
        min - sum_j theta_j
        s.t theta_j <= 0
            |e^{M*theta} - (1 - w)|_2 <= 1
    """
    m, n = M_val.shape
    c = cvxopt.matrix(-1.0, (n,1))

    theta = tensor.row().T
    z = tensor.row().T
    theta_ = theta.flatten()
    z_ = z.flatten()
    M = theano.shared(M_val.astype(theano.config.floatX))
    w = theano.shared(w_val.astype(theano.config.floatX))
    y = (tensor.exp(M.dot(theta_)) - (1 - w)).norm(2) - 1
    y_diff = tensor.grad(y, theta_)
    y_hess = z[0] * theano.gradient.hessian(y, theta_)
    f_x = theano.function([theta], [y, y_diff], allow_input_downcast=True)
    f_xz = theano.function([theta, z], [y, y_diff, y_hess], allow_input_downcast=True)

    def F(x=None, z=None):
        if x is None:
            return 1, cvxopt.matrix(1.0, (n,1))
        elif z is None:
            y, y_diff = f_x(x)
            return cvxopt.matrix(float(y), (1, 1)),\
                    cvxopt.matrix(y_diff.astype("float64")).T
        else:
            y, y_diff, y_hess = f_xz(x, z)
            return cvxopt.matrix(float(y), (1, 1)), \
                    cvxopt.matrix(y_diff.astype("float64")).T, \
                    cvxopt.matrix(y_hess.astype("float64"))

    G = cvxopt.spdiag([1 for i in xrange(n)])
    h = cvxopt.matrix(0.0, (n,1))

    return cvxopt.solvers.cpl(c,F, G, h)['x']


def test():
    """
    unit test
    """
    M_val = np.random.rand(100, 20)
    w_val = np.random.rand(100)
    print l1regls(M_val, w_val)

if __name__=="__main__":
    test()