By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2018-03-22-astroimages/

In the last days I worked on a very little
package: AstroImages.jl.
Written in the Julia programming language and part of
the JuliaAstro organization, this package
provides simple tools to visualize astronomical images saved
in FITS files.

AstroImages.jl is meant as an interface to popular Julia packages
like Images.jl
and Plots.jl, and
uses FITSIO.jl to read FITS files.

Currently the package supports only single-frame images, in gray scale.

AstroImages.jl is licensed under the MIT “Expat” License.

Installation

The package is available for Julia 0.6 but not yet registered, thus if you would
like to try it out use the following command in the Julia REPL:

julia> Pkg.clone("https://github.com/JuliaAstro/AstroImages.jl")

Usage

After installing the package, you can start using it with

julia> using AstroImages

For the time being the package provides only two functions:

• load (extending the load function
  from FileIO.jl), to load an
  extension of a FITS file
• AstroImage, a type constructor to create a new astronomical image object.

Reading extensions from FITS file

You can load and read the the first extension of a FITS file with load:

julia> load("file.fits")
1300×1200 Array{UInt16,2}:
[...]

The second, optional, argument of this load method is the number of the
extension to read. Read the third extension of the file with:

julia> load("file.fits", 3)
1300×1200 Array{UInt16,2}:
[...]

AstroImage type

The package provides a new type, AstroImage to integrate FITS images with
Julia packages for plotting and image processing. The AstroImage function has
the same syntax as load. This command:

julia> img = AstroImage("file.fits")
AstroImages.AstroImage{UInt16,ColorTypes.Gray}[...]

will read the first extension from the file.fits.

If you are working in a Jupyter notebook, an
AstroImage object is automatically rendered as a PNG image:

Plotting an AstroImage

An AstroImage object can be plotted with the Plots.jl package. Just use

julia> using Plots

julia> plot(img)

and the image will be displayed as a heatmap using your favorite backend.

Feedback welcome

If you test the package and would like to suggest improvements, feel free to
file an issue on GitHub or, even better, send a pull request to add new features
yourself!

By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-12-02-julia-assignment/

Today I realized by chance that in
the Julia programming language, like in C/C++,
assignment returns the assigned value.

To be more precise, assignment returns the right-hand side of the assignment
expression. This doesn’t appear to be documented in the current (as of December
2017) stable version of the manual, but
a
FAQ has
been added to its development version. Be aware of this subtlety, because you
can get unexpected results if you rely on the type of the returned value of an
assignment:

julia> let a
         a = 3.0
       end
3.0

julia> let a
         a::Int = 3.0
       end
3.0

In both cases the let blocks return 3.0, even though in the former case a
is really 3.0 and in the latter case it’s been marked to be an Int.

Assignment in the REPL

The fact that assignment returns the assigned value is a very handy feature.
First of all, this is probably the reason why in the Julia’s REPL the assigned
value is printed after an assignment, so that you don’t have to type again the
name of the variable to view its value:

julia> x = sin(2.58) ^ 2 + cos(2.58) ^ 2
1.0

Instead in Python, for example, assignment returns None. Thus, in the
interactive mode of the interpreter you have to type again the name of the
variable to check its value:

>>> from math import *
>>> x = sin(2.58) ** 2 + cos(2.58) ** 2
>>> x
1.0

Using assignment as condition

Another useful consequence of this feature is that you can use assignment as
condition in, e.g., if and while, making them more concise. For example,
consider this brute-force method to compute
the Euler’s number
stopping when the desired precision is reached:

euler = 0.0
i = 0
while true
    tmp = inv(factorial(i))
    euler += tmp
    tmp < eps(euler) && break
    i += 1
end

By using this feature the while can be simplified to:

euler = 0.0
i = 0
while (tmp = inv(factorial(i))) >= eps(euler)
    euler += tmp
    i += 1
end

Can this feature be a problem in Julia?

A common objection against having this feature in Python is that it is
error-prone. Quoting from
a Python tutorial:

Note that in Python, unlike C, assignment cannot occur inside expressions. C
programmers may grumble about this, but it avoids a common class of problems
encountered in C programs: typing = in an expression when == was intended.

This is indeed an issue because in C and Python conditions can be pretty much
anything that can be converted to plain-bit 0 and 1, including, e.g.,
characters. In addition, in Python a variable can freely change the type, so
that

Consider the following C program:

#include <stdio.h>

