% RunMCMCSimForBPHMM
%  Generic harness for running many iterations of MCMC,
%   allows sensible reporting/saving of samples and diagnostics
%USAGE
%  Usually called from more "user-friendly" function "runBPHMM"
%  but if specific data and initial configuration Psi are available,
%  >> RunMCMCForBPHMM( data, Psi, algP, outP )
%INPUT
%  data : SeqData object, defining a collection of sequences to fit model
%  Psi  : initial model state
%           usually generated by a function in the "init" folder
%  algParams : specifies MCMC behavior (# iterations, proposal distribs)
%           see defaults/defaultMCMCParams_BPHMM.m for details
%  outParams : specifies MCMC output behavior
%           how often to save samples, write to disk, etc.
%           see defaults/defaultOutputParams_BPHMM.m for details
%OUTPUT
%  Markov Chain state variables saved at preset frequency to hard drive
%     at filepath location specified in outParams.saveDir

function [ChainHist] = RunTimedMCMCSimForBPHMM( data, Psi, algParams, outParams, model )
tic;

if isfield( Psi, 'F' )
    % Stating chain from scratch
    n = 0;
    logPr = calcJointLogPr_BPHMMState( Psi, data );
    ChainHist = recordMCMCHistory_BPHMM( 0, outParams, [], Psi, logPr  );

    fprintf( 'Initial Config: \n' );
    printTimedMCMCSummary_BPHMM( 0, Psi, logPr, algParams); 
else
    ChainHist = Psi;
    Psi = unpackBPHMMState(  ChainHist.Psi(end), data, model );
    logPr = calcJointLogPr_BPHMMState( Psi, data );
    n = ChainHist.iters.Psi(end );
    fprintf( 'Resumed Config: \n' );
    printTimedMCMCSummary_BPHMM( 0, Psi, logPr, algParams); 
end

fprintf( 'Running MCMC Sampler %d : %d ... \n', outParams.jobID, outParams.taskID );

while toc < algParams.TimeLimit
    n = n + 1;
    Psi.iter = n;

    % Perform 1 iteration of MCMC, moving to next Markov state!
    [Psi, Stats] = BPHMMsample( Psi, data, algParams );
    
    % Diagnose convergence by calculating joint log pr. of all sampled vars
    if n == 1 || rem(n, outParams.logPrEvery)==0
        % NB: not passing "data" as arg means Psi stores all X suff stats
        logPr = calcJointLogPr_BPHMMState( Psi );
    end
    
    %Record current sampler state
    %  NB: internally only records at preset frequency
    ChainHist = recordMCMCHistory_BPHMM( n, outParams, ChainHist, Psi, logPr, Stats );
    
    doSaveToDisk = n==1 || rem(n, outParams.saveEvery)==0 || toc > algParams.TimeLimit;
    if doSaveToDisk
        filename = fullfile( outParams.saveDir,  'SamplerOutput.mat' );
        save(filename, '-struct', 'ChainHist');
    end
    
    if n == 1 || rem(n, outParams.printEvery)==0
       printTimedMCMCSummary_BPHMM( n, Psi, logPr, algParams); 
    end
    
end % loop over sampler iterations

fprintf( '<<<<< --------------------------------------------------- \n');

end % main function

