Skip to content

Commit

Permalink
Simplify and remove redundancy for FSP matrix term formulation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Munsky
Munsky committed Sep 27, 2024
1 parent 3521dd5 commit 50e94d0
Show file tree
Hide file tree
Showing 4 changed files with 36 additions and 100 deletions.
2 changes: 1 addition & 1 deletion Examples/example_EscapeTimes.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,9 @@
Model1.fspOptions.escapeSinks.f = {'rna'};
Model1.fspOptions.verbose = false;
Model1.fspOptions.escapeSinks.b = 50;
Model1 = Model1.formPropensitiesGeneral('Model1');
[fspSoln1,Model1.fspOptions.bounds] = Model1.solve;
Model1.makePlot(fspSoln1,'escapeTimes',[],[],10)
Model1 = Model1.formPropensitiesGeneral('Model1');

%% Example 2 - escape time with time varying transcription rate
% First let's copy and adjust the previous mdoel to add a time varying
Expand Down
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 @@ -147,7 +147,7 @@
outputTimeCount = length(outputTimes);
outputTimes = unique(outputTimes);
if (outputTimeCount ~= length(outputTimes))
disp('Warning: input tspan contains repeated elements. Converted tspan to vector of unique values');
disp('Warning: input tspan contains repeated elements. Truncating tspan to vector of unique values');
outputTimeCount = length(outputTimes);
end

Expand Down
129 changes: 33 additions & 96 deletions src/+ssit/@FspMatrixTerm/FspMatrixTerm.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,15 @@
methods (Static)

function A_fsp = GenerateMatrixTermTimeInvariant(propensity, parameters, state_set, numConstraints, varNames, computeSens, ipar)
% This function generates the individual components of the inf.
% gen. matrix corresponding to each propensity function and
% stoichiometry vector, and also for the state space defined in
% 'state_set'. It also builds the corresponding sink states
% according the current projection space.

% TODO - it seems that this version does not allow for model
% reductions. Not sure if this code is actually reachable, but
% if so it should be updated to allow for reductions.
arguments
propensity
parameters
Expand All @@ -182,70 +191,26 @@
n_states = size(state_set.states, 2);
reachableIndices = state_set.reachableIndices(:, propensity.reactionIndex);

% if propensity.isTimeDependent&&~isempty(propensity.timeDependentFactor)
% % Constants are sometimes lumped into a constant time dependent
% % term for later convenience.
% prop_val = propensity.timeDependentFactor(0).*...
% propensity.evaluateStateFactor(state_set.states,varNames)...
% +0*zeros(1,size(state_set.states,2));
% else
if computeSens
prop_val = propensity.evaluateStateFactor(state_set.states,parameters,varNames,true,ipar)...
+0*zeros(1,size(state_set.states,2));
else
prop_val = propensity.evaluateStateFactor(state_set.states,parameters,varNames)...
+0*zeros(1,size(state_set.states,2));
end
% end

if min(prop_val)<0
warning('PropFun:negProp','Negative propensity function detected. Results may be incorrect.');
warning('OFF','PropFun:negProp')
% prop_val(prop_val<0) = 0;
end

if (size(prop_val, 2) > 1)
prop_val = prop_val';
end

% Filling in values for MATLAB's sparse matrix format (which is coo)
i = zeros(n_states*2, 1);
j = i;
aij = zeros(length(i), 1);

% Coordinates and values for diagonal entries
j(1:n_states) = 1:n_states;
i(1:n_states) = 1:n_states;
aij(1:n_states) = -1.0*prop_val;

% Coordinates and values for normal states
offset = n_states;
j((1:n_states) + offset) = 1:n_states;
i((1:n_states) + offset) = reachableIndices(1:n_states);
aij((1:n_states) + offset) = prop_val;
A_fsp = ssit.FspMatrixTerm.buildAFSP(n_states,prop_val,reachableIndices,propensity,...
state_set,numConstraints,[]);

% Coordinates and values for sink states
ireaction = propensity.reactionIndex;
n_constrs_failed = sum(state_set.outboundTransitions(:, 1 + numConstraints*(ireaction-1):numConstraints*ireaction), 2);
isinks = zeros(nnz(state_set.outboundTransitions(:, 1 + numConstraints*(ireaction-1):numConstraints*ireaction)), 1);
jsinks = isinks;
aijsinks = zeros(length(jsinks), 1);
k = 1;
for c = 1:state_set.numConstraints
ninsert = nnz(state_set.outboundTransitions(:, c + numConstraints*(ireaction-1)));
isinks(k:k+ninsert-1) = n_states + c;
jsinks(k:k+ninsert-1) = find(state_set.outboundTransitions(:,c + numConstraints*(ireaction-1)));
aijsinks(k:k+ninsert-1) = prop_val(jsinks(k:k+ninsert-1))./n_constrs_failed(jsinks(k:k+ninsert-1));
k = k + ninsert;
end

% eliminate out-of-bound coordinates
indices_keep = find( (i > 0) & ( j > 0) );
i = i(indices_keep);
j = j(indices_keep);
aij = aij(indices_keep);