int main()
{
  char a = '\0';

  if (a = 0)
    printf("True\n");
  else
    printf("False\n");

  if (a == 0)
    printf("True\n");
  else
    printf("False\n");

  return 0;
}

which prints

False
True

Mistyping the condition (= instead of == or vice-versa) in this language can
lead to unexpected results.

Instead, this is mostly a non-issue in Julia because conditions can only be
instances of Bool type. One way to make an error in Julia is the following:

julia> a = false
false

julia> if (a = false)
           println("True")
       else
           println("False")
       end
False

julia> if a == false
           println("True")
       else
           println("False")
       end
True

but it is kind of useless to test equality with a Bool as a condition, since
you can use the Bool itself. The only other possible error you can make in
Julia I came up with is the following:

julia> a = false
false

julia> b = false
false

julia> if (a = b)
           println("True")
       else
           println("False")
       end
False

julia> if a == b
           println("True")
       else
           println("False")
       end
True

For this case I don’t have an easy replacement that could prevent you from an
error, but this is also an unlikely situation. The Julia style guide suggests
to
not parenthesize conditions.
Following this recommendation is a good way to catch the error in this case
because the assignment requires parens:

julia> if a = b
           println("True")
       else
           println("False")
       end

ERROR: syntax: unexpected "="

By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-11-21-hexadecimal-pi/

The
Bailey–Borwein–Plouffe formula is
one of the
several
algorithms to compute $\pi$.
Here it is:

$$\pi = \sum_{k = 0}^{\infty}\left[ \frac{1}{16^k} \left( \frac{4}{8k + 1} - \frac{2}{8k + 4} - \frac{1}{8k + 5} - \frac{1}{8k + 6} \right) \right]$$

What makes this formula stand out among other approximations of $\pi$ is that
it allows one to directly extract the $n$-th fractional digit of the
hexadecimal value of $\pi$ without computing the preceding ones.

Image credit: Cormullion, Julia
code
here.

The Wikipedia article about the Bailey–Borwein–Plouffe formula explains that the
$n$-th fractional digit (well, actually it’s the $n+1$-th) $d_{n}$ is
given by

$$d_{n} = 16 \left[ 4 \Sigma(n, 1) - 2 \Sigma(n, 4) - \Sigma(n, 5) - \Sigma(n, 6) \right]$$

where

$$\Sigma(n, j) = \sum_{k = 0}^{n} \frac{16^{n-k} \bmod (8k+j)}{8k+j} + \sum_{k = n+1}^{\infty} \frac{16^{n-k}}{8k+j}$$

Only the fractional part of expression in square brackets on the right side of
$d_{n}$ is relevant, thus, in order to avoid rounding errors, when we compute
each term of the finite sum above we can take only the fractional part. This
allows us to always use ordinary double precision floating-point arithmetic,
without resorting to arbitrary-precision numbers. In addition note that the
terms of the infinite sum get quickly very small, so we can stop the summation
when they become negligible.

Implementation in Julia language

Here is a Julia implementation of the algorithm to
extract the $n$-th fractional digit of $\pi$:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27



# Return the fractional part of x, modulo 1, always positive
fpart(x) = mod(x, one(x))

function Σ(n, j)
    # Compute the finite sum
    s = 0.0
    denom = j
    for k in 0:n
        s = fpart(s + powermod(16, n - k, denom) / denom)
        denom += 8
    end
    # Compute the infinite sum
    num = 1 / 16
    frac = num / denom
    while frac > eps(s)
        s     += frac
        num   /= 16
        denom += 8
        frac   = num / denom
    end
    return fpart(s)
end

pi_digit(n) =
    floor(Int, 16 * fpart(4Σ(n-1, 1) - 2Σ(n-1, 4) - Σ(n-1, 5) - Σ(n-1, 6)))

pi_string(n) = "0x3." * join(hex.(pi_digit.(1:n))) * "p0"

The pi_digit function gives the $n$-th hexadecimal fractional digit of
$\pi$ as a base-10 integer, and the pi_string function returns the first
$n$ hexadecimal digits of $\pi$ as a valid hexadecimal floating-point
literal:

julia> pi_digit(1)
2

julia> pi_digit(6)
10

