By: Julia on μβ

Re-posted from: https://mbesancon.github.io/post/2018-05-23-colgen/

Integer optimization often feels weird (at least to me). Simple reformulations
of a (mixed) integer optimization problem (MIP) can make it way easier to solve.
We’re going to explore one well-known example of such integer problem in two
blog posts. This first part introduces the problem and develops a naive solution.
We’re going to see why it’s complex to solve and why this formulation does not
scale.

Among major reformulations, decomposition techniques leverage special
structures to build an easy-to-solve sub-problem and a “master problem” converging
to the exact solution to the initial problem. That’s what we’re going to see in
a second post.

Integer optimization reminder

An optimization problem takes three components: decisions variables $x$, a set of
constraints telling you if a decision is feasible or not and a cost function
$c(x)$ giving a total cost of a decision. Optimization is a domain of applied
mathematics consisting in finding the best feasible decision for a problem.
Lots of decision problems come with integrality constraints: if $x$ is the
decision, then it can only take integer values 0,1,2… or even only binary
values ${0,1}$. Think of problems involving number of units produced
for a good, yes/no decisions, etc… If a problem has lots of variables, naive
enumerations of feasible solutions becomes impossible: even problems with 50
variables can make your average laptop crash.

The cutting stock problem

The problem is not new and has been given quite some thoughts because of its
different industrial applications, it has been one of the first applications of
the column generation method we are going to use. The key elements of the problems
are: given some large rolls (metal, paper or other), we need to cut smaller
portions of given lengths to satisfy a demand for the different small lengths.
Find more details here.
A small instance might be: given rolls of size $100cm$, we want to cut at least
7 rolls of size $12cm$ and 9 rolls of size $29cm$. The objective is to minimize
the number of big rolls to satisfy this demand.

How do we formulate this mathematically?

Decisions

$Y_i$ is a binary decision indicating if we use the big roll number $i$. $X_{ij}$ is an integer
giving the number of times we cut a small roll $j$ in the big roll $i$.

Constraints

$Y$ are binary variables, $X$ are integer. Now the less trivial constraints:

• Demand satisfaction constraint: the sum over all $i$ big rolls of the cut $j$
  has to satisfy the demand for that cut:
  $$\sum_{i} X_{ij} \geq D_j $$
• Roll size constraint: if a roll is used, we cannot fit more width onto it
  than its total width:
  $$\sum_{j} X_{ij} \cdot W_j \leq L \cdot Y_i $$

A first naive implementation

Let’s first import the necessary packages: we’re using JuMP as a modeling
tool, which is an optimization-specific language embedded in Julia
(compare it to AMPL, GAMS, Pyomo, PuLP).
As I consider it a language, I’ll do a full import into my namespace with using.
We also use Cbc, an open-source solver for integer problems from the Coin-OR
suite.

using JuMP
using Cbc: CbcSolver

function cutting_stock_model(maxwidth, widths, demand, N = sum(demand))
    m = Model(solver = CbcSolver())
    Y = @variable(m, Y[1:N],Bin)
    X = @variable(m, X[i=1:N,j=1:length(widths)],Int)
    demand_satisfac = @constraint(m, [j=1:length(widths)],
        sum(X[i,j] for i in 1:N) >= demand[j]
    )
    roll_size_const = @constraint(m, [i=1:N],
        sum(X[i,j] * widths[j] for j in 1:length(widths)) <= Y[i] * maxwidth
    )
    @objective(m, Min, sum(Y[i] for i in 1:N))
    return (m, X, Y)
end

Here $N$ has to be an upper bound on the number of big rolls to use, otherwise
the problem will be infeasible (not enough big rolls to find a solution
satisfying the demand). An initial naive value for this could be the total
demand, after all one small cut per roll can be considered a worst-case solution.

Let’s see what the model looks like for different instances:

julia> cutting_stock_model(100, [12,10], [85,97], 200)
(Minimization problem with:
 * 602 linear constraints
 * 600 variables: 200 binary, 400 integer
Solver is CbcMathProg,
X[i,j], integer, ∀ i ∈ {1,2,…,199,200}, j ∈ {1,2},
Y[i] ∈ {0,1} ∀ i ∈ {1,2,…,199,200})

