function [ap,ap_set,Sp,OK] = BR_solve_atom_probs(KB,num_trials,thresh,step)
% BR_solve_sentence_probs - solve for atom probabilities given sentences
% On input:
%     KB (KB data structure - vector): CNF sentences
%       (k).clauses (vector): disjunction of literals
%       (k).prob (float): probability of sentence
%     num_trials (int): max number of trials to search for atom values
%     thresh (float): distance of sentence prob vector to atom-based 
%                     sentence probs
%     step (float): step in atom value when searching
% On output:
%     ap (pxn vector): atom probabilities; multiple rows indicate that
%                      some atoms have a range of possible values
%     ap_set (axn vector): 1 if atom prob is deterministic; else 0
%     Sp (1xm vector): sentence probs based on atom probs
%     OK (Boolean): OK if consistent assignment in [0,1] to atoms
% Call:
%     KB2(1).clauses = 1;
%     KB2(2).clauses = [-1,2];
%     KB2(1).prob = 1;
%     KB2(2).prob = 1;
%     [ap,a_set,Sp,OK] = BR_solve_atom_probs(KB2,20000,0.00001,0.001);
% Author:
%     T. Henderson
%     UU
%     Summer 2017
%

ZERO_THRESH = 0.01;
DEL_A = 0.001;

syms x

OK = 1;
num_sentences = length(KB);
vars = BR_vars(KB,[]);
num_vars = length(vars);
var_index = zeros(1,max(vars));
for v = 1:max(vars)
    index = find(vars==v);
    if ~isempty(index)
        var_index(v) = index;
    end
end
ap = 0.5*ones(1,num_vars);
Sp = zeros(1,num_sentences);
ap_set = zeros(1,num_vars);
done = 0;
tch = 0;

while done==0
    done = 1;
    for s = 1:num_sentences
        var_indexes1 = abs(KB(s).clauses);
        var_indexes = var_index(var_indexes1);
        indexes1 = var_indexes(find(ap_set(var_indexes)==0));
        indexes0 = [];
        if ~isempty(indexes1)
            for v = 1:length(indexes1)
                indexes0 = [indexes0,find(var_index==indexes1(v))];
            end
        end
        num_free = length(indexes0);
        if num_free==1
            F = BR_sentence2formula(KB(s).clauses,KB(s).prob);
            G = BR_sub_values_v2(F,ap,indexes0,var_indexes1,var_index);
            ap(indexes1) = solve(G,x);
            ap_set(indexes1) = 1;
            done = 0;
        end
    end
    if min(ap_set)==1
        done = 1;
    end
end
ap_all = zeros(1,max(vars));
for v = 1:max(vars)
    if var_index(v)~=0
        ap_all(v) = ap(var_index(v));
    end
end
if min(ap_set)==1
    for s = 1:num_sentences
        Sp(s) = BR_prob_or_atom_probs(KB(s).clauses,ap_all);
    end
    if min(ap)<0|max(ap)>1
        OK = 0;
    end
end

if min(ap_set)==0
    [ap,Sp] = BR_new_MC_partial(KB,ap_all,ap_set,num_trials,thresh,step);
end
