% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('../DATABASE'));

%% CHAPTER 19: Flux Balance Analysis

%% Linear spaces
% System of linear equations
syms x1 x2 x3
x = [x1 x2 x3]'
A = [1 3 1;4 2 2];
eqs = A*x;
sol = solve(eqs)
sol.x1
sol.x2
x3 = 2
c = subs(x3)
b = subs(sol.x2)
a = subs(sol.x1)
x = double([a b c]/norm([a b c]))'
x = ([a b c]/norm([a b c]))'

A*x
A*2*x
A*(0.0031*x)
null(A)

% System of linear equations with inequalities
syms x1 x2 x3
x = [x1 x2 x3]'
A = [4 2 2];
eqs = A*x;
sol = solve(eqs)

% Assign 1st 2 values to the free variables
x3 = 0.3
x2 = 0.7
c = subs(x3)
b = subs(x2)
a = subs(sol)
x_1 = double([a b c]/norm([a b c]))'
x_1 = ([a b c]/norm([a b c]))'

% Assign 2nd 2 values to the free variables
x3 = -1.3
x2 = 2.5
c = subs(x3)
b = subs(x2)
a = subs(sol)
x_2 = double([a b c]/norm([a b c]))'
x_2 = ([a b c]/norm([a b c]))'

%% The null space of S

S = [1 -1  0  0  0  0;
     0  1 -1 -1  0  0;
     0  0  1  0 -1  0;
     0  0  0  1  0 -1]
S = sym(S)
NS = null(S)

S = [1 -1  0  0  0  0;
     0  1 -1 -1  0  0;
     0  0  1  0 -1  0;
     0  0  0  1  0 -1]

NS = null(S)
NS'*NS

NS = null(S,'r')
NS'*NS

%%
S = [1 -1  -1]

S = sym(S)
NS = null(S)

n1 = NS(:,1)
n2 = NS(:,2)
normal = cross(n1, n2)

% Next, we declare x, y, and z to be a vector of symbolic variables, and we
% multiply this vector by the normal vector to generate the left hand side
% of the equation ax+by+cz = 0.

syms x y z
P = [x,y,z]
planefunction = P*normal

% The equation of our plane is planefunction = 0, and we can solve for z in
% terms of x and y in this equation:

zplane = solve(planefunction, z)

[X,Y] = meshgrid(0:0.02:7,0:0.02:5);
Z = X-Y;
Zneg = Z<0 | Z>6;
Z(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
figure;surf(X,Y,Z,C,'FaceColor','g','EdgeColor','none','FaceAlpha',0.2);

xlabel('X = V2');ylabel('Y = V3');zlabel('Z = V4');
box on

hold on
YminZmin = find(Y==0&Z==0)
YmaxZmin = find(Y==5&Z==0)
YmaxXmax = find(Y==5&X==7)
XmaxZmax = find(X==7&Z==6)
YminZmax = find(Y==0&Z==6)

P1 = [X(YminZmin) Y(YminZmin) Z(YminZmin)] 
P2 = [X(YmaxZmin) Y(YmaxZmin) Z(YmaxZmin)]
P3 = [X(YmaxXmax) Y(YmaxXmax) Z(YmaxXmax)]
P4 = [X(XmaxZmax) Y(XmaxZmax) Z(XmaxZmax)]
P5 = [X(YminZmax) Y(YminZmax) Z(YminZmax)]

Pmat = [P1;P2;P3;P4;P5]
PX = Pmat(:,1);PY = Pmat(:,2);PZ = Pmat(:,3);
patch(PX,PY,PZ,'g','FaceAlpha',0.1);
% alpha(0.1);

t = 0:0.01:1;
plot3(2*t,2*t,0*t,'r','LineWidth',3)
plot3(2*t,0*t,2*t,'r','LineWidth',3)

NS2 = [2 3;1 2;1 1]
plot3(2*t,1*t,1*t,'g','LineWidth',3)
plot3(3*t,2*t,1*t,'g','LineWidth',3)

v = [5 2 3]'
scatter3(v(1),v(2),v(3),150,'ok','MarkerEdgeColor','blue',...
    'LineWidth',2,'MarkerFaceColor','cyan')
NS\v
NS2\v

%% Limiting planes
Xmax = 7;
Ymax = 5;
Zmax = 6;

%% Plane parallel to xz at given y

pointA = [0,4,0];
pointB = [Xmax,4,0];
pointC = [Xmax,4,Zmax];
pointD = [0,4,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Plane parallel to yz at given x

pointA = [6,0,0];
pointB = [6,Ymax,0];
pointC = [6,Ymax,Zmax];
pointD = [6,0,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Plane parallel to xy at given z

pointA = [0,0,5];
pointB = [0,Ymax,5];
pointC = [Xmax,Ymax,5];
pointD = [Xmax,0,5];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Cut polytope
[X,Y] = meshgrid(0:0.02:6,0:0.02:4);
Z = X-Y;
Zneg = Z<0 | Z>5;
Z(Zneg) = NaN;
[nx,ny] = size(X);

YminZmin = find(Y==0&Z==0)
YmaxZmin = find(Y==4&Z==0)
YmaxXmax = find(Y==4&X==6)
XmaxZmax = find(X==6&Z==5)
YminZmax = find(Y==0&Z==5)

P1 = [X(YminZmin) Y(YminZmin) Z(YminZmin)] 
P2 = [X(YmaxZmin) Y(YmaxZmin) Z(YmaxZmin)]
P3 = [X(YmaxXmax) Y(YmaxXmax) Z(YmaxXmax)]
P4 = [X(XmaxZmax) Y(XmaxZmax) Z(XmaxZmax)]
P5 = [X(YminZmax) Y(YminZmax) Z(YminZmax)]

Pmat = [P1;P2;P3;P4;P5]
PX = Pmat(:,1);PY = Pmat(:,2);PZ = Pmat(:,3);
patch(PX,PY,PZ,'r','FaceAlpha',0.1);

%% Optimization: step 1
clear, clc

% Here we set up the constraints:
S = [1 -1 -1];
der = 0;
lb = zeros(3,1); 		% Lower bound
ub = [6 4 5]; 	        % Upper bound
v0 = [4 4 0];			% Initial solution guess

%% Optimization: step 2a
% Here we set up the objective function:
c = -[1 2 1]';

% using dual-simplex:
options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v,fval,exitflag,output,lambda] = linprog(c,[],[],S,der,lb,ub,[],options);
v, fval

%% Optimization: step 2b
% Here we run again with a different objective function
c = -[1 2 3]';

% using dual-simplex:
options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v,fval,exitflag,output,lambda] = linprog(c,[],[],S,der,lb,ub,[],options);
v, fval

%% Optimization: step 2c
% Here we run again with a different objective function
c = -[1 1 1]';

% using dual-simplex:
options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v,fval,exitflag,output,lambda] = linprog(c,[],[],S,der,lb,ub,[],options);
v, fval
  
%% Reduced cost and shadow price
%
% Shadow price
shadow_price = -lambda.eqlin

% Reduced cost
% lambda.upper
% lambda.lower
reduced_costs = lambda.lower - lambda.upper 
s = c-S'*-lambda.eqlin

