function [S] = multigraphical(x,C,alpha)
%
% function [S] = multigraphical(x,C,alpha)
%
% x is a binary matrix of size TxD of T samples from a
% D-dimensional Ising model
%
% C and alpha are scalar regularization values
%
% S contains the graphical model output estimate
%
% Warning, this is slow code. For speed and accuracy, we recommend
% using MOSEKs quadprog instead of the default Matlab quadprog.
%
% Example usage: 
% T=20; D=10; x=round(rand(T,D)); multigraphical(x,10,0.1)
%
% Copyright 2009 by Tony Jebara

maxiters = min(log(1/0.005)/log(min(1+1/alpha,2)),1000);

T=size(x,1);
D=size(x,2);
S=zeros(D,D);
x = round((x - min(min(x))) / (max(max(x))-min(min(x))));

y=2*x-1;
lambdas=zeros(T*D,1);
A=zeros(D,T*D);
for m=1:D
  for t=1:T
    A(m,(m-1)*T+t)=y(t,m);
  end
end
bzeros=zeros(D,1);



% Compute the objective function (made to be zero at zero)
J = sum(sum(lambdas))+D*(D-1)/2*log(alpha+1);
for m=1:D
for d=(m+1):D
  tot1=0.0; tot2=0.0;
  for t=1:T
    tot1 = tot1+lambdas((m-1)*T+t)*y(t,m)*x(t,d);
    tot2 = tot2+lambdas((d-1)*T+t)*y(t,d)*x(t,m);
  end
  stuff = tot1^2+tot2^2;
  J = J - log(alpha*exp(-0.5*stuff)+1)-0.5*stuff;
end
end
fprintf('Iteration=0 J=%e\n',J);
oldJ=J;

going = 1;
iters = 1;
while (going)

  old = lambdas;
  qpchange = 0;
  vlB = zeros(T*D,1);
  vuB = C*ones(T*D,1);
  hess = zeros(T*D,T*D);
  grad = -ones(T*D,1);
 
  for m=1:D
    for d=(m+1):D
      Hmd=zeros(T*D,T*D);
      for n=1:D
        for t=1:T
          for o=1:D
            for u=1:T
              if ((m==n) && (n==o)) Hmd((n-1)*T+t,(o-1)*T+u) = Hmd((n-1)*T+t,(o-1)*T+u) + y(t,m)*y(u,m)*x(t,d)*x(u,d); end
              if ((d==n) && (n==o)) Hmd((n-1)*T+t,(o-1)*T+u) = Hmd((n-1)*T+t,(o-1)*T+u) + y(t,d)*y(u,d)*x(t,m)*x(u,m); end
            end
          end
        end
      end
      scal = 1/(alpha*exp(-0.5*old'*Hmd*old)+1);
Gmd  = tanh(0.5*log(alpha*exp(-0.5*old'*Hmd*old)))/(2*log(alpha*exp(-0.5*old'*Hmd*old)));
grad = grad + scal*Hmd*old - (Hmd + Gmd*Hmd*old*old'*Hmd)*old;
      hess = hess + Hmd + Gmd*Hmd*old*old'*Hmd;
    end
  end
			      hess = 0.5*(hess+hess');
  lambdas = quadprog(hess,grad,[],[],A,bzeros,vlB,vuB,lambdas,[]);
  if (size(lambdas,1)<1) fprintf('Quadprog made no progress.\n'); lambdas=old; going=0; end;

  % Compute the objective function (made to be zero at zero)
  J = sum(sum(lambdas))+D*(D-1)/2*log(alpha+1);
  for m=1:D
    for d=(m+1):D
      tot1=0; tot2=0;
      for t=1:T
        tot1 = tot1+lambdas((m-1)*T+t)*y(t,m)*x(t,d);
        tot2 = tot2+lambdas((d-1)*T+t)*y(t,d)*x(t,m);
      end
      stuff = tot1^2+tot2^2;
      S(m,d) = 1/(1+alpha*exp(-0.5*stuff));
      S(d,m) = 1/(1+alpha*exp(-0.5*stuff));
      J = J - log(alpha*exp(-0.5*stuff)+1)-0.5*stuff;
    end
  end
  fprintf('Iteration=%d J=%e\n',iters,J);
  if (min(J,oldJ)>0) if ((J/oldJ)<1.00000001) going = 0; fprintf('Insufficient progress.\n'); end; end;
  oldJ = J;
  % Check if number of iterations was exceeded
  if (iters>maxiters)
    going = 0;
    fprintf('Exceeded max iterations stoping early.\n');
  end
  iters = iters+1;
end

