Add files via upload

MCMCsampler.m

@@ -0,0 +1,131 @@
+% MCMCsampler.m - performs markov chain monte carlo sampling of parameters
+%
+%     /--------------------------------------------------------------\
+%     | Date:         02/16/2022                                     |
+%     | Author:       Jonathan Cody                                  |
+%     | Affiliation:  Pienaar Computational Systems Pharmacology Lab |
+%     |               Weldon School of Biomedical Engineering        |
+%     |               Purdue University                              |
+%     \--------------------------------------------------------------/
+%
+% TOODOO
+% -> validate with FitBOUNDS as inf
+% -> switch everything to linear scale?
+%
+% This function uses the results file from 'calibrator.m' to setup parallel
+% tempering MCMC. Chains are initialized from 'calibrator.m' parameter
+% results. Log-uniform prior distributions are assumed (i.e. known
+% parameter distributions prior to performing MCMC). These are bounded
+% according to 'setup.FitBOUNDS'.
+%
+% NOTE: Cost function must be negative log likelihood (see 'costFun_NLL').
+%
+% This uses functions from a 3rd-party toolbox (see address for citation)
+% https://github.com/ICB-DCM/PESTO
+%
+% 'performPT.m' is called to perform parallel tempering MCMC
+% https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m
+%
+% 'performPT.m' has function 'updateStatistics' added as subfunction
+% https://github.com/ICB-DCM/PESTO/blob/master/private/updateStatistics.m
+%
+% 'PTOptions.m' is referenced for 'performPT.m' default options
+% https://github.com/ICB-DCM/PESTO/blob/master/%40PTOptions/PTOptions.m
+%
+%% ========================================================================
+%	INPUTS
+%  ========================================================================
+%
+% resultsFile = file name for 'calibrator.m' results file
+%
+% InitialRanks = indices in 'Results' from which to pull initial parameter
+% vectors (each vector will initialize one chain)
+% (EX: InitialRanks = 1:5 will use 5 best parameter sets)
+%
+% numSample = MCMC sample size (aka number of iterations)
+%
+% saveEach = save period for MCMC (used to continue algorithm if stopped)
+% (EX: saveEach = 5000 will save every 5000 iterations)
+% (NOTE: saving too frequently will slow down algorithm)
+%
+%% ========================================================================
+%	OUTPUTS
+%  ========================================================================
+%
+% res.par = matrix of log values of sampled parameters
+%
+% res.logPost = evaluations of cost function ( *-1 )
+%
+% res.opt = settings for MCMC file 'performPT.m' (see file)
+%
+%% ========================================================================
+%	FUNCTION
+%  ========================================================================
+function res = MCMCsampler(resultsFile,InitialRanks,numSamples,saveEach)
+
+% load data files generated by 'calibrator.m' -----------------------------
+
+warning off
+load(resultsFile,'Results')
+setup = Results(1).setup ;
+warning on
+
+% define 'logPostHandle' input for 'performPT.m' --------------------------
+
+% logPostHandle takes in log-scale of sampled parameters, converts them to
+% linear scale, inserts them into the full parameter vector, and evaluates
+% costFunction (similar to calibrator.m). Negative cost is returned.
+logPostHandle = @(LSP) targetFun(LSP,setup)  ;
+
+    function negModelCost = targetFun(LSP,setup)
+
+        try
+            warning off
+            ParameterFit = setup.GuessBOUNDS(1,:)  ;
+            ParameterFit(setup.SampleIndex) = 10.^LSP  ;
+            negModelCost = -setup.costFunction(ParameterFit)  ;
+            warning on
+        catch
+            negModelCost = inf  ;
+        end
+    end
+
+% define 'par' input for 'performPT.m' ------------------------------------
+
+SampleBOUNDS = setup.FitBOUNDS(:,setup.SampleIndex)  ;% parameter bounds
+par.min = log10(SampleBOUNDS(1,:))'     ;% parameter log min
+par.max = log10(SampleBOUNDS(2,:))'     ;% parameter log max
+par.number = length(setup.SampleIndex)  ;% parameter number
+
+% define 'opt' input for 'performPT.m' ------------------------------------
+
+opt.mode = 'text'  ;% only print statements
+opt.debug = false  ;% do not include debug outputs
+opt.saveEach = saveEach  ;% save every 'saveEach' iterations
+opt.saveFileName = [resultsFile '_MCMC'] ;
+
+numberChains = length(InitialRanks)      ;% number of chains
+theta0 = zeros(numberChains,par.number)  ;% initial points
+for n = 1:numberChains
+    theta0(n,:) = Results(InitialRanks(n)).ParameterFit(setup.SampleIndex);
+end
+theta0 = log10(theta0) ;
+
+opt.sigma0      = cov(theta0)   ;% initial covariance matrix
+opt.theta0      = theta0'       ;% initial points
+opt.nIterations = numSamples    ;% number of iterations
+opt.PT.nTemps   = numberChains  ;% number of chains
+
+opt.PT.exponentT      = 4      ;% see 'performPT' for description
+opt.PT.alpha          = 0.51   ;
+opt.PT.temperatureNu  = 1e3    ;
+opt.PT.memoryLength   = 1      ;
+opt.PT.regFactor      = 1e-6   ;
+opt.PT.temperatureEta = 100    ;
+opt.PT.maxT           = 5e4    ;
+
+% run MCMC ----------------------------------------------------------------
+
+res = performPT( logPostHandle, par, opt ) ;
+
+end