%% Test for shadow price
der_vec = [0:.01:0.1];
fval_vec = zeros(length(der_vec),1);
for i = 1:length(der_vec)
[v,fval_vec(i)] = linprog(c,[],[],S,der_vec(i),lb,ub,[],options);
end
shadow_price = (fval_vec(end)-fval_vec(1))/(der_vec(end)-der_vec(1))    
figure;plot(der_vec,fval_vec)

%% Test for reduced costs
lb = zeros(3,1); 		% Lower bound
ub = [6 4 5]; 	        % Upper bound
lb_vec = [0:.01:0.1] + lb(1);
fval_vec = zeros(length(lb_vec),1);
for i = 1:length(lb_vec)
    lb(1) = lb_vec(i)
    [~,fval_vec(i)] = linprog(c,[],[],S,der,lb,ub,[],options);
end

reduced_costs
reduced_costs_v1_lower = (fval_vec(end)-fval_vec(1))/(lb_vec(end)-lb_vec(1))    
figure;plot(lb_vec,fval_vec)

% but
lb = zeros(3,1); 		% Lower bound
ub = [6 4 5]; 	        % Upper bound
ub_vec = [0:.01:0.1] + ub(1);
fval_vec = zeros(length(ub_vec),1);
for i = 1:length(ub_vec)
    ub(1) = ub_vec(i)    
    [~,fval_vec(i)] = linprog(c,[],[],S,der,lb,ub,[],options);
end

reduced_costs
reduced_costs_v1_upper = (fval_vec(end)-fval_vec(1))/(ub_vec(end)-ub_vec(1))    
figure;plot(ub_vec,fval_vec)

%% LP
clear, clc

S = [1 -1  0  0  0  0;
     0  1 -1 -1  0  0;
     0  0  1  0 -1  0;
     0  0  0  1  0 -1];
 
NS = null(S,'r')

dxdt = [0 0 0 0]';
lb = [0 0 0 0 0 0]';    % Lower bound
ub = [3 3 3 3 3 3]';    % Upper bound
c = -[1 1 1 0 1 0]';

% Dual-simplex method is the default
% options = optimoptions(@linprog,'Algorithm','interior-point-legacy');
% options = optimoptions(@linprog,'Algorithm','interior-point');
options = optimoptions(@linprog,'Algorithm','dual-simplex');
% options = optimoptions(@linprog,'Algorithm','simplex');
[v_opt,fval_opt,exitflag,output,lambda] = ...
    linprog(c,[],[],S,dxdt,lb,ub,[],options);

v_opt
fval_opt
lambda_equal = lambda.eqlin;
lambda_upper = lambda.upper
lambda_lower = lambda.lower
shadow_prices = -lambda.eqlin
reduced_costs = c-S'*-lambda.eqlin
reduced_costs = lambda_lower - lambda_upper

%% LP
% clear, clc
% 
% S = [1 -1  0  0  0  0;
%      0  1 -1 -1  0  0;
%      0  0  1  0 -1  0;
%      0  0  0  1  0 -1];
%  
% NS = null(S,'r')
% 
% dxdt1 = [0 0 1 0]';
% lb = [0 0 0 0 0 0]';    % Lower bound
% ub = [3 3 3 3 3 3]';    % Upper bound
% c = -[1 1 1 0 1 0]';
% 
% % Dual-simplex method is the default
% % options = optimoptions(@linprog,'Algorithm','interior-point-legacy');
% % options = optimoptions(@linprog,'Algorithm','interior-point');
% options = optimoptions(@linprog,'Algorithm','dual-simplex');
% % options = optimoptions(@linprog,'Algorithm','simplex');
% [v_opt_1,fval_opt_1,exitflag_1,output_1,lambda_1] = ...
%     linprog(c,[],[],S,dxdt1,lb,ub,[],options);
% 
% v_opt_1
% fval_opt_1
% lambda_equal_1 = lambda_1.eqlin;
% lambda_upper_1 = lambda_1.upper
% lambda_lower_1 = lambda_1.lower
% shadow_prices_1 = -lambda_1.eqlin
% reduced_costs_1 = c-S'*-lambda_1.eqlin
% reduced_costs_1 = lambda_lower_1 - lambda_upper_1

%% Increased lower bound on v6
clear, clc

S = [1 -1  0  0  0  0;
     0  1 -1 -1  0  0;
     0  0  1  0 -1  0;
     0  0  0  1  0 -1];
 
dxdt = [0 0 0 0]';
lb = [0 0 0 0 0 0.1]';    % Lower bound
ub = [3 3 3 3 3 3]';    % Upper bound
c = -[1 1 1 0 1 0]';

% options = optimoptions(@linprog,'Algorithm','interior-point-legacy');
% options = optimoptions(@linprog,'Algorithm','interior-point');
options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v_opt,fval_opt,~,~,lambda] = linprog(c,[],[],S,dxdt,lb,ub,[],options);

v_opt
fval_opt

lambda_equal = lambda.eqlin;
lambda_upper = lambda.upper
lambda_lower = lambda.lower
shadow_prices = -lambda.eqlin
reduced_costs = c-S'*-lambda.eqlin
reduced_costs = lambda_lower - lambda_upper

%% Using the open source glpk package (also dual-simplex algorithm)
% 
ctype = 'SSSS';
vartype = 'CCCCCC';
sense = 1; 
[v_opt_glpk, fval_opt_glpk, status, extra] = glpk (c, S, dxdt, lb, ub, ctype,...
      vartype,sense);
v_opt_glpk 
fval_opt_glpk
shadow_prices_glpk = extra.lambda
reduced_costs_glpk = extra.redcosts
%   
%% The left null space of S
% 
% The rows of L represent metabolites whose concentration sum does not
% change during the reaction: these are the metabolic pools.

