By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/post/log1p/

edit: fixed bogus interaction of MathJax and code highlighting.

edit2: added benchmarks.

This is a follow-up on a question I asked on the Julia forums about calculating

\[
\text{log1p}(x) = \log(1+x)
\]

This calculation is tricky because if \(x \approx 0\),

\[
\log(1+x) \approx x
\]

while as \(x \to \infty\), \(\log(1+x)\) approaches \(\log(x)\), so simply using log(1+x) will not be as accurate as it could be. Numerical analysts have developed specialized methods for these calculations that try to get an accurate answer, and all programming languages serious about numerical calculations have an implementation either in the core language or a library.

Julia’s Base.log1p currently suggests that Base.Math.JuliaLibm.log1p would be preferable, but then I was wondering why isn’t that the default? So I decided to perform a trivial numerical experiment, calculating the error for both, and also benchmark the two methods.

Accuracy

The key question is what to compare the results with. One could compare to an existing "gold standard" implementation, or simply calculate the results using a higher precision representation. Fortunately, Julia has BigFloats available right out of the box.

The graph below shows the base-2 logarithm of the relative error for Base.log1p vs \(\log\_2(1+x)\); horizontal lines are log2(eps()) and log2(eps())+1. This suggests that Base.log1p is very accurate, but not as good as it could be when \(x \approx 0\).

The next plot shows the relative accuracy of the relative error above, comparing Base.Math.JuliaLibm.log1p to Base.log1p (lower values better). In these simulations, Base.Math.JuliaLibm.log1p is never worse, but sometimes significantly better (resulting in an extra binary digit of accuracy). This matters especially when \(x \approx 0\).

The next plot confirms this.

Speed

I also evaluated relative speed. The graph below shows the relative runtimes, obtained using BenchmarkTools.@belapsed. Values below \(1\) mean that Base.Math.JuliaLibm.log1p is faster: indeed, this seems to be the case except for values very close to \(0\), where there is a 10–20% performance penalty. At other values, Base.Math.JuliaLibm.log1p is 30–40% faster.

Conclusion

1. For values near \(0\), Base.Math.JuliaLibm.log1p is more accurate, at a slight performance cost.

2. For values away from \(0\), it is at least as accurate as Base.log1p, and 30—40% faster.

To me, this suggests that Base.Math.JuliaLibm.log1p should be the default method — mostly because the extra accuracy is more important to me than the slight performance cost.

Code is available below.

# consistent random numbers
srand(UInt32[0xfd909253, 0x7859c364, 0x7cd42419, 0x4c06a3b6])

"""
    err(x, [prec])

Return two values, which are the log2 relative errors for calculating
`log1p(x)`, using `Base.log1p` and `Base.Math.JuliaLibm.log1p`.

The errors are calculated by compating to `BigFloat` calculations with the given
precision `prec`.
"""
function err(x, prec = 1024)
    yb = log(1+BigFloat(x, prec))
    e2(y) = Float64(log2(abs(y-yb)/abs(yb)))
    e2(log1p(x)), e2(Base.Math.JuliaLibm.log1p(x))
end

z = exp.(randn(20000)*10)       # z > 0, Lognormal(0, 10)
x = z .- 1                      # x > -1
es = map(err, x)                # errors

using Plots; gr()               # plots

scatter(log2.(z), first.(es), xlab = "log2(x+1)", ylab = "log2 error of Base.log1p",
        legend = false)
hline!(log2(eps())-[0,1])
Plots.svg("Base_log1p_error.svg")
scatter(log2.(z), last.(es) .- first.(es), xlab = "log2(x+1)",
        ylab = "improvement (Base.Math.JuliaLibm.log1p)", legend = false)
Plots.svg("JuliaLibm_improvement.svg")
scatter(log2.(z), last.(es), xlab = "log2(x+1)",
        ylab = "log2 error of Base.Math.JuliaLibm.log1p", legend = false)
hline!(log2(eps())-[0,1])
Plots.svg("JuliaLibm_log1p_error.svg")

######################################################################
# WARNING: these run for a very long time
######################################################################
using BenchmarkTools

z = exp.(vcat(randn(200)*10, rand(200)*0.1)) # z > 0, more values around 
x = z .- 1                                   # x > -1
b1 = [@belapsed log1p($x) for x in x]        # WARNING: takes forever
b2 = [@belapsed Base.Math.JuliaLibm.log1p($x) for x in x] # ditto

scatter(log2.(z), b2 ./ b1, xlab = "log2(x+1)",
        ylab = "time Math.JuliaLibm.log1p / log1p", legend = false, yticks = 0:0.2:1.2)
hline!([1])
Plots.svg("relative_time.svg")

By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/post/emacs-julia-customizations/

I find the following customizations very useful for editing Julia
code in Emacs. Add them to julia-mode-hook, eg

(defun customize-julia-mode ()
  "Customize julia-mode."
  (interactive)
  ;; my customizations go here
  )

(add-hook 'julia-mode-hook 'customize-julia-mode)

Highlight FIXME/TODO/…

When I just want to note something in a comment for future reference,
I prefer to have certain words highlighted. You can use something like
this:

(font-lock-add-keywords nil
                        '(("\\<\\(FIXME\\|TODO\\|QUESTION\\|NOTE\\)"
                        1 font-lock-warning-face t)))

