%  Input:  an nxn matrix A whose non-zero entries correspond to the edge weights
%  of an Ising model, a vector w of singelton weights, and the number of
%  iterations of coordinate ascent to perform, its
%
%  Ouput:  approximate log of the partition function and approximate
%  mean field marginals, q

function [logZ, q] = meanfieldis(A, w, its)

[m,n] = size(A);

%Initialize the q's to a random probability distribution
q = rand(n,2);
for i = 1:n
    zi = q(i,1) + q(i,2);
    q(i, 1) = q(i,1)/zi;
    q(i, 2) = q(i,2)/zi;
end

%main loop
for t = 1:its
    %perform coordinate ascent for q_i(1) and q_i(-1)
    for i = 1:n
        q(i, 1) = exp(-w(i));
        q(i, 2) = exp(w(i));
        for j = 1:n
            if A(i,j) == 1
                q(i, 1) = q(i,1) * exp(-A(i,j)*q(j,2) + A(i,j)*q(j,1));
                q(i, 2) = q(i,2) * exp(A(i,j)*q(j,2) - A(i,j)*q(j,1));
            end
        end
        zi = q(i,1) + q(i,2);
        q(i, 1) = q(i,1)/zi;
        q(i, 2) = q(i,2)/zi;
    end
end

%Compute the estimate of the partition function
logZ = 0;
for i = 1:n
    logZ = logZ - log(q(i,1)^q(i,1)) - log(q(i,2)^q(i,2));
    for j = i+1:n
        if A(i,j) == 1
            logZ = logZ + (q(i,1)*q(j,1) + q(i,2)*q(j,2))*A(i,j) - (q(i,1)*q(j,2) + q(i,2)*q(j,1))*A(i,j);
        end
    end
end