S = [-1 -1 1]'
L = null(S','r')
L'*S
L(:,1) = L(:,1)+L(:,2)
L, L'
a = [2 1]'

%% x vector as intersection of the two planes

x_short = pinv(L')*a

% L_large = [L' ; S']
% a_large = [a;0]
% x_short = L_large\a_large
% x_short'*S

x_long = x_short - S

% psi = [0:.1:1]';
% [psi_mesh,s_mesh] = meshgrid([0:.1:1],S')
% 
% x_param = psi_mesh.*s_mesh + double(x_long(:,ones(1,length(psi))));
% 
% plot3(x_param(1,:),x_param(2,:),x_param(3,:),'--b','lineWidth',3)
% axis equal

%% Elemental matrix
% H-Cl + Na-OH ==> H-OH + Na-Cl
S = [-1 -1 1 1]'
L = null(S','r')
rank(L)
L'*S

E = [1 0 0 1;...
     1 1 2 0;...
     0 1 1 0;...
     0 1 0 1]
rank(E)
E*S

%% Elemental matrix
% H-Cl + Na-OH ==> H-OH + Na-Cl
% Na-Cl ==> Na+ + Cl-
% H-OH ==> H+ + OH-
                 
S = [-1 -1 1 1 0 0 0 0;0 0 0 -1 1 1 0 0;0 0 -1 0 0 0 1 1]'
E = [1 0 0 1 0 1 0 0;...
     1 1 2 0 0 0 1 1;...
     0 1 1 0 0 0 0 1;...
     0 1 0 1 1 0 0 0;...
     0 0 0 0 1 -1 1 -1]
rank(E)
E*S

L = null(S','r')
L(:,1) = L(:,1)+L(:,2); 
rank(L)
L'*S

%% SVD of the stoichiometric matrix
clear, clc
S = [-1 1 ;
     -1 1 ;
     1 -1]
% format rat
% format short
[U,SIG,V] = svd(S)

%% Mixed Integer Linear Programming
%
S = [1 -1  0  0 -1  0 -1  0  0  0  0  0  0  0;
     0  1 -1  0  0  0  0 -1  1 -1  0 -1  0  0;
     0  0  1 -1  0  1  1  0  0  0  0  0  0  1;
     0  0  0  0  1 -1  0  1 -1  0  1  0  0  0;
     0  0  0  0  0  0  0  0  0  0  0  1 -1 -1];
 
[nrows,ncols] = size(S);

NS = null(sym(S));

dxdt = zeros(nrows,1);        
lb = zeros(ncols,1);                    %  Lower bound
ub = ones(ncols,1)*max(abs(S(:)));      %  Upper bound
 
c = -eye(ncols);
v_opt = zeros(ncols,ncols);
fval_opt = zeros(ncols,1);
intcon = [1:ncols];
options = optimoptions('intlinprog','Display','off');

for i = 1:ncols
% For MATLAB versions earlier the 14
% [v_opt(:,i),fval_opt(i),exitflag,output] = ...
%      bintprog(c(i,:),[],[],S,dxdt);
[v_opt(:,i),fval_opt(i),exitflag,outputa] = ...
     intlinprog(c(i,:),intcon,[],[],S,dxdt,lb,ub,options);
end
%v_opt
v_opt_unique = (unique(round(v_opt'),'rows'))'

%% Random sampling of the null space
S = [1 -1  -1]

% S = sym(S)
NS = null(S,'r')

n1 = NS(:,1)
n2 = NS(:,2)
normal = cross(n1, n2)

% Next, we declare x, y, and z to be a vector of symbolic variables, and we
% multiply this vector by the normal vector to generate the left hand side
% of the equation ax+by+cz = 0.

syms x y z
P = [x,y,z]
planefunction = P*normal

% The equation of our plane is planefunction = 0, and we can solve for z in
% terms of x and y in this equation:

zplane = solve(planefunction, z)

[X,Y] = meshgrid(0:0.02:10,0:0.02:5);
Z = X-Y;
Zneg = Z<0 | Z>6;
Z(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
figure;surf(X,Y,Z,C,'FaceColor',[0 1 0],'EdgeColor','none','FaceAlpha',0.3);

xlabel('X = V2');ylabel('Y = V3');zlabel('Z = V4');
box on
hold on

t = 0:0.01:1;
plot3(n1(1)*t,n1(2)*t,n1(3)*t,'-r','LineWidth',2)
scatter3(n1(1),n1(2),n1(3),100,'o','MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n2(1)*t,n2(2)*t,n2(3)*t,'-b','LineWidth',2)
scatter3(n2(1),n2(2),n2(3),100,'o','MarkerEdgeColor','blue',...
    'LineWidth',2,'MarkerFaceColor','yellow')
ax= gca;
ax.BoxStyle = 'full';

%% Limiting planes (use patch or fill3 with identical results)
Xmax = 7;
Ymax = 5;
Zmax = 6;

% Plane parallel to xz at given y

pointA = [0,4,0];
pointB = [Xmax,4,0];
pointC = [Xmax,4,Zmax];
pointD = [0,4,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
% patch(X,Y,Z,'b','FaceAlpha',0.1);
fill3(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to yz at given x

pointA = [6,0,0];
pointB = [6,Ymax,0];
pointC = [6,Ymax,Zmax];
pointD = [6,0,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
% patch(X,Y,Z,'b','FaceAlpha',0.1);
fill3(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to xy at given z

pointA = [0,0,5];
pointB = [0,Ymax,5];
pointC = [Xmax,Ymax,5];
pointD = [Xmax,0,5];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
% patch(X,Y,Z,'b','FaceAlpha',0.1);
fill3(X,Y,Z,'b','FaceAlpha',0.1);

%% Random sampling

n1_rand = rand(1,1000)*4;
n2_rand = rand(1,1000)*6;

N_rand_mat = n1(:,ones(1,1000)).*n1_rand(ones(3,1),:) + ...
    n2(:,ones(1,1000)).*n2_rand(ones(3,1),:);
plot3(N_rand_mat(1,:)',N_rand_mat(2,:),N_rand_mat(3,:),'or','MarkerSize',3,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','m')

x_limit = 6;
y_limit = 4;
z_limit = 5;

x_limited = N_rand_mat(1,:)>x_limit;
y_limited = N_rand_mat(2,:)>y_limit;
z_limited = N_rand_mat(3,:)>z_limit;
xyz_limited = any([x_limited ; y_limited ; z_limited]);
N_rand_mat_limited = N_rand_mat(:,~xyz_limited);
plot3(N_rand_mat_limited(1,:)',N_rand_mat_limited(2,:),...
    N_rand_mat_limited(3,:),'or','MarkerSize',5,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','c')

%% Probability histograms
n1_rand = rand(1,10e6)*4;
n2_rand = rand(1,10e6)*6;

N1 = double(n1);
N2 = double(n2);

N_rand_mat = N1(:,ones(1,10e6)).*n1_rand(ones(3,1),:) + ...
    N2(:,ones(1,10e6)).*n2_rand(ones(3,1),:);

x_limited = N_rand_mat(1,:)>x_limit;
y_limited = N_rand_mat(2,:)>y_limit;
z_limited = N_rand_mat(3,:)>z_limit;
xyz_limited = any([x_limited ; y_limited ; z_limited]);
N_rand_mat_limited = N_rand_mat(:,~xyz_limited);

figure;
[nelems,xcenters] = hist(N_rand_mat_limited(1,:),40)
nelems = nelems/10e6 ;
bar(xcenters,nelems,'BarWidth',1)
xlabel('V2 value');ylabel('V2 probability')
xlim([0 7])

figure;
[nelems,xcenters] = hist(N_rand_mat_limited(2,:),40)
nelems = nelems/10e6;
bar(xcenters,nelems,'BarWidth',1)
xlabel('V3 value');ylabel('V3 probability')
xlim([0 7])

figure;
[nelems,xcenters] = hist(N_rand_mat_limited(3,:),40)
nelems = nelems/10e6;
bar(xcenters,nelems,'BarWidth',1)
xlabel('V4 value');ylabel('V4 probability')
xlim([0 7])

% Alternatively we could determine the pdf
% pdx = fitdist(N_rand_mat_limited(1,:)','Kernel')
% pdf_kernel = pdf(pdx,[1:0.1:8]);
% figure;
% plot([1:0.1:8],pdf_kernel,'Color','b','LineWidth',1);

%% Redundant and dominant constraints
S = [1 -1  -1]

% S = sym(S)
NS = null(S,'r')

n1 = NS(:,1)
n2 = NS(:,2)
normal = cross(n1, n2)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)

[X,Y] = meshgrid(0:0.02:11,0:0.02:6);
Z = X-Y;
Zneg = Z<0 | Z>8;
Z(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
figure;surf(X,Y,Z,C,'EdgeColor','none','FaceAlpha',0.4);

xlabel('X = V2');ylabel('Y = V3');zlabel('Z = V4');
box on
hold on

t = 0:0.01:1;
plot3(n1(1)*t,n1(2)*t,n1(3)*t,'-r','LineWidth',2)
plot3(n1(1),n1(2),n1(3),'-or','MarkerSize',10,'MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n2(1)*t,n2(2)*t,n2(3)*t,'-b','LineWidth',2)
plot3(n2(1),n2(2),n2(3),'-or','MarkerSize',10,'MarkerEdgeColor','blue',...
    'LineWidth',2,'MarkerFaceColor','yellow')

%% 1st case: Limiting planes

% Initial values
% Xmax = 7;
% Ymax = 5;
% Zmax = 6;

Xmax = 7;
Ymax = 5;
Zmax = 7;

% Plane parallel to xz at given y

pointA = [0,4,0];
pointB = [Xmax,4,0];
pointC = [Xmax,4,Zmax+1];
pointD = [0,4,Zmax+1];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to yz at given x

pointA = [6,0,0];
pointB = [6,Ymax,0];
pointC = [6,Ymax,Zmax+1];
pointD = [6,0,Zmax+1];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to xy at given z

pointA = [0,0,Zmax];
pointB = [0,Ymax,Zmax];
pointC = [Xmax,Ymax,Zmax];
pointD = [Xmax,0,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Random sampling

n1_rand = rand(1,1000)*4;
n2_rand = rand(1,1000)*6;

N_rand_mat = N1(:,ones(1,1000)).*n1_rand(ones(3,1),:) + ...
    N2(:,ones(1,1000)).*n2_rand(ones(3,1),:);
plot3(N_rand_mat(1,:)',N_rand_mat(2,:),N_rand_mat(3,:),'or','MarkerSize',3,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','m')

x_limit = 6;
y_limit = 4;
z_limit = 7;

x_limited = N_rand_mat(1,:)>x_limit;
y_limited = N_rand_mat(2,:)>y_limit;
z_limited = N_rand_mat(3,:)>z_limit;
xyz_limited = any([x_limited ; y_limited ; z_limited]);
N_rand_mat_limited = N_rand_mat(:,~xyz_limited);
plot3(N_rand_mat_limited(1,:)',N_rand_mat_limited(2,:),...
    N_rand_mat_limited(3,:),'or','MarkerSize',5,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','c')

hold off

%% 2nd case: Limiting planes
S = [1 -1  -1]

% S = sym(S)
NS = null(S,'r')

n1 = NS(:,1)
n2 = NS(:,2)
normal = cross(n1, n2)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)

[X,Y] = meshgrid(0:0.02:11,0:0.02:6);
Z = X-Y;
Zneg = Z<0 | Z>8;
Z(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
figure;surf(X,Y,Z,C,'EdgeColor','none','FaceAlpha',0.4);

xlabel('X = V2');ylabel('Y = V3');zlabel('Z = V4');
box on
hold on

t = 0:0.01:1;
plot3(n1(1)*t,n1(2)*t,n1(3)*t,'-r','LineWidth',2)
plot3(n1(1),n1(2),n1(3),'-or','MarkerSize',10,'MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n2(1)*t,n2(2)*t,n2(3)*t,'-b','LineWidth',2)
plot3(n2(1),n2(2),n2(3),'-or','MarkerSize',10,'MarkerEdgeColor','blue',...
    'LineWidth',2,'MarkerFaceColor','yellow')

%%
% Initial values
% Xmax = 7;
% Ymax = 5;
% Zmax = 6;

Xmax = 7;
Ymax = 5;
Zmax = 7;

% Plane parallel to xz at given y

pointA = [0,4,0];
pointB = [Xmax,4,0];
pointC = [Xmax,4,Zmax+1];
pointD = [0,4,Zmax+1];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to yz at given x

pointA = [6,0,0];
pointB = [6,Ymax,0];
pointC = [6,Ymax,Zmax+1];
pointD = [6,0,Zmax+1];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to xy at given z

pointA = [0,0,2];
pointB = [0,Ymax,2];
pointC = [Xmax,Ymax,2];
pointD = [Xmax,0,2];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Random sampling

n1_rand = rand(1,1000)*4;
n2_rand = rand(1,1000)*6;

N_rand_mat = N1(:,ones(1,1000)).*n1_rand(ones(3,1),:) + ...
    N2(:,ones(1,1000)).*n2_rand(ones(3,1),:);
plot3(N_rand_mat(1,:)',N_rand_mat(2,:),N_rand_mat(3,:),'or','MarkerSize',3,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','m')

x_limit = 6;
y_limit = 4;
z_limit = 2;

x_limited = N_rand_mat(1,:)>x_limit;
y_limited = N_rand_mat(2,:)>y_limit;
z_limited = N_rand_mat(3,:)>z_limit;
xyz_limited = any([x_limited ; y_limited ; z_limited]);
N_rand_mat_limited = N_rand_mat(:,~xyz_limited);
plot3(N_rand_mat_limited(1,:)',N_rand_mat_limited(2,:),...
    N_rand_mat_limited(3,:),'or','MarkerSize',5,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','c')

hold off
%% 3rd case: Limiting planes

S = [1 -1  -1]

% S = sym(S)
NS = null(S,'r')

n1 = NS(:,1)
n2 = NS(:,2)
normal = cross(n1, n2)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)

[X,Y] = meshgrid(0:0.02:11,0:0.02:6);
Z = X-Y;
Zneg = Z<0 | Z>8;
Z(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
figure;surf(X,Y,Z,C,'EdgeColor','none','FaceAlpha',0.4);

xlabel('X = V2');ylabel('Y = V3');zlabel('Z = V4');
box on
hold on

t = 0:0.01:1;
plot3(n1(1)*t,n1(2)*t,n1(3)*t,'-r','LineWidth',2)
plot3(n1(1),n1(2),n1(3),'-or','MarkerSize',10,'MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n2(1)*t,n2(2)*t,n2(3)*t,'-b','LineWidth',2)
plot3(n2(1),n2(2),n2(3),'-or','MarkerSize',10,'MarkerEdgeColor','blue',...
    'LineWidth',2,'MarkerFaceColor','yellow')


Xmax = 3;
Ymax = 5;
Zmax = 6;

% Plane parallel to xz at given y

pointA = [0,4,0];
pointB = [Xmax,4,0];
pointC = [Xmax,4,Zmax+2];
pointD = [0,4,Zmax+2];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to yz at given x

pointA = [3,0,0];
pointB = [3,Ymax,0];
pointC = [3,Ymax,Zmax+2];
pointD = [3,0,Zmax+2];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

% Plane parallel to xy at given z

pointA = [0,0,5];
pointB = [0,Ymax,5];
pointC = [7,Ymax,5];
pointD = [7,0,5];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.1);

%% Random sampling

n1_rand = rand(1,1000)*4;
n2_rand = rand(1,1000)*6;

N_rand_mat = N1(:,ones(1,1000)).*n1_rand(ones(3,1),:) + ...
    N2(:,ones(1,1000)).*n2_rand(ones(3,1),:);
plot3(N_rand_mat(1,:)',N_rand_mat(2,:),N_rand_mat(3,:),'or','MarkerSize',3,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','m')

x_limit = 3;
y_limit = 4;
z_limit = 5;

x_limited = N_rand_mat(1,:)>x_limit;
y_limited = N_rand_mat(2,:)>y_limit;
z_limited = N_rand_mat(3,:)>z_limit;
xyz_limited = any([x_limited ; y_limited ; z_limited]);
N_rand_mat_limited = N_rand_mat(:,~xyz_limited);
plot3(N_rand_mat_limited(1,:)',N_rand_mat_limited(2,:),...
    N_rand_mat_limited(3,:),'or','MarkerSize',5,'MarkerEdgeColor','blue',...
    'LineWidth',1,'MarkerFaceColor','c')


%% Feasible flux solution space (polytope)

clear, clc

NS = [1 1 0.5;1 0.5 1;0.5 1 1]

n1 = NS(:,1)
n2 = NS(:,2)
n3 = NS(:,3)

% Constraints
Xmax = 7;
Ymax = 7;
Zmax = 7;

% Basis
Flux_solution_space = figure;
t = [0 10];
plot3(n1(1)*t,n1(2)*t,n1(3)*t,'-r','LineWidth',2); 
hold on
scatter3(n1(1)*Xmax,n1(2)*Ymax,n1(3)*Zmax,100,'o','MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n2(1)*t,n2(2)*t,n2(3)*t,'-b','LineWidth',2)
scatter3(n2(1)*Xmax,n2(2)*Ymax,n2(3)*Zmax,100,'o','MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')
plot3(n3(1)*t,n3(2)*t,n3(3)*t,'-g','LineWidth',2)
scatter3(n3(1)*Xmax,n3(2)*Ymax,n3(3)*Zmax,100,'o','MarkerEdgeColor','red',...
    'LineWidth',2,'MarkerFaceColor','yellow')

scatter3(Xmax,Ymax,Zmax,100,'o','MarkerEdgeColor','magenta',...
    'LineWidth',2,'MarkerFaceColor','yellow')

xlim([-3 10]);ylim([-3 10]);zlim([-3 10]);
box on
grid on
xlabel('V1');ylabel('V2');zlabel('V3');

Plot3AxisAtOrigin

%%
% Plane 1
normal = cross(n1, n2)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)

[X,Y] = meshgrid(-3:0.02:10,-3:0.02:10);
Z1 = (3*X)/2 - Y;
Zneg = Z1<-3 | Z1>10;
Z1(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
surf(X,Y,Z1,C,'FaceColor',[0 1 0],'EdgeColor','none','FaceAlpha',0.2);

% Plane 2
normal = cross(n1, n3)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)
[X,Y] = meshgrid(-3:0.02:10,-3:0.02:10);
Z2 = (3*Y)/2 - X;
Zneg = Z2<-3 | Z2>10;
Z2(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
surf(X,Y,Z2,C,'FaceColor',[1 0 0],'EdgeColor','none','FaceAlpha',0.2);

% Plane 3
normal = cross(n2, n3)

syms x y z
P = [x,y,z]
planefunction = P*normal

zplane = solve(planefunction, z)

[X,Y] = meshgrid(-3:0.02:10,-3:0.02:10);
Z3 = (2*X)/3 + (2*Y)/3;
Zneg = Z3<-3 | Z3>10;
Z3(Zneg) = NaN;
[nx,ny] = size(X);
C = ones(nx,ny);
surf(X,Y,Z3,C,'FaceColor',[0 1 1],'EdgeColor','none','FaceAlpha',0.2);

%% Constraints
Xmax = 7;
Ymax = 7;
Zmax = 7;

% Plane parallel to xz at given y

pointA = [0,Ymax,0];
pointB = [10,Ymax,0];
pointC = [10,Ymax,10];
pointD = [0,Ymax,10];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.05);

% Plane parallel to yz at given x

pointA = [Xmax,0,0];
pointB = [Xmax,10,0];
pointC = [Xmax,10,10];
pointD = [Xmax,0,10];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.05);

% Plane parallel to xy at given z

pointA = [0,0,Zmax];
pointB = [0,10,Zmax];
pointC = [10,10,Zmax];
pointD = [10,0,Zmax];

point_mat = [pointA;pointB;pointC;pointD]
X = point_mat(:,1);
Y = point_mat(:,2);
Z = point_mat(:,3);
patch(X,Y,Z,'b','FaceAlpha',0.05);
       
%% polytope open
XYZ = [0 0 0;...
       n1(1)*10 n1(2)*10 n1(3)*10;...
       n2(1)*10 n2(2)*10 n2(3)*10;...
       n3(1)*10 n3(2)*10 n3(3)*10;...
       ]

XYZ1 = XYZ([1 2 4]',:)
patch(XYZ1(:,1),XYZ1(:,2),XYZ1(:,3),'y','FaceAlpha',0.5);
XYZ2 = XYZ([1 2 3]',:)
patch(XYZ2(:,1),XYZ2(:,2),XYZ2(:,3),'y','FaceAlpha',0.5);
XYZ3 = XYZ([1 3 4]',:)
patch(XYZ3(:,1),XYZ3(:,2),XYZ3(:,3),'y','FaceAlpha',0.5);


%% polytope closed
XYZ = [0 0 0;...
       n1(1)*Xmax n1(2)*Ymax n1(3)*Zmax;...
       n2(1)*Xmax n2(2)*Ymax n2(3)*Zmax;...
       n3(1)*Xmax n3(2)*Ymax n3(3)*Zmax;...
       n1(1)*10 n1(2)*10 n1(3)*10;...
       n2(1)*10 n2(2)*10 n2(3)*10;...
       n3(1)*10 n3(2)*10 n3(3)*10;...
       
       Xmax,Ymax,Zmax]

XYZ1 = XYZ([1 2 4]',:)
patch(XYZ1(:,1),XYZ1(:,2),XYZ1(:,3),'y','FaceAlpha',0.5);
XYZ2 = XYZ([1 2 3]',:)
patch(XYZ2(:,1),XYZ2(:,2),XYZ2(:,3),'y','FaceAlpha',0.5);
XYZ3 = XYZ([1 3 4]',:)
patch(XYZ3(:,1),XYZ3(:,2),XYZ3(:,3),'y','FaceAlpha',0.5);
XYZ4 = XYZ([2 3 5]',:)
patch(XYZ4(:,1),XYZ4(:,2),XYZ4(:,3),'y','FaceAlpha',0.5);
XYZ5 = XYZ([2 4 5]',:)
patch(XYZ5(:,1),XYZ5(:,2),XYZ5(:,3),'y','FaceAlpha',0.5);
XYZ6 = XYZ([3 4 5]',:)
patch(XYZ6(:,1),XYZ6(:,2),XYZ6(:,3),'y','FaceAlpha',0.5);

%% polytope ecluded
XYZ = [0 0 0;...                            % 1
       n1(1)*Xmax n1(2)*Ymax n1(3)*Zmax;... % 2
       n2(1)*Xmax n2(2)*Ymax n2(3)*Zmax;... % 3
       n3(1)*Xmax n3(2)*Ymax n3(3)*Zmax;... % 4
       n1(1)*10 n1(2)*10 n1(3)*10;...       % 5
       n2(1)*10 n2(2)*10 n2(3)*10;...       % 6
       n3(1)*10 n3(2)*10 n3(3)*10;...       % 7
       ]

XYZ1 = XYZ([2 5 7 4]',:)
patch(XYZ1(:,1),XYZ1(:,2),XYZ1(:,3),'y','FaceAlpha',0.1);
XYZ2 = XYZ([2 5 6 3]',:)
patch(XYZ2(:,1),XYZ2(:,2),XYZ2(:,3),'y','FaceAlpha',0.1);
XYZ3 = XYZ([3 6 7 4]',:)
patch(XYZ3(:,1),XYZ3(:,2),XYZ3(:,3),'y','FaceAlpha',0.1);

%% resized box
xlim([-3 10]);ylim([-3 10]);zlim([-3 10]);
xlim([-3 11]);ylim([-3 11]);zlim([-3 11]);

ax = gca;
polytope_view = ax.View

view(39.7, 25.2)

%% FBA of E. coli CORE model (only tested up to R2016a)

% First, we need to install and initialize the COBRA toolbox, which can be
% downloaded following the link:
% https://www.mathworks.com/matlabcentral/linkexchange/links/3050. A
% convenient location for the toolbox is our TOOLBOX directory, usually
% under the name cobratoolbox:

cd ../TOOLBOXES/cobratoolbox
initCobraToolbox()
cd ../../CODE

addpath(genpath('../TOOLBOXES/cobratoolbox'));
clear, clc

% Gurobi activation (these directories are likely to be different in your
% installation)
cd /Library/gurobi702/mac64/matlab
gurobi_setup
cd /Users/dgatti/Documents/FOUNDATIONS_COMPUTATIONAL_BIOLOGY/MATHWORKS_COURSE/FOUNDATIONS_COMPUTATIONAL_BIOLOGY/CODE

%% Load the simplified E. coli model.
model = readCbModel('../TOOLBOXES/cobratoolbox/test/additionalTests/testMaps/ecoli_core_model.mat');
default_model = model;
% default_model = modelReg;

% Load the simplified E. coli map.
map = readCbMap('../TOOLBOXES/cobratoolbox/test/additionalTests/testMaps/ecoli_core_map.txt');

% If readCbMap is called with no arguments, a dialog box will prompt the
% user to select a map file. After the map has been read into MATLAB, it
% can be viewed as a MATLAB figure or a scalable vector graphic (svg). To
% view a map as a MATLAB figure the following commands are used

changeCbMapOutput('matlab');
drawCbMap(map);
colorbar('off')
ax = gca;
ax.Visible = 'off';

%% Stoichiometric matrix
S = model.S;
S_full = full(model.S);

figure
spy(model.S);

Sbin = zeros(size(model.S));
Sbin(find(model.S)) = 1;
metConnectivity = sum(Sbin,2)
rxnParticipation = sum(Sbin)'

%% Biomass calculation
S_full = full(model.S);
biomass_met_ind = find(S_full(:,13));
biomass_met_n = length(biomass_met_ind);
biomass_met_names = model.metNames(biomass_met_ind);
biomass_met_stoic = S_full(biomass_met_ind,13);
biomass_met_cell = cell(biomass_met_n,2)
biomass_met_cell(:,1) = biomass_met_names(:)
for i = 1:biomass_met_n
biomass_met_cell(i,2) = {biomass_met_stoic(i)};
end
biomass_met_cell

%%
% To set the maximum glucose uptake rate to 18.5 mmol gDW-1 hr-1
% (millimoles per gram dry cell weight per hour, the default flux units
% used in the COBRA Toolbox), enter into MATLAB:

model = changeRxnBounds(model,'EX_glc(e)',-18.5,'l');

% This changes the lower bound ('l') of the glucose exchange reaction to
% -18.5, a biologically realistic uptake rate. By convention, exchange
% reactions are written as export reactions (e.g. glc[e] <==>), so import
% of a metabolite is a negative flux. To allow unlimited oxygen uptake,
% enter:

model = changeRxnBounds(model,'EX_o2(e)',-1000,'l');

% By setting the lower bound of the oxygen uptake reaction to such a large
% number, it is practically unbounded. Next, to ensure that the biomass
% reaction is set as the objective function, enter:

model = changeObjective(model,'Biomass_Ecoli_core_N(w/GAM)-Nmet2');
% modelReg = changeObjective(modelReg,'Biomass_Ecoli_core_w_GAM');

% To perform FBA with maximization of the biomass reaction as the
% objective, enter:

FBAsolution_aer = optimizeCbModel(model,'max')
% printFluxVector(model,FBAsolution_aer.x)

for i = 1:95
    FluxVector{i,1} = model.rxnNames(i);
    FluxVector{i,2} = FBAsolution_aer.x(i);
end

% FBAsolution_aer.f then gives the value of the objective function (Z). 

%% Display result on the reaction map

% changeCbMapOutput('svg');
ecoli_met_map_aer = figure;
set(gcf,'Unit','Normalized','Position',[0.1 0.1 0.5 0.8])
changeCbMapOutput('matlab');
drawCbMap(map);
drawflux_options = struct('lb',0,'colorScale',cool,...
    'zeroFluxWidth',0.2,'zeroFluxColor',[0.9,0.9,0.9]);
drawFlux(map,model,FBAsolution_aer.x,drawflux_options)

colorbar('off')
ax = gca;
ax.Visible = 'off';
% colorbar

%% FBA solution with linprog

c = zeros(95,1); c(13) = -1
lb = model.lb;
ub = model.ub;
dxdt = zeros(72,1);

options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v_opt,fval_opt,exitflag,output,lambda] = ...
    linprog(c,[],[],S_full,dxdt,lb,ub,[],options);

v_opt
fval_opt
lambda_equal = lambda.eqlin;
lambda_upper = lambda.upper;
lambda_lower = lambda.lower;
shadow_prices = -lambda.eqlin
reduced_costs = c-S_full'*-lambda.eqlin
reduced_costs = lambda_lower - lambda_upper
 
%% FBA solution with GUROBI via MATLAB interface

c = zeros(95,1); c(13) = -1;
lb = model.lb;
ub = model.ub;
dxdt = zeros(72,1);

Aeq = sparse(S_full);
gu_model.A = Aeq;
gu_model.obj = c;
gu_model.rhs = dxdt;
gu_model.sense = repmat('=', size(Aeq,1), 1);
gu_model.lb = lb;
gu_model.ub = ub;
gu_model.modelsense = 'min';


params.outputflag = 0;
result = gurobi(gu_model, params);

v_opt_gu = result.x;
fval_opt_gu = result.objval;
shadow_prices_gu = result.pi;
% reduced_costs_gu = c-S_full'*shadow_prices_gu
reduced_costs_gu = result.rc;

%% FBA solution with GLPK

c = zeros(95,1); c(13) = -1;
ctype = repmat('S',1,72);
vartype = repmat('C',1,95);
sense = 1; 
[v_opt_glpk, fval_opt_glpk, status, extra] = glpk (c, S_full, dxdt, lb, ub, ctype,...
      vartype,sense);
v_opt_glpk 
fval_opt_glpk
shadow_prices_glpk = extra.lambda;
reduced_costs_glpk = extra.redcosts;

%% Compare linprog, gurobi, glpk, COBRA
v_opt_comp = [v_opt v_opt_gu v_opt_glpk FBAsolution_aer.x ];
shadow_prices_comp = [shadow_prices shadow_prices_gu shadow_prices_glpk FBAsolution_aer.y];
reduced_costs_comp = [reduced_costs reduced_costs_gu reduced_costs_glpk FBAsolution_aer.w ];

%%
% The same simulation is performed under anaerobic conditions. With
% the same model, enter:

model = changeRxnBounds(model,'EX_o2(e)',0,'l');

FBAsolution_anaer = optimizeCbModel(model,'max')
% printFluxVector(model,FBAsolution_anaer.x)

for i = 1:95
    FluxVector{i,1} = model.rxnNames(i);
    FluxVector{i,2} = FBAsolution_aer.x(i);
    FluxVector{i,3} = FBAsolution_anaer.x(i);
end

% The lower bound of the oxygen exchange reaction is now 0, so oxygen may
% not enter the system. When optimizeCbModel is used as before, the
% resulting growth rate is now much lower, 0.4706 hr-1. The flux
% distribution shows that oxidative phosphorylation is not used in these
% conditions, and that acetate, formate, and ethanol are produced by
% fermentation pathways.

map = readCbMap('../DATABASE/ecoli_core_map');
% changeCbMapOutput('svg');
ecoli_met_map_anaer = figure;
set(gcf,'Unit','Normalized','Position',[0.1 0.1 0.5 0.8])
changeCbMapOutput('matlab');
drawCbMap(map);
drawflux_options = struct('lb',0,'colorScale',cool,...
    'zeroFluxWidth',0.2,'zeroFluxColor',[0.9,0.9,0.9]);
drawFlux(map,model,FBAsolution_anaer.x,drawflux_options)

colorbar('off')
ax = gca;
ax.Visible = 'off';
% colorbar

%% Calculation of maximal ATP yield
% FBA can also be used to determine the maximum yields of important
% cofactors and biosynthetic precursors from glucose and other substrates.
% In this example, the maximum yields of the cofactors ATP, under aerobic
% conditions is calculated. 

model = default_model;
model = changeRxnBounds(model,'EX_o2(e)',-1000,'l');
model = changeRxnBounds(model,'EX_glc(e)',-1,'b');
model = changeObjective(model,'ATPM');
model = changeRxnBounds(model,'ATPM',0,'l');
model = changeRxnBounds(model,'ATPM',1000,'u');
FBAsolution_ATP = optimizeCbModel(model,'max')

for i = 1:95
    FluxVector{i,1} = model.rxnNames(i);
    FluxVector{i,2} = FBAsolution_ATP.x(i);
end

for i = 1:72
ShadowPriceVector{i,1} = model.metNames(i);
ShadowPriceVector{i,2} = FBAsolution_ATP.y(i);
end

%% FBA solution with MATLAB linprog

c = zeros(95,1); c(11) = -1
lb = model.lb;
ub = model.ub;
dxdt = zeros(72,1);

options = optimoptions(@linprog,'Algorithm','dual-simplex');
[v_opt_ATP,fval_opt_ATP,exitflag,output,lambda_ATP] = ...
    linprog(c,[],[],S_full,dxdt,lb,ub,[],options);

v_opt_ATP
fval_opt_ATP
lambda_equal_ATP = lambda_ATP.eqlin;
lambda_upper_ATP = lambda_ATP.upper;
lambda_lower_ATP = lambda_ATP.lower;
shadow_prices_ATP = -lambda_ATP.eqlin
reduced_costs_ATP = c-S_full'*-lambda_ATP.eqlin
reduced_costs_ATP = lambda_lower_ATP - lambda_upper_ATP

%% FBA solution with GUROBI via MATLAB interface

c = zeros(95,1); c(11) = -1;
lb = model.lb;
ub = model.ub;
dxdt = zeros(72,1);

Aeq = sparse(S_full);
gu_model.A = Aeq;
gu_model.obj = c;
gu_model.rhs = dxdt;
gu_model.sense = repmat('=', size(Aeq,1), 1);
gu_model.lb = lb;
gu_model.ub = ub;
gu_model.modelsense = 'min';

params.outputflag = 0;
result = gurobi(gu_model, params);

v_opt_ATP_gu = result.x
fval_opt_ATP_gu = result.objval;
shadow_prices_ATP_gu = result.pi;
% reduced_costs_gu = c-S_full'*shadow_prices_gu
reduced_costs_ATP_gu = result.rc;

%% FBA solution with GLPK

c = zeros(95,1); c(11) = -1
ctype = repmat('S',1,72);
vartype = repmat('C',1,95);
sense = 1; 
[v_opt_ATP_glpk, fval_opt_ATP_glpk, status, extra] = glpk (c, S_full, dxdt, lb, ub, ctype,...
      vartype,sense);
v_opt_ATP_glpk 
fval_opt_ATP_glpk
shadow_prices_ATP_glpk = extra.lambda
reduced_costs_ATP_glpk = extra.redcosts

%% Compare linprog, gurobi-MATLAB, glpk, COBRA-gurobi 

v_opt_comp = [v_opt_ATP v_opt_ATP_gu v_opt_ATP_glpk FBAsolution_ATP.x];
for i = 1:72
shadow_prices_comp{i,1} = model.mets(i) 
shadow_prices_comp{i,2} = model.metNames(i)
shadow_prices_comp{i,3} = shadow_prices_ATP(i)
shadow_prices_comp{i,4} = shadow_prices_ATP_gu(i) 
shadow_prices_comp{i,5} = shadow_prices_ATP_glpk(i) 
shadow_prices_comp{i,6} = FBAsolution_ATP.y(i);
end

for i = 1:95
reduced_costs_comp{i,1} = model.rxns(i) 
reduced_costs_comp{i,2} = model.rxnNames(i) 
reduced_costs_comp{i,3} = reduced_costs_ATP(i) 
reduced_costs_comp{i,4} = reduced_costs_ATP_gu(i) 
reduced_costs_comp{i,5} = reduced_costs_ATP_glpk(i) 
reduced_costs_comp{i,6} = FBAsolution_ATP.w(i)
end

%% P/O ratio
PO_ratio = FBAsolution_ATP.x(11)/FBAsolution_ATP.x(36)

%% Robustness analysis
% Another method that uses FBA to analyze network properties is robustness
% analysis. In this method, the flux through one reaction is varied and the
% optimal objective value is calculated as a function of this flux. This
% reveals how sensitive the objective is to a particular reaction. To determine the effect of
% varying glucose uptake on growth, first use changeRxnBounds to set the
% oxygen uptake rate (EX_o2(e)) to 17 mmol gDW-1 hr-1, a realistic uptake
% rate. Then, use a for loop to set both the upper and lower bounds of the
% glucose exchange reaction to values between 0 and -20 mmol gDW-1 hr-1,
% and use optimizeCbModel to perform FBA with each uptake rate:

model = default_model;
model = changeRxnBounds(model,'EX_o2(e)',-17,'b');
growthRates = zeros(21,1);
for i = 0:20
model = changeRxnBounds(model,'EX_glc(e)',-i,'b');
FBAsolution = optimizeCbModel(model,'max');
growthRates(i+1) = FBAsolution.f;
end
plot([0:20]',growthRates)
xlabel('Glucose uptake rate')
ylabel('Growth rate')

% The oxygen uptake rate can also be varied with the glucose uptake rate
% held constant. 

model = default_model;
model = changeRxnBounds(model,'EX_glc(e)',-10,'b');
growthRates = zeros(21,1);
for i = 0:30
model = changeRxnBounds(model,'EX_o2(e)',-i,'b');
FBAsolution = optimizeCbModel(model,'max');
growthRates(i+1) = FBAsolution.f;
end
plot([0:30]',growthRates)
xlabel('O_2 uptake rate')
ylabel('Growth rate')

%% Simulating gene knockouts
% Just as growth in different environments can be simulated with FBA, gene
% knockouts can also be simulated by changing reaction bounds. To simulate
% the knockout of any gene, its associated reaction or reactions can simply
% be constrained not to have any flux. 
% FBA may be used to predict the model phenotype with gene knockouts. For
% example, growth can be predicted for E. coli growing aerobically on
% glucose with the gene b1852 (zwf), corresponding to the reaction
% glucose-6-phosphate dehydrogenase (G6PDH2r), knocked out.

model = default_model;
changeCbMapOutput('matlab');
drawCbMap(map);
colorbar('off')
ax = gca;
ax.Visible = 'off';
S = model.S;
S_full = full(model.S);

%%
model = default_model;
[del_model,hasEffect,constrRxnNames,deletedGenes] = ...
 deleteModelGenes(model,model.genes(58))

del_model = changeRxnBounds(del_model,'EX_glc(e)',-18.5,'l');
del_model = changeRxnBounds(del_model,'EX_o2(e)',-1000,'l');
del_model = changeObjective(del_model,'Biomass_Ecoli_core_N(w/GAM)-Nmet2');

FBAsolution = optimizeCbModel(del_model,'max')
% printFluxVector(del_model,FBAsolution.x)

% Display result on the reaction map
% changeCbMapOutput('svg');
changeCbMapOutput('matlab');
drawCbMap(map);
drawflux_options = struct('lb',0,...
    'zeroFluxWidth',0.2,'zeroFluxColor',[0.9,0.9,0.9],'colorScale',spring);
drawFlux(map,del_model,FBAsolution.x,drawflux_options)

colorbar('off')
ax = gca;
ax.Visible = 'off';
colormap('spring')
colorbar

%% Genetic Engineering

%% OPTKNOCK algorithm

clear, clc, close all
model = readCbModel('../TOOLBOXES/cobratoolbox/test/additionalTests/testMaps/ecoli_core_model.mat');
default_model = model;
clear model

% Here we load a model containing information on reactions reversibility.
global CBTDIR

% load model
% model = readCbModel('../TOOLBOXES/cobratoolbox/test/models/ecoli_core_model.mat');
load([CBTDIR, filesep, 'test' filesep 'models' filesep 'ecoli_core_model.mat'], 'model');
% Here we copy the reversibility information to the default model.
default_model.rev = model.rev;
% Here we update the active model from the modified default model.
model = default_model;

% set the tolerance
tol = 1e-3;

model = changeRxnBounds(model, 'EX_o2(e)', 0, 'l'); % anaerobic
model = changeRxnBounds(model, 'EX_glc(e)', -20, 'l'); % set glucose uptake to 20

% selectedRxns (don't want to knockout biomass or exchange rxns)
selectedRxns = {model.rxns{ [1, 3:5, 7:8, 10, 12, 15:16, 18, 40:41, 44,...
    46, 48:49, 51, 53:55, 57, 59:62, 64:68, 71:77, 79:83, 85:86, 89:95]}};

%set OptKnock condition and target
minGrowthRate = 0.05;

% default OptKnock settings
options.numDelSense = 'L';
options.vMax = 1000;
options.solveOptKnock = true;

% Set up the lower limit on growth rate and ATPM flux
% biomassRxnName = 'Biomass_Ecoli_core_w_GAM';
biomassRxnName = 'Biomass_Ecoli_core_N(w/GAM)-Nmet2';
constrOpt.sense = 'GE'; % 'G' greater,'E' equal,'L' less
constrOpt.values = [minGrowthRate 8.39];
constrOpt.rxnList = {biomassRxnName, 'ATPM'};

% Previous solutions that should not be repeated again
previousSolution = [];

% Run OptKnock (change this for different simulations)
% options.targetRxn= 'EX_succ(e)';
options.targetRxn= 'EX_lac-D(e)';
options.numDel = 5;
substrateRxn = 'EX_glc(e)';

% change or set a time limit
maxTime = 3600; %time in seconds
changeCobraSolverParams('MILP', 'timeLimit', maxTime);

% define the solver package to be used
solverPkgs = {'gurobi6'};
changeCobraSolver(solverPkgs, 'MILP', 0);

% the function that actually runs optknock
[optKnockSol,bilevelMILPproblem] = OptKnock(model, selectedRxns, options, constrOpt, previousSolution, 0, 'optknocksol');

% tag for the solution
optKnockSol.substrateRxn = substrateRxn;
optKnockSol.targetRxn = options.targetRxn;
optKnockSol.numDel = options.numDel;

%get biomass reaction number
% biomassRxnID = find(strcmp(biomassRxnName, bilevelMILPproblem.model.rxns));

%check the result from OptKnock
[growthRate, minProd, maxProd] = testOptKnockSol(model, optKnockSol.targetRxn, optKnockSol.rxnList);

% check if valid_sln or unsound_sln
% assert(abs(optKnockSol.obj - maxProd) / maxProd < tol);
% assert((optKnockSol.full(biomassRxnID) - growthRate) / growthRate < tol);

%result display
fprintf('\n\nSubstrate = %s  \nTarget reaction= %s \n', optKnockSol.substrateRxn , options.targetRxn);
fprintf('Optknock solution is: %f \n\n', optKnockSol.obj );
optKnockSol.rxnList
growthRate

%% Genetic Design Local Search (GDLS) algorithm

[gdlsSolution, bilevelMILPproblem, gdlsSolutionStructs] = ...
    GDLS(model, 'EX_lac-D(e)', 'minGrowth', 0.05, 'selectedRxns', ...
    selectedRxns, 'maxKO', 5, 'nbhdsz', 3);

