By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2019-05-03-julia-get-pointer-value/

Yesterday a question in the Julia’s Slack
channel (request an invite
here!) prompted a lively discussion about
how to deal with pointers in Julia. As the history of the channel is ephemeral,
I wrote this post to keep note of some interesting points.

In the C world, the operation of getting the value of a pointer is called
dereferencing. When
dealing with C libraries in Julia, we may need to get the value of a pointer
into a Julia object. The Julia manual is quite detailed about accessing data
in
pointers
and I warmly recommend reading that first.

Note that, usually, in Julia it’s not possible to “dereference” a pointer in the
C-sense, because most of the operations we’re going to see will copy the data
into Julia structures.

unsafe_load is
the generic function to access the value of a pointer and copy it to a Julia
object. However, depending on the specific data types involved, there may be
different and more efficient solutions. For example, you can obtain the value
of the pointer to an array with
unsafe_wrap,
while for strings there is
unsafe_string.

Example 1: get the current time

Let’s have a look at a simple time-related example.

julia> result = Ref{Int64}(0)
Base.RefValue{Int64}(0)

julia> ccall((:time, "libc.so.6"), Int64, (Ptr{Int64},), result)
1556838715

julia> result
Base.RefValue{Int64}(1556838715)

By calling the time function from the C standard library, we have saved in
result the current number of seconds since the Unix epoch. We now want to
turn this number into an instance of the C struct
tm. To do this, we can use the
localtime function. However, if we want to later use the tm struct in Julia
we need to define a Julia mutable structure with the same layout as the C one
and define the show method to have a fancy output:

julia> mutable struct Ctm
           sec::Cint
           min::Cint
           hour::Cint
           mday::Cint
           mon::Cint
           year::Cint
           wday::Cint
           yday::Cint
           isdst::Cint
       end

julia> function Base.show(io::IO, t::Ctm)
           print(io, t.year + 1900, "-", lpad(t.mon, 2, "0"), "-",
                 lpad(t.mday, 2, "0"), "T", lpad(t.hour, 2, "0"), ":",
                 lpad(t.min, 2, "0"), ":", lpad(t.sec, 2, "0"))
       end

We are now going to use unsafe_load to copy the result of the call to
localtime into an instance of the Ctm structure:

julia> localtime = ccall((:localtime, "libc.so.6"), Ptr{Ctm}, (Ptr{Int64},), result)
Ptr{Ctm} @0x00007f730fe9d300

julia> unsafe_load(localtime)
2019-04-03T00:11:55

julia> dump(unsafe_load(localtime))
Ctm
  sec: Int32 55
  min: Int32 11
  hour: Int32 0
  mday: Int32 3
  mon: Int32 4
  year: Int32 119
  wday: Int32 5
  yday: Int32 122
  isdst: Int32 1

Of course, if you want to deal with dates you can use the
Dates standard library
instead of playing with system calls. Have also a look at the other
time-related functions in the Base.Libc
module of Julia.

Example 2: copy from a data type to another one

Now that we’ve seen how to use unsafe_load, we can define the following
function to “dereference” a generic pointer to a Julia object (with the caveat
that we’re copying it’s value to the new object!):

julia> dereference(ptr::Ptr) = unsafe_load(ptr)
dereference (generic function with 1 method)

julia> dereference(T::DataType, ptr::Ptr) = unsafe_load(Ptr{T}(ptr))
dereference (generic function with 2 methods)

This dereference function has two arguments: the type that will wrap the data,
and the input pointer. The first argument is optional, and defaults to the type
of the data pointed by ptr. The first method is effectively an alias of the
plain unsafe_load. The second method is interesting because it allows us to
copy the value of a pointer of a certain data type to an object of another type
with the same memory layout.

For example, let’s define a simple Julia mutable structure called Bar, get
an instance of it, and its pointer:

julia> mutable struct Bar
           x::UInt64
       end

julia> b = Bar(rand(UInt64))
Bar(0x55139e3fd61e43a4)

