N = 4; 
dag = zeros(N,N);
C = 1; S = 2; R = 3; W = 4;
dag(C,[R S]) = 1;
dag(R,W) = 1;
dag(S,W)=1;

discrete_nodes = 1:N;
node_sizes = 2*ones(1,N); 

bnet = mk_bnet(dag, node_sizes, 'discrete', discrete_nodes);
bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);


clear engine;

engine{1} = jtree_inf_engine(bnet);
engine{2} = var_elim_inf_engine(bnet);
engine{3} = belprop_inf_engine(bnet);

maxIter = 100;
numVars = 5;

evidence = cell(1,numVars);
E = length(engine);
ll = zeros(maxIter,3);
totalTimes = zeros(1,E);
times = zeros(maxIter,3);

marginalTimes = zeros(maxIter, 3);
totalMargTimes = zeros(1,E);

for i=1:maxIter
    for j=1:numVars
        evidence{j} = floor(rand*2)+1;
    end
    

    tic; [engine{1}, ll(i,1)] = enter_evidence(engine{1}, evidence); 
    times(i,1) = times(i,1) + toc; totalTimes(1) = totalTimes(1) + times(i,1);
    
    tic; [engine{2}, ll(i,2)] = enter_evidence(engine{2}, evidence); 
    times(i,2) = times(i,2) + toc; totalTimes(2) = totalTimes(2) + times(i,2);
    
    tic; [engine{3}, ll(i,3), niter(3)] = enter_evidence(engine{3}, evidence);
    times(i,3) = times(i,3) + toc; totalTimes(3) = totalTimes(3) + times(i,3);

    for j=1:4
        tic; marginal_nodes(engine{1}, j); 
        marginalTimes(i,1) = marginalTimes(i,1) + toc; totalMargTimes(1) = totalMargTimes(1) + marginalTimes(j,1);

        tic; marginal_nodes(engine{2}, j); 
        marginalTimes(i,2) = marginalTimes(i,2) + toc; totalMargTimes(2) = totalMargTimes(2) + marginalTimes(j,2);

        tic; marginal_nodes(engine{3}, j);
        marginalTimes(i,3) = marginalTimes(i,3) + toc; totalMargTimes(3) = totalMargTimes(3) + marginalTimes(j,3);
    end
    
end

