% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: 4/24/2011
% Version: associated with production of the figures for Birch et al. 2012
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: Calls glpk, with default solver options, for input FBA problem.
% 
%
% Inputs:	stoichMatrix - is stoichiometric matrix S, for problem Sv = b
%           fluxBoundsLower - lower bounds for fluxes (cols of S)
%           fluxBoundsUpper - upper bounds for fluxes (cols of S) (may equal lower if
%                             setting exact value)
%           objectiveCoefficientVector - vector f to maximize value of max(fv)
%           rhsCell - {[Set(other than zero)],[Free], [Upper], [Lower], [Upper & Lower]}
%
% Outputs:  
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

function [fluxes,status,dual] = FBA_EB_glpkSolve(stoichMatrix,fluxBoundsLower,fluxBoundsUpper,objectiveCoefficientVector,rhsCell)

% add part here to check variables and output helpful errors

% create rhs vector, start assuming all constrained exactly to 0
glpkCtype = repmat('S', size(stoichMatrix,1), 1);
rhsVector = zeros(size(stoichMatrix,1),1);
% fill rhs vector acording to inputs
% Set, other than zero
if ~isempty(rhsCell{1})
    % set values, no need to change type
    rhsVector( rhsCell{1}(:,1) ) = rhsCell{1}(:,2);
end % if any Set bounds passed
% Free
if ~isempty(rhsCell{2})
    % set type for free variables, no value change needed
    glpkCtype( rhsCell{2}(:,1) ) = repmat( 'F', size(rhsCell{2},1), 1 );
end % if any Free bounds passed
% Upper
if ~isempty(rhsCell{3})
    % set values
    rhsVector( rhsCell{3}(:,1) ) = rhsCell{3}(:,2);
    % change type
    glpkCtype( rhsCell{3}(:,1) ) = repmat( 'U', size(rhsCell{3},1), 1 );
end % if any Upper bounds passed
% lower
if ~isempty(rhsCell{4})
    % set values
    rhsVector( rhsCell{4}(:,1) ) = rhsCell{4}(:,2);
    % change type
    glpkCtype( rhsCell{4}(:,1) ) = repmat( 'L', size(rhsCell{4},1), 1 );
end % if any Lower bounds passed
% Double
if ~isempty(rhsCell{5})
    % set values
    rhsVector( rhsCell{5}(:,1) ) = rhsCell{5}(:,2);
    % change type
    glpkCtype( rhsCell{5}(:,1) ) = repmat( 'D', size(rhsCell{5},1), 1 );
end % if any Double bounds passed

% these are the settings currently for FBA, including mostly default
glpkSense = -1 ; % maximization
glpkVartype = repmat('C', size(stoichMatrix,2), 1); % all variables continuous
glpkParam = struct([]);

[fluxes, ~, status, extra] = ...
                              glpk (objectiveCoefficientVector, ... % c
                                    stoichMatrix, ... % a
                                    rhsVector, ... % b
                                    fluxBoundsLower, ... % lb
                                    fluxBoundsUpper, ... % ub
                                    glpkCtype, ... % ctype
                                    glpkVartype, ... % vartype
                                    glpkSense, ... % sense
                                    glpkParam ) ; % param
                                
dual = extra.lambda;