function [normpf, normdf, X,Z,theta] = pcthetaimp( n,indsE,indsEc);
%        [normpf, normdf, X,Z,theta] = pcthetaimp( n,indsE,indsEc);

%%   THETA solves the SDP version of the Lovacz-theta function problem
%%      max trace EX    s.t.  trace(X)=1, X_ij=0 for all ij in EE, X psd
%%  where E is the matrix of all ones, EE is the set of edges in the graph
%%   graph is  G=(V,EE) with n=|V| and m=|EE| and mc=(n(n-1)/2)-m;
%  input: n=|V|
%         indsE ... edge set
%         indsEc ... edge set complement
% output: theta ... weighted theta function of graph
%         X   ... primal variables
%         y   ... dual variables
% solves: max trace(EX) s.t. X psd, X(i,j) =0, (ij) nonedge, trace(X)=1
%         min y0 s.t. y0*eye(n) + sum_ij y(ij)A(ij) - E psd
%         E(ij)= (e_i*e_j' + e_je_i')/sqrt(2);  (ij nonedge)


[iE,jE] = ind2sub([n,n],indsE);
[iEc,jEc] = ind2sub([n,n],indsEc);

mc = length(indsEc);
m = length(indsE);
index = 1:m;
t0 = clock;                     % take the time
w = ones(n,1);
digits = 8;                     % stop if required accuracy reached

% starting triple (d,x,y,Z) is feasible and in interior
X = eye(n)/n;
x = zeros(mc,1);
d = zeros(n-1,1);
z0 = 2*n;
%z0=1.1*w'*w;
y = zeros(m,1);
Z = z0*speye(n)-ones(n); 

phi = z0;               % initial dual cost
psi = w'* X * w;        % intial primal cost
delta = phi-psi;        % delta controls stopping condition
mu = Z(:)' * X(:) /(2*n);

iter = 0;
text = ['      iter    alphap    alphad    log(gap)  lower      upper'];
disp( text); disp([ iter 1 1 log10(delta) psi  phi ]);

% main iteration
while delta > (abs(phi)+1) * 10^(-digits)  % while duality gap is too large

        iter = iter + 1;        % start a new iteration
        Zi = inv(Z);
        Zi = (Zi + Zi')/2;      % Zi is symmetric

% form the matrix of the system! 
        R=chol(Z);
        iZX = Zi*X;
        iZXs=(iZX+iZX')/2;
        alf2 = trace(iZX);
	A = []; 
        A1=[];
       A2 = [];
    for cnt=1:m, 
       iZMatEX =( Zi(:,iE(cnt))*X(jE(cnt),:)+Zi(:,jE(cnt))*X(iE(cnt),:) )/sqrt(2); 
       iZMatEX  = (iZMatEX+iZMatEX')/2;
       alf1 = trace(iZMatEX)/alf2;
       ls = - alf1*iZXs + iZMatEX; % left hand side without projection
       lhs = sqrt(2)*ls(indsE);   
       A = [A lhs];
       ls1 = - alf1*iZXs;
       lhs1 = sqrt(2)*ls1(indsE);   
       A1 = [A1 lhs1];
       ls2 = iZMatEX;
       lhs2 = sqrt(2)*ls2(indsE);   
       A2 = [A2 lhs2];
    end;
    % temp. test
    if norm(A-A')>1e-10,
       ['error in A not symmetric  ' num2str(norm(A-A')) ]
       A=(A+A')/2;
       if min(eig(A))<0,
          'error A not p.d.'
       end
     else
       A=(A+A')/2;
       if min(eig(A))<0,
          'error A not p.d.'
       end
     end

        muiZmX = mu*Zi-X;
        trmuiZmX = trace(muiZmX);
        rs = muiZmX - trmuiZmX*iZX/alf2;
        rhs = sqrt(2)*rs(indsE);

        lowtri = chol(A); 
        %clear A;
	lowtri = sparse(lowtri);
        uptri = lowtri'; 
        ytmp = uptri\rhs; 
	dy = lowtri\ytmp;   % dy !!

        usMatEdy = usMatE(dy,n,m,indsE);
        iZusMatX = Zi*(usMatEdy*X);

        dz0 = (- trace(iZusMatX) + trmuiZmX )/alf2;


        % calculating dx:
        dz0iZX = dz0*iZX;
        iZusMatX= (iZusMatX+ iZusMatX')/2;
        rhsx = muiZmX - dz0iZX - iZusMatX;
        %rhsx = ( rhsx + rhsx')/2;       % symmetrize
        dx = sqrt(2)*rhsx(indsEc);  
        usMatEcdx = usMatEc(dx,n,mc,indsEc);

        % calculating dd:
        rhsdd = muiZmX - dz0iZX -iZusMatX; %change this because it is dag(Vdd) 
        ddiag = diag(rhsdd); 
        dd = ddiag(2:n); 

        dX = usMatEcdx + diag([-sum(dd); dd]); 
        dZ = dz0*eye(n) + usMatEdy;


% find steplengths alphap and alphad
        alphap =  1;

        [dummy,posdef] = chol( X + alphap * dX);
        while posdef > 0,
                alphap = alphap * .8;
		[dummy,posdef] = chol( X + alphap * dX);
               end;
        if alphap < 1, alphap = alphap * .95; end;

        alphad = 1;
        [dummy,posdef] = chol( Z + alphad * dZ);
        while posdef > 0;
                alphad = alphad * .8;
		[dummy,posdef] = chol( Z + alphad * dZ);
                end;
        if alphad < 1, alphad = alphad * .95; end;

% update
        X = X + alphap * dX;
        x = x +  alphap *dx;
        z0 = z0 + alphad *dz0;
        y = y + alphad * dy;
        Z = Z + alphad * dZ;
        mu = .6 * X( :)' * Z( :) / n;
        if alphap + alphad > 1.5, mu = mu * .6; end;
        if alphap + alphad > 1.7, mu = mu * .4; end;

        phi = z0;
        psi = sum(sum(X));
        delta = phi-psi;
%'after delta'
%keyboard

% display current iteration
        disp([ iter alphap alphad log10(delta) psi  phi]);
        end;            % end of main loop

theta = phi;
disp('elapsed time:'); disp(etime(clock, t0));

% evaluate feasibilities
normpf=norm(X(indsEc))+abs(trace(X)-1));
%normdf=norm((diag(y)-L)-Z);
disp('at end of thetaimp - using improved HKM direction')
disp(['norm(X(indsEc))+abs(trace(X)-1)) is: ',num2str(normpf)]);

