addLoopLawConstraints

PURPOSE ^

addLoopLawConstraints adds loop law constraints to LP problem or MILP problem.

SYNOPSIS ^

function [MILPproblem] = addLoopLawConstraints(LPproblem, model, rxnIndex)

DESCRIPTION ^

addLoopLawConstraints adds loop law constraints to LP problem or MILP problem.
INPUT
 LPproblem Structure containing the following fields
  A      LHS matrix
  b      RHS vector
  c      Objective coeff vector
  lb     Lower bound vector
  ub     Upper bound vector
  osense Objective sense (-1 max, +1 min)
  csense Constraint senses, a string containting the constraint sense for
         each row in A ('E', equality, 'G' greater than, 'L' less than).
  F (optional)  If *QP problem
  vartype (optional) if MI*P problem
 model     The model for which the loops should be removed

OPTIONAL INPUT
 rxnIndex The index of variables in LPproblem corresponding to fluxes. 
     default = [1:n]


OUTPUT
 Problem structure containing the following fields describing an MILP problem
 A, b, c, lb, ub - same as before but longer
 vartype - variable type of the MILP problem ('C', and 'B')
 x0 = [] Needed for solveMILPproblem

 Jan Schellenberger Sep 27, 2009

 different ways of doing it.  I'm still playing with this.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [MILPproblem] = addLoopLawConstraints(LPproblem, model, rxnIndex)
