Add files via upload

Analyzer/MCMCsampler.m

@@ -1,7 +1,7 @@
 % MCMCsampler.m - performs markov chain monte carlo sampling of parameters
 %
 %     /--------------------------------------------------------------\
-%     | Date:         02/16/2022                                     |
+%     | Date:         03/07/2023                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Pienaar Computational Systems Pharmacology Lab |
 %     |               Weldon School of Biomedical Engineering        |
@@ -102,7 +102,7 @@
 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'] ;
+opt.saveFileName = [resultsFile '_MCMC_fullsave'] ;
 
 numberChains = length(InitialRanks)      ;% number of chains
 theta0 = zeros(numberChains,par.number)  ;% initial points
@@ -124,8 +124,20 @@
 opt.PT.temperatureEta = 100    ;
 opt.PT.maxT           = 5e4    ;
 
-% run MCMC ----------------------------------------------------------------
+% run MCMC and save results -----------------------------------------------
 
 res = performPT( logPostHandle, par, opt ) ;
 
-end
+opt.par = par  ;% adding parameter space to options structure
+res.opt = opt  ;% adding options structure to results structure
+
+% Thinning results by every 100th sample
+thinID = 1:100:numSamples  ;
+res.par     = res.par    (:,thinID)  ;
+res.logPost = res.logPost(  thinID)  ;
+
+save([resultsFile '_MCMC.mat'],'res')
+disp([resultsFile ' done!'])
+disp(' ')
+
+end

Analyzer/calibrator.m

@@ -174,4 +174,4 @@
         extraOutput = [] ;
     end
 end
-end
+end

Analyzer/calibrator_plot.m

@@ -1,7 +1,7 @@
 %% calibrator_plot.m = plots results from calibrator.m
 %
 %     /--------------------------------------------------------------\
-%     | Date:         08/09/2022                                     |
+%     | Date:         02/04/2023                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Purdue University                              |
 %     |               Weldon School of Biomedical Engineering        |
@@ -24,8 +24,13 @@
 %     EXAMPLE: X{i} = [column vector] and Y{i} = [column vector] 
 %
 %     modelFunction = handle for calling model function file
-%     NOTE: Must return a matrix where each column is a different output
-%     EXAMPLE: [ data type 1 column vector , data type 2 column vector ]
+%     NOTE: Must return a cell array where each element is a column vector.
+%     (output){i} must be same size as Xeval{i} 
+%
+%     Xeval = cell arrays for model plot x-values
+%     EXAMPLE: Xeval{i} = [column vector]
+%
+% Fields below are optional:
 %
 %     CensorID = logical array showing censored data (same size as Ydata)
 %                (use 'prepper.m' to generate)
@@ -72,6 +77,7 @@ function calibrator_plot(Results,plotSetup,figName)
 % Load required variables from input structures
 Xdata = plotSetup.Xdata ;
 Ydata = plotSetup.Ydata ;
+Xeval = plotSetup.Xeval ;
 modelFunction = plotSetup.modelFunction ;
 GuessBOUNDS = Results(1).setup.GuessBOUNDS ;
 SampleIndex = Results(1).setup.SampleIndex ;
@@ -88,7 +94,6 @@ function calibrator_plot(Results,plotSetup,figName)
 YLabels   = []  ;% {'y1-axis label','y2-axis label',etc}
 XRANGE    = []  ;% x-axis ranges [ x1min x1max ; x2min x2max ; etc ]
 YRANGE    = []  ;% y-axis ranges [ y1min y1max ; y2min y2max ; etc ]
-xPrecise  = .1  ;% space between x values for plotting
 
 if isfield(plotSetup,'CensorID')  ; CensorID  = plotSetup.CensorID  ; end
 if isfield(plotSetup,'Fits2Plot') ; Fits2Plot = plotSetup.Fits2Plot ; end
@@ -101,7 +106,6 @@ function calibrator_plot(Results,plotSetup,figName)
 if isfield(plotSetup,'YLabels')   ; YLabels   = plotSetup.YLabels   ; end
 if isfield(plotSetup,'XRANGE')    ; XRANGE    = plotSetup.XRANGE    ; end
 if isfield(plotSetup,'YRANGE')    ; YRANGE    = plotSetup.YRANGE    ; end
