wolffd@0: function [strategy, MEU, niter] = solve_limid(engine, varargin)
wolffd@0: % SOLVE_LIMID Find the (locally) optimal strategy for a LIMID
wolffd@0: % [strategy, MEU, niter] = solve_limid(inf_engine, ...)
wolffd@0: %
wolffd@0: % strategy{d} = stochastic policy for node d (a decision node)
wolffd@0: % MEU = maximum expected utility
wolffd@0: % niter = num iterations used
wolffd@0: %
wolffd@0: % The following optional arguments can be specified in the form of name/value pairs:
wolffd@0: % [default in brackets]
wolffd@0: %
wolffd@0: % max_iter - max. num. iterations [ 1 ]
wolffd@0: % tol - tolerance required of consecutive MEU values, used to assess convergence [1e-3]
wolffd@0: % order - order in which decision nodes are optimized [ reverse numerical order ]
wolffd@0: %
wolffd@0: % e.g., solve_limid(engine, 'tol', 1e-2, 'max_iter', 10)
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: 
wolffd@0: % default values
wolffd@0: max_iter = 1;
wolffd@0: tol = 1e-3;
wolffd@0: D = bnet.decision_nodes;
wolffd@0: order = D(end:-1:1);
wolffd@0: 
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'max_iter', max_iter  = args{i+1}; 
wolffd@0:    case 'tol',      tol = args{i+1}; 
wolffd@0:    case 'order',    order = args{i+1}; 
wolffd@0:    otherwise,  
wolffd@0:     error(['invalid argument name ' args{i}]);       
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: CPDs = bnet.CPD;
wolffd@0: ns = bnet.node_sizes;
wolffd@0: N = length(ns);
wolffd@0: evidence = cell(1,N);
wolffd@0: strategy = cell(1, N);
wolffd@0: 
wolffd@0: iter = 1;
wolffd@0: converged = 0;
wolffd@0: oldMEU = 0;
wolffd@0: while ~converged & (iter <= max_iter)
wolffd@0:   for d=order(:)'
wolffd@0:     engine = enter_evidence(engine, evidence, 'exclude', d);
wolffd@0:     [m, pot] = marginal_family(engine, d);
wolffd@0:     %pot = marginal_family_pot(engine, d);
wolffd@0:     [policy, score] = upot_to_opt_policy(pot);    
wolffd@0:     e = bnet.equiv_class(d);
wolffd@0:     CPDs{e} = set_fields(CPDs{e}, 'policy', policy);
wolffd@0:     engine = update_engine(engine, CPDs);
wolffd@0:     strategy{d} = policy;
wolffd@0:   end  
wolffd@0:   engine = enter_evidence(engine, evidence);
wolffd@0:   [m, pot] = marginal_nodes(engine, []);
wolffd@0:   %pot = marginal_family_pot(engine, []);
wolffd@0:   [dummy, MEU] = upot_to_opt_policy(pot);    
wolffd@0:   if approxeq(MEU, oldMEU, tol)
wolffd@0:     converged = 1;
wolffd@0:   end
wolffd@0:   oldMEU = MEU;
wolffd@0:   iter = iter + 1;
wolffd@0: end
wolffd@0: niter = iter - 1;