function multistable = isMultistable(m)
% function multistable = isMultistable(m)
% This function implements the Feinberg criterion elucidated in refs [1-2].
% 
% It checks to see if there are no possible parameter values for which a
% biochemical network can be multistable.  This check is done by computing the
% determinant of 
% 
% Note that because this 
% 
% [1] Craciun G, Feinberg M. SIAM J Appl Math 65(5):1526-1546 (2005)
% [2] Craciun G et al. PNAS 103(23):8697-8702 (2006)

% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.

% get stoichiometry matrix
S = m.S;
% get matrix of complexes
C = -m.S.*(m.S < 0);

inds = find(any(C));
C = C(:,inds);
S = S(:,inds);

% check whether C contains at least one first-order degradation reaction for each species
hasDeg = zeros(m.nX, 1);
for i = 1:size(S, 2)
    if length(find(S(:,i))) == 1 && S(find(S(:,i)), i) == -1
        hasDeg(find(S(:,i))) = 1;
    end
end
assert(all(hasDeg), 'kronecker:isMultistable', 'Species %d does not have a first-order degradation reaction; multistability analysis does not apply.', find(hasDeg == 0, 1));

% check determinants of all n x n submatrices of the stoichiometry-complex matrix
list = nchoosek(1:size(C, 2), size(C, 1));
for k = 1:size(list,1)
    val(k) = det(C(:,list(k,:)))*det(S(:,list(k,:)));
end
multistable = any(val < 0);

if multistable
    disp(sprintf('\n\tFor the model %s, multistability is not precluded.\n', m.name));
else
    disp(sprintf('\n\tFor the model %s, there are no parameters for which the system can be multistable.\n', m.name));
end