calib2_runMultHistSims.m

% Parameter estimation

% Calculates the negative log likelihood of a sample of potential parameter sets
% Accepts:
% 1) paramSetIdx (index in larger sample matrix to start selection of sub-set)
% Saves:
% 1) File: negSumLogL_calib_[date].dat (negative log likelihood for each parameter set in sub-set)

function calib2_runMultHistSims(paramSetIdx , tstep_abc , date_abc , username)

t_curr = tstep_abc;
date = date_abc;

%% Cluster information -- Running on UW cluster
% pc = parcluster('local');    % create a local cluster object
% pc.JobStorageLocation = strcat('/gscratch/csde/' , username , '/' , getenv('SLURM_JOB_ID'))    % explicitly set the JobStorageLocation to the temp directory that was created in the sbatch script
% numCPUperNode = str2num(getenv('SLURM_CPUS_ON_NODE'))
% parpool(pc , numCPUperNode)    % start the pool with max number workers

%% Cluster information -- Erisone
pc = parcluster('local'); 
pc.JobStorageLocation = getenv('TMPDIR'); % how to pull job id?
numWorkers = 9; 
% numCPUperNode = str2double(getenv('LSB_DJOB_NUMPROC')); 
% numCPUperNode = 28; % how to pull CPUs on node? set to 8 as an initial test.
% parpool(pc, numCPUperNode)
parpool(pc, numWorkers)

%%
nPrlSets = 25;

%% Load all particles
paramDir = [pwd , '/Params/'];
% masterSetMatrix = load(string(strjoin([paramDir , 'masterSets_calib_' , date , '_' , num2str(t_curr) , '.dat'],''))); % load most recent parameter sample
masterSetMatrix = load([paramDir , 'masterSets_calib_' , date , '_' , num2str(t_curr) , '.dat']); 
% pIdx = load(string(strjoin([paramDir , 'pIdx_calib_' , date , '_0' , '.dat'],''))); % load parameter indices
pIdx = load([paramDir , 'pIdx_calib_' , date , '_0' , '.dat']);
% orderedLL = load(string(strjoin([paramDir , 'orderedLL_calib_' , date , '_' , num2str(t_curr) , '.dat'],''))); % load most recent ordered log-likelihoods
orderedLL = load([paramDir , 'orderedLL_calib_' , date , '_' , num2str(t_curr) , '.dat']);

%% Get indices and parameter values of numBestFits*2  sets
numBestFits = 25;
numSets50 = size(orderedLL , 1)*((numBestFits*2)/5600);
specIndsList = [1,2,3,6,8,9,11,12,13,15,20,21,22,26,27,32,34,35,38,39,40,41,42,45,47]; 
top50Inds = orderedLL(specIndsList,1); %orderedLL(1:numSets50,1);
top50Params = masterSetMatrix(:,top50Inds);

%% If on Phase 2 of calibration, uncomment the following to plot resampled Ph1 parameters + Ph2 parameters
% pIdx = load(string(strjoin([paramDir,'pIdx_calib_' , date , '_0_wPh1Resample.dat'],'')));
pIdx = load([paramDir,'pIdx_calib_' , date , '_0_wPh1Resample.dat']);
% masterResampleSubsetMatrix = load(string(strjoin([paramDir , 'masterResampleSubsetMatrix_calib_' , date , '_' , num2str(t_curr) , '.dat'],''))); % load most recent Ph1 resampled parameters
masterResampleSubsetMatrix = load([paramDir , 'masterResampleSubsetMatrix_calib_' , date , '_' , num2str(t_curr) , '.dat']);
masterCombinedPhaseMatrix = [masterResampleSubsetMatrix ; masterSetMatrix];
top50Params = masterCombinedPhaseMatrix(:,top50Inds);

%% Set up paramsSub for indexing into paramSet matrix
[paramsAll] = genParamStruct();
paramsSub = cell(length(pIdx),1);
startIdx = 1;
for s = 1 : length(pIdx)
    paramsSub{s}.length = paramsAll{pIdx(s)}.length;
    paramsSub{s}.inds = (startIdx : (startIdx + paramsSub{s}.length - 1));
    startIdx = startIdx + paramsSub{s}.length;
end

%% Obtain model output for each set of sampled parameters
for m = paramSetIdx : nPrlSets : (numBestFits+paramSetIdx-1)
    subMatrixInds = [m : (m + nPrlSets - 1)];
    parfor n = 1 : nPrlSets
        paramSet = top50Params(:,subMatrixInds(n));
        historicalSim(1 , pIdx , paramsSub , paramSet ,specIndsList(m + n - 1) , tstep_abc , date_abc , n);
    end
end