function G = calculate_G_policies(A, B, C, qs_current, policies, epsilon)
% def calculate_G_policies(A, B, C, qs_current, policies):
% 
%   G = np.zeros(len(policies)) # initialize the vector of expected free energies, one per policy
G = zeros(size(policies,2),1);
%   H_A = entropy(A)            # can calculate the entropy of the A matrix beforehand, since it'll be the same for all policies
H_A = -sum(A.*log(A+epsilon),1);
% 
%   for policy_id, policy in enumerate(policies): # loop over policies - policy_id will be the linear index of the policy (0, 1, 2, ...) and `policy` will be a column vector where `policy[t,0]` indexes the action entailed by that policy at time `t`

for policy_id = 1:size(policies,2)% 
%     t_horizon = policy.shape[0] # temporal depth of the policy
t_horizon = size(policies,1);
% 
%     G_pi = 0.0 # initialize expected free energy for this policy
G_pi = 0;
% 
%     for t in range(t_horizon): # loop over temporal depth of the policy
for t = 1:t_horizon
% 
%       action = policy[t,0] # action entailed by this particular policy, at time `t`
action_i = policies(t,policy_id);
% 
%       # get the past predictive posterior - which is either your current posterior at the current time (not the policy time) or the predictive posterior entailed by this policy, one timstep ago (in policy time)
%       if t == 0:
%         qs_prev = qs_current 
%       else:
%         qs_prev = qs_pi_t

if t==1
    qs_prev = qs_current;
else
    qs_prev = qs_pi_t;
end




% 
%       qs_pi_t = get_expected_states(B, qs_prev, action) # expected states, under the action entailed by the policy at this particular time
qs_pi_t = get_expected_states(B, qs_prev, action_i);

%       qo_pi_t = get_expected_observations(A, qs_pi_t)   # expected observations, under the action entailed by the policy at this particular time
qo_pi_t = get_expected_observations(A, qs_pi_t);
% 
%       kld = kl_divergence(qo_pi_t, C) # Kullback-Leibler divergence between expected observations and the prior preferences C
kld = kl_divergence(qo_pi_t, C, epsilon);
% 
%       G_pi_t = H_A.dot(qs_pi_t) + kld # predicted uncertainty + predicted divergence, for this policy & timepoint
    G_pi_t = H_A*qs_pi_t + kld;
    

% 
%       G_pi += G_pi_t # accumulate the expected free energy for each timepoint into the overall EFE for the policy
G_pi = G_pi + G_pi_t;

end% 
%     G[policy_id] += G_pi
G(policy_id) = G_pi;
end

% 
%   return G

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