This is what it looks like:

            chklapackerror(info[])
            if any(ifail .!= 0)
                # TODO: better error message / type
                error("failed to converge eigenvectors:\n$(nonzeros(ifail))")
            end

Highlight symbols

After

(require 'highlight-symbol)

add a hook for

(local-set-key [(control ?c) ?s] 'highlight-symbol-at-point)
(local-set-key [(control ?c) ?n] 'highlight-symbol-next)
(local-set-key [(control ?c) ?p] 'highlight-symbol-prev)

This highlights symbols with C-c s:

function issymmetric(A::AbstractMatrix)
    indsm, indsn = indices(A)
    if indsm != indsn
        return false
    end
    for i = first(indsn):last(indsn), j = (i):last(indsn)
        if A[i,j] != transpose(A[j,i])
            return false
        end
    end
    return true
end

Fill docstrings

This is useful if you want to use M-q on docstrings.

(defun julia-fill-string ()
  "Fill a docstring, preserving newlines before and after triple quotation marks."
  (interactive)
  (if (and transient-mark-mode mark-active)
      (fill-region (region-beginning) (region-end) nil t)
    (cl-flet ((fill-if-string ()
                              (when (or (looking-at (rx "\"\"\""
                                                        (group
                                                         (*? (or (not (any "\\"))
                                                                 (seq "\\" anything))))
                                                        "\"\"\""))
                                        (looking-at (rx "\""
                                                        (group
                                                         (*? (or (not (any "\\"))
                                                                 (seq "\\" anything))))
                                                        "\"")))
                                (let ((start (match-beginning 1))
                                      (end (match-end 1)))
                                  ;; (ess-blink-region start end)
                                  (fill-region start end nil nil nil)))))
      (save-excursion
        (let ((s (syntax-ppss)))
          (when (fourth s) (goto-char (ninth s))))
        (fill-if-string)))))

Add

(local-set-key (kbd "M-q") 'julia-fill-string)

to the mode hook.

Highlight things after column 80

I add this to the mode hook:

(set-fill-column 80)

I also use whitespace globally:

(require 'whitespace)
(setq whitespace-style '(face empty tabs lines-tail trailing))
(global-whitespace-mode t)

This is what it looks like:

    QR{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ)
end
QR(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QR{T,typeof(factors)}(factors, τ)

Hungry delete-mode

Add this to the mode hook:

(hungry-delete-mode)

and then backspace and delete will remove all whitespace near the point in the relevant direction.

In case you are wondering, the theme is alect-light.

By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/post/julia-local-test/

I run Pkg.update() fairly often to stay up to date and benefit from
the latest improvements of various packages. I rarely ever pin to a
specific package version, but I occasionally checkout master for
some packages, especially if I am contributing.

Despite updating regularly, I found that the documentation subtly diverged from what I was experiencing for some packages. After looking into the issue, I learned that I was 2–3 minor versions behind despite updating regularly. For example, when I would

Pkg.update("ForwardDiff")

I would be told that there is a new version, and to get it I should update ReverseDiff and Optim. But

Pkg.update("ForwardDiff", "ReverseDiff", "Optim")

would just run quietly without updating.

I could not figure out the cause for this and did not want to get sidetracked debugging it, so I decided to wipe the package directory and start over. However, in order to do this, I had to make sure that no code is lost, especially for local packages. First, I moved my local packages into a separate directory, and added that to LOAD_PATH in .juliarc.jl:

push!(LOAD_PATH, expanduser("~/src/julia-local-packages/"))

Then I ran multi-git-status to make sure that there were no unpushed changes. Finally, I deleted the package directory and reinstalled everything. Surprisingly, Pkg.add ran much faster than before.

In case I have to do this again, I decided to keep my local packages separate — the only drawback is that Pkg.test now can’t find them. A workaround is below, using some code from Base.Pkg:

"""
    local_test(pkgname, [coverage])

Find and test a package in `LOAD_PATH`. Useful when the package is outside
`Pkg.dir()`.

Assumes the usual directory structure: package has the same name as the module,
the main file is in `src/Pkgname.jl`, while tests are in `test/runtests.jl`.
"""
function local_test(pkgname; coverage::Bool=false)
    module_path = Base.find_in_path(pkgname, nothing)
    src_dir, module_file = splitdir(module_path)
    dir = normpath(src_dir, "..")
    test_path = joinpath(dir, "test", "runtests.jl")
    @assert isfile(test_path) "Could not find $(test_path)"
    Base.cd(dir) do
        try
            color = Base.have_color? "--color=yes" : "--color=no"
            codecov = coverage? ["--code-coverage=user"] : ["--code-coverage=none"]
            compilecache = "--compilecache=" * (Bool(Base.JLOptions().use_compilecache) ? "yes" : "no")
            julia_exe = Base.julia_cmd()
            run(`$julia_exe --check-bounds=yes $codecov $color $compilecache $test_path`)
            info("$module_file tests passed")
        catch err
            Base.Pkg.Entry.warnbanner(err, label="[ ERROR: $module_file ]")
        end
    end
end

Compared to simply include("wherever-it-is/runtests.jl"), this has the advantage of running a separate Julia process, so your workspace does not contaminate the test environment and in case of segfaults, the parent process won’t be affected.

Hopefully, the code above will be obsolete once Pkg3 is released, but until then it is a useful workaround.

edit: function above was corrupted during copy-paste, corrected.