function x = comb_gibbs_step(x, theta)

x = comb_gibbs_helper(x, theta);
temp = comb_gibbs_helper(x(end:-1:1, end:-1:1), theta);
x = temp(end:-1:1, end:-1:1);



function x = comb_gibbs_helper(x, theta)

[M, N] = size(x);
phi_pos = zeros(M, N);

% compute unary potentials in the tree graph
for j = 1:2:N
  phi_pos(2:end-2, j) = phi_pos(2:end-2, j) + theta * x(2:end-2, j+1);
end
for j = 3:2:N
  phi_pos(2:end-2, j) = phi_pos(2:end-2, j) + theta * x(2:end-2, j-1);
end
for j = 2:2:N
  phi_pos(1, j) = phi_pos(1, j) + theta * x(2, j);
end
for j = 1:2:N
  phi_pos(M-1, j) = phi_pos(M-1, j) + theta * x(M, j);
end


% sum product messages
[m_up_neg, m_up_pos, m_down_neg, m_down_pos, m_left_neg, m_left_pos, m_right_neg, m_right_pos] = ...
    comb_sum_product(-phi_pos, phi_pos, theta);

% sample up the leftmost column
for i = M-1:-1:2
  odds = 2 * phi_pos(i, 1) + m_down_pos(i, 1) - m_down_neg(i, 1);
  if i < M-1
    odds = odds + 2 * theta * x(i+1, 1);
  end
  x(i, 1) = sample_from_odds(odds);
end

% sample right in the top row
for j = 1:N
  odds = 2 * phi_pos(1, j) + m_left_pos(1, j) - m_left_neg(1, j);
  
  if j == 1
    odds = odds + 2 * theta * x(2, 1);
  elseif mod(j, 2) == 1
    odds = odds + m_up_pos(1, j) - m_up_neg(1, j);
  end
  
  if j > 1
    odds = odds + 2 * theta * x(1, j-1);
  end
  
  x(1, j) = sample_from_odds(odds);
end

% sample down the remaining columns
for j = 3:2:N
  for i = 2:M-1
    odds = 2 * phi_pos(i, j) + m_up_pos(i, j) - m_up_neg(i, j) + 2 * theta * x(i-1, j);
    x(i, j) = sample_from_odds(odds);
  end
end



function x = sample_from_odds(odds)

prob = 1 / (1 + exp(-odds));
if rand < prob
  x = 1;
else
  x = -1;
end

