function [W, rate]= HBOptimalWeight(A,H)

% % % % A : adjacency matrix
% % % % % % H: Hessian at the optimal point
% % % Add path
% root = 'C:\Program Files\MATLAB\R2010b\Yalmip';
% addpath(strcat((root),'\yalmip'));
% addpath(strcat((root),'\yalmip\extras'));
% addpath(strcat((root),'\yalmip\demos'));
% addpath(strcat((root),'\yalmip\solvers'));
% addpath(strcat((root),'\yalmip\modules'));
% addpath(strcat((root),'\yalmip\modules\parametric'));
% addpath(strcat((root),'\yalmip\modules\moment'));
% addpath(strcat((root),'\yalmip\modules\global'));
% addpath(strcat((root),'\yalmip\modules\sos'));
% addpath(strcat((root),'\yalmip\operators'));
% 
% % % % % 
% % % % SEDUMI 
% % % % 
% root_sedumi = 'C:\Program Files\MATLAB\R2010b\seDumi';
% 
% addpath(strcat((root_sedumi),'\SeDuMi_1_3'));


netSize = size(A,1);
% 

% % here we make sure that A has diagonal 1 element
A= max(A,eye(netSize,netSize));


% constructing P matrix which is a (n-1*n) orthonormal matrix all the
% columns are orthogonal to v1 vector. i.e P*v1=0; and v1 is given by
% H^{-1/2}* 1_vec

v1 =  H^(-.5)*ones(netSize,1);

% construct P_n as linearly indpendent (non-singular) with first column of
% ones
% You can try using MATLAB's qr (or orth) to perform the orthogonalization.

P = v1*ones(1,netSize);
for k=2:netSize,
    P(k-1,k) = P(k-1,k)+1;
end;
P = GramShmidt(P);
P = P(:,2:netSize);


%cast the problem like this min t s.t lambda2 I <= P' Q P <= lambda2 * t* I

% defining the variables
 t = sdpvar (1,1);

Q = sdpvar(netSize,netSize);
% define Q (laplacian type) as a positive semidefinite(stmmetric) and Q1=0; also impose sparsity patterns for given network 


F =  [ t>= 0];
%sparsity pattern

for i=1:netSize,
    for j=i:netSize,
        if A(i,j)==0
            F = F + set(Q(i,j)==0)+ set(Q(j,i)==0);
        end;
    end;
end;
%
% F = F + [Q==Laplacian];
F = F+ [Q >= 0] + [Q*v1 == zeros(netSize,1)];
% constrain corresponds to maximum eigenvale
%  F = F+[ Q <= eye(netSize)];

F = F+[ t*eye(netSize-1)<= P'*Q*P <= eye(netSize-1)];
%  diagnostics = solvesdp(F,t);
diagnostics = solvesdp(F, -t);
Omega = double(Q);
W = H^(-1/2)*Omega*H^(-1/2);
rate = double(t);


