% Study 1
clc, clear
clear global
rng(0,'twister');
global soft_out
options = optimset;% optimization specs
options = optimset(options,'Display', 'off','PlotFcns', {@optimplotx @optimplotfval}, ...
   'Algorithm', 'active-set', 'FinDiffType', 'central');
options.PlotFcns = [];
D=dir('*.csv');% array of subject files to loop through
subjStay=[];
subjProbe=[];
AIC=[];
for s=1:length({D.name})
    %% Data organizing
    format longg
    opts = detectImportOptions(D(s).name);
    opts = setvartype(opts, 'condition', 'char');
    % Load subject-specific data
    dat = readtable(D(s).name, opts);
    % Color-to-condition assignment

    % Shape-to-condition assignment
    dat.Properties.VariableNames(15)={'shapes'};
    shapesCB = erase(split(dat.shapes(1))', ["[", "]", "'", ","]);
    choice=dat.choice(~isnan(dat.choice));   
    trialType=dat.ProbeType(~isnan(dat.choice)); 
    rewLev=dat.rewardLevel(~isnan(dat.choice));
    cons=dat.payoffCons(~isnan(dat.choice));% consensus level on each trial   
    chosUnchosOPT=zeros(height(choice),3);% array to plug in chosen & unchosen trial option + outcome condition
    conds={'1', '2','3', '4', '1_partial', '2_partial'};% consensus & win/lose conditions   
    for cb=1:6
        chosUnchosOPT(strcmp(dat.imageLeft, shapesCB{cb}),1)=cb;% add IDs of available shape options
        chosUnchosOPT(strcmp(dat.imageRight, shapesCB{cb}),2)=cb;
        chosUnchosOPT(strcmp(dat.condition, conds(cb)),3)=cb+6;
    end
    colprob={'42', '24', '13', '31'};% probe color options
    for cb=1:4
        chosUnchosOPT(strcmp(dat.condition, colprob(cb)),1)=str2double(colprob{cb}(1))+6;
        chosUnchosOPT(strcmp(dat.condition, colprob(cb)),2)=str2double(colprob{cb}(2))+6;
    end
    chosUnchosREW=[dat.rewardLeft(~isnan(dat.choice)), dat.rewardRight(~isnan(dat.choice))];% payoff on trial for each option
    % Convert to Chosen & Unchosen options and rewards
    for i = 1:height(choice)
        if choice(i)==39
            chosUnchosOPT(i,1:2)=[chosUnchosOPT(i,2),chosUnchosOPT(i,1)]; 
            chosUnchosREW(i,1:2)=[chosUnchosREW(i,2),chosUnchosREW(i,1)]; 
        end
        if chosUnchosREW(i,1)<chosUnchosREW(i,2)
            chosUnchosREW(i,1)=0;% set lower payoff to $0
        else
            chosUnchosREW(i,2)=0;
        end
    end
%% Behavioral Analyses
    % Stay probabilities
    numRewProbe = zeros(1,12);% mean reward for each shape/color
    rewStay = zeros(2,6);% cumulative reward for Stay/stay trials
    bigStay = zeros(2,6);% stay counter
    numTrials = zeros(2,6);% trial counter
    consWinLev=zeros(height(trialType), 3);% chos and unchos options on potential Stay trial; col 3 is 1 (stayed) or 0 (switched)
    for i = 1:height(trialType)-1 
       if trialType(i)==0% if feedback trial
          if cons(i)==2% condition order is CW,CL,DW,DL,PW,PL
            cond=1;
          elseif cons(i)==1
            cond=5;
          elseif cons(i)==0
            cond=3;
          end
          if ~isnan(cond) && chosUnchosREW(i,1)<chosUnchosREW(i,2)
                cond=cond+1;
          end
          if chosUnchosOPT(i,1)==chosUnchosOPT(i+1,1) || chosUnchosOPT(i,1)==chosUnchosOPT(i+1,2)% if chosen option is repeated on the next trial
             numTrials(rewLev(i)+1, cond)=numTrials(rewLev(i)+1, cond)+1;% counter      
             consWinLev(i+1,rewLev(i)+1)=cond;% index conditions on previous trial
             rewStay(rewLev(i)+1, cond)=rewStay(rewLev(i)+1, cond)+chosUnchosREW(i,1); 
             if chosUnchosOPT(i+1,1)==chosUnchosOPT(i,1) %if stayed
                bigStay(rewLev(i)+1, cond)=bigStay(rewLev(i)+1, cond)+1;
                consWinLev(i+1,3)=1;
             end  
          end 
       end
    end
    preRating = dat.Var17(dat.Var17>-1)'; postRating = dat.Var19(dat.Var19>-1)'; 
    behavStay=bigStay./numTrials;
    for c=1:8 
        choice_probe(c)=length(find(chosUnchosOPT(:,1)==c & trialType(:,1)~=0))/length(find((chosUnchosOPT(:,1)==c | chosUnchosOPT(:,2)==c) & trialType(:,1)~=0));
        rewProbe(c)=nanmean(chosUnchosREW(chosUnchosOPT(:, 1)==c,1));
    end 
    subjSpecs{s}={trialType, chosUnchosOPT, chosUnchosREW, preRating, cons, consWinLev};% Store subject-specific variable arrays
    % Cell array organizing data for plotting means, participants are rows, measures are columns 
    subjStay(s,:)=nanmean(behavStay);% Stay probabilities in each condition, order is CW, CL, DW, DL, PW, PL, for high and low reward distributions 
    subjProbe(s,:)=[nanmean(choice_probe(:, 1:2)), nanmean(choice_probe(:, 4:5)), choice_probe(:,3), ...
        choice_probe(:,6:8)];% Probe choice props for shapes by rew and freq levels, (HH, LH, HL, LL), followed by consensus pref
    totRew(s,:)=[rewStay(1,:), rewStay(2,:), [mean(rewProbe(1:2)),mean(rewProbe(:,4:5)), mean(rewProbe(:,3)), mean(rewProbe(:,6))]];% Rewards received in each Stay & probe condition          
    totTrials(s,:)=[numTrials(1,:), numTrials(2,:)];% Trials in each Stay condition      
    subjRat{s,5}=postRating-preRating;
%% Fit & Compare Models   
    k=[4,3,3,2];% # of parameters for AIC
    for mods=1:4
        paramRec=0;
        numPar=repmat([0.5; 0.01; 0.99], 1, k(mods)-2);
        [params, nll, ~, ~]=fminsearchbnd(@(fp) fitparams(fp, subjSpecs{s}, mods, paramRec),[numPar(1,:),0.5,5],[numPar(2,:),0.01,1],[numPar(3,:), 0.99,10],options);
        subjParams{s,mods}=params; 
        AIC(s,mods)=2*k(mods)+2*(nll);      
    end  
end
save behavSpecs subj* AIC
%% Parameter Recovery 
paramRec=input('Enter 1 for Parameter Recovery: ');
iter=1000;
if paramRec == 1
    % Randomly generated parameter values
    qlr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for money
    ulr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for action imitation
    vlr=(0.99-0.01).*rand(iter,1) + 0.01;% learning rate for social reward
    dn=(10-1).*rand(iter,1) + 1; % Decision noice 
    theta=[qlr, ulr, vlr, dn]; 
    modParam=[];
    m=[1 5 9 13];
    for pr=1:iter
        %  generate experiment events using randomly drawn participant
        r=randi([1,length({D.name})]);
        specs=subjSpecs{r};        
        rep=[1 2 3];
        for x=1:length(rep)% increase training trials
            specs{rep(x)}=repmat(specs{rep(x)}, 10,1);
        end
        specs{5}=randi(2, 1500,1)+randi(2, 1500,1)-2;% generate choices by 2 other agents              
        for mods=1:length(m)
            % simulate weights & decisions given parameters & model
            fitparams(theta(pr,:),specs, mods, paramRec);  
            simChoice{pr,mods}=soft_out; 
            % recover parameters 
            paramRec=-1;
            [params, nll, ~, ~]=fminsearchbnd(@(fp) fitparams(fp, simChoice{pr,mods}, mods, paramRec),[0.3,0.3,0.3,5],[0.01,0.01,0.01,1],[0.99,0.99,0.99,10],options);        
            modParam(pr, m(mods):m(mods)+3)=params;
            paramRec=1;
        end  
    end
    for mods=1:length(m)
        h=modParam(:,m(mods):m(mods)+3);
        for pr=1:width(theta)
            [r,p]=corr(h(:, pr), theta(:,pr));
            paramCorrR(mods,pr)=r;  
            paramCorrP(mods,pr)=p; 
        end
    end
   paramRecovCorr=array2table(paramCorrP, 'RowNames', {'hybrid', 'imitation', 'social', 'baseline'},'VariableNames', {'$lr', 'imitationlr', 'sociallr', 'decision noise'});
   save paramRecov simChoice paramCorr* theta modParam paramRecovCorr  
    %% Model Predictions
    probeTrials=[];
    stayProbs=[];
    mPt=[1 9 17 25];
    mSp=[1 7 13 19];
    for x=1:length(mPt)
        for y=1:height(simChoice)
           % Probe trials  
           hold=simChoice{y,x};
           hold=hold{2}(hold{1}~=0, 1:2);
           pT=zeros(2,8);
           for z=1:8
            pT(1,z)=size(find(hold(:,1)==z),1); 
            pT(2,z)=size(find(hold(:,1:2)==z),1);
           end
           hold=simChoice{y,x};
           hold=hold{2}(hold{1}==0, 1:3);
           hold(:,3)=hold(:,3)-6;
           stayTrials=zeros(2,6);
           for i=1:height(hold)-1
               if hold(i+1,1) == hold(i,1) || hold(i+1,2) == hold(i,1) % if chosen option is repeated on the next trial
                  stayTrials(1, hold(i,3))=stayTrials(1, hold(i,3))+1;% counter       
                   if hold(i+1,1) == hold(i,1) %if stayed
                      stayTrials(2, hold(i,3))=stayTrials(2, hold(i,3))+1;% counter 
                   end  
               end 
           end
           probeTrials(y,mPt(x):mPt(x)+7)=pT(1,:)./pT(2,:);
           stayProbs(y,mSp(x):mSp(x)+5)=stayTrials(2,:)./stayTrials(1,:);
        end    
    end
end
%% Models
function nll=fitparams(fp, specs, mods, paramRec) 
    global soft_out
    % specs: 1. trialType, 2. chosUnchosOPT, 3. chosUnchosREW, 4. rating, 5. cons, 6. consWinLev
    freqVals = ones(1,12);
    if mods<3 % hybrid and imitation
        actVals = zeros(1,12);
        payVals = zeros(1,12)+0.5;
        sumVals=actVals+payVals;
    else
        sumVals = zeros(1,12)+0.5;
    end % social and monetary
    initVals=sumVals;
    M=zeros(height(specs{1}));
    modelprob=zeros(size(M,1),1);
    specsO=specs;
    for i = 1:height(M)
        M(i,1)=sumVals(specs{2}(i,1))+sqrt(2*log(i)/freqVals(specs{2}(i,1)));
        M(i,2)=sumVals(specs{2}(i,2))+sqrt(2*log(i)/freqVals(specs{2}(i,2)));
        modelprob(i)=exp(fp(length(fp))*M(i,1))/(exp(fp(length(fp))*M(i,1))+exp(fp(length(fp))*M(i,2)));% softmax with noise parameter 
        if paramRec==1    
            chosResp = randsample([1 0], 1, true, [modelprob(i) 1-modelprob(i)]);% generate response based on model probs
            if chosResp==0% chosen and unchosen response contingencies 
                specs{2}(i,1)=specsO{2}(i,2); 
                specs{3}(i,1)=specsO{3}(i,2); 
                specs{2}(i,2)=specsO{2}(i,1); 
                specs{3}(i,2)=specsO{3}(i,1);
                specs{5}(i)=2-specsO{5}(i);
            end
            if specs{5}(i)==2; specs{2}(i,3)=7; 
            elseif specs{5}(i)==0; specs{2}(i,3)=9; 
            elseif specs{5}(i)==1; specs{2}(i,3)=11; 
            end
            if specs{3}(i,1)<specs{3}(i,2); specs{2}(i,3)=specs{2}(i,3)+1; end
        end
        if specs{1}(i)==0% filter out non-feedback trials before updating 
            freqVals(specs{2}(i,1))=freqVals(specs{2}(i,1))+1;% increase n for chosen shape on feedback trials
            freqVals(specs{2}(i,3))=freqVals(specs{2}(i,3))+1;% increase n for color on feedback trials
            if mods==1% Hybrid model
               actVals(specs{2}(i,1))=actVals(specs{2}(i,1))+fp(2)*(specs{5}(i)-actVals(specs{2}(i,1)));% action imitation
               payVals(specs{2}(i,1))=payVals(specs{2}(i,1))+fp(1)*(specs{3}(i,1)-(payVals(specs{2}(i,1))))+fp(3)*(specs{5}(i)-(payVals(specs{2}(i,1))));% action value based on money & conformity
               payVals(specs{2}(i,3))=payVals(specs{2}(i,3))+fp(1)*(specs{3}(i,1)-(payVals(specs{2}(i,3))))+fp(3)*(specs{5}(i)-(payVals(specs{2}(i,3))));% state value based on money & conformity
               sumVals = actVals+payVals;
            elseif mods==2% Imitation model
               actVals(specs{2}(i,1))=actVals(specs{2}(i,1))+fp(2)*(specs{5}(i)-actVals(specs{2}(i,1)));% action imitation
               payVals(specs{2}(i,1))=payVals(specs{2}(i,1))+fp(1)*(specs{3}(i, 1)-payVals(specs{2}(i,1)));% action value based on money
               payVals(specs{2}(i,3))=payVals(specs{2}(i,3))+fp(1)*(specs{3}(i, 1)-payVals(specs{2}(i,3)));% state value based on money
               sumVals = actVals+payVals;
            elseif mods==3% Social Reward model
               sumVals(specs{2}(i,1))=sumVals(specs{2}(i,1))+fp(1)*(specs{3}(i,1)-(sumVals(specs{2}(i,1))))+fp(2)*(specs{5}(i)-(sumVals(specs{2}(i,1))));% action value based on money & conformity
               sumVals(specs{2}(i,3))=sumVals(specs{2}(i,3))+fp(1)*(specs{3}(i,1)-(sumVals(specs{2}(i,3))))+fp(2)*(specs{5}(i)-(sumVals(specs{2}(i,3))));% state value based on money & conformity
            elseif mods==4% Money Only model
               sumVals(specs{2}(i,1))=sumVals(specs{2}(i,1))+fp(1)*(specs{3}(i,1)-(sumVals(specs{2}(i,1))));% action value based on money             
               sumVals(specs{2}(i,3))=sumVals(specs{2}(i,3))+fp(1)*(specs{3}(i,1)-(sumVals(specs{2}(i,3))));% state value based on money            
            end
        end
    end
    specs{4}=sumVals;
    nll=sum(log(modelprob))*-1; 
    if paramRec == 1
        soft_out=specs;
    end      
end




