% Copyright (C) Daphne Koller, Stanford University, 2012

function [MEU OptimalDecisionRule] = OptimizeLinearExpectations( I )
  % Inputs: An influence diagram I with a single decision node and one or more utility nodes.
  %         I.RandomFactors = list of factors for each random variable.  These are CPDs, with
  %              the child variable = D.var(1)
  %         I.DecisionFactors = factor for the decision node.
  %         I.UtilityFactors = list of factors representing conditional utilities.
  % Return value: the maximum expected utility of I and an optimal decision rule
  % (represented again as a factor) that yields that expected utility.
  % You may assume that there is a unique optimal decision.
  %
  % This is similar to OptimizeMEU except that we will have to account for
  % multiple utility factors.  We will do this by calculating the expected
  % utility factors and combining them, then optimizing with respect to that
  % combined expected utility factor.
  #MEU = [];
  #OptimalDecisionRule = [];
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % YOUR CODE HERE
  %
  % A decision rule for D assigns, for each joint assignment to D's parents,
  % probability 1 to the best option from the EUF for that joint assignment
  % to D's parents, and 0 otherwise.  Note that when D has no parents, it is
  % a degenerate case we can handle separately for convenience.
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  D = I.DecisionFactors(1);
  EUFcomb = struct('var', [], 'card', [], 'val', []);
  I1.RandomFactors = I.RandomFactors;
  I1.DecisionFactors = I.DecisionFactors;
  for i = 1:numel(I.UtilityFactors)
    I1.UtilityFactors = I.UtilityFactors(i);
    EUFcomb = FactorSum(EUFcomb, CalculateExpectedUtilityFactor(I1));
  endfor

  OptimalDecisionRule = struct('var', EUFcomb.var, 'card', EUFcomb.card, 'val', []);

  [randVar, indxRandVar] = setdiff(EUFcomb.var, D.var(1), "stable");
  [Dvar, indxDvar] = intersect(EUFcomb.var, D.var(1), "stable");
  decision_number = prod(EUFcomb.card(indxRandVar));

  if isempty(indxRandVar)
    [MEU, idxMEU] = max(EUFcomb.val);
    lista = zeros(1, EUFcomb.card);
    lista(idxMEU) = 1;
    OptimalDecisionRule.val = lista;
    return;
  endif

  max_list = -Inf * ones(1, decision_number);
  for i = 1:prod(EUFcomb.card)
    as = AssignmentToIndex( IndexToAssignment(i, EUFcomb.card)(indxRandVar), EUFcomb.card(indxRandVar) );
    max_list(as) = max(max_list(as), EUFcomb.val(i));
  endfor

  MEU = sum(max_list);

  for i = 1:prod(EUFcomb.card)
    as = AssignmentToIndex( IndexToAssignment(i, EUFcomb.card)(indxRandVar), EUFcomb.card(indxRandVar) );
    if (EUFcomb.val(i) == max_list(as))
      OptimalDecisionRule.val(i) = 1;
    else
      OptimalDecisionRule.val(i) = 0;
    endif
  endfor

end
