From 4f41258e50c6e0178428912c861a121091692838 Mon Sep 17 00:00:00 2001 From: "Cody, Jonathan William" Date: Mon, 9 Jan 2023 10:55:59 -0600 Subject: [PATCH] Add files via upload --- N803_collector_2.m | 87 +++++++++++++++++++++++++++++++--------------- 1 file changed, 59 insertions(+), 28 deletions(-) diff --git a/N803_collector_2.m b/N803_collector_2.m index dcc14e8..69def17 100644 --- a/N803_collector_2.m +++ b/N803_collector_2.m @@ -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]' ; @@ -105,6 +121,22 @@ 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 % ======================================================================== @@ -112,10 +144,11 @@ % 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') @@ -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,:)) ; @@ -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 %% ========================================================================