Author Archives: Julia on μβ

Building our own graph type in Julia

By: Julia on μβ

Re-posted from: https://matbesancon.xyz/post/2018-08-17-abstract_graph/


This is an adapted post on the talk we gave with James
at JuliaCon 2018 in London. You can see the
original slides,
the video still requires a bit of post-processing.

Last week JuliaCon in London was a great and very condensed experience.
The two talks on LightGraphs.jl
received a lot of positive feedback and more than that, we saw
how people are using the library for a variety of use cases which is a great
signal for the work on the JuliaGraphs ecosystem
(see the lightning talk).

I wanted to re-build the same graph for people who prefer a post version to
my clumsy live explanations on a laptop not handling dual-screen well
(those who prefer the latter are invited to see the live-stream of the talk).

Why abstractions?

The LightGraphs library is built to contain as few elements as possible to get
anyone going with graphs. This includes:

  • The interface a graph type has to comply with to be used
  • Essential algorithms implemented by any graph respecting that interface
  • A simple, battery-included implementation based on adjacency lists

The thing is, if you design an abstraction which in fact has just one
implementation, you’re doing abstraction wrong. This talks was also a
reality-check for LightGraphs, are we as composable, extensible as we promised?

The reason for abstraction is also that minimalism has its price.
The package was designed as the least amount of complexity required to get
graphs working. When people started to use it, obviously they needed more
features, some of which they could code themselves, some other required
extensions built within LightGraphs. By getting the core abstractions right,
you guarantee people will be able to use it and to build on top with minimal
friction, while keeping it simple to read and contribute to.

Our matrix graph type

Let’s recall that a graph is a collection of nodes and a collection of
edges between these nodes. To keep it simple, for a graph of $n$ edges,
we will consider they are numbered from 1 to n. An edge connects a node $i$
to a node $j$, therefore all the information of a graph can be kept as an
adjacency matrix $M_{ij}$ of size $n \times n$:

$$M_{ij} = \begin{cases} 1, & \mbox{if edge (i $\rightarrow$ j) exists} \\ 0 & \mbox{otherwise}\end{cases}$$

We don’t know what the use cases for our type will be, and therefore,
we will parametrize the graph type over the matrix type:

import LightGraphs; const lg = LightGraphs
mutable struct MatrixDiGraph{MT <: AbstractMatrix{Bool}} <: lg.AbstractGraph{Int}
  matrix::MT
end

The edges are simply mapping an entry (i,j) to a boolean (whether there is an
edge from i to j). Even though creating a graph type that can be directed
or undirected depending on the situation is possible, we are creating a type
that will be directed by default.

Implementing the core interface

We can now implement the core LightGraphs interface for this type, starting
with methods defined over the type itself, of the form function(g::MyType)

I’m not going to re-define each function here, their meaning can be found
by checking the help in a Julia REPL: ?LightGraphs.vertices or on the
documentation page.

lg.is_directed(::MatrixDiGraph) = true
lg.edgetype(::MatrixDiGraph) = lg.SimpleGraphs.SimpleEdge{Int}
lg.ne(g::MatrixDiGraph) = sum(g.m)
lg.nv(g::MatrixDiGraph) = size(g.m)[1]

lg.vertices(g::MatrixDiGraph) = 1:nv(g)

function lg.edges(g::MatrixDiGraph)
    n = lg.nv(g)
    return (lg.SimpleGraphs.SimpleEdge(i,j) for i in 1:n for j in 1:n if g.m[i,j])
end

Note the last function edges, for which the documentation specifies that we
need to return an iterator over edges. We don’t need to collect the comprehension
in a Vector, returning a lazy generator.

Some operations have to be defined on both the graph and a node, of the form
function(g::MyType, node).

lg.outneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[node, v]]
lg.inneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[v, node]]
lg.has_vertex(g::MatrixDiGraph, v::Integer) = v <= lg.nv(g) && v > 0

Out MatrixDiGraph type is pretty straight-forward to work with and all
required methods are easy to relate to the way information is stored in the
adjacency matrix.

The last step is implementing methods on both the graph and an edge of the
form function(g::MatrixDiGraph,e). The only one we need here is:

lg.has_edge(g::MatrixDiGraph,i,j) = g.m[i,j]

Optional mutability

Mutating methods were removed from the core interace some time ago,
as they are not required to describe a graph-like behavior.
The general behavior for operations mutating a graph is to return whether
the operation succeded. They consist in adding or removing elements from
either the edges or nodes.