calibrator.m

@@ -0,0 +1,177 @@
+%% calibrator.m = calibrates model parameters via multi-start local-search
+%
+%     /--------------------------------------------------------------\
+%     | Date:         10/28/2021                                     |
+%     | Author:       Jonathan Cody                                  |
+%     | Affiliation:  Purdue University                              |
+%     |               Weldon School of Biomedical Engineering        |
+%     |               Pienaar Computational Systems Pharmacology Lab |
+%     \--------------------------------------------------------------/
+% 
+%% ========================================================================
+%	INPUTS
+%  ========================================================================
+% 
+% setup = structure variable with fields given below
+% 
+% setup.GuessBOUNDS = array for boundaries of initial parameter guesses ---
+%     EXAMPLE: GuessBOUNDS(:,p) = [minimum maximum]
+%     p = index for parameter
+%     NOTE: Enter equal minimum and maximum to fix a parameter
+% 
+% setup.FitBOUNDS = array for boundaries of parameter fitting -------------
+%     EXAMPLE: FitBOUNDS(:,p) = [minimum maximum]
+%     p = index for parameter
+%     NOTE: Enter -inf or inf for unbounded fitting of that parameter
+% 
+% setup.costFunction = handle for cost function file ----------------------
+%     Must take parameter vector and output [ model cost , extra outputs ]
+%     Extra outputs are optional (can be set to [])
+%     EXAMPLE: see 'costFun_NLL.m'
+% 
+% setup.timeLimit = number of minutes allowed for each local search
+% 
+% filePrefix = string to prefix to file names -----------------------------
+%     EXAMPLE: filePrefix = 'myprefix'
+% 
+% numberGuesses = number of initial parameter guesses ---------------------
+% 
+% numberWorkers = number of parallel processes ----------------------------
+%     NOTE: Enter [] or 0 for serial
+% 
+%% ========================================================================
+%	OUTPUTS
+%  ========================================================================
+% 
+% Results = structure array with fields given below -----------------------
+%     Orded with lowest cost result first.
+%     EXAMPLE: Results(n).ParameterFit = nth lowest cost parameter set
+%     NOTE: Results is saved in file 'prefix_Results.mat'
+% 
+%     ParameterGuess = parameter guess used to initialize optimization
+% 
+%     ParameterFit = fitted parameters obtained from optimization
+% 
+%     modelCost = model cost when using fitted parameters
+% 
+%     extraOutput = extra model outputs (as defined by costFunction)
+% 
+%     modelCalls = number of times optimization algorithm called model
+% 
+%     loopMinutes = minutes to run optimization algorithm
+% 
+%     errorMessage = error message output by optimization algorithm
+% 
+%     setup = copy of 'setup' structure (Inputs) with one addition...
+%
+%     setup.SampleIndex = index in GuessBOUNDS that references parameters
+%     that were sampled (i.e. GuessBOUND(1,p) ~= GuessBOUND(2,p) )
+% 
+%% ========================================================================
+%	FUNCTION
+%  ========================================================================
+function Results = calibrator(setup,filePrefix,numberGuesses,numberWorkers)
+
+% Declare the index in 'GuessBOUNDS' for sampled parameters
+setup.SampleIndex = find(setup.GuessBOUNDS(1,:)~=setup.GuessBOUNDS(2,:)) ;
+% Declare 'setupFunction' for 'looper.m' (see 'calibrator_setup' below)
+setupFunction = @() calibrator_setup(numberGuesses,setup.GuessBOUNDS) ;
+% Declare 'loopFunction' for 'looper.m' (see 'calibrator_loop' below)
+loopFunction  = @(result) calibrator_loop(result,setup) ;
+
+% Execute calibration using 'looper.m'
+Results = looper(filePrefix,setupFunction,loopFunction,...
+                 numberWorkers,'modelCost') ;
+
+% Save 'setup' to 'Results' file
+Results(1).setup = setup ;
+save([filePrefix '_Results.mat'],'Results')
+
+end
+%% ========================================================================
+%	'calibrator_setup' function
+%  ========================================================================
+function Results = calibrator_setup(numberGuesses,GuessBOUNDS)
+
+% Call 'sampler' to generate cell vector of parameter sets and vector
+% of sample parameter indices (within 'GuessBOUNDS')
+ParameterGuess = sampler(GuessBOUNDS,numberGuesses) ;
+
+% Initialize 'Results' structure array with variable field 'ParameterGuess'
+Results = struct('ParameterGuess',ParameterGuess) ;
+
+% Include elements in 'Results' to store calibration results
+[Results(:).ParameterFit] = deal([]) ;
+[Results(:).modelCost]    = deal([]) ;
+[Results(:).extraOutput]  = deal([]) ;
+[Results(:).modelCalls]   = deal([]) ;
+
+end
+%% ========================================================================
+%	'calibrator_loop' function
+%  ========================================================================
+function result = calibrator_loop(result,setup)
+
+% Load variables from 'result' and 'setup' structures
+ParameterGuess = result.ParameterGuess ;
+GuessBOUNDS    = setup.GuessBOUNDS ;
+FitBOUNDS      = setup.FitBOUNDS ;
+costFunction   = setup.costFunction ;
+timeLimit      = setup.timeLimit ;
+SampleIndex    = setup.SampleIndex ;
+
+% Extract sampled parameters from 'ParameterGuess' (leaving fixed)
+% and convert to log-scale
+LogGuess = log10( ParameterGuess(SampleIndex) ) ;
+
+% Extract sampled parameter fit boundaries from 'FitBOUNDS' (leaving fixed)
+% and convert to log-scale (avoiding -Inf entries)
+LogBOUNDS = FitBOUNDS(:,SampleIndex) ;
+LogBOUNDS(LogBOUNDS~=(-Inf)) = log10( LogBOUNDS(LogBOUNDS~=(-Inf)) ) ;
+
+% Initialize number of model calls and start calibration timer
+modelCalls = 0 ;
+calibrationTime = tic ;
+
+% Fit parameters using optimization function 'fmincon' (see MATLAB help)
+% based on target function 'targetfun' (see below)
+LogFit = fmincon( @(P)targetFun(P),LogGuess,[],[],[],[],...
+                  LogBOUNDS(1,:),LogBOUNDS(2,:),[],...
+                  optimoptions(@fmincon,'Display','none')  ) ;
+
+% Retrieve model cost, parameter values, costs for each data type, and
+% model output for fitted parameter set
+[modelCost,ParameterFit,extraOutput] = targetFun(LogFit) ;
+
+% Store results in structure 'result'
+result.ParameterFit = ParameterFit ;
+result.modelCost    = modelCost ;
+result.extraOutput  = extraOutput ;
+result.modelCalls   = modelCalls ;
+
+%  ------------------------------------------------------------------------
+%   'targetfun' function
+%  ------------------------------------------------------------------------
+function [modelCost,ParameterFit,extraOutput] = targetFun(P) 
+
+    % If time limit is exceeded, throw error
+    if toc(calibrationTime) >= 60*timeLimit
+        error('Optimization algorithm timed out.')
+    end
+
+    % Increment number of model calls
+    modelCalls = modelCalls + 1 ;
+
+    % Declare parameter vector
+    ParameterFit = GuessBOUNDS(1,:) ;
+    ParameterFit(SampleIndex) = 10.^(P) ;
+
+    % Use a try-catch statement to continue if an error is thrown
+    try
+        [modelCost,extraOutput] = costFunction(ParameterFit) ;
+    catch
+        modelCost   = nan ;
+        extraOutput = [] ;
+    end
+end
+end