function G = calculate_G(A, B, C, qs_current, actions, epsilon)

  G = zeros(length(actions),1); % vector of expected free energies, one per action

  H_A = -sum(A.*log(A+epsilon),1); % entropy(A) entropy of the observation model, P(o|s)

  for action_i = 1:length(actions)
    
    qs_u = get_expected_states(B, qs_current, action_i); % expected states, under the action we're currently looping over
    qo_u = get_expected_observations(A, qs_u);           % expected observations, under the action we're currently looping over

    pred_uncertainty = H_A*qs_u; % predicted uncertainty, i.e. expected entropy of the A matrix
    pred_div = kl_divergence(qo_u, C, epsilon); % predicted divergence

    G(action_i) = pred_uncertainty + pred_div; % sum them together to get expected free energy
  end
  

function qs_u = get_expected_states(B, qs_current, action_i)
qs_u = B(:,:,action_i)*qs_current;
end

function qo_u = get_expected_observations(A, qs_u)
qo_u = A*qs_u;
end

function pred_div = kl_divergence(qo_u, C, epsilon)
pred_div = (log(qo_u+epsilon)' - log(C+epsilon))*qo_u;
end

end

