addConstraint(constraintsMatrix , bVector , pijmw_constraintMatrix , fijsd_constraintMatrix , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD)

%  addConstraint
% 
%>> Usage: function [constraintsMatrix , bVector] = addConstraint ...
%           (constraintsMatrix , bVector , pijmw_constraintMatrix ,
%               fijsd_constraintMatrix , bValue , constraintType ,
%                      NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD)
%
%>> Abstract: This function adds a new constraint to the constraint matrix
%   "A" and to the vector "B" of upper bounds on constraints. 
%
%>> Arguments:
% o  In: 
%   constraintsMatrix: Linear constraint matrix, dense or sparse m x n
%    matrix. 
%
%   bVector: Upper bounds on linear constraints, 1 x m.
%
%   pijmw_constraintMatrix: N-by-N-by-M_W matrix with the new constraint
%    cofficients correspondig to all the "lightpath wavelength link
%    indicators" (MILP decision variables 'pijmw' for each lightpath source
%    node 'i', for each lightpath end node 'j' and for each wavelength channel
%    'mw' in the network (wavelength 'w' on  physical link 'm')
%
%   fijsd_constraintMatrix: N-by-N-by-N-by-N matrix with the new constraint
%    cofficients correspondig to all the "components of traffic flow carried
%    by a given lightpath <i,j> due to a given node pair <s,d>" (MILP
%    decision variables 'f_ijsd')
%
%   bValue: Value of the bound on the new linear constraint
%
%   constraintType: Integer number
%            -1, constraint <= 
%             0, constraint == 
%            +1, constraint >=
%
%   NUMBERVARIABLES_PIJMW: Number of decision variables of the type
%   'pijmw': "lightpath wavelength link indicators" .
%
%   NUMBERVARIABLES_FIJSD: Number of decision variables of the type
%   'f_ijsd': "components of traffic flow carried by a given lightpath
%   <i,j> due to a given node pair <s,d>" .
%
% o Out:
%   constraintsMatrix: Linear constraint matrix, dense or sparse m x n
%    matrix, after adding the new constraint.
%
%   bVector: Upper bounds on linear constraints, 1 x m, after adding the
%    new constraint.
%
%
function [constraintsMatrix , bVector] = addConstraint(constraintsMatrix , bVector , pijmw_constraintMatrix , fijsd_constraintMatrix , bValue , constraintType , NUMBERVARIABLES_PIJMW , NUMBERVARIABLES_FIJSD)

if numel(pijmw_constraintMatrix)==0
    pijmw_constraint = sparse (1 ,NUMBERVARIABLES_PIJMW );
else
    % the array pijmw_constraintMatrix(NxNxMW)is converted into a row vector (1xN^2*MW) 
    pijmw_constraint  = sparse (array2vect(pijmw_constraintMatrix));    
end

if numel(fijsd_constraintMatrix)== 0
    fijsd_constraint = sparse (1 , NUMBERVARIABLES_FIJSD);
else
    % the array fijsd_constraintMatrix(NxNxNxN)is converted into a row vector (1xN^4) 
    fijsd_constraint  = sparse (array2vect(fijsd_constraintMatrix));
end

switch constraintType,  % -1 <= ; 0 == ; +1 >=
    case -1  % a <= b 
        %We add the constraints to the total matrix A and to the toal vector B
        constraintsMatrix = [constraintsMatrix; pijmw_constraint  fijsd_constraint];
        bVector = [bVector; bValue];
        
    case 0  % a = b ==> a >= b 
            %          -a <= -b
        constraintsMatrix = [constraintsMatrix; pijmw_constraint  fijsd_constraint ; -pijmw_constraint  -fijsd_constraint];
        bVector = [bVector; bValue ; -bValue];
        
    case 1  % a >= b ==> -a <= -b
        constraintsMatrix = [constraintsMatrix; -pijmw_constraint  -fijsd_constraint];
        bVector = [bVector; -bValue];        
end