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

NUM_INITIAL = num_initial;
NEAR = near;
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);
err = zeros(max_samples,1);
samps = rand(NUM_INITIAL,n)<0.5;
for s = 1:NUM_INITIAL
    for c = 1:m
        clause = KB(c).clauses;
        if length(clause)==1
            samps(s,abs(clause)) = ~(clause<0);
        end
    end
end
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;
err(1) = norm(p-q);
tr = q';
num_samples = NUM_INITIAL;
total_samps = num_samples;
wb = waitbar(0,'MC');
while num_samples<max_samples&norm(q-p)>NEAR
    waitbar(num_samples/max_samples);
    total_samps = total_samps + 1;
    FOUND = 0;
    count = 0;
    while FOUND==0
        count = count + 1;
        if mod(count,100)==0
            tch = 0;
        end
        samp = rand(1,n)<0.5*ones(1,n);
        for c = 1:m
            clause = KB(c).clauses;
            if length(clause)==1
                samp(abs(clause)) = ~(clause<0);
            end
        end
        c = 0;
        qt = 0*q;
        t_counts = true_counts;
        for c = 1:m
            t_counts(c) = t_counts(c) + BR_sat(KB(c).clauses,samp);
        end
        qt = t_counts/(num_samples+1);
        if norm(qt-p)<norm(q-p)
            FOUND = 1;
        end
        if count==100
            FOUND = 1;
        end
    end
    num_samples = num_samples + 1;
    samps(num_samples,:) = samp;
    [num_samples,max_samples]
    count = 0;
    for c = 1:m
        true_counts(c) = true_counts(c)...
            + BR_sat(KB(c).clauses,samp);
    end
    for k = 1:n
        true_query(k) = true_query(k) + BR_sat([k],samp);
    end
    q = true_counts/num_samples;
    err(num_samples) = norm(q-1);
%    q'
    tr = [tr;q'];
end
close(wb);
prob = true_query/num_samples;
tch = 0;
