%% This is the code for the algorithm - simulator system

% Basic definitions %
m = 10;                     % Number of Features
n = 4;                      % Number of products (J)

k = 0.7;                            % Error parameter alpha'
J = n;                              % Number of choices for a question
kk = (J*k - 1)/(J-1);                % Error paramter converted to weight parameter i.e. alpha

V = [];
P = [];
ACC = zeros(m,1);
M = 2;                      % Number of uniform distributions which approximate the gaussian
wm = [];
wm(1) = 0.3;
wm(2) = 0.4;
wm(3) = 0.3;

Cm = [];                    % The coefficients for every uniform component
Cm(1) = 0.05;
Cm(2) = 0.10;
Cm(3) = 0.15;

utrue = 10*rand(m,1);          % Freeze the true partworth vector

PF = randi(10,m,n);            % Generate the feature-product maxtrix

Q = 10;                     % Limit on number of questions
xmax = zeros(m,Q);

%% Zero'th stage %%



for q=1:Q

    % Get the product feature vector which has the maximum utility
    dummy = max_utility(utrue, PF,m,n);
    xmax(:,q) = dummy ;

    % Generate the inequalities arising from deterministic utility
    % Format is (u')A <= b

    for y=1:n
        A(:,y) = -(xmax(:,q) - PF(:,y)) ;
        b(y) = 0 ;
    end

    % Generate the gaussian approxation inequalities
    % Append the same to the existing inequality matrix

    % I think we can directly incorporate them as constraints in the
    % analytical center optimization problem

    %% Generate all inequalities which define polytope %%
    % Store information up to the qth question to leverage the subsets %

    % Generate subsets

    Setq = [];
    Size = 0;
    Subset = {};

    % Creating subsets

    for i=1:q
        Setq(i)=i;
        Sq{i} = combntns(1:q,i);
        Size = Size + size(Sq{i},1);
    end

    SS = [];

    for i=1:q
        SS = Sq{i};
        for j=1:size(Sq{i},1)
            Subset = [Subset ; SS(j,:)];
        end
    end


    %% Analytical center and ellipsoids %%

    for k=1:M

        for i=1:Size

            s = size(Subset{i},2);
            ws(i) = (kk^s)*((1-kk)^(q-s));

            % Generating the respondent's matrix %
            % Get into familiar notation i.e A*x1 <= b %
            bb = b';
            AA = A';
            AA = ones(n,m) + AA;

%             AA = ((-1)^m)*AA;


            lb = (utrue - ((Cm(k)*(1*(ones(m,1))))));
            ub = (utrue + ((Cm(k)*(1*(ones(m,1))))));

            % Analytic center - My Code (V2) %
                cvx_begin
                    variable x1(m);
                    minimize -sum(log(eps + (AA*x1 - bb)));
                    subject to
                        sum(x1) <= 100;
                        x1 >= ones(m,1);
%                         for t=1:m;
%                               lb(t)  <= x1(t) <=  ub(t) ;
%                         end

%                         (ones(1,m))*x1 == 10;

                cvx_end

            x = x1;

            % Ellipsoid - Using Hessian computation %

            AC = (wm(k)*ws(i)).*(x) ;

%             H = g1(AC, AA, bb) ;

            syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
            % syms x1 x2 x3
            Vv = [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10];
            % V = [x1;x2;x3];
            f = -sum(log(AA*Vv-bb));
            % HH = jacobian(gradient(f));
            HH = hessian(f,Vv);
             H = (subs(HH, [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10], AC ));


            Ht = inv(H);
%             Ht = H;
            U = sqrtm(Ht);

            uhat = x1;
            [EV,D] = eig(H);
            EV1 = EV(:,n);
            VSM = EV1';
            V = [V;VSM];

            P = [P; wm(k)*ws(i)];

            ACC = ACC + AC ;

        end


    end

    %% Evaluating New Data Points %%

    Pi = diag(P);
    ES = (V')*(Pi)*(V) ;

    [EVect, EVal] = eig (ES) ;

    % for i=1:J/2
    i=0;
        AX = EVect(:,n-i);
        b0 = [b; dot(AX,x1)];
        A0 = [AA; AX];
        XJ = linsolve(A0,b0);

        JJ = horzcat(JJ, XJ);


    % Remind to update PF %

end