julia> cutting_stock_model(100, [12,10,25], [85,97,52], 300)
(Minimization problem with:
 * 1203 linear constraints
 * 1200 variables: 300 binary, 900 integer
Solver is CbcMathProg,
X[i,j], integer,∀ i ∈ {1,2,…,299,300}, j ∈ {1,2,3},
Y[i] ∈ {0,1} ∀ i ∈ {1,2,…,299,300})

julia> cutting_stock_model(100, [12,10,25,40,30,41], [85,97,52,63,77,31], 500)
(Minimization problem with:
 * 3506 linear constraints
 * 3500 variables: 500 binary, 3000 integer
Solver is CbcMathProg,
X[i,j], integer, ∀ i ∈ {1,2,…,499,500}, j ∈ {1,2,…,5,6},
Y[i] ∈ {0,1} ∀ i ∈ {1,2,…,499,500})

We see the number of variables and constraints explode as we add more possible
cut sizes, the model also becomes more and more difficult to solve. Without
going into details on the solving process, two things make the problem difficult
to solve:

1. Symmetry: if we place cuts on a roll $Y_1$ and leave another $Y_2$ unused,
   the resulting solution is concretely the same as using $Y_2$ and leaving $Y_1$
   unused.
2. Bad relaxation: integer solvers mostly work by solving a “relaxed” version
   of the problem without the integrality constraint, and then iteratively
   restricting the problem to find the best integer solution. If the relaxed
   version of the problem yields solutions far away from an integer one, the solver
   will have more work to get there.

Difficulty (1) is pretty intuitive, but we could get some insight on (2).
Let’s define our relaxed problem. We’re going to use the Clp solver, which
will solve the same problem, but without the Int restriction for $X$
nor the Bin restriction for $Y$:

function relaxed_cutting_stock(maxwidth, widths, demand, N = sum(demand))
   m = Model(solver = ClpSolver())
   Y = @variable(m, 0 <= Y[1:N] <= 1)
   X = @variable(m, X[1:N,1:length(widths)] >= 0)
   demand_satisfac = @constraint(m, [j=1:length(widths)], sum(X[i,j] for i in 1:N) >= demand[j])
   roll_size_const = @constraint(m, [i=1:N], sum(X[i,j] * widths[j] for j in 1:length(widths)) <= Y[i] * maxwidth)
   @objective(m, Min, sum(Y[i] for i in 1:N))
   return (m,Y,X)
end

Let’s see the results:

julia> res = [(i,getvalue(Y[i])) for i in 1:N if getvalue(Y[i]) ≉ 0]
33-element Array{Tuple{Int64,Float64},1}:
 (1, 1.0)
 (2, 1.0)
 (3, 1.0)
 (4, 1.0)
 (5, 1.0)
 (6, 1.0)
 (7, 1.0)
 (8, 1.0)
 (9, 1.0)
 (10, 1.0)
 (11, 1.0)
 (12, 1.0)
 (13, 1.0)
 (14, 1.0)
 (15, 1.0)
 (16, 1.0)
 (17, 1.0)
 (18, 1.0)
 (19, 1.0)
 (20, 1.0)
 (21, 1.0)
 (22, 1.0)
 (23, 1.0)
 (24, 1.0)
 (25, 1.0)
 (26, 1.0)
 (27, 1.0)
 (28, 1.0)
 (29, 1.0)
 (30, 1.0)
 (31, 1.0)
 (32, 0.9)
 (84, 1.0)

idxs = [i for (i,_ ) in res]
julia> [getvalue(X)[i,:] for i in idxs]
33-element Array{Array{Float64,1},1}:
 [0.0, 7.0, 1.2]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 0.0, 4.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [0.0, 0.0, 4.0]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [8.0, 0.0, 0.16]
 [5.8, 0.0, 1.216]
 [7.2, 0.0, 0.144]
 [8.0, 0.0, 0.16]