-if isfield(plotSetup,'xPrecise')  ; xPrecise  = plotSetup.xPrecise; end
 
 % Reformat data variables (if needed)
 if ~iscell(Xdata) ; Xdata = {{Xdata}} ; end
@@ -186,28 +190,13 @@ function calibrator_plot(Results,plotSetup,figName)
                     ./ NormGuessRange ;
 end
 
-% Generate Xeval using minimum and maximum of data x-values
-xMin = min(Xdata{1}{1}) ;
-xMax = max(Xdata{1}{1}) ;
-for n = 1:length(Xdata)
-    for m = 1:length(Xdata{n})
-        xMin = min([ min(Xdata{n}{m}) xMin ]) ;
-        xMax = max([ max(Xdata{n}{m}) xMax ]) ;
-    end
-end
-if ~isempty(XRANGE)
-    xMin = min(XRANGE(:,1)) ;
-    xMax = max(XRANGE(:,2)) ;
-end
-Xeval = xMin:xPrecise:xMax ;
-
 % Evaluate model at Xeval
 YPlots = cell(nPlots,1) ;
 ErrorSkips = zeros(nPlots,1) ;
 for m = 1:nPlots
     try
         Pars = Results(FitIndex(Fits2Plot(m))).ParameterFit ;
-        YPlots{m} = modelFunction(Xeval,Pars) ;
+        YPlots{m} = modelFunction(Pars) ;
     catch
         ErrorSkips(m) = 1 ;
     end
@@ -256,13 +245,13 @@ function calibrator_plot(Results,plotSetup,figName)
             disp(['Skipped plot number ' num2str(m) ' due to error.'])
             continue
         end
-        Yplot = YPlots{m}(:,n) ;
+        Yplot = YPlots{m}{n} ;
         lineIndex = wrapindex(m,length(ModLines)) ;
         if length(Fits2Plot) < nFits2Sum
-            plot(Xeval,Yplot,ModLines{lineIndex},'LineWidth',1.5)
+            plot(Xeval{n},Yplot,ModLines{lineIndex},'LineWidth',1.5)
         else
             % Otherwise, use grayscale lines
-            plot(Xeval,Yplot,'Color',GrayScale{m},'LineWidth',1.5)
+            plot(Xeval{n},Yplot,'Color',GrayScale{m},'LineWidth',1.5)
         end
     end
 
@@ -348,4 +337,4 @@ function calibrator_plot(Results,plotSetup,figName)
 
 %#ok<*AGROW>
 
-end
+end

Analyzer/costFun_NLL.m

@@ -1,7 +1,7 @@
 %% costFun_NLL = model cost calculator (NLL,SSE,MSE w/additional options)
 %
 %     /--------------------------------------------------------------\
-%     | Date:         06/25/2022                                     |
+%     | Date:         05/06/2023                                     |
 %     | Author:       Jonathan Cody                                  |
 %     | Affiliation:  Purdue University                              |
 %     |               Weldon School of Biomedical Engineering        |
@@ -15,80 +15,75 @@
 % Parameters = vector of parameters for 'modelFunction' -------------------
 %
 % modelFunction = handle for calling model function -----------------------
-%     EXAMPLE: modelFunction = @(Time,Parameters) ...
+%     EXAMPLE: modelFunction = @(Parameters) ...
 %              yourFunctionFilename(Time,Parameters,extraInputs)
-%     NOTE: Should return a matrix where each column is a different output
-%     EXAMPLE: [ data type 1 column vector , data type 2 column vector ]
 %
-% Xdata/Ydata = cell arrays for data x-values and y-values ----------------
-%     EXAMPLE: X{i} = [column vector] and Y{i} = [column vector] 
-%     i = index for data type (e.g. different cells or molecules)
+% Ydata = vector for data y-values ----------------------------------------
+%
+% NOTE: For multiple data types, Y should be a cell vector, where each
+% element is a data type. 'modelFunction' should return an analgous
+% cell vector.
 %
 %% ========================================================================
 %	OPTIONS
 %  ========================================================================
 % NOTE: 'varargin' is used to add optional inputs
-%     EXAMPLE: costFun_NLL(P,modFun,X,Y,'Input Name',[Input Value])
+%     EXAMPLE: costFun_NLL(P,modFun,Y,'Input Name',[Input Value])
 %
 % censorSide = string inidicating direction of censoring ------------------
 %     (details pending, see code)
 %
 % CensorIndex = logical index indicating data that is censored ------------