0002 %addLoopLawConstraints adds loop law constraints to LP problem or MILP problem.
0003 %INPUT
0004 % LPproblem Structure containing the following fields
0005 %  A      LHS matrix
0006 %  b      RHS vector
0007 %  c      Objective coeff vector
0008 %  lb     Lower bound vector
0009 %  ub     Upper bound vector
0010 %  osense Objective sense (-1 max, +1 min)
0011 %  csense Constraint senses, a string containting the constraint sense for
0012 %         each row in A ('E', equality, 'G' greater than, 'L' less than).
0013 %  F (optional)  If *QP problem
0014 %  vartype (optional) if MI*P problem
0015 % model     The model for which the loops should be removed
0016 %
0017 %OPTIONAL INPUT
0018 % rxnIndex The index of variables in LPproblem corresponding to fluxes.
0019 %     default = [1:n]
0020 %
0021 %
0022 %OUTPUT
0023 % Problem structure containing the following fields describing an MILP problem
0024 % A, b, c, lb, ub - same as before but longer
0025 % vartype - variable type of the MILP problem ('C', and 'B')
0026 % x0 = [] Needed for solveMILPproblem
0027 %
0028 % Jan Schellenberger Sep 27, 2009
0029 %
0030 % different ways of doing it.  I'm still playing with this.
0031 method = 2; % methd = 1 - separete af,ar;  method = 2 - only af;  method 3 - same as method 2 except use b_L, b_U instad of b and csense;
0032 reduce_vars = 1; % eliminates additional integer variables.  Should be faster in all cases but in practice may not be for some weird reason.
0033 combine_vars = 0; % combines flux coupled reactions into one variable.  Should be faster in all cases but in practice may not be.
0034 
0035 if nargin < 3
0036    if size(LPproblem.A,2) == size(model.S,2); % if the number of variables matches the number of model reactions
0037        rxnIndex = 1:size(model.S,2);
0038    elseif size(LPproblem.A,2) > size(model.S,2)
0039        display('warning:  extra variables in LPproblem.  will assume first n correspond to v')
0040        rxnIndex = 1:size(model.S,2);
0041    else
0042        display('LPproblem must have at least as many variables as model has reactions');
0043        return;
0044    end
0045 elseif length(find(rxnIndex)) ~= size(model.S,2)
0046     display('rxnIndex must contain exactly n entries');
0047     return;
0048 end
0049 if any(rxnIndex > size(LPproblem.A,2))
0050     display('rxnIndex out of bounds');
0051     return;
0052 end
0053 
0054 MILPproblem = LPproblem;
0055 
0056 S = model.S;
0057 [m,n] = size(LPproblem.A);
0058 nontransport = (sum(S ~= 0) > 1)'; %reactions which are not transport reactions.
0059 %nnz(nontransport)
0060 nontransport = (nontransport | (model.lb ==0 & model.ub == 0));
0061 %nnz(nontransport)
0062 %pause;
0063 if reduce_vars == 1
0064     active1 = ~(model.lb ==0 & model.ub == 0);
0065     S2 = S(:,active1); % exclude rxns with ub/lb ==0
0066     
0067     N2 = sparseNull(sparse(S2));
0068     N = zeros(length(active1), size(N2,2));
0069     N(active1,:) = N2;
0070     %size(N)
0071     active = any(abs(N) > 1e-6, 2); % exclude rxns not in null space
0072     %size(active)
0073     %size(nontransport)
0074     nontransport = nontransport & active;
0075 end
0076 
0077 Sn = S(:,nontransport);
0078 
0079 Ninternal = sparseNull(sparse(Sn));
0080 %max(max(abs(Ninternal)))
0081 %pause
0082 linternal = size(Ninternal,2);
0083 
0084 nint = length(find(nontransport));
0085 temp = sparse(nint, n);
0086 temp(:, rxnIndex(nontransport)) = speye(nint);
0087 
0088 
0089 if method == 1 % two variables (ar, af)
0090     MILPproblem.A = [LPproblem.A, sparse(m,3*nint);   % Ax = b (from original LPproblem)
0091         temp, -10000*speye(nint), sparse(nint, 2*nint); % v < 10000*af
0092         temp, sparse(nint, nint), 10000*speye(nint), sparse(nint, nint); % v > -10000ar
0093         sparse(nint, n), speye(nint), speye(nint), sparse(nint, nint);  % ar + af <= 1
0094         sparse(nint, n), -100*speye(nint), 1*speye(nint), speye(nint);  % E < 100 af - ar
0095         sparse(nint, n), -1*speye(nint), 100*speye(nint), speye(nint);  % E > af - 100 ar
0096         sparse(linternal, n+2*nint), Ninternal']; % N*E = 0
0097 
0098     MILPproblem.b = [LPproblem.b;
0099         zeros(2*nint,1);
0100         ones(nint,1);
0101         zeros(2*nint + linternal,1);];
0102 
0103     MILPproblem.c = [LPproblem.c;
0104         zeros(3*nint,1)];
0105 
0106     MILPproblem.csense = LPproblem.csense;
0107     for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end   % v < 1000*af
0108     for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end  % v > -1000ar
0109     for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end  % ar + af < 1
0110     for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end  % E <
0111     for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end  % E >
0112     for i = 1:linternal, MILPproblem.csense(end+1,1) = 'E';end % N*E = 0
0113 
0114     MILPproblem.vartype = [];
0115     if isfield(LPproblem, 'vartype')
0116         MILPproblem.vartype = LPproblem.vartype;  % keep variables same as previously.
0117     else
0118         for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end; %otherwise define as continuous (used for all LP problems)
0119     end
0120     for i = 1:2*nint, MILPproblem.vartype(end+1,1) = 'B';end;
0121     for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end;
0122 
0123     if isfield(LPproblem, 'F') % used in QP problems
0124         MILPproblem.F = sparse(size(MILPproblem.A,2),   size(MILPproblem.A,2));
0125         MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0126     end
0127 
0128 
0129     MILPproblem.lb = [LPproblem.lb; 
0130         zeros(nint*2,1);
0131         -1000*ones(nint,1);];
0132     MILPproblem.ub = [LPproblem.ub; 
0133         ones(nint*2,1);
0134         1000*ones(nint,1);];
0135 
0136     MILPproblem.x0 = [];
0137     
0138 elseif method == 2 % One variables (a)
0139     MILPproblem.A = [LPproblem.A, sparse(m,2*nint);   % Ax = b (from original LPproblem)
0140         temp, -10000*speye(nint), sparse(nint, nint); % v < 10000*af
0141         temp, -10000*speye(nint), sparse(nint, nint); % v > -10000 + 10000*af
0142         sparse(nint, n), -101*speye(nint), speye(nint);  % E < 100 af - ar
0143         sparse(nint, n), -101*speye(nint), speye(nint);  % E > af - 100 ar
0144         sparse(linternal, n + nint), Ninternal']; % N*E = 0
0145 
0146     MILPproblem.b = [LPproblem.b; % Ax = b (from original problem)
0147         zeros(nint,1); % v < 10000*af
0148         -10000*ones(nint, 1); % v > -10000 + 10000*af
0149         -ones(nint,1); % e<
0150         -100*ones(nint, 1); % e>
0151         zeros(linternal,1);];
0152 
0153     MILPproblem.c = [LPproblem.c;
0154         zeros(2*nint,1)];
0155 
0156     MILPproblem.csense = LPproblem.csense;
0157     for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end   % v < 1000*af
0158     for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end  % v > -1000ar
0159     for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end  % E <
0160     for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end  % E >
0161     for i = 1:linternal, MILPproblem.csense(end+1,1) = 'E';end % N*E = 0
0162 
0163     MILPproblem.vartype = '';
0164     if isfield(LPproblem, 'vartype')
0165         MILPproblem.vartype = LPproblem.vartype;  % keep variables same as previously.
0166     else
0167         for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end; %otherwise define as continuous (used for all LP problems)
0168     end
0169     for i = 1:nint, MILPproblem.vartype(end+1,1) = 'B';end; % a variables
0170     for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end; % G variables
0171 
0172     if isfield(LPproblem, 'F') % used in QP problems
0173         MILPproblem.F = sparse(size(MILPproblem.A,2),   size(MILPproblem.A,2));
0174         MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0175     end
0176 
0177 
0178     MILPproblem.lb = [LPproblem.lb; 
0179         zeros(nint,1);
0180         -1000*ones(nint,1);];
0181     MILPproblem.ub = [LPproblem.ub; 
0182         ones(nint,1);
0183         1000*ones(nint,1);];
0184 
0185     MILPproblem.x0 = [];
0186 elseif method == 3 % like method 3 except reduced constraints.
0187         MILPproblem.A = [LPproblem.A, sparse(m,2*nint);   % Ax = b (from original LPproblem)
0188         temp, -10000*speye(nint), sparse(nint, nint); % -10000 < v -10000*af < 0
0189         %temp, -10000*speye(nint), sparse(nint, nint); % v > -10000 + 10000*af
0190         sparse(nint, n), -101*speye(nint), speye(nint);  %  -100 < E - 101 af < -1
0191         %sparse(nint, n), -101*speye(nint), speye(nint);  % E > af - 100 ar
0192         sparse(linternal, n + nint), Ninternal']; % N*E = 0
0193 
0194     MILPproblem.b_L = [LPproblem.b; % Ax = b (from original problem)
0195         %zeros(nint,1); % v < 10000*af
0196         -10000*ones(nint, 1); % v > -10000 + 10000*af
0197         %-ones(nint,1); % e<
0198         -100*ones(nint, 1); % e>
0199         zeros(linternal,1);];
0200     MILPproblem.b_U = [LPproblem.b; % Ax = b (from original problem)
0201         zeros(nint,1); % v < 10000*af
0202         %-10000*ones(nint, 1); % v > -10000 + 10000*af
0203         -ones(nint,1); % e<
0204         %-100*ones(nint, 1); % e>
0205         zeros(linternal,1);];
0206 
0207     MILPproblem.b_L(find(LPproblem.csense == 'E')) = LPproblem.b(LPproblem.csense == 'E');
0208     MILPproblem.b_U(find(LPproblem.csense == 'E')) = LPproblem.b(LPproblem.csense == 'E');
0209     MILPproblem.b_L(find(LPproblem.csense == 'G')) = LPproblem.b(LPproblem.csense == 'G');
0210     MILPproblem.b_U(find(LPproblem.csense == 'G')) = inf;
0211     MILPproblem.b_L(find(LPproblem.csense == 'L')) = -inf;
0212     MILPproblem.b_U(find(LPproblem.csense == 'L')) = LPproblem.b(LPproblem.csense == 'L');
0213     
0214     MILPproblem.c = [LPproblem.c;
0215         zeros(2*nint,1)];
0216         
0217     MILPproblem.csense = [];
0218     
0219     MILPproblem.vartype = [];
0220     if isfield(LPproblem, 'vartype')
0221         MILPproblem.vartype = LPproblem.vartype;  % keep variables same as previously.
0222     else
0223         for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end; %otherwise define as continuous (used for all LP problems)
0224     end
0225     for i = 1:nint, MILPproblem.vartype(end+1,1) = 'B';end; % a variables
0226     for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end; % G variables
0227 
0228     if isfield(LPproblem, 'F') % used in QP problems
0229         MILPproblem.F = sparse(size(MILPproblem.A,2),   size(MILPproblem.A,2));
0230         MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0231     end
0232 
0233 
0234     MILPproblem.lb = [LPproblem.lb; 
0235         zeros(nint,1);
0236         -1000*ones(nint,1);];
0237     MILPproblem.ub = [LPproblem.ub; 
0238         ones(nint,1);
0239         1000*ones(nint,1);];
0240 
0241     MILPproblem.x0 = [];
0242 else
0243     display('method not found')
0244     method
0245     pause;
0246 end
0247 
0248 if combine_vars && method == 2
0249 %    MILPproblem
0250     %pause;
0251     Ns = N(nontransport,:);
0252     %full(Ns)
0253     %pause;
0254     %Ns = sparseNull(S(:,nontransport));
0255     %size(Ns)
0256     Ns2 = Ns;
0257     for i = 1:size(Ns2,1)
0258         m = sqrt(Ns2(i,:)*Ns2(i,:)');
0259         Ns2(i,:) = Ns2(i,:)/m;
0260     end
0261     %min(m)
0262     t = Ns2 * Ns2';
0263 %     size(t)
0264      %spy(t> .99995 | t < -.99995);
0265     %full(t)
0266      %pause;
0267      %t = corrcoef([Ns, sparse(size(Ns,1),1)]');
0268      %full(t)
0269 %     size(t)
0270      %spy(t> .99995 | t < -.99995);
0271      %pause;
0272     cutoff = .9999999;
0273     %[m1, m2] = find(t>.99 & t < .999);
0274     %for i = 1:length(m1)
0275 %         t(m1(i), m2(i))
0276 %         [m1(i), m2(i)]
0277 %         [Ns(m1(i),:); Ns(m2(i),:)]
0278 %         corr(Ns(m1(i),:)', Ns(m2(i),:)')
0279 %         pause;
0280     %end
0281     %pause;
0282     t2 = sparse(size(t,1), size(t, 2));
0283     t2(abs(t) > cutoff) = t(abs(t) > cutoff);
0284     t = t2;
0285     checkedbefore = zeros(nint,1);    
0286 
0287     for i = 2:nint
0288         x = find(t(i,1:i-1)>cutoff);
0289         if ~isempty(x)
0290             checkedbefore(i) = x(1);
0291         end
0292         y = find(t(i,1:i-1)<-cutoff);
0293         if ~isempty(y)
0294             checkedbefore(i) = -y(1);
0295         end
0296         if ~isempty(x) && ~isempty(y);
0297             if x(1) < y(1)
0298                 checkedbefore(i) = x(1);
0299             else
0300                 checkedbefore(i) = -y(1);
0301             end
0302         end
0303     end
0304     %sum(checkedbefore ~= 0)
0305     %pause;
0306     %[find(nontransport)', (1:length(checkedbefore))', checkedbefore]
0307     %nint
0308     %checkedbefore
0309     %checkedbefore(55)
0310     %    t(55,29)
0311 
0312     %pause;
0313     %checkedbefore(56:end) = 0;
0314     offset = size(LPproblem.A, 2);
0315     for i = 1:length(checkedbefore)
0316         if checkedbefore(i) ==0
0317             continue;
0318         end
0319         pretarget = abs(checkedbefore(i)); % variable that this one points to.
0320  %       [pretarget,i]
0321         if checkedbefore(i) > 0
0322             if any(MILPproblem.A(:,offset+pretarget).*MILPproblem.A(:,offset+i))
0323                 display('trouble combining vars'),pause;
0324             end
0325             MILPproblem.A(:,offset+pretarget) = MILPproblem.A(:,offset+pretarget) + MILPproblem.A(:,offset+i);
0326         else
0327             MILPproblem.A(:,offset+pretarget) = MILPproblem.A(:,offset+pretarget) - MILPproblem.A(:,offset+i);
0328             MILPproblem.b = MILPproblem.b - MILPproblem.A(:,offset+i);
0329         end
0330     end
0331     %markedfordeath = offset + find(checkedbefore > .5);
0332     markedforlife = true(size(MILPproblem.A,2), 1);
0333     markedforlife(offset + find(checkedbefore > .5)) = false;
0334 %    size(markedforlife)
0335     MILPproblem.markedforlife = markedforlife;
0336     MILPproblem.A = MILPproblem.A(:,markedforlife);
0337     MILPproblem.c = MILPproblem.c(markedforlife);
0338     MILPproblem.vartype = MILPproblem.vartype(markedforlife);
0339     MILPproblem.lb = MILPproblem.lb(markedforlife);
0340     MILPproblem.ub = MILPproblem.ub(markedforlife);
0341 %    MILPproblem.nontransport = full(double(nontransport))';
0342 %    MILPproblem.energies = zeros(size(MILPproblem.A,2), 1);
0343 %    MILPproblem.energies((end-nint+1):end) = 1;
0344 %    MILPproblem.checkedbefore = checkedbefore;
0345 %     MILPproblem.as = zeros(size(MILPproblem.A,2), 1);
0346 %     MILPproblem.as((offset+1):(offset+nint)) = 1;
0347     %pause;
0348 end

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003