function lp = cobra_struct2coo(model,filename)
%cobra_struct2coo 
% Changes COBRA Toolbox model structure to LP standard form 
% and exports the data into three data files used by DFBSIM.
%
% lp = cobra_struct2coo(model,filename)
%
%INPUTS
% model             COBRA Toolbox model structure
% filename          Name of data files in which the LP is saved
%
%OUTPUT
% lp                Structure containing the LP in standard form
%                   lp.A, lp,b, lp.c
%
% Kai Hoeffner 8/15/12

c = model.c ;
S = model.S ;
b = model.b ;
lb = model.lb ;
ub = model.ub ;

%% Transform into LP standard form
n = length(lb);
m = length(b);
A1=zeros(m,2*n);
A2=zeros(1,2*n);
b1=zeros(m,1);
b2=0;

c1=zeros(1,2*n);
rhslb=zeros(m,1);
rhsub=zeros(m,1);

k=0;
% There are several different possibilities
for i=1:n
    c1(i) = c(i);
    A1(:,i) = S(:,i);    
    if (lb(i)==-Inf)
        if (ub(i)==Inf) 
            % variable is unconstraint
            A1(:,n+i) = -S(:,i);
            c1(n+i) = -c(i); 
        else
            % change sign of the variable and shift upper bound to zero            
            A1(:,i) = -S(:,i);
            c1(i) = -c(i);       
            rhsub = rhsub + S(:,i)*ub(i);
            ub(i) = -(lb(i)-ub(i));
        end
    else % lb~=-Inf
        % shift finite lower bound to zero
        if (lb(i)==ub(i))
           k = k + 1;
           A2(k,i) = 1;
           b2(k) = ub(i);
        end
        if (ub(i)==Inf)   
            rhslb = rhslb + S(:,i)*lb(i);  
            lb(i) = 0;
        elseif (ub(i)>lb(i))
           rhslb = rhslb + S(:,i)*lb(i);
           ub(i) = ub(i) - lb(i);
           lb(i) = 0;
           % add upper bound constraint
           k = k + 1;
           A2(k,i) = 1;
           A2(k,n+i) = 1;
           b2(k) = ub(i);           
        end
    end
end

A = [A1;A2];
b1 = b - rhslb - rhsub;
b = [b1;b2'];
c = c1;

lp.A = A ;
lp.b = b ;
lp.c = c ;
str=fprintf('Size of A : %d %d %d \n',[size(A),nnz(A)]);

%% Write LP data to files
% A
fid = fopen(strcat(filename,'.A'),'w');
fprintf(fid,'%d , %d , %d\n',[size(A),nnz(A)]);

[i,j,val] = find(A);
coo = [i,j,val];
for i=1:nnz(A)
fprintf(fid,'%i , %i , %f\n',coo(i,1),coo(i,2),coo(i,3));
end
fclose(fid);

%b
fid = fopen(strcat(filename,'.b'),'w');
fprintf(fid,'%i\n',size(b,1));
fprintf(fid,'%f\n',b);
fclose(fid);

%c
fid = fopen(strcat(filename,'.c'),'w');
fprintf(fid,'%i\n',size(c,1));
fprintf(fid,'%f\n',c);
fclose(fid);