julia> pointer_from_objref(b)
Ptr{Nothing} @0x00007f1f72eb7850

julia> ptr = Ptr{Bar}(pointer_from_objref(b))
Ptr{Bar} @0x00007f1f72eb7850

pointer_from_objref
gave us a pointer to Nothing (that is, the C void), we have then casted it
to a pointer to Bar and assigned this pointer to ptr. We now define another
data structure with the same memory layout as Bar and use
dereference to “dereference” ptr as an instance of the new structure:

julia> mutable struct Foo
           a::UInt16
           b::UInt16
           c::UInt32
       end

julia> f = dereference(Foo, ptr)
Foo(0x43a4, 0xd61e, 0x55139e3f)

Note that

julia> UInt64(f.a) | UInt64(f.b) << 16 | UInt64(f.c) << 32
0x55139e3fd61e43a4

julia> b.x
0x55139e3fd61e43a4

as expected.

Remember that the term “dereference” is a stretch, because data is copied from
the pointer. In fact:

julia> dereference(Foo, ptr) |> pointer_from_objref
Ptr{Nothing} @0x00007f1f72a62440

julia> dereference(Foo, ptr) |> pointer_from_objref
Ptr{Nothing} @0x00007f1f72fb3840

julia> dereference(Bar, ptr) |> pointer_from_objref
Ptr{Nothing} @0x00007f1f730054e0

julia> dereference(Bar, ptr) |> pointer_from_objref
Ptr{Nothing} @0x00007f1f73007340

Conclusions

What we’ve seen here is nothing new for programmers already familiar with the C
language, but at the same time shows how easy and natural can be to communicate
from Julia with C programs.

Special bonus

If you feel really bold and creative, you can even abuse the * operator to do
very weird things like:

julia> Base.:*(ptr::Ptr) = dereference(ptr)

julia> Base.:*(T::DataType, ptr::Ptr) = dereference(T, ptr)

julia> *(ptr)
Bar(0x55139e3fd61e43a4)

julia> Foo *ptr
Foo(0x43a4, 0xd61e, 0x55139e3f)

julia> Bar *ptr
Bar(0x55139e3fd61e43a4)

However, the * operator shouldn’t be really used to mock C syntax is an
awkward way (we’re not actually dereferencing and the syntax Foo *ptr is for
declaration of a pointer rather than actual dereferencing), so don’t tell people
that I told you to use this!

By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2019-02-21-julia-do-vs-python-with/

I’ve recently started coding in Python more than I’ve used to, and it’s
interesting for me to have a closer look to its features and compare them with
similar Julia’s ones.

The do-block in the Julia programming language and the with statement in
Python can appear to be similar because of some common uses (most notably,
opening a stream and automatically closing it at the end), but are fundamentally
different. In this post I will review them.

with statement in Python

The with
statement has
been introduced in Python with PEP
343 to simplify the control flow of
try/finally
statements.

A good explanation of this feature is given at
http://effbot.org/zone/python-with-statement.htm. In general, you can use
try/finally to handle unmanaged resources like streams:

def controlled_execution(callback):
    set things up
    try:
        callback(thing)
    finally:
        tear things down

def my_function(thing):
    do something with thing

controlled_execution(my_function)

For an object of a class, it is possible to reorganize the above control flow by
using the with statement:

class controlled_execution:
    def __enter__(self):
        set things up
        return thing
    def __exit__(self, type, value, traceback):
        tear things down

with controlled_execution() as thing:
    do something with thing

The __enter__ method prepares the object for later use and returns thing,
which can be then consumed in the body of with. The __exit__ method will be
called in any case at the end of the body of with and can be used to close the
previously opened resources and handle any exception raised within the with
body, if necessary.

Opening a file with with

A common application of the with statement is the open
function:

with open(filename, mode='w') as f:
    do something with f

Just to practically see how the __enter__ and __exit__ methods can be
defined, we implement a very simple class, MyFile, to replicate the behaviour
of the builtin open function when opening a file:

class MyFile:
    def __init__(self, file):
        self.file = file
    def __enter__(self):
        return self.file
    def __exit__(self, type, value, traceback):
        return self.file.close()
    def read(self):
        return self.file.read()

def MyOpen(name, mode='r'):
    return MyFile(open(name, mode))

with MyOpen('file.txt') as f:
    f.read()

do-block in Julia

The
do-block
in Julia addresses a different problem: write in a different way the passing of
a function as argument to another function.

In a simple way, any function func accepting a function as first argument:

function func(f::Function, x, y, z...)
    # do something with arguments x, y, z..., and call f at some point
end

can be called as

func(x, y, z...) do args
    # put here the body of function f(args) that will
    # be passed as first argument to func(f, x, y, z...)
end

For example, the
filter
enables filtering a collection by removing those elements for which the provided
function returns false. We can select the odd numbers between 1 and 100
divisible by 3 with

julia> filter(x -> isinteger(x/3) && isodd(x), 1:100)
17-element Array{Int64,1}:
  3
  9
 15
 21
 27
 33
 39
 45
 51
 57
 63
 69
 75
 81
 87
 93
 99

Here we passed an anonymous
function
as first argument to filter in order to select the desired elements of the
collection passed as second argument. We can rewrite this by using the
do-block in the following way:

julia> filter(1:100) do x
           isinteger(x/3) && isodd(x)
       end
17-element Array{Int64,1}:
  3
  9
 15
 21
 27
 33
 39
 45
 51
 57
 63
 69
 75
 81
 87
 93
 99

The code in the body of the do-block is used to create an anonymous function,
whose arguments are specified after do, that is passed as first argument to
the function before the do.

Opening a file with do

A common use of the do-block is to open a file and automatically close it when
you’re done, just like the with statement in Python. For example:

open("outfile", "w") do io
    write(io, data)
end

In Julia’s Base code, this is achieved by extending the open
method in the
following way:

