Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@cody2
cody2 authored Jan 9, 2023
1 parent ec6f24d commit 4f41258
Showing 1 changed file with 59 additions and 28 deletions.
87 changes: 59 additions & 28 deletions N803_collector_2.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,9 +38,25 @@
Dose{5} = 0 ;
Dose{6} = 7*[0 1 2 7] ;

% Save the parameter sample or initial condition sample
SavePars{1} = true ; SaveYics{1} = true ;
SavePars{2} = false ; SaveYics{2} = true ;
SavePars{3} = false ; SaveYics{3} = false ;
SavePars{4} = false ; SaveYics{4} = false ;
SavePars{5} = false ; SaveYics{5} = false ;
SavePars{6} = false ; SaveYics{6} = false ;

% Save the following outputs (see list below and in 'N803_model_2.m')
SaveOutput{1} = [1 2 4 27 29 39 49 51 53 58 60] ;
SaveOutput{2} = [1 2 4 27 29 39 49 51 53 58 60] ;
SaveOutput{3} = [1 2] ;
SaveOutput{4} = [1 2] ;
SaveOutput{5} = 21 ;
SaveOutput{6} = [3 21] ;

% Descriptions of additional model outputs (see 'N803_model_2.m')
Output{01} = 'virions [log fold change]' ;
Output{02} = 'CD8+ T cells [fold change]' ;
Output{02} = 'CD8+ T cells [fold change]' ;
Output{03} = 'V [log wrt C1 pre]' ;
Output{04} = 'S_0 [log wrt C1 pre]' ;
Output{05} = 'S_1 [log wrt C1 pre]' ;
Expand Down Expand Up @@ -105,17 +121,34 @@
Output{64} = 'R generation by aN [log wrt C1 pre]' ;
Output{65} = 'R gen by aN fraction' ;

% Add figure labels to the following outputs (for Scenarios 1-6)
Figure = cell(size(Output)) ;
Figure{01} = {'2A','2D','S3A','S3D', [], []} ;
Figure{02} = {'2B','2E','S3B','S3E', [], []} ;
Figure{03} = { [], [], [], [], [],'S4C'} ;
Figure{04} = {'4D','4D', [], [], [], []} ;
Figure{21} = { [], [], [], [],'S4A','S4D'} ;
Figure{27} = {'3A','3C', [], [], [], []} ;
Figure{29} = {'4B','4B', [], [], [], []} ;
Figure{39} = {'4F','4F', [], [], [], []} ;
Figure{49} = {'4I','4I', [], [], [], []} ;
Figure{51} = {'4H','4H', [], [], [], []} ;
Figure{53} = {'4G','4G', [], [], [], []} ;
Figure{58} = {'4A','4A', [], [], [], []} ;
Figure{60} = {'4C','4C', [], [], [], []} ;

%% ========================================================================
% BODY
% ========================================================================

% Percolate results structure
N = length(Dose) ;
C = struct('Scenario',Scenario) ;
[C(:).X] = deal(X{:}) ;% x-values (days) corresponding to outputs
[C(:).X] = deal(X{:}) ;% x-values (days) corresponding to outputs
[C(:).Dose] = deal(Dose{:}) ;% x-values (days) for doses given
[C(:).YSummary] = deal([]) ;% output summary statistics
[C(:).ParSample] = deal([]) ;% parameter samples
[C(:).ParSample] = deal([]) ;% parameter samples from Bayesian MCMC
[C(:).YicSample] = deal([]) ;% initial condition samples
[C(:).YSummary] = deal([]) ;% output summary statistics (95% CI)

% load(name,'C')

Expand Down Expand Up @@ -144,22 +177,25 @@
'fullOut',true,'ivDosing',880*[.06 .06 .06 1 0]/.03) ;
end

% Load model parameters, etc., from MSLS and MCMC data files
SamID = nThin:nThin:nIter ;% thin sample index (within MCMC)
SamID(end) = nIter-1 ;% using 2nd-to-last because of an MCMC glitch
nSam = length(SamID) ;% number of thinned MCMC samples

