clear all; close all; clc;
changeCobraSolver('glpk','LP');
model = readCbModel('Ec_iJR904_GlcMM.xml');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added reactions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model = addReaction(model,'Acyl-ACP_Reductase',{'ACP[c]','nadph[c]',...
    'Aldehyde','nadp[c]'},[-1,-1,1,1],false,-1000,1000);
model = addReaction(model,'Aldehyde_Decarbonylase',{'Aldehyde',...
    'nadph[c]', 'Alkane','nadp[c]','for[c]'},[-1,-1,1,1,1],false,0,1000);
model = addReaction(model,'Alkane_Exchange',{'Alkane','Alkane[e]'},...
    [-1,1],false,0,1000);

%%%%%%%%%%%%%%%%%%%%%%%%%% Glucose exchange rate %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% http://www.pnas.org/content/101/8/2235/T3.expansion.html %%%%%%%%%
model = changeRxnBounds(model, 'EX_glc(e)', -11.9, 'l');

%%%%%%%%%%%%%%%%%%%%%%%%% Thiamine exchange rate %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Assumption: exchange rate is proportional to that of glucose %%%%%%
model = changeRxnBounds(model, 'THMabc', -0.0004,'l');

% model = changeRxnBounds(model, 'GLYCt', -4, 'l');
% model = changeRxnBounds(model, 'SULabc', -2,'l'); 

% model = changeObjective(model,{'ACPS1','Alkane_Exchange'});
model = changeObjective(model,{'BiomassEcoli'});

%%%%%%%%%%%%%%%%%%%%%% Removing all existing bounds %%%%%%%%%%%%%%%%%%%%%%%
% model = changeRxnBounds(model, 'EX_glc(e)', -1000, 'l');
% model = changeRxnBounds(model, 'EX_o2(e)', -1000, 'l');
% model = changeRxnBounds(model, 'ATPM', -1000, 'l');
% model = changeRxnBounds(model, 'ATPM', 1000, 'u');

%%%%%%%%%%%%%%% Making the synthesis rate of A-ACP nonzero %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Manual limiting of other A-ACP reactions %%%%%%%%%%%%%%%%%
% model = changeRxnBounds(model, 'ACPS1', 10, 'l');
% model = changeRxnBounds(model, 'AACPS1', 10, 'u');
% model = changeRxnBounds(model, 'AACPS2', 10, 'u');
% model = changeRxnBounds(model, 'AACPS3', 10, 'u');
% model = changeRxnBounds(model, 'AACPS4', 10, 'u');
% model = changeRxnBounds(model, 'AACPS5', 10, 'u');

%%%%%%%%%%%%%%%% Attempt to step through all permutations %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% Memory and time restrictions %%%%%%%%%%%%%%%%%%%%%%
% possibleCombinations = nchoosek(model.rxns,10);
% fluxesOfInterest = [];
% for i = 1:length(possibleCombinations)
%     randBound = -rand(1,10)*100;
%     for j = 1:10
%         model = changeRxnBounds(model, possibleCombinations(i,j),...
%     randBound(j), 'l');
%         fluxesOfInterest(i,3) = randBound(j);
%     end
%     solution = optimizeCbModel(model,'max');
%     fluxesOfInterest(i,1) = solution.x(end-1); 
%     fluxesOfInterest(i,2) = solution.x(end);
% end

%%%%%%%%%%%%%% Attempt to generate generxnList for optGene.m %%%%%%%%%%%%%%
% numPos = ceil(rand(1)*100);
% for i = 1:numPos
%     posList = length(model.rxns)+1;
%     while(posList > length(model.rxns))
%         posList = ceil(rand(1,numPos)*1000);
%     end
%     rxnList(i,1) = model.rxns(posList(i));
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% detectDeadEnds.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
changeCobraSolver('glpk','MILP');
temp = detectDeadEnds(model);
deadEndMets = model.mets(temp);
% gapFind(model)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% optGene.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Requires fuctioning gurobi linear programming solver %%%%%%%%%%
% [x,poluation,scores,optGeneSol]=optGene(model,...
%    'Aldehyde_Decarbonylase','BiomassEcoli', model.rxns(1:end-3),5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GDLS.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs]=GDLS(model,...
    {'ACPS1','BiomassEcoli'});
model2=removeRxns(model,gdlsSolution.KOs);

%%%%%%%%%%%%%%%%%%% Optimize model with biomass function %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Print solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solution = optimizeCbModel(model,'max'); % Original model
printFluxVector(model, solution.x,true);
solution2 = optimizeCbModel(model2,'max'); % Model with reaction knock outs
printFluxVector(model2, solution2.x,true);
