function [prob,q,tr] = BR_MC_query_MH(KB,num_initial,near,...
    max_samples)
%

NUM_INITIAL = num_initial;
NEAR = near;

%figure(1);
%clf
%hold on

m = length(KB);
p = zeros(m,1);
for c = 1:m
    p(c) = KB(c).prob;
end

vars = BR_vars(KB,[]);
n = length(vars);
prob = zeros(n,1);
samps = rand(NUM_INITIAL,n)<0.5;
true_counts = zeros(m,1);
true_query = zeros(n,1);
q = zeros(m,1);
for s = 1:NUM_INITIAL
    for c = 1:m
        true_counts(c) = true_counts(c) + BR_sat(KB(c).clauses,samps(s,:));
    end
    for k = 1:n
        true_query(k) = true_query(k) + BR_sat([k],samps(s,:));
    end
end
q = true_counts/NUM_INITIAL;
tr = q';
num_samples = NUM_INITIAL;
total_samps = num_samples;
wb = waitbar(0,'MC');
samp = samps(end,:);
while num_samples<max_samples&norm(q-p)>NEAR
    waitbar(num_samples/max_samples);
    total_samps = total_samps + 1;
    k = randi(n);
    proposed = samp;
    proposed(k) = mod(proposed(k)+1,2);
    tc = true_counts;
    ts = true_counts;
    for c = 1:m
        tc(c) = tc(c) + BR_sat(KB(c).clauses,proposed);
        ts(c) = ts(c) + BR_sat(KB(c).clauses,samp);
    end
    qc = tc/(num_samples+1);
    qs = ts/(num_samples+1);
    alpha = min(1,-norm(p-qc)/-norm(p-qs));
    if rand<alpha
        samp = proposed;
        true_counts = tc;
    else
        true_counts = ts;
    end
    samps(end+1,:) = samp;
    num_samples = num_samples + 1;
    for k = 1:n
        true_query(k) = true_query(k) + BR_sat([k],samp);
    end
    q = true_counts/num_samples;
    q'
    tr = [tr;q'];
    %    plot(1,num_samples,'ro');
%    plot(2,total_samps,'g*');
%    drawnow
end
close(wb);
prob = true_query/num_samples;
tch = 0;