% Load parameters for initial condition and parameter calculation
% functions: 'N803_shared.m', 'N803_single.m'
warning off
load(saveFile,'Results')
load([saveFile '_MCMC'],'res')
warning on
SamID = nThin:nThin:nIter ;% thin sample index (within MCMC)
SamID(end) = nIter-1 ;% using 2nd-to-last because of an MCMC glitch !!
nSam = length(SamID) ;% thinned number of MCMC samples
PARS = repmat(Results(1).ParameterFit,nSam,1) ;% sample of parameters
ParID = Results(1).setup.SampleIndex ;% index within par array
PARS(:,ParID) = 10.^res.par(:,SamID)' ;% sample of parameters
ParID = Results(1).setup.SampleIndex ;
PARS(:,ParID) = 10.^res.par(:,SamID)' ;

% Evaluate model and collect terms
% Evaluate model and collect outputs, along with parameter and initial
% conditions (ICs) for 'N803_model_2.m' (called by 'shared' & 'single')
Ysam = cell(nSam,1) ;% sample of model outputs
Pcoh = cell(nSam,1) ;% sample of cohort parameters
passID = true(nSam,1) ;% index of models that did not fail
Pcoh = cell(nSam,1) ;% sample of cohort [ parameters , ICs ]
passID = true(nSam,1) ;% index of models that did not fail
parfor (sam = 1:nSam, nWork)
try
[Ysam{sam},Pcoh{sam}] = modelFun(PARS(sam,:)) ;
Expand All @@ -173,42 +209,37 @@
nSam = sum(passID) ;
Ysam = Ysam(passID) ;
Pcoh = Pcoh(passID) ;
PCOH = cat(1,Pcoh{:}) ;% sample of calculated parameters

C(n).NLL = - res.logPost( SamID(passID) ) ;
PCOH = cat(1,Pcoh{:}) ;
if SavePars{n} ; C(n).ParSample = PCOH(:, 1:31 ) ; end
if SaveYics{n} ; C(n).YicSample = PCOH(:,32:end) ; end

% Store confidence regions (OR samples)
OUTPUTS = zeros(length(X{n}),nSam) ;% storing output sample (temp)
nOuts = size(Ysam{1},2) ;% number of outputs
Ystats = struct('Output',Output) ;
[Ystats(:).min] = deal([]) ;% minimum values across samples
[Ystats(:).lo95] = deal([]) ;% 2.5th percentile
[Ystats(:).med] = deal([]) ;% median
[Ystats(:).hi95] = deal([]) ;% 97.5th percentile
[Ystats(:).max] = deal([]) ;% maximum
[Ystats(:).Figure] = deal([]) ;% figure label
[Ystats(:).lo95] = deal([]) ;% 2.5th percentile
[Ystats(:).hi95] = deal([]) ;% 97.5th percentile

if n == 1
NORMFACTS = ones(nOuts,nSam) ;% matrix of normalization factors
end
for out = 1:nOuts
for out = SaveOutput{n}
for sam = 1:nSam
OUTPUTS(:,sam) = Ysam{sam}(:,out) ;
end
if any( n == 1:2 ) && any( out == [ 3:15 , 18:19 , 29:61 , 63:64 ] )
if any( n == 1:2 ) && any( out == [ 3:15, 18:19, 29:61, 63:64 ] )
if n == 1
NORMFACTS(out,:) = OUTPUTS(1,:) ;
end
OUTPUTS = OUTPUTS - NORMFACTS(out,:) ;
end
Ystats(out).min = prctile(OUTPUTS', 0 )' ;
Ystats(out).Figure = Figure{out}{n} ;
Ystats(out).lo95 = prctile(OUTPUTS', 2.5)' ;
Ystats(out).med = prctile(OUTPUTS', 50 )' ;
Ystats(out).hi95 = prctile(OUTPUTS', 97.5)' ;
Ystats(out).max = prctile(OUTPUTS',100 )' ;
end

C(n).YSummary = Ystats ;% output summary statistics
C(n).ParSample = PCOH ;% sample of calculated parameters
C(n).YSummary = Ystats ;% output summary statistics
end

%% ========================================================================
Expand Down

0 comments on commit 4f41258

Please sign in to comment.