Skip to content

Commit

Permalink
Merge branch 'AnalyticalSensitivityCalculations'
Browse files Browse the repository at this point in the history
  • Loading branch information
@Munsky
Munsky committed Nov 3, 2023
2 parents 1332e11 + affa412 commit 39acae1
Show file tree
Hide file tree
Showing 27 changed files with 18,608 additions and 16,514 deletions.
83 changes: 51 additions & 32 deletions CommandLine/SSIT.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,11 @@
end


function obj = formPropensitiesGeneral(obj,prefixName)
function obj = formPropensitiesGeneral(obj,prefixName,computeSens)
arguments
obj
prefixName = [];
computeSens = true;
end
% This function starts the process to write m-file for each
% propensity function.
Expand All @@ -103,11 +104,14 @@
end

if obj.useHybrid
PropensitiesGeneral = ssit.Propensity.createAsHybridVec(sm, obj.stoichiometry,...
obj.parameters, obj.species, obj.hybridOptions.upstreamODEs, logicTerms, prefixName);
PropensitiesGeneral = ...
ssit.Propensity.createAsHybridVec(sm, obj.stoichiometry,...
obj.parameters, obj.species, obj.hybridOptions.upstreamODEs,...
logicTerms, prefixName, computeSens);
else
PropensitiesGeneral = ssit.Propensity.createAsHybridVec(sm, obj.stoichiometry,...
obj.parameters, obj.species, [], logicTerms, prefixName);
PropensitiesGeneral = ...
ssit.Propensity.createAsHybridVec(sm, obj.stoichiometry,...
obj.parameters, obj.species, [], logicTerms, prefixName, computeSens);
end

obj.propensitiesGeneral = PropensitiesGeneral;
Expand Down Expand Up @@ -855,28 +859,25 @@ function summarizeModel(obj)
pause;
end
case 'fspSens'
if strcmp(obj.sensOptions.solutionMethod,'forward')&&isempty(obj.propensitiesGeneral{1}.sensTimeFactorVec)
obj = formPropensitiesGeneral(obj,'Sensitivities',true);
end
if ~isempty(obj.parameters)
model = ssit.SrnModel(obj.stoichiometry,...
obj.propensityFunctions,...
obj.parameters(:,1),...
obj.inputExpressions);
app.ReactionsTabOutputs.parameters = obj.parameters(:,1);
else
model = ssit.SrnModel(obj.stoichiometry,...
obj.propensityFunctions,...
[],...
obj.inputExpressions);
app.ReactionsTabOutputs.parameters = [];
end
app.ReactionsTabOutputs.varNames = obj.species;

