%% Binary Integer Programming
% This demo shows how to formulate a binary integer programming problem
% called the office assignment problem. It then demonstrates how to solve
% the problem using BINTPROG.

%   Copyright 1990-2004 The MathWorks, Inc. 
%   $Revision: 1.1.4.2 $  $Date: 2004/12/06 16:36:17 $

%% Office assignment problem
% We want to assign 6 people, Marcelo, Rakesh, Peter, Tom, Marjorie, and
% Mary Ann, to 7 offices.  Each office can have no more than one person,
% and each person gets exactly one office.  The people can give preferences
% for the offices, and their preferences are considered based on their
% seniority (the longer they have been at the MathWorks, the higher the
% seniority). Some offices have windows, some do not, and some windows are
% bigger than others. Additionally, Peter and Tom often work together, so
% we would like them to be no more than one office away, and Marcelo and
% Rakesh often work together so they also should be no more than one office
% away. 

%% Office layout
% Offices 1, 2, 3, and 4 are inside offices (no windows). Offices 5, 6 and 7
% have windows, but the window in office 5 is smaller than the other two.
% Here is how the offices are arranged.
text(0.1, .73, 'office1');
text(.35, .73, 'office2');
text(.60, .73, 'office3');
text(.82, .73, 'office4');

text(.35, .42, 'office5');
text(.60, .42, 'office6');
text(.82, .42, 'office7');
title('Office layout: window offices are in the bottom row'); 
axis off
set(gcf,'color','w');

%% Problem formulation
% We need to select our problem formulation. The first step is to choose
% what each element of our variable |x| that we are solving for represents in
% our problem. Since this is a binary integer problem, a good choice is
% that each element represents a person assigned to an office. If the
% person is assigned to the office, the variable has value 1. If the person
% is not assigned to the office, the variable has value 0. We'll consider
% the people always in this order: 

%%
%  1. Mary Ann
%  2. Marjorie
%  3. Tom
%  4. Peter
%  5. Marcelo
%  6. Rakesh

%%
% We need |x| to be a vector. Then |x(1)| to |x(7)| is Mary Ann assigned to
% office 1, office 2, etc., to office 7. The next 7 elements are Marjorie
% assigned to the seven offices, etc. In all, our |x| vector has 42 elements
% since we have six people to assign to seven offices.

%% Seniority
% We will want to weight the preferences based on seniority so that the
% longer you have been at the MathWorks, the more your preferences count.
% The seniority is as follows: Mary Ann 9 years, Marjorie 10 years, Tom 5
% years, Peter 3 years, Marcelo 1.5 years, and Rakesh 2 years. We can
% create a normalized weight vector based on seniority.

seniority = [9 10 5 3 1.5 2];
weightvector = seniority/sum(seniority);

%% Preferences
% We set up a preference matrix where the rows correspond to offices and
% the columns correspond to people. We ask each person to give values for
% each office so that the sum of all their choices, i.e., their column,
% sums to 100.  A higher number means the person prefers the office.  Tom
% is represented by column three, so Tom's desire to be in office six is
% "40", that is |prefmatrix(6,3)| is 40 because the sixth row is for the
% sixth office, and Tom is person 3 (column 3).

prefmatrix = [ 0  0  0  1  3 10 ;
               0  0  0  3  4 10 ;
               0  0  0  3  1 10 ;
               0  0  0  3  2 10 ;  
              10 20 30 10 10 20 ;
              40 40 40 40 40 20 ;
              50 40 30 40 40 20 ];

%%
% We want to weight the preferences matrix by the |weightvector| to scale the
% columns by seniority.  Also, it will be more convenient if we reshape this
% matrix as a vector in column-order so that it corresponds to our |x| vector.  
PM = prefmatrix * diag(weightvector);
c = PM(:);

%% Objective function
% Our objective is to maximize the satisfaction of the preferences weighted
% by seniority. This is a linear objective function
%%
%       max c'*x
%%
% or equivalently
%%
%       min -c'*x

