function [prob,q,tr] = BR_MCMC_query(KB,query,num_initial,near,max_samples)
%

NUM_INITIAL = num_initial;
NEAR = near;

prob = 0;

%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);
query_clause = query(1).clauses;
samps = rand(NUM_INITIAL,n)<0.5;
true_counts = zeros(m,1);
true_query = 0;
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
    true_query = true_query + BR_sat(query_clause,samps(s,:));
end
q = true_counts/NUM_INITIAL;
tr = q';
num_samples = NUM_INITIAL;
total_samps = num_samples;
while num_samples<max_samples
    total_samps = total_samps + 1;
    below = q<p;
    diffs = abs(q-p);
    [vals,indexes] = sort(diffs,'descend');
    max_dif = vals(1);
    max_index = indexes(1);
    done = 0;
    while done==0
        samp = rand(1,n)<0.5*ones(1,n);
        s_val = BR_sat(KB(max_index).clauses,samp);
        if s_val==0&q(max_index)>p(max_index)...
                |s_val==1&q(max_index)<p(max_index)
            done = 1;
            FOUND = 1;
        end
    end
%    FOUND = 1;
%    c = 0;
%    while FOUND==1&c<m
%        c = c + 1;
%        c_val = BR_sat(KB(c).clauses,samp);
%        if (c_val==1)&(q(c)<p(c))...
%                |abs(q(c)-p(c))<NEAR
%            FOUND = FOUND*1;
%        elseif (c_val==0)&(q(c)>p(c))...
%                |abs(q(c)-p(c))<NEAR
%            FOUND = FOUND*1;
%        else
%            FOUND = FOUND*0;
%        end
%    end
    if FOUND==1
        num_samples = num_samples + 1;
        for c = 1:m
            true_counts(c) = true_counts(c)...
                + BR_sat(KB(c).clauses,samp);
        end
        true_query = true_query + BR_sat(query_clause,samp);
        q = true_counts/num_samples;
        q'
        tr = [tr;q'];
    end
%    plot(1,num_samples,'ro');
%    plot(2,total_samps,'g*');
%    drawnow
end
prob = true_query/num_samples;
tch = 0;