julia> pi_string(1000)
"0x3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b54709179216d5d98979fb1bd1310ba698dfb5ac2ffd72dbd01adfb7b8e1afed6a267e96ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc16636920d871574e69a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59b59c30d5392af26013c5d1b023286085f0ca417918b8db38ef8e79dcb0603a180e6c9e0e8bb01e8a3ed71577c1bd314b2778af2fda55605c60e65525f3aa55ab945748986263e8144055ca396a2aab10b6b4cc5c341141e8cea15486af7c72e993b3ee1411636fbc2a2ba9c55d741831f6ce5c3e169b87931eafd6ba336c24cf5c7a325381289586773b8f48986b4bb9afc4bfe81b6628219361d809ccfb21a991487cac605dec8032ef845d5de98575b1dc262302eb651b8823893e81d396acc50f6d6ff383f442392e0b4482a484200469c8f04a9e1f9b5e21c66842f6e96c9a670c9c61abd388f06a51a0d2d8542f68960fa728ab5133a36eef0b6c137a3be4ba3bf0507efb2a98a1f1651d39af017666ca593e82430e888cee8619456f9fb47d84a5c33b8b5ebee06f75d885c12073401a449f56c16aa64ed3aa62363f77061bfedf72429b023d37d0d724d00a1248db0fead3p0"

While I was preparing this post I found an unregistered
package PiBBP.jl that implements the
Bailey–Borwein–Plouffe formula. This is faster than my code above, mostly
thanks to a function
for
modular exponentiation
more efficient than that available in Julia standard library.

Test results

Let’s check if the function is working correctly. We can use
the
parse function
to convert the string to a decimal floating point
number. IEEE 754 double precision
floating-point numbers have a 53-bit mantissa, amounting to $53 / \log_{2}(16) \approx 13$ hexadecimal digits:

julia> pi_string(13)
"0x3.243f6a8885a30p0"

julia> parse(Float64, pi_string(13))
3.141592653589793

julia> Float64(π) == parse(Float64, pi_string(13))
true

Generator expressions allow
us to obtain the decimal value of the number in a very simple way, without using
the hexadecimal string:

julia> 3 + sum(pi_digit(n)/16^n for n in 1:13)
3.141592653589793

We can use
the
arbitrary-precision BigFloat
to check the correctness of the result for even more digits. By default,
BigFloat numbers in Julia have a 256-bit mantissa:

julia> precision(BigFloat)
256

The result is correct for the first $256 / \log_{2}(16) = 64$ hexadecimal
digits:

julia> pi_string(64)
"0x3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89p0"

julia> BigFloat(π) == parse(BigFloat, pi_string(64))
true

julia> 3 + sum(pi_digit(n)/big(16)^n for n in 1:64)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

It’s possible to increase the precision of BigFloat numbers, to further test
the accuracy of the Bailey–Borwein–Plouffe formula:

julia> setprecision(BigFloat, 4000) do
           BigFloat(π) == parse(BigFloat, pi_string(1000))
       end
true

Multi-threading

Since the Bailey–Borwein–Plouffe formula extracts the $n$-th digit of $\pi$
without computing the other ones, we can write a multi-threaded version of
pi_string, taking advantage of native support
for
multi-threading
in Julia:

1
2
3
4
5
6
7



function pi_string_threaded(N)
    digits = Vector{Int}(N)
    Threads.@threads for n in eachindex(digits)
        digits[n] = pi_digit(n)
    end
    return "0x3." * join(hex.(digits)) * "p0"
end

For example, running Julia with 4 threads gives a 2× speed-up:

julia> Threads.nthreads()
4

julia> using BenchmarkTools

julia> pi_string(1000) == pi_string_threaded(1000)
true

julia> @benchmark pi_string(1000)
BenchmarkTools.Trial:
  memory estimate:  105.28 KiB
  allocs estimate:  2016
  --------------
  minimum time:     556.228 ms (0.00% GC)
  median time:      559.198 ms (0.00% GC)
  mean time:        559.579 ms (0.00% GC)
  maximum time:     564.502 ms (0.00% GC)
  --------------
  samples:          9
  evals/sample:     1

julia> @benchmark pi_string_threaded(1000)
BenchmarkTools.Trial:
  memory estimate:  113.25 KiB
  allocs estimate:  2018
  --------------
  minimum time:     270.577 ms (0.00% GC)
  median time:      271.075 ms (0.00% GC)
  mean time:        271.598 ms (0.00% GC)
  maximum time:     278.350 ms (0.00% GC)
  --------------
  samples:          19
  evals/sample:     1