-%     NOTE: Must be the same shape as Xdata/Ydata
+%     NOTE: Must be the same shape as Ydata
 %
 % costFunction = string inidicating form of cost function -----------------
 %     DEFAULT: negative log-likelihood
 %
-% TypeWeights = vector of weights to be applied to each data subject cost
-%
-% Xevaluate = vector of x-values to evaluate model function
-%     DEFAULT: all unique values in Xdata
+% CostWeights = vector of weights to be applied to each data subject cost
 %
-% Xpost = vector of x-values corresponding to post-data region
-%     NOTE: Xpost is a cell array (one element for each data type)
-%     DEFAULT: Xpost is empty
+% ScalingFacts = vector of scaling factors to be applied to each subject
+%     model/data BEFORE cost function is calculated
 %
-% YpostFun = special cost function handles applied to post-data region
-%     NOTE: YpostFun is a cell array (one element for each data type)
-%     DEFAULT: YpostFun is empty cell array
+% XtraCostFuns = handle(s) for extra cost function(s) (additive)
+%     NOTE: Two inputs are passed: ( Parameters, modelFunction output )
+%     NOTE: Use a cell array for multiple extra cost functions
+%     DEFAULT: XtraCostFuns is empty cell array
 %
 %% ========================================================================
 %	OUTPUTS
 %  ========================================================================
 %
 % cost = scalar model cost ------------------------------------------------
 %
-% TypeCosts = vector of costs for each data type (with TypeWeights applied)
+% TypeCosts = vector of costs for each data type (with CostWeights applied)
+%     NOTE: Costs from XtraCostFuns are appended, if used
 %
 %% ========================================================================
 %	FUNCTION
 %  ========================================================================
 function [cost,TypeCosts] = ...
-    costFun_NLL(Parameters,modelFunction,Xdata,Ydata,varargin)
+    costFun_NLL(Parameters,modelFunction,Ydata,varargin)
 
 % Reformat data variables (if needed)
-if ~iscell(Xdata) ; Xdata = {Xdata} ; end
 if ~iscell(Ydata) ; Ydata = {Ydata} ; end
 
 nTypes = length(Ydata) ;% number of data types
 
 % Declare optional input default values (then overwrite, if specified)
 censorSide   = 'left'            ;% direction of censoring
-CensorIndex  = cell(size(Xdata)) ;% indices in 'Xdata' for censored data
+CensorIndex  = cell(size(Ydata)) ;% indices in 'Xdata' for censored data
 costFunction = 'NLL'             ;% type of cost function
-TypeWeights  = ones(1,nTypes)    ;% cost weights for different data types
-Xevaluate = unique(cat(1,Xdata{:})) ;% x-values to evaluate model function
-Xpost    = []                ;% x-values for post-data behavior
-YpostFun = cell(size(Xdata)) ;% allowable range for post-data behavior
+CostWeights  = ones(1,nTypes)    ;% cost weights for each data type
+ScalingFacts = ones(1,nTypes)    ;% model/data scaling factors
+XtraCostFuns = []                ;% additional user-defined cost functions
 
 for n = 1:length(varargin)
 if strcmp('censorSide',  varargin(n)) ; censorSide   = varargin{n+1} ; end
 if strcmp('CensorIndex', varargin(n)) ; CensorIndex  = varargin{n+1} ; end
 if strcmp('costFunction',varargin(n)) ; costFunction = varargin{n+1} ; end
-if strcmp('TypeWeights', varargin(n)) ; TypeWeights  = varargin{n+1} ; end
-if strcmp('Xevaluate',   varargin(n)) ; Xevaluate    = varargin{n+1} ; end
-if strcmp('Xpost',       varargin(n)) ; Xpost        = varargin{n+1} ; end
-if strcmp('YpostFun',    varargin(n)) ; YpostFun     = varargin{n+1} ; end
+if strcmp('CostWeights', varargin(n)) ; CostWeights  = varargin{n+1} ; end
+if strcmp('ScalingFacts',varargin(n)) ; ScalingFacts = varargin{n+1} ; end
+if strcmp('XtraCostFuns',varargin(n)) ; XtraCostFuns = varargin{n+1} ; end
 end
 
 % Declare cost function handle
