function [H,loss,HGy,nmes,ct]= miqpex_HFull_indicator(Gy,Gdy,Juu,Jud,Wd,Wn,nm)
% Enter and optimize a mixed integer quadratic programming problem
%
% This function fills in the data structures for the mixed integer
% quadratic program:
%
%       Maximize
%        obj: x1 + 2 x2 + 3 x3 + x4
%             - 0.5 ( 33x1*x1 + 22*x2*x2 + 11*x3*x3
%                    -  12*x1*x2 - 23*x2*x3 )
%       Subject To
%        c1: - x1 + x2 + x3 + 10x4  <= 20
%        c2: x1 - 3 x2 + x3         <= 30
%        c3:       x2       - 3.5x4  = 0
%       Bounds
%        0 <= x1 <= 40
%        0 <= x2
%        0 <= x3
%        0 <= x4 <= 3
%       Integers
%         x4
%       End

%  ------------------------------------------------------------------------
%  File: examples/src/matlab/miqpex1.m
%  Version 12.0
%  ------------------------------------------------------------------------
%   (C) Copyright ILOG S.A.S. 2008, 2009
%   All Rights Reserved.
%   Permission is expressly granted to use this example in the
%   course of developing applications that use ILOG products.
%  ------------------------------------------------------------------------

try
    % Initialize the CPLEX object
    cplex = Cplex('miqpex_HFull_indicator');
    cplex.DisplayFunc = []; %@disp;
    cplex.Param.clocktype.Cur = 0;
    cplex.Param.timelimit.Cur = 5400;
    % Fill in the data for the problem using populatebyrow
    
    populatebyrow(cplex,Gy,Gdy,Juu,Jud,Wd,Wn,nm);
    
    % Optimize the problem
    cplex.solve();
    
    % Write the solution
    fprintf('\nSolution status = %s \n',cplex.Solution.statusstring);
    fprintf('Solution value = %f \n',cplex.Solution.objval);
    disp('Values =');
    %     disp(cplex.Solution.x);
    disp('Slacks =');
    %     disp(cplex.Model.rhs-cplex.Solution.ax);
    xn = cplex.Solution.x;
    ct = cplex.Solution.time;
    ny = size(Gy,1);
    nd = size(Gdy,2);
    H = reshape(xn(1:ny*nu),ny,nu)';
    C = inv(Juu)*Jud;
    F = [(Gy*C - Gdy)*Wd Wn];
    loss = norm(sqrtm(Juu)*inv(H*Gy)*(H*F),'fro')^2/2;
    HGy = H*Gy;
    no = (1:ny)';
    nsy = xn(nu*ny+1:end).*no;
    nsy(find(abs(nsy) < 1e-3)) = 0;
    nmes = nonzeros(nsy);
catch m
    disp(m.message);
end
    function populatebyrow(cplex,Gy,Gdy,Juu,Jud,Wd,Wn,nmeas)
        C = inv(Juu)*Jud;
        B = [(Gy*C - Gdy)*Wd Wn];
        Juuhalf = sqrtm(Juu);
        [ny,nu]= size(Gy);
        nd = size(Gdy,2);
        Jn = []; Fn = []; GnyT = [];
        for i = 1:nu
            GnyT = daug(GnyT,Gy');
            Jn = [Jn;Juuhalf(i,:)'];
            Fn = daug(Fn,B);
        end
        Jn = Jn;
        Fxb = Fn*Fn';
        c = zeros(size(Fxb,1),1);
        F    = 2*daug(Fxb,zeros(ny,ny));
        c    = zeros(size(F,1),1);
        GnewyT = [GnyT zeros(nu*nu,ny)];
        LconB = [zeros(1,ny*nu) ones(1,ny)];
%         A = [A1;GnewyT;LconB];
        x_L  = [-Inf*ones(nu*ny,1);zeros(ny,1)];
        x_U  = [Inf*ones(nu*ny,1);ones(ny,1)];
        for xy = 1:ny*nu
            str(1,xy)='C';
        end
        for xy = 1:ny
            str(1,ny*nu+xy)='B';
        end
        zz = 1;
        
        A = [GnewyT;LconB];
        b_L  = [Jn;nmeas];
        b_U  = [Jn;nmeas];
        cplex.Model.sense = 'minimize';
        cplex.addCols(c, [], x_L,x_U,char(str));
        cplex.Model.ctype = char(str);
        for i = 1:ny
            cplex.addIndicators(2*ny+i,1,[zeros(1,i-1),1,zeros(1,2*ny-i),zeros(1,ny)]','E',0);
            cplex.addIndicators(2*ny+i,1,[zeros(1,ny+i-1),1,zeros(1,ny-i),zeros(1,ny)]','E',0);
        end
        cplex.addRows(b_L, A, b_U);
        cplex.Model.Q = F;
%                 cplex.writeModel('test.lp')
%                 type test.lp
    end
end


% cplex = Cplex()
% cplex.addCols(c, [], x_L,x_U,char(str));
% cplex.Model.ctype = char(str)
% for i = 1
%     cplex.addIndicators(ny+1,1,[zeros(1,i-1),1,zeros(1,2*ny-i),zeros(1,ny)]','E',0)
%     cplex.addIndicators(ny+1,1,[zeros(1,ny+i-1),1,zeros(1,ny-i),zeros(1,ny)]','E',0)
% end
% cplex.writeModel('test.lp')
% type test.lp
%
%
% cplex = Cplex()
% cplex.addCols([1;2;3;4])
% cplex.Model.ctype = 'CCBB'
% cplex.addIndicators(3,1,[1,0,0,0]','E',0)
% cplex.addIndicators(4,0,[0,1,zeros(1,2)]','E',0)
% cplex.writeModel('test.lp')
% type test.lp
