Category Archives: Julia

Finding the maximum cardinality matching in a bipartite graph

This is kind of part 14 of the series about: How to build a Constraint Programming solver?
But it is also interesting for you if you don’t care about constraint programming.
In this post you’ll learn about bipartite graphs and matchings.

Anyway these are the other posts so far in the ongoing series:

All of that can only happen because of the generous support of some people.
I’m very happy that I can cover my server expenses now because of my Patrons.
To keep this going please consider supporting this blog on Patreon. If you do I will create a better/faster blog experience soon and you will get posts like this earlier!

If you want to just keep reading -> Here you go 😉

Last time I realized that the library I’m using doesn’t implement the easy algorithm I need which is bipartite maximum cardinality matching instead of bipartite maximum matching. This means they are solving a harder problem which I don’t have 😀

Additionally I like coding things from scratch so using a library is not the best choice but I am of course also limited in the amount of time I have and if people built it before it is often a good practice to not reinvent the wheel. Nevertheless this is a learning project where it is often good to know how a wheel works from the ground up.

The problem

We want to know whether a setting of numbers and possibilities in an all_different constraint is satisfiable or not.

For example we have 3 cells we want to fill with the numbers 1-3 and each should be used exactly once.

In the following examples our values are always at the top and the variables at the bottom. The goal is to match our variables to the values such that the all different constraint is satisfied. In this example there is no such matching:

Bipartite maximum cardinality smaller than #values

In a more precise definition a matching is a set of edges such that two edges don’t share the same vertex: Each vertex has at most one matching edge connected to it.

We are always interested in a matching which uses the most edges (maximum cardinality matching).

This is for example a matching (blue edges):

Small matching

We need this maximum cardinality matching to check whether it is possible to fulfill the constraint and later to also remove edges which would make it impossible when used (explained in the all different section).

How can we find those matchings?

Let’s have a look at a slightly bigger example:

Bipartite graph

I should mention that these graphs are called bipartite graphs as they have two…

Learning algorithmic techniques: dynamic programming

By: Julia on EPH

Re-posted from: https://ericphanson.com/blog/2019/learning-algorithmic-techniques-dynamic-programming/

Three nice techniques In the past months, I’ve found myself really appreciating some nice techniques for constructing algorithms. These are probably quite familiar to those with a computer science background, but are new to me. There are three in particular that I have in mind; I’ll just highlight the first two here, and then discuss the third in detail.
The first such technique I learned about was solving the traveling salesman problem via mixed-integer programming, using lazy constraints.

Evaluating Arbitrary Precision Integer Expressions in Julia using Metaprogramming

While watching the Mathologer masterclass on power sums

I came across a challenge to evaluate the following sum
1^{10}+2^{10}+\cdots+1000^{10}

This can be easily evaluated via brute force in Julia to yield 

function power_cal(n)
  a=big(0)
  for i=1:n
    a+=big(i)^10
  end
  a 
end
 
julia> power_cal(1000)
91409924241424243424241924242500

Note the I had to use big to makes sure we are using BigInt in the computation. Without that, we would be quickly running in into an overflow issue and we will get a very wrong number.

In the comment section of the video I found a very elegant solution to the above sum, expressed as

(1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 – 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000 = 91409924241424243424241924242500

If I try to plug this into the Julia, I get 

julia> (1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + r1000^5-(1/2) * 1000^3 + (5/66) * 1000
-6.740310541071357e18

This negative answer is not surprising at all, because we obviously ran into an overflow. We can, of course, go through that expression and modify all instances of Int64 with BigInt by wrapping it in the big function. But that would be cumbersome to do by hand.

The power of Metaprogramming

In Julia, metaprogramming allows you to write code that creates code, the idea here to manipulate the abstract syntax tree (AST) of a Julia expression. We start to by “quoting” our original mathematical expressing into a Julia expression. In the at form it is not evaluated yet, however we can always evaluate it via eval. 

julia> ex1=:((1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000)
:((((((1 / 11) * 1000 ^ 11 + (1 / 2) * 1000 ^ 10 + (5 / 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 / 2) * 1000 ^ 3) + (5 / 66) * 1000)
 
julia> dump(ex1)
Expr
  head: Symbol call
  args: Array{Any}((3,))
    1: Symbol +
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol -
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol +
            2: Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol -
                2: Expr
                3: Expr
            3: Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol ^
                2: Int64 1000
                3: Int64 5
        3: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol *
            2: Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol /
                2: Int64 1
                3: Int64 2
            3: Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol ^
                2: Int64 1000
                3: Int64 3
    3: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol *
        2: Expr
          head: Symbol call
          args: Array{Any}((3,))
            1: Symbol /
            2: Int64 5
            3: Int64 66
        3: Int64 1000
 
julia> eval(ex1)
-6.740310541071357e18

The output of dump show the follow AST in all its glory (…well almost the depth is a bit truncated). Notice that here all our numbers are interpreted as Int64.
Now we walk through the is AST and change all occurrences of Int64 with BigInt by using the big function.

function makeIntBig!(ex::Expr)
  args=ex.args
  for i in eachindex(args)
       if args[i] isa Int64
          args[i]=big(args[i])
       end
       if args[i] isa Expr
          makeIntBig!(args[i])
       end
   end
end
 
julia> ex2=copy(ex1)
:((((((1 / 11) * 1000 ^ 11 + (1 / 2) * 1000 ^ 10 + (5 / 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 / 2) * 1000 ^ 3) + (5 / 66) * 1000)
 
julia> makeIntBig!(ex2)
 
julia> eval(ex2)
9.14099242414242434242419242425000000000000000000000000000000000000000000000014e+31

We see an improvement here, but the results are not very satisfactory. The divisions yield BigFloat results, which had a tiny bit of floating point errors. Can we do better?
Julia has support for Rational expressions baked in. We can use that improve the results. We just need to search for call expressions the / symbol and replace it by the // symbol. For safety we just have to makes sure the operands are as subtype of Integer. 

function makeIntBig!(ex::Expr)
    args=ex.args
    if ex.head == :call && args[1]==:/ && 
              length(args)==3 && 
              all(x->typeof(x) <: Integer,args[2:end]) 
        args[1]=://
        args[2]=big(args[2])
        args[3]=big(args[3])
    else 
        for i in eachindex(args)
           if args[i] isa Int64
              args[i]=big(args[i])
           end
           if args[i] isa Expr
              makeIntBig!(args[i])
           end
        end
    end
end
 
julia> ex2=copy(ex1);
 
julia> makeIntBig!(ex2)
 
julia> eval(ex2)
91409924241424243424241924242500//1

Now that is much better! We have not lost any precision and we ended us with a Rational expression.

Finally, we can build a macro so the if we run into such expressions in the future and we want to evaluate them, we could just conveniently call it. 

macro eval_bigInt(ex)
   makeIntBig!(ex)
   ex
end

and we can now simply evaluate our original expression as 

julia> @eval_bigInt (1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000
91409924241424243424241924242500//1