%% Constraints
% The first set of constraints requires that each person gets exactly one
% office, that is for each person, the sum of the |x| values corresponding to
% that person is exactly one. For example, since Marjorie is the second
% person, this means that |sum(x(8:14))=1|. We can represent these linear
% constraints in an equality matrix Aeq and right hand side vector beq,
% where |Aeq*x = beq|, by building the appropriate matrices.  The matrix
% |Aeq| consists of ones and zeros. For example, the second row of |Aeq|
% will correspond to Marjorie getting exactly one office, so the row will 
% look like
%%
%       0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0
%%
% that is, seven ones in columns 8 through 14 and zeros elsewhere. Then
% |Aeq(2,:)*x = 1| is equivalent to |sum(x(8:14)) = 1|.
numOffices = 7;
numPeople = 6;
numDim = numOffices * numPeople;
onesvector = ones(1,numOffices);
% Each row of Aeq corresponds to one person. 
Aeq = blkdiag(onesvector,onesvector,onesvector,onesvector, ...
    onesvector,onesvector);
beq = ones(numPeople,1);
% View the structure of Aeq, that is, where there are nonzeros (ones)
figure;
spy(Aeq)
set(gcf,'color','w');
title('Equality constraints: each person gets exactly one office')

%%
% The second set of constraints are inequalities. These constraints specify
% that each office has no more than one person in it, i.e., each office 
% has one person in it, or is empty. We build the matrix |A| and the vector
% |b| such that |A*x <= b| to capture these constraints. Each row of |A|
% and |b| corresponds to an office and so row 1 corresponds to people
% assigned to office 1.  This time the rows have this type of pattern, for
% row 1:
%%
%      1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0  ... 1 0 0 0 0 0 0
%%
% Each row after this is similar but shifted (circularly) to the right by
% one spot. For example, row 3 corresponds to office 3 and says that
% |A(3,:)*x <= 1|, i.e., office 3 cannot have more than one person in it.
e = eye(1,numOffices);
A = repmat(e,numOffices,numPeople);
for i = 2:7
    A(i,:) =  circshift(A(i-1,:),[0,1]);
end
b = ones(numOffices,1);
% View the structure of Aeq, that is, where there are nonzeros (ones)
figure; 
spy(A)
set(gcf,'color','w');
title('Inequality constraints: no more than one person per office')

%%
% The next set of constraints are also inequalities, so we will add them to
% the matrix |A| and vector |b|, which already contain the inequalities
% from above. We want Tom and Peter no more than one office away, and the
% same with Marcelo and Rakesh. First we'll build the distance matrix of
% the offices based on their physical locations. This is a symmetric matrix.
D = zeros(numOffices);
%   Set up the top right half of the matrix
D(1,2:end) = [1 2 3 2 3 4];
D(2,3:end) = [1 2 1 2 3];
D(3,4:end) = [1 2 1 2];
D(4,5:end) = [3 2 1];
D(5,6:end) = [1 2];
D(6,end)   = 1;
% The lower left half is the same as the upper right
D = triu(D)' + D;

% We find the offices that are more than one away.
[officeA,officeB] = find(D > 1);
numPairs = length(officeA)

%%
% This finds |numPairs| pairs of offices. For these |numPairs|, if Tom
% occupies one office in the pair, then Peter cannot occupy the other
% office in the pair, otherwise they will be more than one office away. The
% same for Marcelo and Rakesh. This gives |2*numPairs| more inequality
% constraints which we will add to |A| and |b|.

% Add enough rows to A to accommodate these constraints.
numrows = 2*numPairs + numOffices; 
A((numOffices+1):numrows, 1:numDim) = zeros(2*numPairs,numDim);