A_fsp = sparse([i;isinks], [j;jsinks], [aij;aijsinks], n_states + numConstraints, n_states + numConstraints);
end

function A_fsp = generateTimeVaryingMatrixTerm(t, propensity, state_set, parameters, numConstraints, modRedTransformMatrices)
Expand All @@ -270,55 +235,8 @@
prop_val = prop_val*ones(size(state_set.states,2),1);
end

% Filling in values for MATLAB's sparse matrix format (which is coo)
i = zeros(n_states*2, 1);
j = i;
aij = zeros(length(i), 1);

% Coordinates and values for diagonal entries
j(1:n_states) = 1:n_states;
i(1:n_states) = 1:n_states;
aij(1:n_states) = -1.0*prop_val;

offset = n_states;
j((1:n_states) + offset) = 1:n_states;
i((1:n_states) + offset) = reachableIndices(1:n_states);
aij((1:n_states) + offset) = prop_val(:);

% Coordinates and values for sink states
ireaction = propensity.reactionIndex;
n_constrs_failed = sum(state_set.outboundTransitions(:, 1 + numConstraints*(ireaction-1):numConstraints*ireaction), 2);
isinks = zeros(nnz(state_set.outboundTransitions(:, 1 + state_set.numConstraints*(ireaction-1):state_set.numConstraints*ireaction-1)), 1);
jsinks = isinks;
aijsinks = zeros(length(jsinks), 1);
k = 1;
for c = 1:state_set.numConstraints
ninsert = nnz(state_set.outboundTransitions(:, c + state_set.numConstraints*(ireaction-1)));
isinks(k:k+ninsert-1) = n_states + c;
jsinks(k:k+ninsert-1) = find(state_set.outboundTransitions(:,c + state_set.numConstraints*(ireaction-1)));
aijsinks(k:k+ninsert-1) = prop_val(jsinks(k:k+ninsert-1))./n_constrs_failed(jsinks(k:k+ninsert-1));

k = k + ninsert;
end

% eliminate out-of-bound coordinates
indices_keep = find( (i > 0) & ( j > 0) );
i = i(indices_keep);
j = j(indices_keep);
aij = aij(indices_keep);

A_fsp = sparse([i;isinks], [j;jsinks], [aij;aijsinks], n_states + numConstraints, n_states + numConstraints);
if ~isempty(modRedTransformMatrices)
A_fsp = modRedTransformMatrices.phi_inv*...
A_fsp(1:end-numConstraints,1:end-numConstraints)*...
modRedTransformMatrices.phi;
if ~isempty(modRedTransformMatrices.phiScale)
A_fsp = A_fsp - diag(diag(A_fsp));
A_fsp = A_fsp.*modRedTransformMatrices.phiScale;
A_fsp = A_fsp - diag(sum(A_fsp));
end
end

A_fsp = ssit.FspMatrixTerm.buildAFSP(n_states,prop_val,reachableIndices,propensity,...
state_set,numConstraints,modRedTransformMatrices);

end

Expand All @@ -345,6 +263,23 @@
prop_val = prop_val*ones(size(state_set.states,2),1);
end

A_fsp = ssit.FspMatrixTerm.buildAFSP(n_states,prop_val,reachableIndices,propensity,...
state_set,numConstraints,modRedTransformMatrices);

end

function A_fsp = buildAFSP(n_states,prop_val,reachableIndices,propensity,...
state_set,numConstraints,modRedTransformMatrices)
arguments
n_states
prop_val
reachableIndices
propensity
state_set
numConstraints
modRedTransformMatrices = []
end

% Filling in values for MATLAB's sparse matrix format (which is coo)
i = zeros(n_states*2, 1);
j = i;
Expand All @@ -361,7 +296,7 @@
aij((1:n_states) + offset) = prop_val(:);

% Coordinates and values for sink states
ireaction = hybridPropensity.reactionIndex;
ireaction = propensity.reactionIndex;
n_constrs_failed = sum(state_set.outboundTransitions(:, 1 + numConstraints*(ireaction-1):numConstraints*ireaction), 2);
isinks = zeros(nnz(state_set.outboundTransitions(:, 1 + state_set.numConstraints*(ireaction-1):state_set.numConstraints*ireaction-1)), 1);
jsinks = isinks;
Expand Down Expand Up @@ -396,5 +331,7 @@

end



end
end
3 changes: 1 addition & 2 deletions src/CommandLine/SSIT.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -657,7 +657,6 @@ function exportToSBML(obj,modelName)
obj.initialProbs = max(0,real(fspVector.p.data.vals));

end


function [pdo] = generatePDO(obj,pdoOptions,paramsPDO,fspSoln,variablePDO,maxSize)
arguments
Expand Down Expand Up @@ -865,7 +864,7 @@ function summarizeModel(obj)
switch obj.solutionScheme
case 'FSP'
if ~isempty(stateSpace)&&size(stateSpace.states,2)~=stateSpace.state2indMap.Count
error('HERE')
error('Mismatch in statespace definition.')
end

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

0 comments on commit 50e94d0

Please sign in to comment.