import LightGraphs: rem_edge!, rem_vertex!, add_edge!, add_vertex!

function add_edge!(g::MatrixDiGraph, e)
    has_edge(g,e) && return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = true
end

function rem_edge!(g::MatrixDiGraph,e)
    has_edge(g,e) || return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = false
    return true
end

function add_vertex!(g::MatrixDiGraph)
    n = nv(g)
    m = zeros(Bool,n+1,n+1)
    m[1:n,1:n] .= g.m
    g.m = m
    return true
end

Testing our graph type on real data

We will use the graph type to compute the PageRank of

import SNAPDatasets
data = SNAPDatasets.loadsnap(:ego_twitter_d)
twitter_graph = MatrixDiGraph(lg.adjacency_matrix(data)[1:10,1:10].==1);
ranks = lg.pagerank(twitter_graph)

Note the broadcast check .==1, adjacency_matrix is specified to yield a
matrix of Int, so we use this to cast the entries to boolean values.

I took only the first 10 nodes of the graph, but feel free to do the same with
500, 1000 or more nodes, depending on what your machine can stand 🙈

Overloading non-mandatory functions

Some methods are already implemented for free by implementing the core interface.
That does not mean it should be kept as-is in every case. Depending on your
graph type, some functions might have smarter implementations, let’s see one
example. What MatrixDiGraph is already an adjacency_matrix, so we know
there should be no computation required to return it (it’s almost a no-op).

using BenchmarkTools: @btime

@btime adjacency_matrix(bigger_twitter)
println("why did that take so long?")
lg.adjacency_matrix(g::MatrixDiGraph) = Int.(g.m)
@btime A = lg.adjacency_matrix(bigger_twitter)
println("that's better.")

This should yield roughly:

13.077 ms (5222 allocations: 682.03 KiB)
why did that take so long?
82.077 μs (6 allocations: 201.77 KiB)
that's better.

You can fall down to a no-op by storing the matrix entries as Int directly,
but the type ends up being a bit heavier in memory, your type, your trade-off.

Conclusion

We’ve implemented a graph type suited to our need in a couple lines of Julia,
guided by the LightGraphs interface specifying how to think about our
graph instead of getting in the way of what to store. A lighter version
of this post can be read as slides.

As usual, ping me on Twitter for any
question or comment.

Bonus

If you read this and want to try building your own graph type, here are two
implementations you can try out, put them out in a public repo and show them off
afterwards:

  1. We created a type just for directed graphs, why bothering so much? You can create your own type which can be directed or not,
    either by storing the information in the struct or by parametrizing the type
    and getting the compiler to do the work for you.
  2. We store the entries as an AbstractMatrix{Bool}, if your graph is dense
    enough (how dense? No idea), it might be interesting to store entries as as
    BitArray.

Image source: GraphPlot.jl

Building our own graph type in Julia

By: Julia on μβ

Re-posted from: https://matbesancon.github.io/post/2018-08-17-abstract_graph/


This is an adapted post on the talk we gave with James
at JuliaCon 2018 in London. You can see the
original slides,
the video still requires a bit of post-processing.

Last week JuliaCon in London was a great and very condensed experience.
The two talks on LightGraphs.jl
received a lot of positive feedback and more than that, we saw
how people are using the library for a variety of use cases which is a great
signal for the work on the JuliaGraphs ecosystem
(see the lightning talk).

I wanted to re-build the same graph for people who prefer a post version to
my clumsy live explanations on a laptop not handling dual-screen well
(those who prefer the latter are invited to see the live-stream of the talk).

Why abstractions?

The LightGraphs library is built to contain as few elements as possible to get
anyone going with graphs. This includes:

  • The interface a graph type has to comply with to be used
  • Essential algorithms implemented by any graph respecting that interface
  • A simple, battery-included implementation based on adjacency lists

The thing is, if you design an abstraction which in fact has just one
implementation, you’re doing abstraction wrong. This talks was also a
reality-check for LightGraphs, are we as composable, extensible as we promised?

The reason for abstraction is also that minimalism has its price.
The package was designed as the least amount of complexity required to get
graphs working. When people started to use it, obviously they needed more
features, some of which they could code themselves, some other required
extensions built within LightGraphs. By getting the core abstractions right,
you guarantee people will be able to use it and to build on top with minimal
friction, while keeping it simple to read and contribute to.

Our matrix graph type