%%
% For each pair of offices in numPairs, for the |x(i)| that corresponds to
% Tom in |officeA| and for the |x(j)| that corresponds to Peter in
% |OfficeB|, 
%%
%   x(i) + x(j) <= 1
%%
% that is, either Tom or Peter can occupy one of these offices, but they
% both cannot.
for i = 1:numPairs
    %  Tom is person 3 and Peter is person 4
    tom = 3;
    peter = 4;
    tomInOfficeA = sub2ind([numOffices numPeople],officeA(i),tom);
    peterInOfficeB = sub2ind([numOffices numPeople],officeB(i),peter);
    A(i+numOffices, [tomInOfficeA, peterInOfficeB]) = 1;
   
    % Repeat for Marcelo and Rakesh, adding more rows to A and b.
    % Marcelo is person 5 and Rakesh is person 6
    marcelo = 5;
    rakesh = 6;
    marceloInOfficeA = sub2ind([numOffices numPeople],officeA(i),marcelo);
    rakeshInOfficeB = sub2ind([numOffices numPeople],officeB(i),rakesh);
    A(i+numPairs+numOffices, [marceloInOfficeA, rakeshInOfficeB]) = 1;
end
b(numOffices+1:numOffices+2*numPairs) = ones(2*numPairs,1);
% View the structure of the newly added constraints in A, that is, 
%   where there are nonzeros (ones)
figure;
spy( A((numOffices+1):numrows,:) )
set(gcf,'color','w');
title('Inequality constraints: Tom & Peter nearby; Marcelo & Rakesh nearby')

%% Solving using BINTPROG
% Our problem formulation consists of a linear objective function
%%
%     min -c'*x
%%
% subject to the linear constraints
%%
%     Aeq*x = beq
%     A*x <= b
%%
% As this is the format BINTPROG expects, we can pass these matrices to
% BINTPROG.  

% Let BINTPROG choose the start point.
x0 = [];
f = -c;
% Show the iterative output for each node displayed in the branch and 
% bound algorithm.
options = optimset('Display','iter','NodeDisplayInterval',1);
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
% View the solution to see who got what office.
printofficeassign(x);
title('Solution for default BranchStrategy and NodeSearchStrategy');

%%
% For this problem, the satisfaction of the preferences by seniority is
% maximized to the value of |-fval|.  A positive |exitflag| tells us that
% BINTPROG converged.  The output structure gives us information about
% how many nodes were explored, how long the computation took, and how many
% cumulative iterations were used in the solving LP relaxations
% subproblems.

%% Changing the options
% To reduce the number of nodes explored, the time, or number of iterations
% taken, you can try varying a couple of options.  BINTPROG uses a branch
% and bound algorithm and we can use the options to adjust that algorithm.
% For example, the default branching strategy is |'maxinfeas'|, which says
% to choose the variable with the maximum integer infeasibility to branch
% on next, that is, the variable whose value is closest to 0.5. We can run
% our problem again with the branching strategy set to |'mininfeas'|, which
% says to choose the variable with the minimum integer infeasibility to
% branch on next, that is, the variable whose value is closest to 0 or 1
% but not equal to 0 or 1.

% Try BranchStrategy = mininfeas
options = optimset(options,'BranchStrategy','mininfeas');
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
printofficeassign(x);
title('Solution for BranchStrategy=mininfeas and default NodeSearchStrategy');

%%
% For this problem, the number of nodes and iterations are reduced with the
% alternate branching strategy, but the same solution is found as before.

%% 
% Finally, we can also choose a different node search strategy, i.e., the
% strategy the branch-and-bound algorithm uses to select the next node to
% search in the search tree. By default, the search is done using the
% |'bf'| or best node first search strategy, which chooses the node with
% lowest bound on the objective function to search next. We can change this
% to use the |'df'| or depth first search strategy. At each node in the
% search tree, if there is a child node one level down in the tree that has
% not already been explored, the algorithm chooses one such child to
% search. Otherwise, the algorithm moves to the node one level up in the
% tree and chooses a child node one level down from that node. 

% Try NodeSearchStrategy = df
options = optimset(options,'NodeSearchStrategy','df');
[x,fval,exitflag,output] = bintprog(f,A,b,Aeq,beq,x0,options);
fval
exitflag
output
printofficeassign(x);
title('Solution for BranchStrategy=mininfeas and NodeSearchStrategy=df');

%%
% For this problem, the number of nodes, iterations, and time are reduced
% with the alternate node search strategy. The same solution is found as
% before.
