clear all
betavals = (0.1:0.1:3);
l = 10;
E = 4*l*(l-1);
Tmaxvals = [15 25];
delta = zeros(length(betavals), length(Tmaxvals));
bhorfixed = (rand(l, l-1)-0.5);
bvertfixed = (rand(l-1, l)-0.5);
bnodesfixed = (rand(l, l)-0.5);
bhorfixed = (rand(l, l-1));
bvertfixed = (rand(l-1, l));
bnodesfixed = (rand(l, l));
for iter = 1:length(betavals)
	Tmax = max(Tmaxvals); 
	beta = betavals(iter); 
    
	% uniform random initialization of theta's in [0,beta]
	bhor = bhorfixed*beta;
	bvert = bvertfixed*beta;
	bnodes = bnodesfixed*beta;

	% The following four matrices contain the messages (in log-likelihood format) 
    % such that hhorright(i,j) = (1/2)log(nu_{(i,j)->(i,j+1)}(+1)/nu_{(i,j)->(i,j+1)}(-1))
	hhorright = zeros(l, l-1);  % hhorright(i,j) is the message from (i,j) to (i,j+1)
	hhorleft = zeros(l, l-1);   % hhorleft(i,j) is the message from (i,j+1) to (i,j)
	hvertup = zeros(l-1, l);    % hvertup(i,j) is the message from (i+1,j) to (i,j)
	hvertdown = zeros(l-1, l);  % hvertdown(i,j) is the message from (i,j) to (i+1,j)

	% The following matrices store the updated messages, 
	hhorrightnew = zeros(l, l-1);
	hhorleftnew = zeros(l, l-1);
	hvertupnew = zeros(l-1, l);
	hvertdownnew = zeros(l-1, l);
    
    % for each instance of Tmax
	for t = 1:Tmax

		%looping over edges with 3 inputs
		for i = 2:(l-1)
			for j = 2:(l-1)
				hhorrightnew(i, j) = bnodes(i, j) + sum(atanh([ tanh(bhor(i,j-1))*tanh(hhorright(i, j-1)) tanh(bvert(i-1,j))*tanh(hvertdown(i-1, j)) tanh(bvert(i, j))*tanh(hvertup(i, j))]));
				hvertdownnew(i, j) = bnodes(i, j) + sum(atanh( [ tanh(bvert(i-1,j))*tanh(hvertdown(i-1, j)) tanh(bhor(i, j-1))*tanh(hhorright(i, j-1)) tanh(bhor(i, j))*tanh(hhorleft(i,j))]));
			end
        end
		for i = 2:(l-1)
			for j = 1:(l-2)
				hhorleftnew(i, j) = bnodes(i, j+1) + sum(atanh([tanh(bhor(i, j+1))*tanh(hhorleft(i, j+1)) tanh(bvert(i-1, j+1))*tanh(hvertdown(i-1, j+1)) tanh(bvert(i, j+1))*tanh(hvertup(i, j+1))]));
				hvertupnew(j, i) = bnodes(j+1, i) + sum(atanh([tanh(bvert(j+1, i))*tanh(hvertup(j+1, i)) tanh(bhor(j+1, i-1))*tanh(hhorright(j+1, i-1)) tanh(bhor(j+1, i))*tanh(hhorleft(j+1, i))]));
			end
		end

		%looping over edges with 2 inputs
		for i = 2:(l-1)
			hvertupnew(l-1, i) = bnodes(l, i) + sum(atanh([tanh(bhor(l, i-1))*tanh(hhorright(l, i-1)) tanh(bhor(l, i))*tanh(hhorleft(l, i))]));
			hvertdownnew(1, i) = bnodes(1, i) + sum(atanh([tanh(bhor(1, i-1))*tanh(hhorright(1, i-1)) tanh(bhor(1, i))*tanh(hhorleft(1, i))]));
            		hhorrightnew(i, 1) = bnodes(i, 1) + sum(atanh([tanh(bvert(i-1, 1))*tanh(hvertdown(i-1, 1)) tanh(bvert(i, 1))*tanh(hvertup(i, 1))]));
			hhorleftnew(i, l-1) = bnodes(i, l) + sum(atanh([tanh(bvert(i-1,l))*tanh(hvertdown(i-1, l)) tanh(bvert(i, l))*tanh(hvertup(i, l))]));

        end
		for i = 2:(l-1)
			hvertdownnew(i, 1) = bnodes(i, 1) + sum(atanh([tanh(bvert(i-1, 1))*tanh(hvertdown(i-1, 1)) tanh(bhor(i, 1))*tanh(hhorleft(i, 1))]));
			hvertdownnew(i, l) = bnodes(i, l) + sum(atanh([tanh(bvert(i-1, l))*tanh(hvertdown(i-1, l)) tanh(bhor(i, l-1))*tanh(hhorright(i-1, l-1))]));
			hhorrightnew(1, i) = bnodes(1, i) + sum(atanh([tanh(bhor(1, i-1))*tanh(hhorright(1, i-1)) tanh(bvert(1, i))*tanh(hvertup(1, i))]));
			hhorrightnew(l, i) = bnodes(l, i) + sum(atanh([tanh(bhor(l, i-1))*tanh(hhorright(l, i-1)) tanh(bvert(l-1,i))*tanh(hvertdown(l-1, i))]));
        end
		for i = 1:(l-2)
			hvertupnew(i, 1) = bnodes(i+1, 1) + sum(atanh([tanh(bvert(i+1,1))*tanh(hvertup(i+1, 1)) tanh(bhor(i+1, 1))*tanh(hhorleft(i+1, 1))]));
			hvertupnew(i, l) = bnodes(i+1, l) + sum(atanh([tanh(bvert(i+1, l))*tanh(hvertup(i+1,l)) tanh(bhor(i+1, l-1))*tanh(hhorright(i+1, l-1))]));
			hhorleftnew(1, i) = bnodes(1, i+1) + sum(atanh([tanh(bhor(1, i+1))*tanh(hhorleft(1, i+1)) tanh(bvert(1, i+1))*tanh(hvertup(1, i+1))]));
			hhorleftnew(l, i) = bnodes(l, i+1) + sum(atanh([tanh(bhor(l, i+1))*tanh(hhorleft(l, i+1)) tanh(bvert(l-1, i+1))*tanh(hvertdown(l-1, i+1))]));
		end

		%Edges with single input
		hvertdownnew(1, 1) = bnodes(1, 1) + atanh(tanh(bhor(1, 1))*tanh(hhorleft(1, 1)));
		hvertdownnew(1, l) = bnodes(1, l) + atanh(tanh(bhor(l-1,1)));
        hvertupnew(l-1, 1) = bnodes(l, 1) + atanh(tanh(bhor(l, 1))*tanh(hhorleft(l, 1)));
        hvertupnew(l-1, l) = bnodes(l, l) + atanh(tanh(bhor(l, l-1))*tanh(hhorright(l, l-1)));

        hhorrightnew(1, 1) = bnodes(1, 1) + atanh(tanh(bvert(1, 1))*tanh(hvertup(1, 1)));
        hhorrightnew(l, 1) = bnodes(l, 1) + atanh(tanh(bvert(l-1, 1))*tanh(hvertdown(l-1, l)));
        hhorleftnew(1, l-1) = bnodes(1, l) + atanh(tanh(bvert(1, l-1))*tanh(hvertup(1, l-1)));
        hhorleftnew(l, l-1) = bnodes(l, l) + atanh(tanh(bvert(l-1, l))*tanh(hvertdown(l-1, l)));
    
        % convert from log-likelihoods to probabilities
		nuvertdownnew = exp(2*hvertdownnew)./(1+exp(2*hvertdownnew));
		nuvertupnew   = exp(2*hvertupnew)./(1+exp(2*hvertupnew));
		nuhorrightnew = exp(2*hhorrightnew)./(1+exp(2*hhorrightnew));
		nuhorleftnew  = exp(2*hhorleftnew)./(1+exp(2*hhorleftnew));

		nuvertdown = exp(2*hvertdown)./(1+exp(2*hvertdown));
		nuvertup   = exp(2*hvertup)./(1+exp(2*hvertup));
		nuhorright = exp(2*hhorright)./(1+exp(2*hhorright));
		nuhorleft  = exp(2*hhorleft)./(1+exp(2*hhorleft));
        
        % store the new messages as current messages
		hvertdown = hvertdownnew;
		hvertup   = hvertupnew;
		hhorleft  = hhorleftnew;
		hhorright = hhorrightnew;
		
        % compute the difference delta
		if(sum(t==Tmaxvals))
			iter2 = find(t==Tmaxvals);
			delta(iter, iter2) = (sum(sum(abs(nuvertdownnew-nuvertdown))) + sum(sum(abs(nuvertupnew-nuvertup))) + sum(sum(abs(nuhorrightnew-nuhorright))) + sum(sum(abs(nuhorleftnew-nuhorleft))))/E;
		end
	end
end
% plot resulting delta
semilogy(betavals, delta)
grid on
