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&gt; 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&gt; (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&gt; 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&gt; 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&gt; 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&gt; 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&gt; makeIntBig!(ex2)
 
julia&gt; 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

How to profile julia code?

This is kind of part 13 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.
It’s just about: How to make your code faster.

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

In the last post the benchmark looked good for classic sudokus but not so much for killer sudokus and okayish for graph coloring I would say. As always one can try to improve the performance by rethinking the code but today I mostly want to share how one can check which parts of the code are slow and maybe improve on that i.e by reducing the memory consumption.

First step in improving the performance is always to know what is slow at the moment. We don’t want to take hours in improving one line of code which is only responsible for 1% of the code run time.

An important step is to reduce the memory used.

Memory usage

Let’s find out how much memory we use in the the normal sudoku case at first because we call less functions there:

> julia
julia> using Revise
julia> include("benchmark/sudoku/cs.jl")
julia> main()
julia> @time main(;benchmark=true)

quite a normal setup.
The first main() call is for compiling and checking that we solve every sudoku. The @time macro gives us the time but also the memory usage.

1.009810 seconds (3.69 M allocations: 478.797 MiB, 18.51% gc time)

Attention: This includes all kind of timing not just the solve and is run on my Laptop so don’t compare that to the benchmark results from last time.

We use ~500MiB of memory for solving 95 sudokus. Aha and now?
Yeah that doesn’t give us much detail. Therefore we make memory measurements per line now. Close the current session first.

> julia --track-allocation=user
julia> include("benchmark/sudoku/cs.jl")
julia> main(;benchmark=true)

Close the session.

In the file explorer of your IDE you can now see a couple of *.mem files.
I suppose the most memory is used in the constraint all_different.jl

i.e at the bottom of the file you see:

           - function all_different(com::CoM, constraint::Constraint, value::Int, index::Int)
 12140913     indices = filter(i->i!=index, constraint.indices)
  7885120     return !any(v->issetto(v,value), com.search_space[indices])
           - end

7885120 means that ~7.5MiB are used in that line. (This is a cumulative measurement).

The function itself checks whether we can set value at index or whether that isn’t possible (violates the constraint).

We do a filtering step first here which is creates a new array in the first line and for the second line we again create an array and fill it with true or false and then iterate over it to check…

Coloring in Scientific Publications

By: Hendrik Ranocha -- Julia blog

Re-posted from: https://ranocha.de/blog/colors/

In this blog post, I want to share my current default choice of plotting options for articles
and presentations. After some motivation, you will find the code at the end.

When you have been in science a couple of years ago and have used MATLAB or matplotlib in Python,
you probably know the old default color map jet (or rainbow). I don’t want to reiterate the
arguments against this specific color map, since there are lots of nice ressources about it,
such as the
2007 IEEE article “Rainbow Color Map (Still) Considered Harmful” of Borland and Russel.
Today, you probably won’t use jet anymore, since it has been replaced by more sensible default
color maps in many visualization tools. If you can spare twenty minutes, I recommend watching
the talk describing the color maps in matplotlib 2
or reading matplotlib’s color map tutorial.

However, there are also other visualizations besides pseudocolor plots.
In particular, line and scatter plots are very common, e.g. for convergence studies or to show
the time evolution of a quantity of interest in a numerical simulation. And there are often several
of these quantities which shall be visualized in the same figure.
While good visualization tools vary the appearance of subsequent plots in the same figure,
just changing the color of the lines is often not sufficient, at least in my opinion. Firstly,
relying solely on colors can make the life of people affected by some sort of color blindness
really hard. Secondly, if the figure shall appear in a paper, it will probably be printed in
grayscale, since many publishers charge a lot for color figures in the printed version.
While most people will use the online
version of articles, many of them prefer to print an article to read it. And printing
something in color still costs a lot more than just printing everything in grayscale. If you
had some color figures in you PhD thesis and needed to print it several times for examination
or needed to submit it to a publisher, you will probably know what I mean. (I had only two
color pages in my PhD thesis, which is
totally okay.)

To sum up, I think it is worth investing some time in creating figures for publications.
Since I don’t want to rely exclusively on colors for accessible plots, I prefer to vary
the general shape of the plots additionally. For line plots, the lines can be solid, dotted, dashed,
dash-dotted, and there are lots of variations of the latter. To keep things simple, I prefer
to use only these four standard choices of linestyles in most figures.
For scatter plots, there are lots of markers which can be used. But I don’t want to specify
the linestyle and marker type for every plot, which can be annoying for plots
generated in a loop. Luckily, matplotlib provides a nice interface to set the style of plots.
In Python, I use

import matplotlib.pyplot as plt

# line cyclers adapted to colourblind people
from cycler import cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

If you use only line plots via plt.plot, you can set the default cycler via

plt.rc("axes", prop_cycle=line_cycler)

If you want to choose the default cycler for each new axes, you can use

fig, ax = plt.subplots(1, 1)
ax.set_prop_cycle(line_cycler)

The corresponding Julia code is

import PyPlot; plt = PyPlot
using PyCall
using LaTeXStrings

# line cyclers adapted to colourblind people
cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

# matplotlib's standard cycler
standard_cycler = cycler("color", ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"])

plt.rc("axes", prop_cycle=line_cycler)

At the time of writing, it is not possible to use different markers with plt.scatter plots,
as described in this issue on GitHub.
Nevertheless, you can get scatter plots by using the marker_cycler defined above.

I also use some additional customizations that can be found in the source files linked at the end of this blog.
More information can be found in the customizing matplotlib tutorial.

Let’s now have a look at the results. Using the line_cycler defined above, a plot looks like
line_plot_python
The default cycler of matplotlib yields
line_plot_python_standard
These lines are also easily distinguishable (unless you’re red blind; but even in this case, there should be
some differences of the brightness, although it is not as easy as for other people).
However, if printed in grayscale, it becomes hardly possible to identify the lines correctly using the
standard cycler of matplotlib.


`line_cycler`	`standard_cycler`

The results for scatter plots are similar: Using different marker shapes really helps to identify the plots correctly.

`marker_cycler`	`standard_cycler`

In my opinion, it is worth investing some time in creating accessible plots. That’s why I like to use the
color blindness simulator
to check how figures look like for people affected by some sort of color blindness and in grayscale (monochromatic view).
The Python code and
Julia code used to produce
the figures in this post contain also additional plotting options I like to use. As all code published
in this blog, these scripts are released under the terms of the MIT license.

Python Code

#!/usr/bin/python3

import matplotlib.pyplot as plt

# line cyclers adapted to colourblind people
from cycler import cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

# matplotlib's standard cycler
standard_cycler = cycler("color", ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"])

plt.rc("axes", prop_cycle=line_cycler)

plt.rc("text", usetex=True)
plt.rc("text.latex", preamble=r"\usepackage{newpxtext}\usepackage{newpxmath}\usepackage{commath}\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=200)
plt.rc("legend", loc="best", fontsize="medium", fancybox=True, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=10, markeredgewidth=2.5)


import os
directory = os.path.dirname(os.path.abspath(__file__))


# plots
plt.close("all")

import numpy as np
x = np.linspace(0., 1., 200)

plt.rc("axes", prop_cycle=line_cycler)
plt.figure()
plt.plot(x,   np.sin(np.pi*x), label="$\sin(\pi x)$")
plt.plot(x,   np.cos(np.pi*x), label="$\cos(\pi x)$")
plt.plot(x, 2*np.sin(np.pi*x), label="$2 \sin(\pi x)$")
plt.plot(x, 2*np.cos(np.pi*x), label="$2 \cos(\pi x)$")
plt.legend()
plt.xlabel("$x$")
plt.xlim(x[0], x[-1])
plt.ylabel("Function Values")
plt.savefig(os.path.join(directory, "line_plot_python.png"), bbox_inches="tight")

plt.rc("axes", prop_cycle=standard_cycler)
plt.figure()
plt.plot(x,   np.sin(np.pi*x), label="$\sin(\pi x)$")
plt.plot(x,   np.cos(np.pi*x), label="$\cos(\pi x)$")
plt.plot(x, 2*np.sin(np.pi*x), label="$2 \sin(\pi x)$")
plt.plot(x, 2*np.cos(np.pi*x), label="$2 \cos(\pi x)$")
plt.legend()
plt.xlabel("$x$")
plt.xlim(x[0], x[-1])
plt.ylabel("Function Values")
plt.savefig(os.path.join(directory, "line_plot_python_standard.png"), bbox_inches="tight")


import numpy as np
x = np.linspace(0., 1., 20)

plt.rc("axes", prop_cycle=marker_cycler)
plt.figure()
plt.plot(x,   np.sin(np.pi*x), label="$\sin(\pi x)$")
plt.plot(x,   np.cos(np.pi*x), label="$\cos(\pi x)$")
plt.plot(x, 2*np.sin(np.pi*x), label="$2 \sin(\pi x)$")
plt.plot(x, 2*np.cos(np.pi*x), label="$2 \cos(\pi x)$")
plt.legend()
plt.xlabel("$x$")
plt.xlim(x[0], x[-1])
plt.ylabel("Function Values")
plt.savefig(os.path.join(directory, "marker_plot_python.png"), bbox_inches="tight")

plt.rc("axes", prop_cycle=standard_cycler)
plt.rc("lines", linewidth=1.5, markersize=6, markeredgewidth=1)
plt.figure()
plt.scatter(x,   np.sin(np.pi*x), label="$\sin(\pi x)$")
plt.scatter(x,   np.cos(np.pi*x), label="$\cos(\pi x)$")
plt.scatter(x, 2*np.sin(np.pi*x), label="$2 \sin(\pi x)$")
plt.scatter(x, 2*np.cos(np.pi*x), label="$2 \cos(\pi x)$")
plt.legend()
plt.xlabel("$x$")
plt.xlim(x[0], x[-1])
plt.ylabel("Function Values")
plt.savefig(os.path.join(directory, "marker_plot_python_standard.png"), bbox_inches="tight")

Julia Code

import PyPlot; plt = PyPlot
using PyCall
using LaTeXStrings

# line cyclers adapted to colourblind people
cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

# matplotlib's standard cycler
standard_cycler = cycler("color", ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"])

plt.rc("axes", prop_cycle=line_cycler)

plt.rc("text", usetex=true)
plt.rc("text.latex", preamble="\\usepackage{newpxtext}\\usepackage{newpxmath}\\usepackage{commath}\\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=200)
plt.rc("legend", loc="best", fontsize="medium", fancybox=true, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=10, markeredgewidth=2.5)


directory = @__DIR__


# plots
plt.close("all")

x = range(0., 1., length=200)

plt.rc("axes", prop_cycle=line_cycler)
plt.figure()
plt.plot(x,   sinpi.(x), label=L"$\sin(\pi x)$")
plt.plot(x,   cospi.(x), label=L"$\cos(\pi x)$")
plt.plot(x, 2*sinpi.(x), label=L"$2 \sin(\pi x)$")
plt.plot(x, 2*cospi.(x), label=L"$2 \cos(\pi x)$")
plt.legend()
plt.xlabel(L"$x$")
plt.xlim(x[1], x[end])
plt.ylabel("Function Values")
plt.savefig(joinpath(directory, "line_plot_julia.png"), bbox_inches="tight")

plt.rc("axes", prop_cycle=standard_cycler)
plt.figure()
plt.plot(x,   sinpi.(x), label=L"$\sin(\pi x)$")
plt.plot(x,   cospi.(x), label=L"$\cos(\pi x)$")
plt.plot(x, 2*sinpi.(x), label=L"$2 \sin(\pi x)$")
plt.plot(x, 2*cospi.(x), label=L"$2 \cos(\pi x)$")
plt.legend()
plt.xlabel(L"$x$")
plt.xlim(x[1], x[end])
plt.ylabel("Function Values")
plt.savefig(joinpath(directory, "line_plot_julia_standard.png"), bbox_inches="tight")


x = range(0., 1., length=20)

plt.rc("axes", prop_cycle=marker_cycler)
plt.figure()
plt.plot(x,   sinpi.(x), label=L"$\sin(\pi x)$")
plt.plot(x,   cospi.(x), label=L"$\cos(\pi x)$")
plt.plot(x, 2*sinpi.(x), label=L"$2 \sin(\pi x)$")
plt.plot(x, 2*cospi.(x), label=L"$2 \cos(\pi x)$")
plt.legend()
plt.xlabel(L"$x$")
plt.xlim(x[1], x[end])
plt.ylabel("Function Values")
plt.savefig(joinpath(directory, "marker_plot_julia.png"), bbox_inches="tight")

plt.rc("axes", prop_cycle=standard_cycler)
plt.rc("lines", linewidth=1.5, markersize=6, markeredgewidth=1)
plt.figure()
plt.scatter(x,   sinpi.(x), label=L"$\sin(\pi x)$")
plt.scatter(x,   cospi.(x), label=L"$\cos(\pi x)$")
plt.scatter(x, 2*sinpi.(x), label=L"$2 \sin(\pi x)$")
plt.scatter(x, 2*cospi.(x), label=L"$2 \cos(\pi x)$")
plt.legend()
plt.xlabel(L"$x$")
plt.xlim(x[1], x[end])
plt.ylabel("Function Values")
plt.savefig(joinpath(directory, "marker_plot_julia_standard.png"), bbox_inches="tight")

juliabloggers.com

A Julia Language Blog Aggregator