Hi Stéphane, nice post! I have a number of comments and suggestions which
you may find useful. I accompanied these comments with some demo code, you
can find here <https://gist.github.com/carlobaldassi/10312215>.
A) Generic to Julia
A.1) as Ivar said, you can use Unicode characters if you want;
A.2) sarting from the next Julia version, an expression like "1 - eye(n)"
will raise a warning and urge you to use the dotted-minus operator, i.e. "1
.- eye(n)". Better start getting in the habit of doing that! Also, in Julia
v0.3 you can write the same expression as "ones(n,n) - I", which is kind of
nicer maybe;
A.3) nested loops like
for i = 1:n
for j = 1:n
# ...
end
end
can be written concisely as
for i = 1:n, j=1:n
# ...
end
B) Specific to GLPK
B.1) it is not advisable to do "using GLPK", it's better to do "import GLPK",
because of the large number of common names exported by GLPK, which makes
name collisions likely
B.2) one easy way to get the "ia, ja, ar" vector for the sparse
representation is to call the function "findnz" on a sparse matrix, i.e.
you just do
ia,ja,ar = findnz(sparse(A))
Furthermore, the GLPK.jl package does that for you (it's one of the few
improvements over the native GLPK), so that you can simply call
GLPK.load_matrix(lp, sparse(A))
B.3) Indeed, using GNU MP, GLPK.exact on this problem is as fast as
GLPK.simplex. The solutiion found is "0.10714285714285714" (the best
approximation of 3/28 in the 64 bit precision) instead of
"0.10714285714285715" (which is the next floating point value, i.e. there's
a 1 ulp difference. It's likely that setting the tolerance will fill this
gap too. Unfortunately, it's not possible to get the rational value from
GLPK.exact, and there is no linear programming solver which accepts
rational arguments in Julia at the moment).
C) On linear programming in Julia
C.1) Instead of using GLPK directly, it's way easier to use a higher level
Julia package, such as
MathProgBase<https://github.com/JuliaOpt/MathProgBase.jl>,
which uses GLPK internally (via a glue
package<https://github.com/JuliaOpt/GLPKMathProgInterface.jl>)
but does most of the low-level settings for you. It also has the advantage
that the same program can use different solvers (GLPK, Clp, Gurobi, CPLEX,
Mosek, ...) by simply changing a function argument, which I think is pretty
cool.
In the gist which I linked
before<https://gist.github.com/carlobaldassi/10312215>there is an example which
demonstrates how easy and concise and more
readable the program becomes.
C.2) Going up yet another level, you can use
JuMP<https://github.com/JuliaOpt/JuMP.jl>,
which uses MathProgBase itself, but introduces special syntax aimed at
allowing to express these kind of optimization problems much more
naturally. Again, it lets you choose which solver you want to use under the
hood.
Hope this helps, cheers.
On Wednesday, April 9, 2014 7:08:19 PM UTC+2, Stéphane Laurent wrote:
>
> Hello guys,
>
> I hope you'll enjoy this article on my
> blog<http://stla.github.io/stlapblog/posts/KantorovichWithJulia.html>
> .
>
> If you're able to use GNU MP on your machine, would you be able to find
> *3/28* ?
>
> Any other comment is welcomed !
>
>