Let’s recall that a graph is a collection of nodes and a collection of
edges between these nodes. To keep it simple, for a graph of $n$ edges,
we will consider they are numbered from 1 to n. An edge connects a node $i$
to a node $j$, therefore all the information of a graph can be kept as an
adjacency matrix $M_{ij}$ of size $n \times n$:

$$M_{ij} = \begin{cases} 1, & \mbox{if edge (i $\rightarrow$ j) exists} \\ 0 & \mbox{otherwise}\end{cases}$$

We don’t know what the use cases for our type will be, and therefore,
we will parametrize the graph type over the matrix type:

import LightGraphs; const lg = LightGraphs
mutable struct MatrixDiGraph{MT <: AbstractMatrix{Bool}} <: lg.AbstractGraph{Int}
  matrix::MT
end

The edges are simply mapping an entry (i,j) to a boolean (whether there is an
edge from i to j). Even though creating a graph type that can be directed
or undirected depending on the situation is possible, we are creating a type
that will be directed by default.

Implementing the core interface

We can now implement the core LightGraphs interface for this type, starting
with methods defined over the type itself, of the form function(g::MyType)

I’m not going to re-define each function here, their meaning can be found
by checking the help in a Julia REPL: ?LightGraphs.vertices or on the
documentation page.

lg.is_directed(::MatrixDiGraph) = true
lg.edgetype(::MatrixDiGraph) = lg.SimpleGraphs.SimpleEdge{Int}
lg.ne(g::MatrixDiGraph) = sum(g.m)
lg.nv(g::MatrixDiGraph) = size(g.m)[1]

lg.vertices(g::MatrixDiGraph) = 1:nv(g)

function lg.edges(g::MatrixDiGraph)
    n = lg.nv(g)
    return (lg.SimpleGraphs.SimpleEdge(i,j) for i in 1:n for j in 1:n if g.m[i,j])
end

Note the last function edges, for which the documentation specifies that we
need to return an iterator over edges. We don’t need to collect the comprehension
in a Vector, returning a lazy generator.

Some operations have to be defined on both the graph and a node, of the form
function(g::MyType, node).

lg.outneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[node, v]]
lg.inneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[v, node]]
lg.has_vertex(g::MatrixDiGraph, v::Integer) = v <= lg.nv(g) && v > 0

Out MatrixDiGraph type is pretty straight-forward to work with and all
required methods are easy to relate to the way information is stored in the
adjacency matrix.

The last step is implementing methods on both the graph and an edge of the
form function(g::MatrixDiGraph,e). The only one we need here is:

lg.has_edge(g::MatrixDiGraph,i,j) = g.m[i,j]

Optional mutability

Mutating methods were removed from the core interace some time ago,
as they are not required to describe a graph-like behavior.
The general behavior for operations mutating a graph is to return whether
the operation succeded. They consist in adding or removing elements from
either the edges or nodes.

import LightGraphs: rem_edge!, rem_vertex!, add_edge!, add_vertex!

function add_edge!(g::MatrixDiGraph, e)
    has_edge(g,e) && return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = true
end

function rem_edge!(g::MatrixDiGraph,e)
    has_edge(g,e) || return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = false
    return true
end

function add_vertex!(g::MatrixDiGraph)
    n = nv(g)
    m = zeros(Bool,n+1,n+1)
    m[1:n,1:n] .= g.m
    g.m = m
    return true
end

Testing our graph type on real data

We will use the graph type to compute the PageRank of

import SNAPDatasets
data = SNAPDatasets.loadsnap(:ego_twitter_d)
twitter_graph = MatrixDiGraph(lg.adjacency_matrix(data)[1:10,1:10].==1);
ranks = lg.pagerank(twitter_graph)

Note the broadcast check .==1, adjacency_matrix is specified to yield a
matrix of Int, so we use this to cast the entries to boolean values.

I took only the first 10 nodes of the graph, but feel free to do the same with
500, 1000 or more nodes, depending on what your machine can stand 🙈

Overloading non-mandatory functions

Some methods are already implemented for free by implementing the core interface.
That does not mean it should be kept as-is in every case. Depending on your
graph type, some functions might have smarter implementations, let’s see one
example. What MatrixDiGraph is already an adjacency_matrix, so we know
there should be no computation required to return it (it’s almost a no-op).

using BenchmarkTools: @btime

@btime adjacency_matrix(bigger_twitter)
println("why did that take so long?")
lg.adjacency_matrix(g::MatrixDiGraph) = Int.(g.m)
@btime A = lg.adjacency_matrix(bigger_twitter)
println("that's better.")

