Skip to content

Commit

Permalink
Merge pull request #2 from MunskyGroup/main
Browse files Browse the repository at this point in the history 
Catchup merge from main
  • Loading branch information
@dmitrisvetlov
dmitrisvetlov authored Nov 4, 2023
2 parents be4487f + 39acae1 commit 8022a40
Show file tree
Hide file tree
Showing 28 changed files with 18,718 additions and 16,513 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 8022a40

Please sign in to comment.