We notice the $Y$ variables are overall pretty saturated and almost integer,
but the $X$ variables are highly fractional: the linear cuts are divided such
that they fit perfectly the big rolls. This will make the variable hard to
get to an integer solution.

Conclusion

This was a quick intro to the cutting stock problem to get a grasp of its
structure and difficulty, the goal was not to get too technical and keep a
broad target audience.

Hope you enjoyed it, if that’s the case, I’ll see you on the next article,
we’ll implement a column generation algorithm from scratch to solve it.
If you have any question/remarks, feel free to get in touch.

Image source: https://www.flickr.com/photos/30478819@N08/38272827564

By: Julia on μβ

Re-posted from: https://mbesancon.github.io/post/2017-09-11-graph-theory-expenses-management/

With the end-of-year celebrations, we all had some expenses to manage,
some of them shared with friends, and we all have this eternal problem
of splitting them fairly.

Les bons comptes font les bons amis.
French wisdom

Applications like Tricount or
Splitwise became famous precisely by
solving this problem for you: just enter the expenses one by one, with who
owes whom and you’ll get the simplest transactions to balance the amounts at
the end.

In this post, we’ll model the expense balancing problem from a graph
perspective and see how to come up with a solution using Julia and the
JuliaGraphs ecosystem [1].

Table of Contents

The expenses model

Say that we have $n$ users involved in the expenses. An expense
$\delta$ is defined by an amount spent $\sigma$, the user who paid the
expense $p$ and a non-empty set of users who are accountable for
this expense $a$.

$\delta = (\sigma, p, a)$

The total of all expenses $\Sigma$ can be though of as: for any two users $u_i$ and $u_j$,
the total amount that $u_i$ spent for $u_j$. So the expenses are a vector of
triplets (paid by, paid for, amount).

As an example, if I went out for
pizza with Joe and paid 8GPHC for the two of us, the expense is modeled as:

$\delta = (\sigma: 8GPHC, p: Mathieu, a: [Mathieu, Joe])$.

Now considering I don’t keep track of money I owe myself, the sum of all expenses
is the vector composed of one triplet:

$\Sigma = [(Mathieu, Joe, \frac{8}{2} = 4)]$

In Julia, the expense information can be translated to a structure:

User = Int
GraphCoin = Float16
struct Expense
    payer::User
    amount::GraphCoin
    users::Set{User}
end

Reducing expenses

Now that we have a full representation of the expenses,
the purpose of balancing is to find a vector of transactions which cancels out
the expenses. A naive approach would be to use the transposed expense matrix
as a transaction matrix. If $u_i$ paid $\Sigma_{i,j}$ for $u_j$,
then $u_j$ paying back that exact amount to $u_i$ will solve the problem.
So we need in the worst case as many transactions after the trip as
$|u| \cdot (|u| – 1)$. For 5 users, that’s already 20 transactions,
how can we improve it?

Breaking strongly connected components

Suppose that I paid the pizza slice to Joe for 4GPHC, but he bought me an ice
cream for 2GPHC the day after. In the naive models, we would have two
transactions after the trip: he give me 4GPHC and I would give him 2GPHC. That
does not make any sense, he should simply pay the difference between what he
owes me and what I owe him. For any pair of users, there should only be
at most one transaction from the most in debt to the other, this result in the
worst case of $\frac{|u| \cdot (|u| – 1)}{2}$ transactions, so 10 transactions
for 5 people.

Now imagine I still paid 4GPHC for Joe, who paid 2GPHC for Marie, who paid 4GPHC
for me. In graph terminology, this is called a
strongly connected component.
The point here is that transactions will flow from one user to the next one,
and back to the first.

If there is a cycle, we can find the minimal due sum within it. In our 3-people
case, it is 2GPHC. That’s the amount which is just moving from hand to hand and
back at the origin: it can be forgotten. This yields a new net debt:
I paid 2GPHC for Joe, Marie paid 2GPHC for me. We reduced the number of
transactions and the amount due thanks to this cycle reduction.

Expenses as a flow problem

To simplify the problem, we can notice we don’t actually care about who paid
whom for what, a fair reimbursement plan only requires two conditions:

