function B = FindAssignment (A, papers_per_reviewer, reviewers_per_paper)
%% Find the optimal assignment between papers and reviewers that satisfies
%% the load constraints given an affinity matrix A.
%% Inputs:
%% A : A sparse affinity matrix where the rows correspond to papers and the
%% columns to reviewers
%% papers_per_reviewer : Number of reviewers to be assigned to each paper -
%% this can be a scalar
%% reviewers_per_paper : The maximum number of papers to be assigned to
%% each reviewer - this can be a scalar
%%
%% Output:
%% B : A sparse binary matrix indicating the optimal assignment between
%% papers and reviewers.

%% Number of papers
npapers = size(A,1);

%% Number of reviewers
nreviewers = size(A,2);

%% Find the number of assignments in the matrix
nedges = nnz(A);

[I J V] = find(A);

%% Generate the node edge assignment matrix associated with A - This is a
%% binary matrix where the rows correspond to the nodes and the columns to
%% the edges.

NE = spalloc(npapers+nreviewers, nedges, 2*nedges);

%% Fill out the array
NE(sub2ind(size(NE),   I, [1:nedges]')) = 1;
NE(sub2ind(size(NE), J+npapers, [1:nedges]')) = 1;

%% Fill in the adjacency constraints

d = zeros(1,npapers+nreviewers);
d(1:npapers) = reviewers_per_paper;
d(npapers+1:end) = papers_per_reviewer;

%% Here we test to see whether all of the weights in the affinity matrix
%% are integral. If so we add a random preturbation to help ensure that the
%% final solution to the optimization problem is integral.
if (all(V == round(V)))
    rho = 0.5 / (npapers * max(reviewers_per_paper));
    V = V + (rho * rand(size(V)));
end

%% Note that the negative sign is used on the weight vector to reflect the
%% fact that we are interested in maximizing the affinity while linprog is
%% set up to minimize the inner product.

x = linprog(-V, NE, d', [], [], zeros(nedges,1), ones(nedges,1));

%% N.B since the linear constraints are totally unimodular we can have
%% confidence that the final result produced here should be binary 
%% modulo rounding errors. We apply a threshold here to get the final
%% selection vector t. 

t = x > 0.5;

%% Final assignment matrix we get this result by using the logical vector t
%% to select the edges in the original matrix that we will finally assign
B = spalloc(npapers,nreviewers, nnz(t));

B(sub2ind(size(B), I(t), J(t))) = 1;

%% You can check the final result by summing the rows of the assignment
%% matrix in search of papers with incomplete assignments - All of
%% the reviewers should have an acceptable number of assignments - we can
%% check this by summing the columns. N.B. we could assign different
%% reviewers different capacities without an issue - they will all stay
%% integral constraints.

%% You can use spy(B) to visualize the results of the assignment procedure