function [SimNLL, BestParameters] = fitDDMwSim_1gradualValShiftFull(subject)

%pathtofile = mfilename('fullpath');
fileName = mfilename;
pathtofile = which(fileName, '-ALL');
if length(pathtofile) > 1
    pathtofile = pathtofile{2};
else
    pathtofile = pathtofile{1};
end

path = pathtofile(1:(regexp(pathtofile,'DDMAnalysis') - 1));
model = 'M1_gradualValShiftFull';

%% load behavioral data
fid = fopen(fullfile(path,'SubjectData','AllSubjectsData.txt'));
t = textscan(fid,'%s');
t = reshape(t{1},13,length(t{1})/13)';
t(searchcell(t(:,9),'NULL'),10) = {'NaN'}; % replace non-responses with NaNs
t(searchcell(t(:,9),'NULL'),9) = {'NaN'}; % replace non-responses with NaNs
t = t(searchcell(t(:,1),subject),:); % get subject-specific data

% wrap up the data in a convenient data structure
clear Data;
Data.resp = cellfun(@str2num,t(:,9));
Data.resp = +(Data.resp > 2);
%NonResp = strcmp(Resp,'NULL');
Data.rt = cellfun(@str2num,t(:,10));
%Resp(NonResp) = {NaN};
%RT(NonResp) = {NaN};
TrialType = cellfun(@str2num,t(:,3));
Data.dim1 = cellfun(@str2num,t(:,5)) - 2.5; % taste rating
Data.dim2 = cellfun(@str2num,t(:,7)) - 2.5; % health rating

% load up Self condition weights from model assuming no shifts over time (used to set initial
% "anchoring" value)
SelfFile = fullfile(path,'DDMAnalysis','FittedResults','M1_gradualValShiftReduced', ['DDM_fittedResults_M1_gradualValShiftReduced_' subject '_Self.mat']);
selfData = load(SelfFile);
conds =  {'Self','Sim','Diff'};
% estimate anchoring-and-adjustment model for similar and dissimilar others
for c = [2, 3] 
    selectedTrials = TrialType == c;
    
    tempData.resp = Data.resp(selectedTrials);
    tempData.rt = Data.rt(selectedTrials);
    tempData.attribVals = [Data.dim1(selectedTrials) ...
                           Data.dim2(selectedTrials)];
    
    %% course grid search for parameters
    combo.preweight1 = selfData.BestParameters(3); % assume Ss start with self weight
    combo.preweight2 = selfData.BestParameters(4); % assume Ss start with self weight
    combo.weight1 = [-.02:.01:-.01, -.005, 0, .005, .01:.01:.14];
    combo.weight2 = [-.02:.01:-.01, -.005, 0, .005, .01:.01:.14];
    combo.barrier = (8:2:14)/100;
    combo.timeShift = 0:.1:.5; 
    combo.rateShift = [1000 2.56, 1.28 .64 .32];
    combo.nonDec = .2:.1:.7;
%     combo.startBias = -.3:.1:.3;

    % create all combinations of desired variables
    comboVars = fieldnames(combo);
    parameterCombos = combo.(comboVars{1})';
    for v = 2:length(comboVars)
        parameterCombos = [repmat(parameterCombos, length(combo.(comboVars{v})), 1), ...
            sort(repmat(combo.(comboVars{v})',size(parameterCombos,1), 1))];
    end
    
    SimNLL = zeros(length(parameterCombos),1);
    nLoops = numel(SimNLL);
   
    ddmFile = fullfile(path,'DDMAnalysis','FittedResults',model, ['DDM_fittedResults_m' model '_' subject '_' conds{c} '.mat']);
    if exist(ddmFile, 'file')
        % delete(ddmFile) % for use if you want to quickly delete old results
        load(ddmFile)
    end
    
    if any(SimNLL(:) == 0) % run simulations if not complete yet
        
        for loop = 1:nLoops
            % check to see if a kill command has been
            % issued (useful for gracefully exiting when using a cluster)
            if exist(fullfile(path,['kill' model '.mat']),'file')
                break;
            end
            
            if SimNLL(loop) == 0
                
                % get negative log likelihood of the data for each combination of
                % parameter values
                SimNLL(loop) = NLLSimData_gradualValShiftSim(tempData,parameterCombos(loop,:));
                
                if mod(loop,100) == 0 || (loop > nLoops - 101)
                    save(ddmFile,'SimNLL','loop');
                end
            end
            clc;
            fprintf('Loop %d completed: %.1f percent finished\n\n',loop,loop/nLoops * 100)
        end
    end
    
    % find the best-fitting parameters (i.e., lowest negative log-likelihood)
    i = find(SimNLL == min(SimNLL));
    if length(i) > 1
        i = i(randperm(length(i))); % pick a random set of parameters if more than one fits
        i = i(1);
    end
    
    BestParameters = parameterCombos(i, :);
    
    % simulate set of choices on full dataset using best parameters
    SimData.resp = [];
    SimData.rt = [];
    SimData.drift = [];
    SimData.drift2 = [];
    for t = 1:length(tempData.attribVals(:,1))
        params.initDrift = tempData.attribVals(t,1:2)*BestParameters(1:2)';
        params.finalDrift = tempData.attribVals(t,1:2)*BestParameters(3:4)';
        params.barrier = BestParameters(5);
        params.timeShift = BestParameters(6);
        params.rateShift = BestParameters(7);
        params.nonDec = BestParameters(8);
        params.nSim = 1000;
        tmp = simulDDM_valShiftNormApprox(params);
        
        SimData.resp = [SimData.resp; tmp.simData.resp];
        SimData.rt = [SimData.rt; tmp.simData.rt];
        SimData.drift = [SimData.drift; repmat(tempData.attribVals(t,1),1000,1)];
        SimData.drift2 = [SimData.drift2; repmat(tempData.attribVals(t,2),1000,1)];
    end
    
    
    
    % save all parameters used in simulations
    params.candidate_driftTaste = combo.weight1;
    params.candidate_driftHealth = combo.weight2;
    params.candidate_nondec = combo.nonDec;
    params.candidate_barrier = combo.barrier;
    params.candidate_timeshift = combo.timeShift;
    params.candidate_speedshift = combo.rateShift;
    params.simCode = fileread(fullfile(path,'DDMAnalysis','simulDDM_valShiftNormApprox.m')); % save the DDM used to simulate the data, just in case

    save(ddmFile,'SimNLL','BestParameters','params','SimData');

end