This should yield roughly:

13.077 ms (5222 allocations: 682.03 KiB)
why did that take so long?
82.077 μs (6 allocations: 201.77 KiB)
that's better.

You can fall down to a no-op by storing the matrix entries as Int directly,
but the type ends up being a bit heavier in memory, your type, your trade-off.

Conclusion

We’ve implemented a graph type suited to our need in a couple lines of Julia,
guided by the LightGraphs interface specifying how to think about our
graph instead of getting in the way of what to store. A lighter version
of this post can be read as slides.

As usual, ping me on Twitter for any
question or comment.

Bonus

If you read this and want to try building your own graph type, here are two
implementations you can try out, put them out in a public repo and show them off
afterwards:
1. We created a type just for directed graphs, why bothering so much? You can create your own type which can be directed or not,
either by storing the information in the struct or by parametrizing the type
and getting the compiler to do the work for you.
2. We store the entries as an AbstractMatrix{Bool}, if your graph is dense
enough (how dense? No idea), it might be interesting to store entries as as
BitArray.


Image source: GraphPlot.jl

The cutting stock problem: part 2, solving with column generation

By: Julia on μβ

Re-posted from: https://mbesancon.github.io/post/2018-05-25-colgen2/


In the previous post, we explored a well-known integer optimization situation
in manufacturing, the cutting stock problem. After some details on the
decisions, constraints and objectives, we implemented a naive model in JuMP.

One key thing to notice is the explosion of number of variables and constraints
and the fact that relaxed solutions (without constraining variables to be
integers) are very far from actual feasible solutions.

We will now use an other way of formulating the problem, using a problem
decomposition and an associated solution method (column generation).

Re-stating the cutting stock problem

In the last post, we used two decisions: $Y_i$ stating if the big roll $i$ is
used and $X_{ij}$ expressing the number of cuts $j$ made in the roll $i$.
To minimize the number of rolls, it makes sense to put as many small cuts
as possible on a big roll. We could therefore identify saturating patterns,
that is, a combination of small cuts fitting on a big roll, such that no
additional cut can be placed, and then find the smallest combination of the
pattern satisfying the demand.

One problem remains: it is impossible to compute, or even to store in memory all
patterns, their number is exponentially big with the number of cuts, so we will
try to find the best patterns and re-solve the problem, using the fact that not
all possible patterns will be necessary.

This is exactly what the Dantzig-Wolfe decomposition does, it splits the problem
into a Master Problem MP and a sub-problem SP.

  • The Master Problem, provided a set of patterns, will find the best combination
    satisfying the demand.
  • The sub-problem, given an “importance” of each cut provided by the master
    problem, will find the best cuts to put on a new pattern.

This is an iterative process, we can start with some naive patterns we can think
of, compute an initial solution for the master problem, which will be feasible
but not optimal, move on to the sub-problem to try to find a new pattern
(or column in the optimization jargon, hence the term of column generation).

How do we define the “importance” of a cut $j$? The value of the dual variable
associated with this constraint will tell us that. This is not a lecture in
duality theory, math-eager readers can check out further documentation on the
cutting stock problem and duality in linear optimization.

Moreover, we are going to add one element to our model: excess cuts can be sold
at a price $P_j$, so that we can optimize by minimizing the net cost (production
cost of the big rolls minus the revenue from excess cuts).

New formulation

As in the previous post, we’re going to formulate the decisions and constraints
for the new version of the problem.

Decisions

At the master problem level, given a pattern $p$, the decision will be
$\theta_p$ (theta, yes Greek letters are awesome), the number of big rolls which
will be used with this pattern. $\theta_p$ is a positive integer.

The decision at the sub-problem level will be to find how many of each cut $j$
to fit onto one big roll, $a_j$.

For a pattern $p$, the number of times a cut $j$ appears is given by $a_{jp}$.

Constraints

The big roll size constraint is kept in the sub-problem, a pattern built
has to respect this constraint:
$$ \sum_j a_{j} \cdot W_j \leq L $$

The demand $D_j$ is met with all rolls of each pattern so it is kept at the master
level. The number of cuts of type $j$ produced is the sum of the number of this
cut on each patterns times the number of the pattern in a solution:

$$ NumCuts_j = \sum_p a_{jp} \cdot \theta_p \geq D_j$$

Objective formulation

At the master problem, we minimize the number of rolls, which is simply:
$$ \sum_{p} \theta_p $$