1. All people who are owed some money are given at least that amount
2. People who owe money don’t pay more than the net amount they ought to pay

We can define a directed flow network with users split in two sets of vertices,
depending on whether they owe or are owed money. We call these two sets $V_1$
and $V_2$ respectively.

• There is a directed edge from any node from $V_1$ to $V_2$.
• We define a source noted $s$ connected to all vertices in $V_1$, the edge
  from $s$ to any node of $V_1$ has a capacity equal to what they owe.
• There is an edge from any node of $V_1$ to any node of $V_2$.
• We define a sink noted $t$ to which all vertices in $V_2$ connect, with
  infinite capacity and a demand (the minimal flow that has to pass through) equal
  to what they are owed.

With this model, GraphCoins will flow from user owing money to users who are
owed money, see Wikipedia description of the flow problem.

Computing net owed amount per user

Given a vector of expenses, we should be able to build the matrix holding what
is owed in net from a user to another:

"""
    Builds the matrix of net owed GraphCoins
"""
function compute_net_owing(expenses::Vector{Expense}, nusers::Int)
    owing_matrix = zeros(GraphCoin, nusers, nusers)
    # row owes to column
    for expense in expenses
        for user in expense.users
            if user != expense.payer
                owing_matrix[user,expense.payer] += expense.amount / length(expense.users)
            end
        end
    end
    # compute net owed amount
    net_owing = zeros(GraphCoin, nusers, nusers)    
    for i in 1:nusers-1
        for j in i+1:nusers
            if owing_matrix[i,j] > owing_matrix[j,i]
                net_owing[i,j] = owing_matrix[i,j] - owing_matrix[j,i]
            elseif owing_matrix[i,j] < owing_matrix[j,i]
                net_owing[j,i] = owing_matrix[j,i] - owing_matrix[i,j]
            end
        end
    end
    return net_owing::Matrix{GraphCoin}
end

From that matrix, we should determine the net amount any user owes or is owed:

"""
    What is owed to a given user (negative if user owes money)
"""
function net_owed_user(net_owing::Matrix{GraphCoin})
    return (sum(net_owing,1)' - sum(net_owing,2))[:,1]
end

The sum function used with 1 or 2 sums a matrix over its rows, columns
respectively. This computes a difference between what a user is owed and what
they owe.

Building the graph and the corresponding flow problem

A flow problem is determined by the directed graph (nodes and directed edges),
the minimal flow for any edge, a maximal flow or capacity for any edge and a
cost of having a certain flow going through each edge.

First, we need to import LightGraphs, the core package of the JuliaGraph
ecosystem containing essential types.

import LightGraphs; const lg = LightGraphs

Note that I use explicit package import (not using), an habit I
kept from using Python and that I consider more readable than importing
the whole package into the namespace. lg has become my usual name for the
LightGraphs package.

function build_graph(net_owing::Matrix{GraphCoin})
    nusers = size(net_owing,1)
    g = lg.DiGraph(nusers + 2)
    source = nusers + 1
    sink = nusers + 2
    net_user = net_owed_user(net_owing)
    v1 = [idx for idx in 1:nusers if net_user[idx] < 0]
    v2 = [idx for idx in 1:nusers if net_user[idx] >= 0]
    capacity = zeros(GraphCoin, nusers+2,nusers+2)
    demand = zeros(GraphCoin, nusers+2,nusers+2)
    maxcap = sum(net_owing)
    for u1 in v1
        lg.add_edge!(g,source,u1)
        capacity[source,u1] = -net_user[u1]
        for u2 in v2
            lg.add_edge!(g,u1,u2)
            capacity[u1,u2] = maxcap
        end
    end
    for u2 in v2
        lg.add_edge!(g,u2,sink)
        capacity[u2,sink] = maxcap
        demand[u2,sink] = net_user[u2]
    end
    (g, capacity, demand)
end

This function builds our graph structure and all informations we need attached.

Solving the flow problem

Now that the components are set, we can solve the problem using another
component of the JuliaGraphs ecosystem specialized for flow problems:

using LightGraphsFlows: mincost_flow
using Clp: ClpSolver

We also need a Linear Programming solver to pass to the flow solver, all we
have to do is bundle the pieces together:

function solve_expense(expenses::Vector{Expense}, nusers::Int)
    (g, capacity, demand) = build_graph(compute_net_owing(expenses, nusers))
    flow = mincost_flow(g, capacity, demand, ones(nusers+2,nusers+2), ClpSolver(), nusers+1, nusers+2)
    return flow[1:end-2,1:end-2]
end

We truncate the flow matrix because we are only interested in what users
are paying each other, not in the flows from and to the source and sink.

Trying out our solution

Now that all functions are set, we can use it on any expense problem:

expenses = [
    Expense(1, 10, Set([1,2])),
    Expense(1, 24, Set([1,2,3])),
    Expense(3, 10, Set([2,3]))
]
solve_expense(expenses, 3)

3×3 Array{Float64,2}:
  0.0  0.0  0.0
 18.0  0.0  0.0
  3.0  0.0  0.0

In the result, each row pays to each column and voilà! Our three users don’t
have to feel the tension of unpaid debts anymore.

Conclusion, perspective and note on GPHC

We managed to model our specific problem using LightGraphs.jl and the
associated flow package pretty easily. I have to admit being biased since
I contributed to the JuliaGraphs ecosystem, if your impression is different
or if you have some feedback, don’t hesitate to file an issue on the
corresponding package, some awesome people
will help you figure things out as they helped me.

There is one thing we ignored in our model, it’s the number of transactions
realized. Using this as an objective turns the problem into a
Mixed-Integer Linear Programming one,
which are much harder to solve and cannot use simple flow techniques. However,
I still haven’t found a case where our simple approach does not yield the
smallest number of transactions.

Final word: I started the idea of this article long before the crypto-madness
(September actually), when currencies where still considered as boring,
nerdy or both, sorry about following the (late) hype. I even changed
GraphCoin symbol to GPHC because I found another one with which my initial
name conflicted.

If you have questions or remarks on LightGraphs, LightGraphsFlows, the article
or anything related, don’t hesitate to ping me!

Edits:
Special thanks to Seth Bromberger for the review.

The cover image was created using
GraphPlot.jl.

[1] James Fairbanks Seth Bromberger and other contributors. Juliagraphs/LightGraphs.jl:
Lightgraphs, 2017, https://doi.org/10.5281/zenodo.889971. DOI: 10.5281/zenodo.889971

By: Julia on μβ

Re-posted from: https://matbesancon.xyz/post/2017-09-11-graph-theory-expenses-management/

With the end-of-year celebrations, we all had some expenses to manage,
some of them shared with friends, and we all have this eternal problem
of splitting them fairly.

Les bons comptes font les bons amis.
French wisdom

Applications like Tricount or
Splitwise became famous precisely by
solving this problem for you: just enter the expenses one by one, with who
owes whom and you’ll get the simplest transactions to balance the amounts at
the end.

In this post, we’ll model the expense balancing problem from a graph
perspective and see how to come up with a solution using Julia and the
JuliaGraphs ecosystem [1].

The expenses model

Say that we have $n$ users involved in the expenses. An expense
$\delta$ is defined by an amount spent $\sigma$, the user who paid the
expense $p$ and a non-empty set of users who are accountable for
this expense $a$.

$\delta = (\sigma, p, a)$

The total of all expenses $\Sigma$ can be though of as: for any two users $u_i$ and $u_j$,
the total amount that $u_i$ spent for $u_j$. So the expenses are a vector of
triplets (paid by, paid for, amount).

As an example, if I went out for
pizza with Joe and paid 8GPHC for the two of us, the expense is modeled as:

$\delta = (\sigma: 8GPHC, p: Mathieu, a: [Mathieu, Joe])$.

Now considering I don’t keep track of money I owe myself, the sum of all expenses
is the vector composed of one triplet:

$\Sigma = [(Mathieu, Joe, \frac{8}{2} = 4)]$

In Julia, the expense information can be translated to a structure:

const User = Int
const GraphCoin = Float16
struct Expense
    payer::User
    amount::GraphCoin
    users::Set{User}
end

Reducing expenses

Now that we have a full representation of the expenses,
the purpose of balancing is to find a vector of transactions which cancels out
the expenses. A naive approach would be to use the transposed expense matrix
as a transaction matrix. If $u_i$ paid $\Sigma_{i,j}$ for $u_j$,
then $u_j$ paying back that exact amount to $u_i$ will solve the problem.
So we need in the worst case as many transactions after the trip as
$|u| \cdot (|u| – 1)$. For 5 users, that’s already 20 transactions,
how can we improve it?

Breaking strongly connected components

Suppose that I paid the pizza slice to Joe for 4GPHC, but he bought me an ice
cream for 2GPHC the day after. In the naive models, we would have two
transactions after the trip: he give me 4GPHC and I would give him 2GPHC. That
does not make any sense, he should simply pay the difference between what he
owes me and what I owe him. For any pair of users, there should only be
at most one transaction from the most in debt to the other, this result in the
worst case of $\frac{|u| \cdot (|u| – 1)}{2}$ transactions, so 10 transactions
for 5 people.

Now imagine I still paid 4GPHC for Joe, who paid 2GPHC for Marie, who paid 4GPHC
for me. In graph terminology, this is called a
strongly connected component.
The point here is that transactions will flow from one user to the next one,
and back to the first.

If there is a cycle, we can find the minimal due sum within it. In our 3-people
case, it is 2GPHC. That’s the amount which is just moving from hand to hand and
back at the origin: it can be forgotten. This yields a new net debt:
I paid 2GPHC for Joe, Marie paid 2GPHC for me. We reduced the number of
transactions and the amount due thanks to this cycle reduction.

Expenses as a flow problem

To simplify the problem, we can notice we don’t actually care about who paid
whom for what, a fair reimbursement plan only requires two conditions:

1. All people who are owed some money are given at least that amount
2. People who owe money don’t pay more than the net amount they ought to pay

We can define a directed flow network with users split in two sets of vertices,
depending on whether they owe or are owed money. We call these two sets $V_1$
and $V_2$ respectively.

• There is an edge from any node of $V_1$ to any node of $V_2$.
• We define a source noted $s$ connected to all vertices in $V_1$, the edge
  from $s$ to any node of $V_1$ has a capacity equal to what they owe.
• We define a sink noted $t$ to which all vertices in $V_2$ connect, with
  infinite capacity and a demand (the minimal flow that has to pass through) equal
  to what they are owed.

With this model, GraphCoins will flow from user owing money to users who are
owed money, see Wikipedia description of the flow problem.

Computing net owed amount per user

Given a vector of expenses, we should be able to build the matrix holding what
is owed in net from a user to another:

"""
Builds the matrix of net owed GraphCoins
"""
function compute_net_owing(expenses::Vector{Expense}, nusers::Int)
    owing_matrix = zeros(GraphCoin, nusers, nusers)
    # row owes to column
    for expense in expenses
        for user in expense.users
            if user != expense.payer
                owing_matrix[user,expense.payer] += expense.amount / length(expense.users)
            end
        end
    end
    # compute net owed amount
    net_owing = zeros(GraphCoin, nusers, nusers)    
    for i in 1:nusers-1
        for j in i+1:nusers
            if owing_matrix[i,j] > owing_matrix[j,i]
                net_owing[i,j] = owing_matrix[i,j] - owing_matrix[j,i]
            elseif owing_matrix[i,j] < owing_matrix[j,i]
                net_owing[j,i] = owing_matrix[j,i] - owing_matrix[i,j]
            end
        end
    end
    return net_owing
end

From that matrix, we should determine the net amount any user owes or is owed:

"""
    What is owed to a given user (negative if user owes money)
"""
function net_owed_user(net_owing::Matrix{GraphCoin})
    return (sum(net_owing,1)' - sum(net_owing,2))[:,1]
end

The sum function used with 1 or 2 sums a matrix over its rows, columns
respectively. This computes a difference between what a user is owed and what
they owe.

Building the graph and the corresponding flow problem

A flow problem is determined by the directed graph (nodes and directed edges),
the minimal flow for any edge, a maximal flow or capacity for any edge and a
cost of having a certain flow going through each edge.

First, we need to import LightGraphs, the core package of the JuliaGraph
ecosystem containing essential types.

import LightGraphs; const lg = LightGraphs

Note that I use explicit package import (not using), an habit I
kept from using Python and that I consider more readable than importing
the whole package into the namespace. lg has become my usual name for the
LightGraphs package.

function build_graph(net_owing::Matrix{GraphCoin})
    nusers = size(net_owing,1)
    g = lg.DiGraph(nusers + 2)
    source = nusers + 1
    sink = nusers + 2
    net_user = net_owed_user(net_owing)
    v1 = [idx for idx in 1:nusers if net_user[idx] < 0]
    v2 = [idx for idx in 1:nusers if net_user[idx] >= 0]
    capacity = zeros(GraphCoin, nusers+2,nusers+2)
    demand = zeros(GraphCoin, nusers+2,nusers+2)
    maxcap = sum(net_owing)
    for u1 in v1
        lg.add_edge!(g,source,u1)
        capacity[source,u1] = -net_user[u1]
        for u2 in v2
            lg.add_edge!(g,u1,u2)
            capacity[u1,u2] = maxcap
        end
    end
    for u2 in v2
        lg.add_edge!(g,u2,sink)
        capacity[u2,sink] = maxcap
        demand[u2,sink] = net_user[u2]
    end
    (g, capacity, demand)
end

This function builds our graph structure and all data we need attached.

Solving the flow problem

Now that the components are set, we can solve the problem using another
component of the JuliaGraphs ecosystem specialized for flow problems:

using LightGraphsFlows: mincost_flow
using Clp: ClpSolver

We also need a Linear Programming solver to pass to the flow solver, all we
have to do is bundle the pieces together:

function solve_expense(expenses::Vector{Expense}, nusers::Int)
    (g, capacity, demand) = build_graph(compute_net_owing(expenses, nusers))
    flow = mincost_flow(g, capacity, demand, ones(nusers+2,nusers+2), ClpSolver(), nusers+1, nusers+2)
    return flow[1:end-2,1:end-2]
end

We truncate the flow matrix because we are only interested in what users
are paying each other, not in the flows from and to the source and sink.

Trying out our solution

Now that all functions are set, we can use it on any expense problem:

expenses = [
    Expense(1, 10, Set([1,2])),
    Expense(1, 24, Set([1,2,3])),
    Expense(3, 10, Set([2,3]))
]
solve_expense(expenses, 3)

3×3 Array{Float64,2}:
  0.0  0.0  0.0
 18.0  0.0  0.0
  3.0  0.0  0.0

In the result, each row pays to each column and voilà! Our three users don’t
have to feel the tension of unpaid debts anymore.

Conclusion, perspective and note on GPHC

We managed to model our specific problem using LightGraphs.jl and the
associated flow package pretty easily. I have to admit being biased since
I contributed to the JuliaGraphs ecosystem, if your impression is different
or if you have some feedback, don’t hesitate to file an issue on the
corresponding package, some awesome people
will help you figure things out as they helped me.

There is one thing we ignored in our model, it’s the number of transactions
realized. Using this as an objective turns the problem into a
Mixed-Integer Linear Programming one,
which are much harder to solve and cannot use simple flow techniques. However,
I still haven’t found a case where our simple approach does not yield the
smallest number of transactions.

Final word: I started the idea of this article long before the crypto-madness
(September actually), when currencies where still considered as boring,
nerdy or both, sorry about following the (late) hype. I even changed
GraphCoin symbol to GPHC because I found another one with which my initial
name conflicted.

If you have questions or remarks on LightGraphs, LightGraphsFlows, the article
or anything related, don’t hesitate to ping me!

Edits:
Special thanks to Seth Bromberger for the review.

The cover image was created using
GraphPlot.jl.

[1] James Fairbanks Seth Bromberger and other contributors. Juliagraphs/LightGraphs.jl:
Lightgraphs, 2017, https://doi.org/10.5281/zenodo.889971. DOI: 10.5281/zenodo.889971