function open(f::Function, args...; kwargs...)
    io = open(args...; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

This method, which accepts a function as first argument, is what makes it
possible to use the do-block in the example above: after effectively opening
the file, the body of the do-block is executed in the try
statement,
and in the end the finally
clause
takes care of safely closing the file, even in case of errors while executing
the function f.

Summary

The with statement in Python is a syntactic sugar for try/finally. This
construct is bound to classes, this means that the __enter__/__exit__
methods are bound to a specific class and they can be defined only once as they
have a fixed signature, so you can’t have multiple __enter__/__exit__
handlers within the same class, or independently from a class. In order for the
with statement to work, both __enter__ and __exit__ must be defined in any
case, however if you don’t really need the try/finally construct you can
still take advantage of the handy with syntax, for example by making the
__exit__ method dummy.

The do-block in Julia is not bound to any type, the only condition to use it
is to have a function taking a function as first argument. The function may
contain a control flow construct, like try/finally, but this is not
mandatory. You can also have different methods for a function to be used with
the do-block, depending on the types and/or the number of the other arguments.

By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2018-12-12-cosmology-unitful/

Cosmology.jl is the
cosmology calculator for
Julia. After you define a cosmological model, you can use it to perform some
operations like computation of angular diameter distance or comoving radial
distance at given redshift values, or calculate the age of the Universe in a
certain cosmological model at a specific redshift.

The package doesn’t have a full documentation, but the information in the
README.md file should be sufficient to get you started, if you’re already
familiar with the topic.

Cosmology.jl is a registered package, so you can install it in the Julia REPL
with the package manager:

pkg> add Cosmology

New release

A few days ago version
v0.5.0 has
been released. The main new feature is the long-awaited integration with
Unitful.jl and
UnitfulAstro.jl. This means
that the functions defined in this package now return a number with proper
physical units. For example:

julia> using Cosmology

julia> c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
Cosmology.FlatWCDM{Float64}(0.7, 0.0, 0.7, 0.3, 0.0, -0.9, 0.1)

julia> comoving_radial_dist(c, 0.6)
2162.83342244173 Mpc

Previously, instead, the units were attached to the name of the function,
without possibility to easily change units. These methods are now deprecated:

julia> comoving_radial_dist_mpc(c, 0.6)
┌ Warning: `comoving_radial_dist_mpc(c::AbstractCosmology, z; kws...)` is deprecated, use `ustrip(comoving_radial_dist(c::AbstractCosmology, z; kws...))` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
2162.83342244173

In addition, now you can select a different unit for the output by simply
passing the desired output unit as first argument:

julia> using Unitful

julia> comoving_radial_dist(u"ly", c, 0.6)
7.054219146683016e9 ly

Integration with Measurements.jl

Cosmology.jl is one of the several examples of how easy is in Julia to combine
together structures from different and independent packages. Thanks to Julia’s
type system:

• the data structure used in this package to represent a cosmological model and
  the units data structure in Unitful.jl don’t know anything about each other
  (in Cosmology.jl the numbers are just multiplied by the appropriate unit, in
  the most natural way possible),
• units in Unitful and the Measurement structure in
  Measurements.jl don’t
  know anything about each other,
• by now you can already tell that Measurement and the data structure for
  cosmological models don’t know anything about each other,

yet, we can use numbers with uncertainties as parameters of cosmological models
or for the redshift. This enable use to propagate the uncertainties through the
operations performed by Cosmology.jl getting the results with the appropriate
physical units.

As an example, we can define a cosmological model using the parameters
determined by the Planck
collaboration in 2015:

julia> using Cosmology, Measurements

julia> planck2015 = cosmology(h=0.6774±0.0046, OmegaM=0.3089±0.0062, Tcmb=2.718±0.021, Neff=3.04±0.33)
Cosmology.FlatLCDM{Measurement{Float64}}(0.6774 ± 0.0046, 0.6910099044166828 ± 0.006200002032729847, 0.3089 ± 0.0062, 9.009558331733912e-5 ± 5.020543222494152e-6)

Now we calculate the comoving volume at redshift z = 3.62±0.04:

julia> z = 3.62 ± 0.04
3.62 ± 0.04

julia> cv = comoving_volume(planck2015, z)
1471.004984671155 ± 45.25029615715926 Gpc^3

Measurements.jl provides a utility,
Measurements.uncertainty_components,
to determine the contribution to the total uncertainty of a quantity:

julia> Measurements.uncertainty_components(cv.val)
Dict{Tuple{Float64,Float64,UInt64},Float64} with 5 entries:
  (3.04, 0.33, 0x0000000000000005)     => 0.0499315
  (0.3089, 0.0062, 0x0000000000000003) => 27.5208
  (3.62, 0.04, 0x0000000000000006)     => 19.8259
  (0.6774, 0.0046, 0x0000000000000002) => 29.952
  (2.718, 0.021, 0x0000000000000004)   => 0.0348057

This means that the major contributions to the uncertainty of cv comes from
the Hubble parameter (0.6774±0.0046) and the matter density (0.3089±0.0062).

We can also compute, in this cosmological model, the age of the Universe today
and at redshift z = 1.42:

julia> age(planck2015, 0) # Age of the Universe today
13.79748128449975 ± 0.12164948254546123 Gyr

julia> age(planck2015, 1.42) # Age of the Universe at redshift z = 1.42
4.481579852131797 ± 0.05168067053818648 Gyr

lookback_time gives the difference between age at redshift 0 and age at
redshift z:

julia> lookback_time(planck2015, 1.42)
9.315901432385681 ± 0.07270835053850357 Gyr

Note that uncertainties are always propagated taking care of the correlation
between quantities. Indeed, we can check that the sum of the lookback time at
redshift z and the age of the Universe at the same redshift is equal to the age
of the Universe today, up to numerical errors of the integrals involved in the
calculations:

julia> lookback_time(planck2015, 1.42) + age(planck2015, 1.42)
13.797481284517477 ± 0.12164948254345108 Gyr