At the sub-problem, we are trying to maximize the gain associated with the need
for the demand + the residual price of the cuts. If we can find a worth using
producing compared to its production cost, it is added.

Implementation

As before, we will formulate the master and sub-problem using Julia with JuMP.
Just like the previous post, we use the Clp and Cbc open-source solvers.
We read the problem data (prices, sizes, demand) from a JSON file.

using JuMP
using Cbc: CbcSolver
using Clp: ClpSolver
import JSON

const res = open("data0.json", "r") do f
    data = readstring(f)
    JSON.Parser.parse(data)
end

const maxwidth = res["maxwidth"]
const cost = res["cost"]
const prices = Float64.(res["prices"])
const widths = Float64.(res["widths"])
const demand = Float64.(res["demand"])
const nwidths = length(prices)

cost is the production cost of a big roll.

Sub-problem

The subproblem is a function taking reduced costs of each cut and maximizing
the utility of the pattern it creates:

"""
    subproblem tries to find the best feasible pattern
    maximizing reduced cost and respecting max roll width
    corresponding to a multiple-item knapsack
"""
function subproblem(reduced_costs, sizes, maxcapacity)
    submodel = Model(solver = CbcSolver())
    n = length(reduced_costs)
    xs = @variable(submodel, xs[1:n] >= 0, Int)
    @constraint(submodel, sum(xs. * sizes) <= maxcapacity)
    @objective(submodel, Max, sum(xs. * reduced_costs))
    solve(submodel)
    return round.(Int,getvalue(xs)), round(Int,getobjectivevalue(submodel))
end

Initial master problem

We saw that the master problem finds a solution and then requires a new pattern
from the sub-problem. This is therefore preferable to start from an initial
feasible, otherwise we fall into a special case we’re not discussing here.
One initial solution would be to build one pattern per cut, with as many cuts as
we can, which is $floor(L/w_j)$.

function init_master(maxwidth, widths, rollcost, demand, prices)
    n = length(widths)
    ncols = length(widths)
    patterns = spzeros(UInt16,n,ncols)
    for i in 1:n
        patterns[i,i] = min(floor(Int,maxwidth/widths[i]),round(Int,demand[i]))
    end
    m = Model(solver = ClpSolver())
    θ = @variable(m, θ[1:ncols] >= 0)
    @objective(m, Min,
        sum(θ[p] * (rollcost - sum(patterns[j,p] * prices[j] for j=1:n)) for p in 1:ncols)
    )
    @constraint(m, demand_satisfaction[j=1:n], sum(patterns[j,p] * θ[p] for p in 1:ncols)>=demand[j])
    if solve(m) != :Optimal
        warn("No optimal")
    end
    return (m, getvalue(θ), demand_satisfaction, patterns)
end

We can compute the reduced costs from the dual values associated with the
demand and the prices of cuts

# getting the model and values
(m, θ, demand_satisfaction, patterns) = init_master(maxwidth, widths, cost, demand, prices);

# compute reduced costs
reduced_costs = getdual(demand_satisfaction)+prices;

# ask sub-problem for new pattern
newcol, newobj = subproblem(reduced_costs, widths, maxwidth)

Putting it all together

We can now build a column generation function putting all elements together and
performing the main iteration:

function column_generation(maxwidth, widths, rollcost, demand, prices; maxcols = 5000)
    (m, θ, demand_satisfaction, patterns) = init_master(maxwidth, widths, rollcost, demand, prices)
    ncols = nwidths
    while ncols <= maxcols
        reduced_costs = getdual(demand_satisfaction) + prices
        newcol, newobj = subproblem(reduced_costs, widths, maxwidth)
        netcost = cost - sum(newcol[j] * (getdual(demand_satisfaction)[j]+prices[j]) for j in 1:nwidths)
        println("New reduced cost: $netcost")
        if netcost >= 0
            return (:Optimal, patterns, getvalue(θ))
        end
        patterns = hcat(patterns, newcol)
        ncols += 1
        m = Model(solver = ClpSolver())
        θ = @variable(m, θ[1:ncols] >= 0)
        @objective(m, Min,
            sum(θ[p] * (rollcost - sum(patterns[j,p] * prices[j] for j=1:nwidths)) for p in 1:ncols)
        )
        @constraint(m, demand_satisfaction[j=1:nwidths], sum(patterns[j,p] * θ[p] for p in 1:ncols)>=demand[j])
        if solve(m) != :Optimal
            warn("No optimal")
            return (status(m), patterns, getvalue(θ))
        end
    end
    return (:NotFound, patterns, :NoVariable)
