Skip to content

Commit

Permalink
Merge branch 'EricsBranch'
Browse files Browse the repository at this point in the history 
This merge fixes a couple bugs in the plotting routines and improves the handling of logical functions in the propensities.
  • Loading branch information
@Munsky
Munsky committed Feb 7, 2024
1 parent a8add7a commit 6e7a990
Show file tree
Hide file tree
Showing 11 changed files with 339 additions and 209 deletions.
6 changes: 5 additions & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,4 +48,8 @@ WorkSpace/OldProablyTrash/*/*
WorkSpace/OldProablyTrash/*/*/*
*/*/*/tmpPropensityFunctions/*
*/*/tmpPropensityFunctions/*
*/tmpPropensityFunctions/*
*/tmpPropensityFunctions/*
Examples/Iterative*.mat
WorkSpace/Josh/Fake*.csv
Workspace/*/*
Workspace/*
Binary file removed WorkSpace/EricModel/TMPEricMH.mat
View file
Binary file not shown.
Binary file removed WorkSpace/EricModel/TMPEricMHDusp1.mat
View file
Binary file not shown.
Binary file removed WorkSpace/EricModel/TMPEricMHGR.mat
View file
Binary file not shown.
2 changes: 1 addition & 1 deletion src/+ssit/+fsp/adaptiveFspSolve.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -520,7 +520,7 @@
if propensities{i}.isFactorizable
y(i) = wt(i).*propensities{i}.stateDependentFactor(x,parameters);
else
y(i) = propensities{i}.hybridJointFactor(t,[]);
y(i) = propensities{i}.hybridJointFactor(t,x,parameters,v);
end
end
end
Expand Down
13 changes: 10 additions & 3 deletions src/+ssit/@FspMatrix/FspMatrix.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,8 +166,15 @@
w = wt(1)*[obj.terms{1}.matrix*vJ1;...
obj.terms{1}.propensity.ODEstoichVector];
for i = 2:length(obj.terms)
w = w + wt(i)*[obj.terms{i}.matrix*vJ1;...
obj.terms{i}.propensity.ODEstoichVector];
if obj.terms{i}.isFactorizable
w = w + wt(i)*[obj.terms{i}.matrix*vJ1;...
obj.terms{i}.propensity.ODEstoichVector];
else
tmp = ssit.FspMatrixTerm.generateHybridMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, ...
obj.terms{i}.numConstraints, v(end-length(upstreamODEs)+1:end)');
w = w + [tmp*vJ1;...
obj.terms{i}.propensity.ODEstoichVector];
end
end

end
Expand Down Expand Up @@ -310,7 +317,7 @@
% A = A + obj.terms{i}.propensity.hybridFactor(t,v2)*obj.terms{i}.matrix;
A = A + wt(i)*obj.terms{i}.matrix;
else
A = A + ssit.FspMatrixTerm.generateHybridgMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, obj.terms{i}.numConstraints, v2);
A = A + ssit.FspMatrixTerm.generateHybridMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, obj.terms{i}.numConstraints, v2);
end
end

Expand Down
58 changes: 32 additions & 26 deletions src/+ssit/@Propensity/Propensity.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,8 +112,8 @@
else
if ~isempty(obj.jointDependentFactor)
y = obj.jointDependentFactor(t, x, parameters);
elseif ~isempty(obj.hybridFactor)
y = obj.hybridFactor(t, x, parameters, upstreamSpecies);
elseif ~isempty(obj.hybridJointFactor)
y = obj.hybridJointFactor(t, x, parameters, upstreamSpecies);
end
end
end
Expand Down Expand Up @@ -244,7 +244,7 @@
oneSym = str2sym('1');

% change to parfor?
parfor iRxn = 1:n_reactions
for iRxn = 1:n_reactions
prop_vars = symvar(symbolicExpression{iRxn});
hybridFactor =[];
prefixNameLocal = [prefixName,'_',num2str(iRxn)];
Expand Down Expand Up @@ -332,7 +332,7 @@
% are created as much slower anonymous functions. this
% will require keeping track since the variables need
% to be sent individually to the anonymous functions.
if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logT')||isfield(logicTerms{iRxn},'logE'))
if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logT')||isfield(logicTerms{iRxn},'logJ'))
obj{iRxn}.anonymousT = true;
hybridFactor = sym2propfun(expr_t, true, false, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn));
anyLogical(iRxn) = true;
Expand All @@ -346,7 +346,7 @@
% hybridFactor = sym2mFun(expr_t, true, false, nonXTpars(:,1), speciesStoch, varODEs);
end

if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logX')||isfield(logicTerms{iRxn},'logE'))
if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logX')||isfield(logicTerms{iRxn},'logJ'))
obj{iRxn}.anonymousX = true;
anyLogical(iRxn) = true;
% elseif sum(contains(string(symvar(expr_x)),'logT'))
Expand Down Expand Up @@ -420,16 +420,21 @@
for i2=1:length(upstreamODEs)
expr_tx = subs(expr_tx,upstreamODEs{i2},varODEs(i2));
end
obj{iRxn}.(jntFactorName) = sym2propfun(expr_tx, true, true, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn));
[obj{iRxn}.(jntFactorName),expr_dt_vec_dodei] = ...
sym2propfun(expr_tx, true, true, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn), true);
obj{iRxn}.isTimeDependent = true;
expr_t_vec(iRxn) = sym(NaN);
expr_x_vec(iRxn) = sym(NaN);

% if computeSens
% for iPar = 1:n_pars
% expr_tx_vec_sens(iRxn,iPar) = signHybridFactor*diff(expr_tx,nonXTpars{iPar,1});
% end
% end

if ~isempty(expr_dt_vec_dodei)
expr_dt_vec_dode(iRxn,:) = expr_dt_vec_dodei;
end

if computeSens
for iPar = 1:n_pars
expr_tx_vec_sens(iRxn,iPar) = signHybridFactor*diff(expr_tx,nonXTpars{iPar,1});
end
end

end
end
Expand Down Expand Up @@ -560,13 +565,13 @@
n = [0,0,0];
stNew = st;
for i=1:3
J = strfind(st,logTypes{i});
J = strfind(stNew,logTypes{i});
for j = 1:length(J)
K = strfind(st,'(');
K = strfind(stNew,'(');
k1 = max(K(K<J(j)));
K = strfind(st,')');
K = strfind(stNew,')');
k2 = min(K(K>J(j)));
logE = st(k1:k2);
logE = stNew(k1:k2);
if contains(logE,'t')&&max(contains(logE,species))
n(1)=n(1)+1;
logicTerms.logJ{n(1),1} = logE;
Expand Down Expand Up @@ -704,14 +709,6 @@
end
end

for i = 1:length(nonXTpars)
exprStr = strrep(exprStr, nonXTpars{i}, ['Parameters(',num2str(i),')']);
end

for i = length(nonXTpars):-1:1
exprStr = strrep(exprStr, ['parameters',num2str(i)], ['Parameters(',num2str(i),')']);
end

for i=1:length(logicTerms)
if isfield(logicTerms{i},'logT')
for j=1:size(logicTerms{i}.logT,1)
Expand All @@ -724,15 +721,24 @@
exprStr=strrep(exprStr,logicTerms{i}.logX{j,2},logicTerms{i}.logX{j,1});
end
end
if isfield(logicTerms{i},'logE')
if isfield(logicTerms{i},'logJ')
% time_dep = true;
% state_dep = true;
for j=1:size(logicTerms{i}.logX,1)
exprStr=strrep(exprStr,logicTerms.logX{j,2},logicTerms.logX{j,1});
exprStr=strrep(exprStr,logicTerms{i}.logX{j,2},logicTerms{i}.logX{j,1});
end
end
end

for i = 1:length(nonXTpars)
exprStr = strrep(exprStr, nonXTpars{i}, ['Parameters(',num2str(i),')']);
end

for i = length(nonXTpars):-1:1
exprStr = strrep(exprStr, ['parameters',num2str(i)], ['Parameters(',num2str(i),')']);
end


if (time_dep && state_dep)
fhandle_var = ['@(t, x',parStr,')'];
for i = 1:length(varNames)
Expand Down
107 changes: 62 additions & 45 deletions src/CommandLine/SSIT.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1049,23 +1049,31 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
end

function [fimTotal,mleCovEstimate,fimMetrics] = evaluateExperiment(obj,...
fimResults,cellCounts)
fimResults,cellCounts,priorCoVariance)
% This function evaluates the provided experiment design (in
% "cellCounts" and produces an array of FIMs (one for each
% parameter set.
arguments
obj
fimResults
cellCounts
priorCoVariance = []
end

Ns = size(fimResults,2);
Nt = size(fimResults,1);
Np = size(fimResults{1,1},1);
fimTotal = cell(1,Ns);
mleCovEstimate = cell(1,Ns);

if isempty(priorCoVariance)
PriorFIM = zeros(Np);
else
PriorFIM = inv(priorCoVariance);
end

for is=1:Ns
fimTotal{is} = 0*fimResults{1,is};
fimTotal{is} = PriorFIM; %0*fimResults{1,is};

for it=1:Nt
fimTotal{is} = fimTotal{is} + cellCounts(it)*fimResults{it,is};
Expand Down Expand Up @@ -1225,17 +1233,21 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
obj.dataSet.linkedSpecies = linkedSpecies;

Q = contains(obj.dataSet.dataNames,{'time','Time','TIME'});
if sum(Q)==1
obj.dataSet.app.ParEstFitTimesList.Items = {};
if sum(Q)>=1
obj.dataSet.app.ParEstFitTimesList.Value = {};
col_time = find(Q);
obj.dataSet.app.ParEstFitTimesList.Items = {};
if sum(Q)>1
warning('Provided data more than one entry with keyword "time"')
end
col_time = find(Q,1);
obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_time_index = col_time;
obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times = sort(unique(cell2mat(obj.dataSet.DATA(:,col_time))));
for i=1:length(obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times)
obj.dataSet.app.ParEstFitTimesList.Items{i} = num2str(obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times(i));
obj.dataSet.app.ParEstFitTimesList.Value{i} = num2str(obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times(i));
end
% We need to make sure that the fitting times are included in the solution times.

else
error('Provided data set does not have required column named "time"')
end
Expand Down Expand Up @@ -1717,6 +1729,13 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
fitAlgorithm = 'fminsearch';
end

% parse fitting options
allFitOptions.suppressFSPExpansion = true;
fNames = fieldnames(fitOptions);
for i=1:length(fNames)
allFitOptions.(fNames{i}) = fitOptions.(fNames{i});
end

if isempty(obj.propensitiesGeneral)
obj = formPropensitiesGeneral(obj);
end
Expand All @@ -1735,12 +1754,11 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
if strcmp(obj.solutionScheme,'FSP') % Set solution scheme to FSP.
[FSPsoln,bounds] = obj.solve; % Solve the FSP analysis
obj.fspOptions.bounds = bounds;% Save bound for faster analyses
objFun = @(x)-obj.computeLikelihood(exp(x),FSPsoln.stateSpace); % We want to MAXIMIZE the likelihood.

if isfield(fitOptions,'suppressFSPExpansion')&&fitOptions.suppressFSPExpansion
if allFitOptions.suppressFSPExpansion
tmpFSPtol = obj.fspOptions.fspTol;
obj.fspOptions.fspTol = inf;
end
objFun = @(x)-obj.computeLikelihood(exp(x),FSPsoln.stateSpace); % We want to MAXIMIZE the likelihood.
elseif strcmp(obj.solutionScheme,'ode') % Set solution scheme to ode.
objFun = @(x)-obj.computeLikelihoodODE(exp(x)); % We want to MAXIMIZE the likelihood.
end
Expand All @@ -1749,19 +1767,13 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)

switch fitAlgorithm
case 'fminsearch'
allFitOptions.suppressFSPExpansion = true;
fNames = fieldnames(fitOptions);
for i=1:length(fNames)
allFitOptions.(fNames{i}) = fitOptions.(fNames{i});
end

[x0,likelihood,~,otherResults] = fminsearch(objFun,x0,allFitOptions);

case 'fminunc'
obj.fspOptions.fspTol = inf;
objFun = @obj.minusLogL; % We want to MAXIMIZE the likelihood.
x0 = log(parGuess);
[x0,likelihood] = fminunc(objFun,x0,fitOptions,FSPsoln.stateSpace,true);
[x0,likelihood] = fminunc(objFun,x0,allFitOptions,FSPsoln.stateSpace,true);

case 'particleSwarm'
obj.fspOptions.fspTol=inf;
Expand All @@ -1772,25 +1784,27 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)
initSwarm = repmat(x0',fitOptions.SwarmSize-1,1);
initSwarm = [x0';initSwarm.*(1+0.1*randn(size(initSwarm)))];
fitOptions.InitialSwarmMatrix = initSwarm;
[x0,likelihood] = particleswarm(OBJps,length(x0),LB,UB,fitOptions);
[x0,likelihood] = particleswarm(OBJps,length(x0),LB,UB,allFitOptions);

case 'mlSearch'
% Not yet working efficiently.
allFitOptions.maxIter=1000;
allFitOptions.burnIn=30;
allFitOptions.updateRate=10;
allFitOptions.guessRate=1000;
allFitOptions.proposalDistribution=@(x)x+0.01*randn(size(x));
allFitOptions.useFIMforSearch = false;
allFitOptions.CovFIMscale = 0.6;
allFitOptions.suppressFSPExpansion = true;
allFitOptions.logForm = true;
allFitOptions.plotFunVals = false;
allFitOptions.proposalDistributionWide=@(x)x+randn(size(x));

fNames = fieldnames(fitOptions);
defaultFitOptions.maxIter=1000;
defaultFitOptions.burnIn=30;
defaultFitOptions.updateRate=10;
defaultFitOptions.guessRate=1000;
defaultFitOptions.proposalDistribution=@(x)x+0.01*randn(size(x));
defaultFitOptions.useFIMforSearch = false;
defaultFitOptions.CovFIMscale = 0.6;
defaultFitOptions.suppressFSPExpansion = true;
defaultFitOptions.logForm = true;
defaultFitOptions.plotFunVals = false;
defaultFitOptions.proposalDistributionWide=@(x)x+randn(size(x));

fNames = fieldnames(defaultFitOptions);
for i=1:length(fNames)
allFitOptions.(fNames{i}) = fitOptions.(fNames{i});
if ~isfield(allFitOptions,fNames{i})
allFitOptions.(fNames{i}) = defaultFitOptions.(fNames{i});
end
end

if allFitOptions.logForm
Expand All @@ -1805,26 +1819,29 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)

case 'MetropolisHastings'

allFitOptions.isPropDistSymmetric=true;
allFitOptions.thin=1;
allFitOptions.numberOfSamples=1000;
allFitOptions.burnIn=100;
allFitOptions.progress=true;
allFitOptions.proposalDistribution=@(x)x+0.1*randn(size(x));
allFitOptions.numChains = 1;
allFitOptions.useFIMforMetHast = false;
allFitOptions.CovFIMscale = 0.6;
allFitOptions.suppressFSPExpansion = true;
allFitOptions.logForm = true;
defaultFitOptions.isPropDistSymmetric=true;
defaultFitOptions.thin=1;
defaultFitOptions.numberOfSamples=1000;
defaultFitOptions.burnIn=100;
defaultFitOptions.progress=true;
defaultFitOptions.proposalDistribution=@(x)x+0.1*randn(size(x));
defaultFitOptions.numChains = 1;
defaultFitOptions.useFIMforMetHast = false;
defaultFitOptions.CovFIMscale = 0.6;
defaultFitOptions.suppressFSPExpansion = true;
defaultFitOptions.logForm = true;

j=1;
while exist(['TMPmh_',num2str(j),'.mat'],'file')
j=j+1;
end
allFitOptions.saveFile = ['TMPmh_',num2str(j),'.mat'];
fNames = fieldnames(fitOptions);
defaultFitOptions.saveFile = ['TMPmh_',num2str(j),'.mat'];

fNames = fieldnames(defaultFitOptions);
for i=1:length(fNames)
allFitOptions.(fNames{i}) = fitOptions.(fNames{i});
if ~isfield(allFitOptions,fNames{i})
allFitOptions.(fNames{i}) = defaultFitOptions.(fNames{i});
end
end

if allFitOptions.logForm
Expand Down Expand Up @@ -1909,7 +1926,7 @@ function sampleDataFromFSP(obj,fspSoln,saveFile)

pars = exp(x0);

if strcmp(obj.solutionScheme,'FSP')&&isfield(fitOptions,'suppressFSPExpansion')
if strcmp(obj.solutionScheme,'FSP')&&allFitOptions.suppressFSPExpansion
obj.fspOptions.fspTol = tmpFSPtol;
end

Expand Down
Loading

0 comments on commit 6e7a990

Please sign in to comment.