%% simulating based on the same trial structure:
clearvars
OVERWRITE = 1;                                                              % redo the simulations figure 7 
currentFolder = pwd;
addpath(genpath('Matlab_functions')) 
path = currentFolder(1:end-8); 
data_folder = 'data\processed\';
fits_folder = 'data\fits\';
%%
NS = [1:10,12:24];
CW = [-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,1,1,1,-1,-1,-1,1,1,-1];      % -1 no/CCW-CW; 1 yes/CW-CCW
%% Motor primitives
P.a = 0.97;
P.sigma = 30;
P.theta =360;
P.nj = 100;
[~,g] = MotorPrimitives(P);
%% T
%% trialtypes 
%   Forces
TT.ff = [0 1 2] ;                                                           
% 0: no force trials,                
% 1: first force field trials,          
% 2: second force field trials
%   Clamp Trials
TT.chan =[3 4 5 6 7];                                                       
% 5: clamp trials in no force field, 
% 3: clamp trials in first force field,
% 6: clamp trials in second force field,
% 4: spontaneous recovery, 
% 7: post break clamp trials
TT.fit = [3 4 6];
% trials used during fitting
%% amount of trials for simulation:
%   fam  BL  PBT  FF1   FF2     SR
N = 132 +44  +11 +250  +15+     44;
%% load in fitting results
 
for s = NS
    if  exist(sprintf('%s%ssim_same_length%02d.MAT',path,'data\simulations\',s),'file')~=2 || OVERWRITE
        clear D X fval ind_end ind_remove
        %% plan based
        file = sprintf('%s%sfits_%s_%02d.MAT',path,fits_folder,'plan_based',s);
        load(file,'T')
        file = sprintf('%s%s%s%02d.MAT',path,'data\BIC\','BIC',s);
        load(file,'X_MB','X_PB')
        %% plan based model prediction with best parameters:
        %% changing T to equal length per participant:
        % only select familiarization and first block
        ind_end = find(T.block == 2,1)-11; % - post break testing
        
        % making the second force field equally long (15 trials)
        ind_FF2 = find((T.trialtype == 2 | T.trialtype ==6)& T.block <2);
        
        if ind_FF2 > 15
            ind_remove = [ind_FF2(16:end)];
        end
        %% changing T to new trial length
        fields = fieldnames(T);
        for i = 1:length(fields)
            if ~strcmp(fields(i),'nTrials')
                dummy = getfield(T,string(fields(i)));
                dummy(ind_end:end) = [];
                if exist('ind_remove','var')
                    dummy(ind_remove) = [];
                end
                T = setfield(T, string(fields(i)),dummy);
            else
                T = setfield(T, string(fields(i)),N);
            end
        end
        T.Fext = CW(s)*T.Fext;
        
        [~,learningPB,succes(s,1),~,modelpredPB]= GeneralizationModel(T, X_PB,g,TT,0);
        %% motion based
        clear D X fval
        file = sprintf('%s%sfits_%s_%02d.MAT',path,fits_folder,'motion_based',s);
        load(file,'T')
        %% changing T to new trial length
        fields = fieldnames(T);
        for i = 1:length(fields)
            if ~strcmp(fields(i),'nTrials')
                dummy = getfield(T,string(fields(i)));
                dummy(ind_end:end) = [];
                if exist('ind_remove','var')
                    dummy(ind_remove) = [];
                end
                T = setfield(T, string(fields(i)),dummy);
            else
                T = setfield(T, string(fields(i)),N);
            end
        end
        T.Fext = CW(s)*T.Fext;
        %% 
        ang = T.theta_planned;
        ang(ang>180) = ang(ang>180)-360;
        T.UpdateAngle = CW(s)*(T.UpdateAngle-ang'); % mirrored update angle relative to the target
        %% motion based model prediction with best parameters:
        [~,learningMB,succes(s,2),~,modelpredMB]= GeneralizationModel(T, X_MB,g,TT,0);
        %%
        file = sprintf('%s%ssim_same_length%02d.MAT',path,'data\simulations\',s);
        save(file,'learningMB','learningPB','modelpredMB','modelpredPB','T')
    end
end