end

We’ve printed information along the computation to see what’s going on more
clearly, now launching it:

status, patterns, θ = column_generation(maxwidth, widths, cost, demand, prices, maxcols = 500);
New reduced cost: -443.18181818181824
New reduced cost: -375.0
New reduced cost: -264.0
New reduced cost: -250.0
New reduced cost: -187.5
New reduced cost: -150.0
New reduced cost: -150.0
New reduced cost: -107.14285714285711
New reduced cost: -97.5
New reduced cost: -107.14285714285734
New reduced cost: -72.0
New reduced cost: -53.571428571428555
New reduced cost: -53.125
New reduced cost: -50.0
New reduced cost: -43.40625
New reduced cost: -36.0
New reduced cost: -34.625
New reduced cost: -41.5
New reduced cost: -21.8515625
New reduced cost: -22.159090909090878
New reduced cost: -20.625
New reduced cost: -16.304347826086314
New reduced cost: -16.304347826086996
New reduced cost: -20.310344827586277
New reduced cost: -18.0
New reduced cost: -8.837209302325732
New reduced cost: -6.060606060606119
New reduced cost: 0.0

While the cost of a new pattern is negative, we can add it to the master and
keep running. This seems to make sense. Now, one thing to note, we have not
yet specified the integrality constraints, meaning that we don’t have integer
number of patterns. We can see that on the $\theta$ variable:

println(θ)
[0.0, 0.0, 0.0, ... 70.0, 0.0, 0.0, 0.0, 12.56, 46.86, 0.0, 0.0, 0.0, 0.0,
3.98, 0.0, 0.0, 21.5, 5.0, 31.12, 61.12, 33.58, 0.0, 0.0, 32.2, 44.0,
46.88, 19.0, 1.88, 16.42]
println(sum(θ))
446.1000000000001

We saw in the last post that the problem without integrality constraints is
a relaxation and therefore, can only yield a better result. This means that we
cannot have an integer solution using 446 big rolls or less, the minimum will
be 447 rolls. Let’s solve the problem with the same patterns, but adding the
integrality:

# compute initial integer solution:
# take worse case from linear solution, round up
intial_integer = ceil.(Int,θ);


"""
    From patterns built in the column generation phase, find an integer solution
"""function branched_model(patterns, demand, rollcost, prices; npatts = size(patterns)[2], initial_point = zeros(Int,npatts))
    npatts = size(patterns)[2]
    m = Model(solver = CbcSolver())
    θ = @variable(m, θ[p = 1:npatts] >= 0, Int, start = initial_point[p])
    @objective(m, Min,
        sum(θ[p] * (rollcost - sum(patterns[j,p] * prices[j] for j=1:nwidths)) for p in 1:npatts)
    )
    @constraint(m, demand_satisfaction[j=1:nwidths], sum(θ[p] * patterns[j,p] for p in 1:npatts) >= demand[j])
    status = solve(m)
    return (status, round.(Int,(getvalue(θ))))
end

Let’s see what the results look like:

status, θ_final = branched_model(patterns, demand, cost, prices; initial_point = intial_integer)
println(status)
:Optimal
println(sum(θ_final))
447

Given that we cannot do better than 447, we know we have the optimal
number of rolls.

Conclusion

After seeing what a mess integer problems can be in the first part, we used a
powerful technique called Dantzig-Wolfe decomposition, spliting the problem into
master and sub-problem, each handling a subset of the constraints.

Column generation is a technique making this decomposition usable in practice,
by adding only one or few columns (patterns) at each iteration, we avoid
an exponentially growing number of variables. The fact that JuMP is built as
an embedded Domain Specific Language in Julia makes it a lot easier to specify
problems and play around them. Most optimization specific modeling languages
are built around declarative features and get messy very quickly when
introducing some logic (like column generation iterations). Developers
could relate this technique to lazy value computation: we know all values are
there, but we just compute them whenever needed.

Hope you enjoyed reading this second post on the cutting stock problem. A
Jupyter notebook summing up all code snippets can be found at
this repository.
Feel free to ping me for feedback!

Note on performance

The column generation approach we just saw scales well to huge problems, but
this particular implementation can feel a bit slow at first. One recommended
thing is to do in such case is “warm-starting” the solver: give it a good
initial solution to start from. Since we built both the master and subproblem
as stateless functions, the model is being re-built from scratch each time.
The advantage is that any solver can be used, since some of them don’t support
warm starts.


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