function New_Factor_Levels=Factor_Levels_Updater(Current_Complete_Data, CURRENT_PARS, PRIORS, KK)
% function New_Factor_Levels=Factor_Levels_Updater(Current_Complete_Data, CURRENT_PARS, PRIORS, KK)
% UPDATES the factor levels. 

% extract the error variance:
Err_Var=CURRENT_PARS.sigma2;
% extract the reduced-space factor variance:
Fact_Var=eval(['CURRENT_PARS.Red_space_var_', num2str(KK)]);
%extract the transform matrix for this factor:
Q_M=eval(['PRIORS.Q_M_', num2str(KK)]);

% we want to arrange the dimensions of Current_Complete_Data so that the
% dimension corresponding to the factor we are currently updating is FIRST:

N_factors=ndims(Current_Complete_Data);
N_levels=size(Current_Complete_Data);
% rotate it
CCD_rotated=permute(Current_Complete_Data, [KK, setdiff([1:1:N_factors], KK)]);

% to get the Psi parameter, we need the PRODUCT of the number of levels in each
% factor, EXCEPT the factor we are currently updating:

Prod_N_levels_NotKK=prod(N_levels(setdiff([1:1:N_factors], KK)));

Psi_Par=(Prod_N_levels_NotKK/Err_Var + 1/Fact_Var)^(-1);


% to get the V parameter, need to subtract the current estimate of the
% grand mean from CCD_Rotated, and then add along all but the first
% dimension.

% a loop to add over the last dimension of the array, N_factors-1 times.
% because the dimension corrsponding to the current factor is first, the
% result is a summation over all but the one currently being updated

% set the array to be the permuted data array, less the grand mean:
Sum_over_NotKK=CCD_rotated-CURRENT_PARS.GM;
%each iteration of the loop will kill on dimensions:
for JJ=1:1:N_factors-1
    Sum_over_NotKK = sum(Sum_over_NotKK, ndims(Sum_over_NotKK));
end
   
% can now get the V parameter:
V_Par = (1/Err_Var) * Q_M' * Sum_over_NotKK;

%get the posterior mean and variance:
Post_Var = Psi_Par;
Post_Mean = Psi_Par * V_Par;

%make the draw:
New_Factor_Levels_reduced = Post_Mean + sqrt(Post_Var) * randn(N_levels(KK)-1, 1);
New_Factor_Levels = Q_M * New_Factor_Levels_reduced;
