In [ ]:

    

NU = 6


    

In [ ]:

    

import numpy as np
import pandas as pd
import os

dir_sim_ = './simulations'
dir_tr_input_ = './transformed_input'

udl_distrib = pd.read_csv(os.path.join(dir_sim_, 'nu_eq_%i.csv' % NU), header=0, index_col=0)
members_pos = pd.read_csv(os.path.join(dir_tr_input_, 'positions.csv'), header=0, index_col=0)

loss_and_profit = np.dot(members_pos.values, udl_distrib.values).transpose()

print loss_and_profit.shape


    

In [ ]:

    

from lib.msra_loss import MSRALossFunctionAbs


    

In [ ]:

    

class LossFunction3(MSRALossFunctionAbs):
    # \mathbb{E} \left( \ell(X - m) \right)
    def shortfall_risk(self, m=None):
        m = self._check_argument(m)
        x_minus_m = np.subtract(self.x, m)
        
        sum_x_minus_m_plus = np.maximum(x_minus_m, 0.).sum(axis=1)
        sum_x_minus_m_minus = -np.minimum(x_minus_m, 0.).sum(axis=1)

        diff = sum_x_minus_m_plus - 0.5 * sum_x_minus_m_minus
        
        dbl_sum = None
        for i in xrange(self.dim):
            for j in xrange(i + 1, self.dim):
                xi_p_xj = x_minus_m[:, i] + x_minus_m[:, j]
                tmp = np.maximum(xi_p_xj, 0.) - 0.5 * (-np.minimum(xi_p_xj, 0.))
                if dbl_sum is None:
                    dbl_sum = tmp
                else:
                    dbl_sum += tmp
        
        res = diff + dbl_sum
        
        return res.mean()

    # \mathbb{E} \left( \nabla \ell(X - m) \right)
    def shortfall_risk_jac(self, m):
        m = self._check_argument(m)
        x_minus_m = np.subtract(self.x, m)
        
        sgn = np.sign(x_minus_m)
        
        sgn_plus = np.maximum(sgn, 0.)
        sgn_minus = -np.minimum(sgn, 0.)
        
        diff1 = sgn_plus + 0.5 * sgn_minus
        diff1_mean = diff1.mean(axis=0)
        
        for i in xrange(self.dim):
            partial_i = 0.
            x_minus_m_i = x_minus_m[:, i]
            for j in xrange(self.dim):
                if i != j:
                    ind = np.sign(x_minus_m_i + x_minus_m[:, j])
                    pos_ind = np.maximum(ind, 0.)
                    pos_neg = -np.minimum(ind, 0.)
                    
                    partial_i += pos_ind.mean() + 0.5 * pos_neg.mean()
                    
            diff1_mean[i] += partial_i
            
        return diff1_mean


    

In [ ]:

    

loss3 = LossFunction3(loss_and_profit, 0.)


    

In [ ]:

    

RESULTS = pd.DataFrame(index=members_pos.index)


    

In [ ]:

    

from scipy.optimize import minimize

methods = ['SLSQP']
maxiter = 100000


    

In [ ]:

    

guess = np.maximum(np.amax(loss_and_profit, axis=0), 0.)

cont = loss3.ineq_constraint(guess) >= 0.
while cont:
    guess = 0.99 * guess
    tmp = loss3.ineq_constraint(guess)
    cont = tmp >= 0.

print loss3.ineq_constraint(guess)


    

In [ ]:

    

cons = ({'type': 'ineq',
         'fun' : lambda x: loss3.ineq_constraint(x),
         'jac' : lambda x: loss3.ineq_constraint_jac(x)})
    
bounds = [(0, None) for _ in xrange(loss3.dim)]

res = minimize(loss3.objective, guess, 
               jac=loss3.objective_jac, 
               constraints=cons, 
               method='SLSQP',
               bounds=bounds,
               options={'maxiter': maxiter, 'disp': True})

RESULTS[r'$\nu = %i$ bounds' % NU] = res.x
print loss3.ineq_constraint(res.x)


    

In [ ]:

    

dir_ = './results'

if not os.path.exists(dir_):
    os.makedirs(dir_)


    

In [ ]:

    

RESULTS.to_csv(os.path.join(dir_, 'l3_nu_eq_%i.csv' % NU))


    

In [ ]:

    

RESULTS