Daniel@0: % oil wildcatter influence diagram in Cowell et al p172
Daniel@0: 
Daniel@0: % T = test for oil?
Daniel@0: % UT = utility (negative cost) of testing
Daniel@0: % O = amount of oil = Dry, Wet or Soaking
Daniel@0: % R = results of test = NoStrucure, OpenStructure, ClosedStructure or NoResult
Daniel@0: % D = drill?
Daniel@0: % UD = utility of drilling
Daniel@0: 
Daniel@0: % Decision sequence = T R D O
Daniel@0: 
Daniel@0: T = 1; UT = 2; O = 3; R = 4; D = 5; UD = 6;
Daniel@0: N = 6;
Daniel@0: dag = zeros(N);
Daniel@0: dag(T, [UT R D]) = 1;
Daniel@0: dag(O, [R UD]) = 1;
Daniel@0: dag(R, D) = 1;
Daniel@0: dag(D, UD) = 1;
Daniel@0: 
Daniel@0: ns = zeros(1,N);
Daniel@0: ns(O) = 3; ns(R) = 4; ns(T) = 2; ns(D) = 2; ns(UT) = 1; ns(UD) = 1;
Daniel@0: 
Daniel@0: limid = mk_limid(dag, ns, 'chance', [O R], 'decision', [T D], 'utility', [UT UD]);
Daniel@0: 
Daniel@0: limid.CPD{O} = tabular_CPD(limid, O, [0.5 0.3 0.2]);
Daniel@0: tbl = [0.6 0 0.3 0 0.1 0  0.3 0 0.4 0 0.4 0  0.1 0 0.3 0 0.5 0  0 1 0 1 0 1];
Daniel@0: limid.CPD{R} = tabular_CPD(limid, R, tbl);
Daniel@0: 
Daniel@0: limid.CPD{UT} = tabular_utility_node(limid, UT,  [-10 0]);
Daniel@0: limid.CPD{UD} = tabular_utility_node(limid, UD, [-70 50 200  0 0 0]);
Daniel@0: 
Daniel@0: if 1
Daniel@0:   % start with uniform policies
Daniel@0:   limid.CPD{T} = tabular_decision_node(limid, T);
Daniel@0:   limid.CPD{D} = tabular_decision_node(limid, D);
Daniel@0: else
Daniel@0:   % hard code optimal policies
Daniel@0:   limid.CPD{T} = tabular_decision_node(limid, T, [1.0 0.0]);        
Daniel@0:   a = 0.5; b = 1-a; % arbitrary value
Daniel@0:   tbl = myreshape([0 a 1 a 1 a a a  1 b 0 b 0 b b b], ns([T R D]));
Daniel@0:   limid.CPD{D} = tabular_decision_node(limid, D,  tbl);
Daniel@0: end
Daniel@0: 
Daniel@0: %fname = '/home/cs/murphyk/matlab/Misc/loopybel.txt';
Daniel@0: 
Daniel@0: engines = {};
Daniel@0: engines{end+1} = global_joint_inf_engine(limid);
Daniel@0: engines{end+1} = jtree_limid_inf_engine(limid);
Daniel@0: %engines{end+1} = belprop_inf_engine(limid, 'max_iter', 3*N, 'filename', fname);
Daniel@0: 
Daniel@0: exact = [1 2];
Daniel@0: %approx = 3;
Daniel@0: approx = [];
Daniel@0: 
Daniel@0: E = length(engines);
Daniel@0: strategy = cell(1, E);
Daniel@0: MEU = zeros(1, E);
Daniel@0: for e=1:E
Daniel@0:   [strategy{e}, MEU(e)] = solve_limid(engines{e});
Daniel@0:   MEU
Daniel@0: end
Daniel@0: MEU
Daniel@0: 
Daniel@0: for e=exact(:)'
Daniel@0:   assert(approxeq(MEU(e), 22.5))
Daniel@0:   % U(T=yes)  U(T=no)
Daniel@0:   % 1         0
Daniel@0:   assert(argmax(strategy{e}{T}) == 1); % test = yes
Daniel@0:   t = 1; % test = yes
Daniel@0:   % strategy{D} T       R      U(D=yes=1)  U(D=no=2)
Daniel@0:   %             1=yes   1=noS  0           1         Don't drill
Daniel@0:   %             2=no    1=noS  1           0
Daniel@0:   %             1=yes   2=opS  1           0
Daniel@0:   %             2=no    2=opS  1           0
Daniel@0:   %             1=yes   3=clS  1           0
Daniel@0:   %             2=no    3=clS  1           0
Daniel@0:   %             1=yes   4=unk  1           0
Daniel@0:   %             2=no    4=unk  1           0
Daniel@0:   
Daniel@0:   for r=[2 3] % OpS, ClS
Daniel@0:     assert(argmax(squeeze(strategy{e}{D}(t,r,:))) == 1); % drill = yes
Daniel@0:   end
Daniel@0:   r = 1; % noS
Daniel@0:   assert(argmax(squeeze(strategy{e}{D}(t,r,:))) == 2); % drill = no
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: for e=approx(:)'
Daniel@0:   approxeq(strategy{exact(1)}{T}, strategy{e}{T})
Daniel@0:   approxeq(strategy{exact(1)}{D}, strategy{e}{D})
Daniel@0: end