@@ -108,47 +103,44 @@
 end
 
 % Run model
-ModelOUTPUTS = modelFunction(Xevaluate,Parameters) ;
+Ymod = modelFunction(Parameters) ;
+
+% Reformat variables (if needed)
+if ~iscell(CensorIndex) ; CensorIndex = {CensorIndex} ; end
+if ~iscell(Ymod)        ; Ymod        = {Ymod}        ; end
 
 % Calculate cost for each data type
-TypeCosts = TypeWeights ;
+TypeCosts = CostWeights ;
 for n = 1:nTypes
 
-    % Interpolate model outputs for 'Xdata' from outputs for 'Xevaluate'
-    Ymod = interp1(Xevaluate,ModelOUTPUTS(:,n),Xdata{n}) ;
+    if ~all( size(Ymod{n}) == size(Ydata{n}) ) 
+        error('Size of model output does not equal size of data.')
+    end
 
     % Account for censored data (if requested) by equating output to
     % left/right- censored data points when the model is below/above data
     % OR omit censored data from calibration entirely
-    if ~isempty(CensorIndex)
-        if any(CensorIndex{n})
-            if     strcmp(censorSide,'left' ) ; indx = Ymod<=Ydata{n} ;
-            elseif strcmp(censorSide,'right') ; indx = Ymod>=Ydata{n} ;
-            else                              ; indx = Ymod==Ymod ;
-            end
-            indx = indx & CensorIndex{n} ;
-            Ymod(indx) = Ydata{n}(indx) ;
+    if any(CensorIndex{n})
+        if     strcmp(censorSide,'left' ) ; indx = Ymod{n}<=Ydata{n} ;
+        elseif strcmp(censorSide,'right') ; indx = Ymod{n}>=Ydata{n} ;
+        else                              ; indx = Ymod{n}==Ymod{n} ;
         end
+        indx = indx & CensorIndex{n} ;
+        Ymod{n}(indx) = Ydata{n}(indx) ;
     end
 
     % Calculate cost for data type (with weights)
-    TypeCosts(n) = costHandle(Ymod,Ydata{n})*TypeWeights(n) ;
+    TypeCosts(n) = costHandle( Ymod{n} *ScalingFacts(n) ...
+                             , Ydata{n}*ScalingFacts(n) ) *CostWeights(n) ;
 end
 
-% Applying post-data cost function (if requested)
-if ~isempty(Xpost)
-    PostCosts = TypeWeights ;
-    for n = 1:nTypes
-        if isempty(YpostFun{n})
-            PostCosts(n) = 0 ;
-        else
-            Ypost = interp1(Xevaluate,ModelOUTPUTS(:,n),Xpost{n}) ;
-            PostCosts(n) = YpostFun{n}(Ypost) ;
-        end
+% Applying extra cost functions (if requested)
+if ~isempty(XtraCostFuns)
+    for n = 1:length(XtraCostFuns)
+        TypeCosts(nTypes+n) = XtraCostFuns{n}(Parameters,Ymod) ;
     end
-    TypeCosts = [TypeCosts PostCosts] ;
 end
 
 cost = sum(TypeCosts) ;% total cost
 
-end
+end

Analyzer/looper.m

@@ -123,4 +123,4 @@
 
 save( fullfile('TEMP',['loop_' num2str(n) '.mat']) ,'result') ;
 
-end
+end

Analyzer/performPT.m

@@ -1,11 +1,7 @@
 % Source: https://github.com/ICB-DCM/PESTO/blob/master/private/performPT.m
-%
-% Two edits from source made by Jonathan Cody:
-%
-% #1 Added function 'updateStatistics' as subfunction (line 325)
+
+% Added function 'updateStatistics' as subfunction
 % https://github.com/ICB-DCM/PESTO/blob/master/private/updateStatistics.m
-%
-% #2 Progress printing frequency decreased (line 195)
 
 function res = performPT( logPostHandle, par, opt )
 
@@ -348,4 +344,4 @@
 
 end
 
-end
+end

Analyzer/prepper.m

@@ -121,4 +121,4 @@
     CensorIndex{n} = cat(1,CensorIndex{n}{:}) ;
 end
 
-end
+end

Analyzer/sampler.m

@@ -90,4 +90,4 @@
     ParameterSets{n} = ParameterSETS(n,:) ;
 end
 
-end
+end