[Solution.sens, bConstraints] = ...
ssit.sensitivity.computeSensitivity(model,...
ssit.sensitivity.computeSensitivity(...
obj.parameters,...
obj.propensitiesGeneral,...
obj.tSpan,...
obj.fspOptions.fspTol,...
obj.initialCondition,...
1.0,...
obj.initialProbs,...
obj.stoichiometry, ...
obj.fspConstraints.f,...
obj.fspConstraints.b,...
[], obj.fspOptions.verbose, 0,...
Expand All @@ -895,9 +896,6 @@ function summarizeModel(obj)
% Solution.plotable = exportSensResults(app);

case 'ode'

% specificPropensities = SSIT.parameterizePropensities(obj.propensitiesGeneral,[obj.parameters{:,2}]');

[~,Solution.ode] = ssit.moments.solveOde2(obj.initialCondition, obj.tSpan, ...
obj.stoichiometry, obj.propensitiesGeneral, [obj.parameters{:,2}]', obj.fspOptions.initApproxSS);
end
Expand Down Expand Up @@ -939,7 +937,6 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
for it=1:Nt
A.time((it-1)*obj.ssaOptions.nSimsPerExpt+1:it*obj.ssaOptions.nSimsPerExpt) = obj.tSpan(it);
for ie = 1:obj.ssaOptions.Nexp
% for is=1:obj.ssaOptions.nSimsPerExpt
for s = 1:size(Solution.trajs,1)
warning('off')
A.(['exp',num2str(ie),'_s',num2str(s)])((it-1)*obj.ssaOptions.nSimsPerExpt+(1:obj.ssaOptions.nSimsPerExpt)) = ...
Expand All @@ -949,7 +946,6 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
Solution.trajsDistorted(s,it,(ie-1)*obj.ssaOptions.nSimsPerExpt+(1:obj.ssaOptions.nSimsPerExpt));
end
end
% end
end
end
writetable(A,saveFile)
Expand Down Expand Up @@ -2258,35 +2254,54 @@ function makeFitPlot(obj,fitSolution,smoothWindow,fignums)
makeSeparatePlotOfData(fitSolution,smoothWindow,fignums)
end

function plotMHResults(obj,mhResults,FIM)
function plotMHResults(obj,mhResults,FIM,fimScale,mhPlotScale,scatterFig)
arguments
obj
mhResults = [];
FIM =[];
fimScale = 'lin';
mhPlotScale = 'log10';
scatterFig = [];
end

parNames = obj.parameters(obj.fittingOptions.modelVarsToFit,1);

if ~isempty(FIM)
pars = [obj.parameters{obj.fittingOptions.modelVarsToFit,2}];

parsLog = log10(pars);

if strcmp(mhPlotScale,'log10')
parsScaled = log10(pars);
elseif strcmp(mhPlotScale,'log')
parsScaled = log(pars);
elseif strcmp(mhPlotScale,'lin')
parsScaled = pars;
end

if ~iscell(FIM)
FIM = diag(pars)*...
FIM(obj.fittingOptions.modelVarsToFit,obj.fittingOptions.modelVarsToFit)*...
diag(pars);
covFIM{1} = FIM^(-1)/log(10)^2;
else
for i=1:length(FIM)
FIM = {FIM};
end
% FIM = diag(pars)*...
% FIM(obj.fittingOptions.modelVarsToFit,obj.fittingOptions.modelVarsToFit)*...
% diag(pars);
% covFIM{1} = FIM^(-1)/log(10)^2;
% else
for i=1:length(FIM)
if strcmp(fimScale,'lin')
FIMi = diag(pars)*...
FIM{i}(obj.fittingOptions.modelVarsToFit,obj.fittingOptions.modelVarsToFit)*...
diag(pars);
else
FIMi = FIM{i};
end
if strcmp(mhPlotScale,'log10')
covFIM{i} = FIMi^(-1)/log(10)^2;
else
covFIM{i} = FIMi^(-1);
end
end
% end
end


if ~isempty(mhResults)
% Make figures for MH convergence
figure;
Expand All @@ -2306,7 +2321,11 @@ function plotMHResults(obj,mhResults,FIM)
ylabel('Auto-correlation')
title('MH Convergence')

figure
if isempty(scatterFig)
figure
else
figure(scatterFig);
end
[valDoneSorted,J] = sort(mhResults.mhValue);
smplDone = mhResults.mhSamples(J,:);
Np = size(mhResults.mhSamples,2);
Expand All @@ -2325,7 +2344,7 @@ function plotMHResults(obj,mhResults,FIM)
end
if ~isempty(FIM)
for iFIM = 1:length(covFIM)
ssit.parest.ellipse(parsLog([j,i]),icdf('chi2',0.9,2)*covFIM{iFIM}([j,i],[j,i]),fimCols{mod(iFIM,5)+1},'linewidth',2)
ssit.parest.ellipse(parsScaled([j,i]),icdf('chi2',0.9,2)*covFIM{iFIM}([j,i],[j,i]),fimCols{mod(iFIM,5)+1},'linewidth',2)
end
end
if ~isempty(mhResults)
Expand Down
16 changes: 16 additions & 0 deletions Examples/FakeExperiment_1.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
time,exp1_s1,exp1_s2
0,0,0
0,0,0
0,0,0
0,0,0
0,0,0
5,38,19
5,58,12
5,9,81
5,79,16
5,18,65
10,11,81
10,12,54
10,11,57
10,9,80
10,89,9
31 changes: 31 additions & 0 deletions Examples/FakeExperiment_2.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
time,exp1_s1,exp1_s2
0,0,0
0,0,0
0,0,0
0,0,0
0,0,0
5,38,19
5,58,12
5,9,81
5,79,16
5,18,65
10,11,81
10,12,54
10,11,57
10,9,80
10,89,9
0.5,21,26
0.5,18,27
0.5,21,23
0.5,26,22
0.5,31,16
0.5,31,19
0.5,33,14
0.5,20,31
1.5,27,36
10,6,77
10,102,6
10,13,71
10,63,17
10,69,15
10,90,13
46 changes: 46 additions & 0 deletions Examples/FakeExperiment_3.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
time,exp1_s1,exp1_s2
0,0,0
0,0,0
0,0,0
0,0,0
0,0,0
5,38,19
5,58,12
5,9,81
5,79,16
5,18,65
10,11,81
10,12,54
10,11,57
10,9,80
10,89,9
0.5,21,26
0.5,18,27
0.5,21,23
0.5,26,22
0.5,31,16
0.5,31,19
0.5,33,14
0.5,20,31
1.5,27,36
10,6,77
10,102,6
10,13,71
10,63,17
10,69,15
10,90,13
0.5,14,33
0.5,31,19
0.5,28,25
0.5,34,17
0.5,15,33
1.5,15,56
1.5,33,26
10,20,38
10,76,12
10,18,62
10,102,14
10,12,75
10,73,14
10,8,73
10,38,25
61 changes: 61 additions & 0 deletions Examples/FakeExperiment_4.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
time,exp1_s1,exp1_s2
0,0,0
0,0,0
0,0,0
0,0,0
0,0,0
5,38,19
5,58,12
5,9,81
5,79,16
5,18,65
10,11,81
10,12,54
10,11,57
10,9,80
10,89,9
0.5,21,26
0.5,18,27
0.5,21,23
0.5,26,22
0.5,31,16
0.5,31,19
0.5,33,14
0.5,20,31
1.5,27,36
10,6,77
10,102,6
10,13,71
10,63,17
10,69,15
10,90,13
0.5,14,33
0.5,31,19
0.5,28,25
0.5,34,17
0.5,15,33
1.5,15,56
1.5,33,26
10,20,38
10,76,12
10,18,62
10,102,14
10,12,75
10,73,14
10,8,73
10,38,25
0.5,25,26
0.5,20,23
0.5,32,20
2,17,59
2,11,80
2,40,24
2,30,25
10,69,13
10,12,72
10,12,77
10,7,82
10,13,80
10,80,16
10,15,60
10,16,73
Loading

0 comments on commit 39acae1

Please sign in to comment.