text:
# ---
# layout: post
# title: "Automata"
# date: 2018-11-29 00:00:00 +0000
# categories: blog
# mathjax: true
# ---

# > To play along with this post, you'll be able to find the source on [github](https://github.com/cormullion/cormullion.github.io), and you'll need Julia (version 1), and the packages Luxor, Colors, and ColorSchemes. This post was written using Fredrik Ekre's [Literate.jl](https://github.com/fredrikekre/Literate.jl){:target="_blank"}.

# > I don't need to give you any references—you'll be spoilt for choice if you start googling.

# Here in the UK we're usually better at closing railway stations than opening new ones. Dr Beeching famously closed thousands of them (more than 50% of the total) in the 1960s. Occasionally, though, new stations are built, and one recent addition to the rail network is the station at [Cambridge North](https://www.networkrail.co.uk/running-the-railway/our-routes/anglia/cambridge-north-station/){:target="_blank"}, built mainly to serve the Cambridge Science Park. (I used to walk along the railway tracks there at weekends before the old disused line was re-opened, but that's another story...)

# When the station was opened, in 2017, there was a bit of chatter about the decorative panels used for the building.

# ![image label](IMAGEFOLDER/Cmglee_Cambridge_North_front_night.jpg)

# (Image by [cmglee at wikipedia](https://commons.wikimedia.org/wiki/File:Cmglee_Cambridge_North_front_night.jpg){:target="_blank"}, licensed CC-SA)

# Thsee patterns were described by the designers at Atkins (the contractors) like this:

# >The station is wrapped in three equal bands of aluminium panels which have been perforated with a design derived from John Horton Conway’s “Game of Life” theories which he established while at Gonville and Caius College, Cambridge in 1970. These beautiful, delicate panels ensure passive security to ground floor glazed areas, assist with wayfinding while crossing the footbridge and allow the building to transform its appearance between day and night through sensitive backlighting.

# This quote is from the original article ([archived here](https://web.archive.org/web/20170219171551/https://www.atkinsglobal.com/en-GB/projects/cambridge-north-station){:target="_blank"}), but it's since been corrected, most probably in response to the anguished cries of thousands of Cambridge nerds who swiftly pointed out that, in fact, these patterns weren't *Life* as we know it, but were actually one-dimensional cellular automata, and so linked, not so much with Cambridge's John Horton Conway, as with noted Oxford/CalTech alumnus Stephen "Mr Mathematica" Wolfram.

# I think the designers did a great job using these simple graphics to "create a harmonic relationship with the scientific research and industry of the Cambridge Colleges and nearby Science Park".

# ### So what is a one-dimensional cellular automaton?

# A cellular automaton is a mathematical model that creates patterns automatically according to simple rules. In its simplest, one-dimensional form, it's a row of empty squares, each of which can be occupied or empty. The rules for the model determine how sensitive the occupants of each square are to their immediate left and right neighbours, and whether, after some unspecified time, they survive or die. Sometimes empty squares can miraculously become occupied. Because a square has just two immediate neighbours, there are 8 different cases to consider, ranging from all empty ("□□□" or `000` in binary) to all full ("■■■" or `111`).

# The rule determines how one generation changes to the next by specifying the outcome for each of the 8 cases. So, for example, an empty square surrounded by empty neighbours can continue to be empty (false, or 0), or it can produce a new occupant in the next generation, (true, or 1). There are 256 different combinations, since each of the 8 cases can be either 0 or 1.

# The wikipedia has [this nice animation](https://commons.wikimedia.org/wiki/File:One-d-cellular-automate-rule-30.gif){:target="_blank"} showing how the rule produces the next generation of squares.

# ![image label](IMAGEFOLDER/one-d-cellular-automaton-rule-110.gif)

# (Well, it has this nice animation now...)

# ### First steps

# To explore these simple automata, I started ([^1]) by making a Julia structure:

mutable struct CA
    rule::Int64
    cells::BitArray{1}
    colorstops::Array{Float64, 1}
    ruleset::BitArray{1}
    generation::Int64
    function CA(rule, ncells = 100)
        cells                    = falses(ncells)
        colorstops               = zeros(Float64, ncells)
        ruleset                  = binary_to_array(rule)
        cells[length(cells) ÷ 2] = true
        generation               = 1
        new(rule, cells, colorstops, ruleset, generation)
    end
end

# The `cells` array can hold trues or falses. The `colorstops` array is eventually going to hold some color information. The middle cell is seeded with a single starter value.

# The `binary_to_array()` function just converts a binary number to a bit array (I suspect there's a quicker way).

function binary_to_array(n)
    a = BitArray{1}()
    for c in 7:-1:0
        k = n >> c
        push!(a, (k & 1 == 1 ? true : false))
    end
    return a
end

# The `rules()` function takes the values of an individual and its neighbours and applies the rule that determines its state for the next generation:

function rules(ca::CA, a, b, c)
    lng = length(ca.ruleset)
    return ca.ruleset[mod1(lng - (4a + 2b + c), lng)]
end

# And a `nextgeneration()` function applies the rule to all the cells. I decided to make it wrap around, so that the final cell considers the first cell as one of its neighbours.

function nextgeneration(ca::CA)
    l = length(ca.cells)
    nextgen = falses(l)
    for i in 1:l
        left   = ca.cells[mod1(i - 1, l)]
        me     = ca.cells[mod1(i, l)]
        right  = ca.cells[mod1(i + 1, l)]
        nextgen[i] = rules(ca, left, me, right)
    end
    ca.cells = nextgen
    ca.generation += 1
    return ca
end

# We'll also teach Julia how to show an automaton in the terminal:

Base.show(io::IO, ::MIME"text/plain", ca::CA) =
    print(io, join([c ? "■" : " " for c in ca.cells]))

# So now we can create a cellular automaton by providing a rule number (using the default of 100 cells):

ca = CA(30)

# update the automaton like this:

nextgeneration(ca)

# and show a historical diagram of its evolution:

for i in 1:30
    display(nextgeneration(ca))
end

# ![image label](IMAGEFOLDER/terminal-output.png)

# ### Some graphics

# The REPL display is more or less functional, but I want to play with the graphic output, so (you guessed):

using Luxor, Colors

function draw(ca::CA, cellwidth=10)
    lng = length(ca.cells)
    for i in 1:lng
        if ca.cells[i] == true
            pt = Point(-(lng ÷ 2) * cellwidth + i * cellwidth, 0)
            box(pt, cellwidth, cellwidth, :fill)
        end
    end
end

@png begin
    ca = CA(30, 200)
    sidelength = 4
    ## start at the top
    translate(boxtopcenter(BoundingBox()) + sidelength)
    for i in 1:200
        draw(ca, sidelength)
        nextgeneration(ca)
        translate(Point(0, sidelength))
    end
end 800 850 "images/automata/simple-ca.png"

# ![image label](IMAGEFOLDER/simple-ca.png)

# You can call `nextgeneration()` without displaying the results, of course. This lets you jump into the future history of an automaton at warp speed.

@png begin
    ca = CA(110, 200)
    translate(boxtopcenter(BoundingBox()) + sidelength)
    sidelength = 4
    ## into the future
    for _ in 1:200_000
        nextgeneration(ca)
    end
    for _ in 1:195
        draw(ca, sidelength)
        nextgeneration(ca)
        translate(Point(0, sidelength))
    end
end 800 800 "images/automata/simple-ca-future.png"

# ![image label](IMAGEFOLDER/simple-ca-future.png)

# I found that sometimes drawing them from left to right looked better (like videos recorded on phones...?):

@png begin
    ca = CA(30)
    translate(boxmiddleleft(BoundingBox()) + sidelength)
    rotate(-π/2)
    sidelength = 3.5
    for i in 1:320
        draw(ca, sidelength)
        nextgeneration(ca)
        translate(Point(0, sidelength))
    end
end 1000 400 "images/automata/simple-landscape-ca.png"

# ![image label](IMAGEFOLDER/simple-landscape-ca.png)

# ### And now in color

# So far I haven't used the color information that's stored.

# The `nextgeneration()` function can be updated with instructions about how to modify the color of the next generation, based on the current set of three cells.

function nextgeneration(ca::CA)
    l = length(ca.cells)
    nextgen = falses(l)
    for i in 1:l
        left   = ca.cells[mod1(i - 1, l)]
        me     = ca.cells[mod1(i, l)]
        right  = ca.cells[mod1(i + 1, l)]
        nextgen[i] = rules(ca, left, me, right)
        if me == 1
            ca.colorstops[i] = mod(ca.colorstops[i] + 0.1, 1)
        elseif left == 1 && right == 1
            ca.colorstops[i] = mod(ca.colorstops[i] - 0.1, 1)
        else
            ca.colorstops[i] = 0
        end
    end
    ca.cells = nextgen
    ca.generation += 1
    return ca
end

# and the `draw()` function can be adapted to make use of the color information. I decided to avoid tackling RGB color value transformations for a first pass, so the single value between 0 and 1 is used to select a color from a color map.

using ColorSchemes

function drawcolor(ca::CA, cellwidth=10;
        scheme=ColorSchemes.leonardo)
    lng = length(ca.cells)
    for i in 1:lng
        if ca.cells[i] == true
            sethue(get(scheme, ca.colorstops[i]))
            pt = Point(-(lng ÷ 2) * cellwidth + (i * cellwidth), 0)
            box(pt, cellwidth, cellwidth, :fill)
        end
    end
end

@svg begin
    background("darkorchid4")
    ca = CA(135, 65)
    ## randomize start state
    ca.cells = rand(Bool, length(ca.cells))
    sidelength = 5
    translate(boxmiddleleft(BoundingBox()) + sidelength)
    rotate(-π/2)
    for i in 1:195
        drawcolor(ca, sidelength, scheme=ColorSchemes.cubehelix)
        nextgeneration(ca)
        translate(Point(0, sidelength))
    end
end 1000 400 "images/automata/simple-color-ca.svg"

# ![image label](IMAGEFOLDER/simple-color-ca.svg)

# This could lead to hours of entertainment (depending on your definition of fun, of course). I uploaded a few experiments that didn't turn out too badly [on Flickr](https://www.flickr.com/photos/153311384@N03/){:target="_blank"}. The current rules show a kind of literal winning streak, as a cell that remains occupied for many generations ends up being brightly illuminated.

# I think these images look quite good when scaled up. It only takes about a second to draw these, but would take much longer to stick them on the wall:

# ![image label](IMAGEFOLDER/office-layout-1.png)

# ![image label](IMAGEFOLDER/office-layout-2.png)

# The rules for specifying a change in color could do with some kind of systematic definition, perhaps, such that, say, "rule C81" means "increase colorstop by amount if previous parent is 1, decrease it if previous uncle-aunt is 1", and so on. Then you could pass a set of color rules to the drawing function. (*uncle-aunt*? I couldn't find a word for something that is either an uncle or an aunt, but not a parent...)

# Instead of drawing simple squares, it's possible to draw other shapes. I'm quite fond of the squircle — you can change the `rt` parameter to get different shapes:

function drawcolorsquircle(ca::CA, cellwidth=10;
        scheme=ColorSchemes.leonardo)
    lng = length(ca.cells)
    for i in 1:lng
        if ca.cells[i] == true
            sethue(get(scheme, ca.colorstops[i]))
            pt = Point(-(lng ÷ 2) * cellwidth + (i * cellwidth), 0)
            squircle(pt, cellwidth, cellwidth, rt=6.0, :fill)
        end
    end
end

@svg begin
    background("navyblue")
    ca = CA(110, 30)
    ## randomize start state
    ca.cells = rand(Bool, length(ca.cells))
    translate(boxmiddleleft(BoundingBox()) + sidelength)
    rotate(-π/2)
    sidelength = 10
    for i in 1:100
        drawcolorsquircle(ca, sidelength, scheme=ColorSchemes.Dark2_8)
        nextgeneration(ca)
        translate(Point(0, sidelength))
    end
end 1000 400 "images/automata/simple-color-ca-squircle.svg"

# ![image label](IMAGEFOLDER/simple-color-ca-squircle.svg)

# (This SVG is quite big, and won't display in Juno. But it should load in a browser.)

# ### Getting around to it

# It occurred to me that you could take a rectangular array and wrap it into a circle.

function drawsector(ca::CA, cellwidth=10;
        scheme=ColorSchemes.leonardo,
        centralradius = 10)
    lng = length(ca.cells)
    width = lng * cellwidth
    angulargap = 2π/lng
    for i in 1:lng
        sethue(get(scheme, ca.colorstops[i]))
        innerradius = centralradius
        outerradius = centralradius + cellwidth
        startang    = rescale(i, 1, lng, 0, 2π)
        endang      = startang + angulargap
        if ca.cells[i] == true
            sector(O, innerradius, outerradius, startang, endang, :fillstroke)
        end
    end
end

function drawrule(rulenumber, pos)
    @layer begin
        translate(pos)
        sethue("black")
        text(string(rulenumber), halign=:center, valign=:middle)
        rotate(-π/2)
        ca = CA(rulenumber, 50)
        ## randomize start state
        ## ca.cells = rand(Bool, length(ca.cells))
        sidelength = 2
        setline(0.0)
        for n in 5:sidelength:30
            drawsector(ca, sidelength,
            scheme=ColorSchemes.klimt,
            centralradius = n)
            nextgeneration(ca)
        end
    end
end

# This shows all the rules in this circular form.

let
    ## best in SVG or PDF, but PNG is faster
    Drawing(1200, 1200, "images/automata/color-sector-assembly.svg")
    origin()
    background("azure")
    fontsize(8)
    for (pos, n) in Tiler(1200, 1200, 16, 16)
        drawrule(n-1, pos)
    end
    finish()
    preview()
end

# ![image label](IMAGEFOLDER/color-sector-assembly.svg)

# (Again, this image is quite demanding for SVG and can take a while to load, even though it doesn't take long to generate.)

# After playing with this idea, I thought it would make some nice jewellery:

# ![image label](IMAGEFOLDER/ca-jewellery-like.png)

# This is rule 150. I mapped the array to a semicircle and drew it twice.

# ### And it moves

# Many of the rules have limited career paths. Some fizzle out very quickly, others settle down into a stable if repetitive life style. There are a few that continue to make patterns as the number of generations heads off into the thousands.

# The unpredictable high achiever of the 1D cellular automata world is Rule 110. The rule itself is so simple, it could have been described in a Shakespeare play:

# > [Enter HAMLET, stage left.]
# >
# > For Zeros become Ones at all positions,
# >
# > Where the value to the right is One.
# >
# > Yet Ones are changed to Zeros where'er
# >
# > The values to left and right are One.

# And yet this simple rule has been the subject of an astonishing amount of analysis, such as [this in-depth paper by Martinez et al](http://uncomp.uwe.ac.uk/genaro/Papers/Papers_on_CA_files/MARTINEZ.pdf){:target="_blank"} published in the wonderfully-named International Journal of Unconventional Computing. [Matthew Cook's famous paper](https://arxiv.org/pdf/0906.3248v1.pdf){:target="_blank"} proved that Rule 110 is capable of emulating the activity of a Turing machine.

# I understood very little of those wonderful papers, but they did make me want to at least see Rule 110 in action, just in case I could spot all those *cyclic tags*, *meta-gliders*, and *pseudo-solitons*.

function frame(scene, framenumber, ca, cahistory, sidelength;
            smoothscrolling=10)
    background(colorant"navy")
    fontsize(12)
    sethue(colorant"azure")
    text(string(ca.generation), boxtopright(BoundingBox()) + (-30, 20),
        halign=:right)
    setline(0.1)

    ## get in position for the first row
    translate(boxmiddleleft(BoundingBox()))
    rotate(-π/2)

    ## for smooth scrolling
    translate(0, -(mod1(framenumber, smoothscrolling))
        * sidelength/smoothscrolling)

    lng = length(ca.cells)
    for gen in cahistory
        for (n, cell) in enumerate(gen)
            if cell == true
                pt = Point(-(lng / 2) * sidelength + (n * sidelength), 0)
                box(pt, sidelength - 2, sidelength - 2, :fill)
            end
        end
        translate(Point(0, sidelength))
    end
    ## "beautiful buttery scrolling"
    if framenumber % smoothscrolling == 0
        ## drop oldest, add a new generation
        popfirst!(cahistory)
        nextgeneration(ca)
        push!(cahistory, ca.cells)
    end
end

function makeanimation(rule, filename)
    width, height = (1920, 1080)
    sidelength = 6
    cellularmovie = Movie(width, height, "cellularmovie")
    ca = CA(rule, convert(Int, height÷sidelength))
    cahistory = []
    ## initial
    for _ in 1:width÷sidelength
        push!(cahistory, ca.cells)
        nextgeneration(ca)
    end
    animate(cellularmovie,
        [Scene(cellularmovie, (s, f) ->
            frame(s, f, ca, cahistory, sidelength, smoothscrolling=4),
                1:500)],
    pathname="$(filename)")
end

makeanimation(110, "images/automata/animated-cellular-automaton.gif")

# ![image label](IMAGEFOLDER/animated-cellular-automaton.gif)

# If each frame moved the history 'window' by one generation, you'd get a jerky animation. The smooth-scrolling used here shifts the contents by a few pixels in each frame but changes the contents less often. This requires more frames than generations, so 300 frames with a `smoothscrolling` value of 10 shows only 30 new generations after the initial bunch.

# This GIF does have a few problems; the conflicting demands of file size, image size, image quality, and scroll speed work against each other. There's a slight flickering, due I think to the rounding or aliasing in the GIF. You could make the GIF as long as you want, subject of course to the maximum size of animated GIF that you're prepared to handle. (And yes, this page does take a long time to load. Sorry.)

# Videos also have their problems with this kind of content, because the downsampling commonly applied can affect the detail. To see what I mean, try watching [Rule 110: The Movie](https://www.youtube.com/watch?v=iKpIuCmATYE&feature=youtu.be){:target="_blank"} (on YouTube). (This is my submission for YouTube's *Most Boring Video 2018* competition.)

# ### Solid work

# If you want to do more than just look at pictures on the screen—perhaps you're building your own railway station?—here are some tips for exporting the graphics to your favourite 3D modeller.

# First, place a box or shape around the outside of the design, and use the `:path` action rather than `:fill` or `:stroke`. Then, to draw what are now the 'holes', use `newsubpath()` before drawing each hole, and make sure to draw the holes with reversed paths (using the `reversepath` keyword, for example).

function drawaspath(ca::CA;
        cellwidth=10)
    lng = length(ca.cells)
    for i in 1:lng
        if ca.cells[i] == true
            newsubpath()
            pt = Point(-(lng ÷ 2) * cellwidth + (i * cellwidth), 0)
            poly(box(pt, cellwidth - 2, cellwidth - 2, vertices=true),
                :path, reversepath=true, close=true)
        end
    end
end

@svg begin
    background("black")
    ca = CA(110, 40)
    ca.cells = rand(Bool, length(ca.cells))
    cellwidth = 7
    bxw = boxwidth(BoundingBox() * 0.98)/2
    squircle(O, bxw, bxw, rt=0.1, :path)
    translate(boxtopcenter(BoundingBox()) + (-5, 7))
    for i in 1:40
        translate(0, cellwidth)
        drawaspath(ca, cellwidth=cellwidth)
        nextgeneration(ca)
    end
    sethue("grey80")
    fillpath()
end 300 300 "images/automata/ca-as-path.svg"

# ![image label](IMAGEFOLDER/ca-as-path.svg){: .center-image}

# So this is a single path, and the holes are reversed subpaths. When you import this single path into a 3D modelling program, you *should* be able to extrude it into a solid without problems. Here's one modelled as a thin gold reflective sheet.

# ![image label](IMAGEFOLDER/ca-3d-render-1.png)

# It's hard to resist playing with 3D modellers... I added some fog and a rusty metal effect:

# ![image label](IMAGEFOLDER/ca-3d-render-2.png)

# If you omit the surrounding path, it might still work visually. It would be a challenge to hang all these pieces on a gallery wall:

# ![image label](IMAGEFOLDER/ca-3d-render-3.png)

# These little critters are quite interesting, and there are a [fair number of scientific papers](https://arxiv.org/search/?query=%22cellular+automata%22&searchtype=title&abstracts=show&order=-announced_date_first&size=50){:target="_blank"} with "cellular automata" in the title . I think this must be because, although things don't get much simpler than a cellular automaton, perhaps they might—in some strange way—be similar to how the universe itself works when you zoom in close enough.

# [2018-11-29]

# [^1]: At the time of writing this I couldn't find any working Julia packages that did one-dimensional cellular automata. But if I had found one, I probably wouldn't have written any of the above anyway.

# ![cormullion signing off](http://steampiano.net/cormullionknot.gif?automata){: .center-image}

                                                                                #src
using Literate                                                                              #src
# preprocess for notebooks                                                                  #src
function setimagefolder(content)                                                            #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                             #src
    return content                                                                          #src
end                                                                                #src

# for Jupyter notebook, put images in subfolder                                             #src
IMAGEFOLDER = "images/automata"                                                             #src
Literate.notebook("source/automata.jl", "notebooks", preprocess = setimagefolder)           #src

# for Markdown/Jekyll notebook, put images in "/images"                                     #src
IMAGEFOLDER = "/images/automata"                                                            #src
Literate.markdown("source/automata.jl", ".", name="_posts/2018-11-29-automata",             #src
 preprocess = setimagefolder,                                                               #src
 codefence = "{% highlight julia %}" => "{% endhighlight julia %}",                         #src
 documenter=false)                                                                          #src
                                                                                #src


text:
module SampleListTest

using Test
using ReactiveMP
using Random
using Distributions
using LinearAlgebra
using StaticArrays
using StableRNGs

import ReactiveMP: deep_eltype, get_samples, get_weights, sample_list_zero_element
import ReactiveMP: get_meta, is_meta_present, get_unnormalised_weights, get_entropy, get_logproposal, get_logintegrand
import ReactiveMP: call_logproposal, call_logintegrand
import ReactiveMP: transform_samples, transform_weights!
import ReactiveMP: approximate_prod_with_sample_list

@testset "SampleList" begin

    @testset "Internal functions" begin
        @test sample_list_zero_element(SampleList([ 1.0, 1.0 ])) === 0.0
        @test sample_list_zero_element(SampleList([ [ 1.0, 1.0 ], [ 1.0, 1.0 ] ])) == [ 0.0, 0.0 ]
        @test sample_list_zero_element(SampleList([ [ 1.0 1.0; 1.0 1.0 ], [ 1.0 1.0; 1.0 1.0 ] ])) == [ 0.0 0.0; 0.0 0.0 ]
        @test sample_list_zero_element(SampleList([ [ 1.0 1.0 1.0; 1.0 1.0 1.0 ], [ 1.0 1.0 1.0; 1.0 1.0 1.0 ] ])) == [ 0.0 0.0 0.0; 0.0 0.0 0.0 ]
        @test sample_list_zero_element(SampleList([ [ 1.0; 1.0], [ 1.0; 1.0 ] ])) == [ 0.0; 0.0 ]
    end

    @testset "Constructor" begin

        rng = StableRNG(1234)

        for N in [ 5, 10, 100 ], type in [ Float64, Float32, BigFloat ]
            scalar_samples    = rand(rng, type, N)
            scalar_samplelist = SampleList(scalar_samples)

            @test collect(get_samples(scalar_samplelist)) == first.(collect(scalar_samplelist))
            @test collect(get_weights(scalar_samplelist)) == fill(one(type) / N, N)
            @test deep_eltype(scalar_samplelist)              === type
            @test eltype(scalar_samplelist)                   === Tuple{type, type}
            @test eltype(get_weights(scalar_samplelist))      === type
            @test variate_form(scalar_samplelist) === Univariate
            @test is_meta_present(scalar_samplelist) === false

            scalar_weights    = rand(rng, type, N)

            @test_throws AssertionError SampleList(scalar_samples, scalar_weights)
            scalar_weights ./= sum(scalar_weights)
            scalar_samplelist = SampleList(scalar_samples, scalar_weights)

            @test collect(get_samples(scalar_samplelist)) == first.(collect(scalar_samplelist))
            @test collect(get_weights(scalar_samplelist)) == scalar_weights
            @test deep_eltype(scalar_samplelist)              === type
            @test eltype(scalar_samplelist)                   === Tuple{type, type}
            @test eltype(get_weights(scalar_samplelist))      === type
            @test variate_form(scalar_samplelist) === Univariate
            @test is_meta_present(scalar_samplelist) === false

            vector_samples    = [ rand(rng, type, 2) for _ in 1:N ]
            vector_samplelist = SampleList(vector_samples)

            @test collect(get_samples(vector_samplelist)) == first.(collect(vector_samplelist))
            @test collect(get_weights(vector_samplelist)) == fill(one(type) / N, N)
            @test deep_eltype(vector_samplelist)              === type
            @test eltype(vector_samplelist)                   === Tuple{ SVector{ 2, type }, type }
            @test eltype(get_weights(vector_samplelist))      === type
            @test variate_form(vector_samplelist) === Multivariate
            @test is_meta_present(vector_samplelist) === false

            vector_weights    = rand(rng, type, N)
            @test_throws AssertionError SampleList(vector_samples, vector_weights)
            vector_weights ./= sum(vector_weights)
            vector_samplelist = SampleList(vector_samples, vector_weights)

            @test collect(get_samples(vector_samplelist)) == first.(collect(vector_samplelist))
            @test collect(get_weights(vector_samplelist)) == vector_weights
            @test deep_eltype(vector_samplelist)              === type
            @test eltype(vector_samplelist)                   === Tuple{ SVector{ 2, type }, type }
            @test eltype(get_weights(vector_samplelist))      === type
            @test variate_form(vector_samplelist) === Multivariate
            @test is_meta_present(vector_samplelist) === false

            matrix_samples    = [ rand(rng, type, 2, 2) for _ in 1:N ]
            matrix_samplelist = SampleList(matrix_samples)

            @test collect(get_samples(matrix_samplelist)) == first.(collect(matrix_samplelist))
            @test collect(get_weights(matrix_samplelist)) == fill(one(type) / N, N)
            @test deep_eltype(matrix_samplelist)              === type
            @test eltype(matrix_samplelist)                   === Tuple{ SMatrix{ 2, 2, type, 4 }, type }
            @test eltype(get_weights(matrix_samplelist))      === type
            @test variate_form(matrix_samplelist) === Matrixvariate
            @test is_meta_present(matrix_samplelist) === false

            matrix_weights    = rand(rng, type, N)
            @test_throws AssertionError SampleList(matrix_samples, matrix_weights)
            matrix_weights ./= sum(matrix_weights)
            matrix_samplelist = SampleList(matrix_samples, matrix_weights)

            @test collect(get_samples(matrix_samplelist)) == first.(collect(matrix_samplelist))
            @test collect(get_weights(matrix_samplelist)) == matrix_weights
            @test deep_eltype(matrix_samplelist)              === type
            @test eltype(matrix_samplelist)                   === Tuple{ SMatrix{ 2, 2, type, 4 }, type }
            @test eltype(get_weights(matrix_samplelist))      === type
            @test variate_form(matrix_samplelist) === Matrixvariate
            @test is_meta_present(matrix_samplelist) === false
        end

        @test_throws AssertionError SampleList(rand(10), rand(5))
        @test_throws AssertionError SampleList(rand(5), rand(10))
        @test_throws AssertionError SampleList(rand(5), [ -1 for _ in 1:5 ])

        @test_throws AssertionError SampleList([ rand(10) for _ in 1:10 ], rand(5))
        @test_throws AssertionError SampleList([ rand(5) for _ in 1:5 ], rand(10))
        @test_throws AssertionError SampleList([ rand(5) for _ in 1:5 ], [ -1 for _ in 1:5 ])

        @test_throws AssertionError SampleList([ rand(10, 10) for _ in 1:10 ], rand(5))
        @test_throws AssertionError SampleList([ rand(5, 5) for _ in 1:5 ], rand(10))
        @test_throws AssertionError SampleList([ rand(5, 5) for _ in 1:5 ], [ -1 for _ in 1:5 ])

    end

    @testset "Statistics" begin

        rng = StableRNG(42)

        # All
        for N in [ 5, 10, 100 ]

            scalar_samples = rand(rng, N)
            scalar_weights = rand(rng, N)
            scalar_weights ./= sum(scalar_weights)
            scalar_samplelist = SampleList(scalar_samples, scalar_weights)

            # Checking i = 1:2 that cache is not corrupted
            for i in 1:2
                @test mean(scalar_samplelist)    ≈ sum(scalar_weights .* scalar_samples)
                @test logmean(scalar_samplelist) ≈ sum(scalar_weights .* log.(scalar_samples))
                @test meanlogmean(scalar_samplelist) ≈ sum(scalar_weights .* scalar_samples .* log.(scalar_samples))
            end

            vector_samples = [ rand(rng, 2) for _ in 1:N ]
            vector_weights = rand(rng, N)
            vector_weights ./= sum(vector_weights)
            vector_samplelist = SampleList(vector_samples, vector_weights)

            # Checking i = 1:2 that cache is not corrupted
            for i in 1:2
                @test mean(vector_samplelist)    ≈ sum(vector_weights .* vector_samples)
                @test logmean(vector_samplelist) ≈ sum(vector_weights .* map(e -> log.(e), (vector_samples)))
                @test meanlogmean(vector_samplelist) ≈ sum(vector_weights .* map(e -> e .* log.(e), (vector_samples)))
            end

            matrix_samples = [ rand(rng, 2, 2) for _ in 1:N ]
            matrix_weights = rand(rng, N)
            matrix_weights ./= sum(matrix_weights)
            matrix_samplelist = SampleList(matrix_samples, matrix_weights)

            # Checking i = 1:2 that cache is not corrupted
            for i in 1:2
                @test mean(matrix_samplelist) ≈ sum(matrix_weights .* matrix_samples)
                @test logmean(matrix_samplelist) ≈ sum(matrix_weights .* map(e -> log.(e), matrix_samples))
                @test meanlogmean(matrix_samplelist) ≈ sum(matrix_weights .* map(e -> e .* log.(e), matrix_samples))
            end

        end

        # Logmean and meanlogmean
        uni_distribution = Gamma(rand(rng) + 1, rand(rng) + 2)
        uni_samples      = rand(rng, uni_distribution, 20_000)
        uni_sample_list  = SampleList(uni_samples)

        uni_distribution2 = NormalMeanVariance(rand(rng) + 1, rand(rng) + 2)
        uni_samples2      = rand(rng, uni_distribution2, 20_000)
        uni_sample_list2  = SampleList(uni_samples2)

        m = rand(rng, 3)
        r = rand(rng, 3)
        Σ = I + 2r*r'
        mv_distribution = MvNormal(m, Σ)
        mv_samples      = [ rand(rng, mv_distribution) for _ in 1:20_000 ]
        mv_sample_list  = SampleList(mv_samples)

        W1 = rand(rng, 3, 4)
        r2 = rand(rng, 3)
        W2 = I + 2r2*r2'
        r3 = rand(rng, 4)
        W3 = I + 2r3*r3'
        mxv_distribution = MatrixNormal(W1, W2, W3)
        mxv_samples      = [ rand(rng, mxv_distribution) for _ in 1:20_000 ]
        mxv_sample_list  = SampleList(mxv_samples)

        # Checking i = 1:2 that cache is not corrupted
        for i in 1:2
            @test isapprox(mean(uni_sample_list), mean(uni_distribution), atol = 0.1)
            @test isapprox(mean(mv_sample_list), mean(mv_distribution), atol = 0.1)
            @test isapprox(mean(mxv_sample_list), mean(mxv_distribution), atol = 1.0)

            @test isapprox(weightedmean(uni_sample_list2), weightedmean(uni_distribution2), atol = 0.1)
            @test all(isapprox.(mean_var(uni_sample_list2), mean_var(uni_distribution2), atol = 0.1))
            @test all(isapprox.(mean_cov(uni_sample_list2), mean_cov(uni_distribution2), atol = 0.1))
            @test all(isapprox.(mean_precision(uni_sample_list2), mean_precision(uni_distribution2), atol = 0.1))
            @test all(isapprox.(mean_invcov(uni_sample_list2), mean_invcov(uni_distribution2), atol = 0.1))
            @test all(isapprox.(weightedmean_cov(uni_sample_list2), weightedmean_cov(uni_distribution2), atol = 0.1))
            @test all(isapprox.(weightedmean_invcov(uni_sample_list2), weightedmean_invcov(uni_distribution2), atol = 0.1))
            @test all(isapprox.(weightedmean_precision(uni_sample_list2), weightedmean_precision(uni_distribution2), atol = 0.1))

            @test isapprox(var(uni_sample_list), var(uni_distribution), atol = 0.5)
            @test isapprox(cov(uni_sample_list), var(uni_distribution), atol = 0.5)
            @test isapprox(var(mv_sample_list), var(mv_distribution), atol = 0.1)
            @test isapprox(cov(mv_sample_list), cov(mv_distribution), atol = 0.1)
            @test isapprox(var(mxv_sample_list), var(mxv_distribution), atol = 0.1)
            @test isapprox(cov(mxv_sample_list), cov(mxv_distribution), atol = 1.0)

            @test isapprox(std(uni_sample_list), std(uni_distribution), atol = 0.2)
            @test isapprox(std(mv_sample_list), cholsqrt(cov(mv_distribution)), atol = 0.2)

            @test isapprox(invcov(uni_sample_list), inv(var(uni_distribution)), atol = 0.2)
            @test isapprox(invcov(mv_sample_list), cholinv(cov(mv_distribution)), atol = 0.2)
        end

        # TODO logmean for matrix variate distribution?

        mv_distribution = Dirichlet(rand(rng, 3))
        mv_samples      = [ rand(rng, mv_distribution) for _ in 1:20_000 ]
        mv_sample_list  = SampleList(mv_samples)

        @test isapprox(logmean(uni_sample_list), logmean(uni_distribution); atol = 0.02)
        @test isapprox(logmean(mv_sample_list), logmean(mv_distribution); atol = 0.5)

        # TODO meanlogmean for multivariate and matrix variate distribution?

        @test isapprox(meanlogmean(uni_sample_list), meanlogmean(uni_distribution); atol = 0.25)

        r4 = rand(rng, 5)
        W4 = I + 2r4*r4'

        mxv_distribution = Wishart(5, W4)
        mxv_samples      = [ rand(rng, mxv_distribution) for _ in 1:20_000 ]
        mxv_sample_list  = SampleList(mxv_samples)

        # Checking i = 1:2 that cache is not corrupted
        for i in 1:2
            @test isapprox(mean(mxv_sample_list), mean(mxv_distribution), atol = 1.0)
            @test isapprox(var(mxv_sample_list), var(mxv_distribution), atol = 1.0)
            @test isapprox(cov(mxv_sample_list), cov(mxv_distribution), atol = 5.0)
        end
    end

    @testset "vague" begin

        @test variate_form(vague(SampleList))    === Univariate
        @test variate_form(vague(SampleList, 2)) === Multivariate
        @test variate_form(vague(SampleList, 2, 2)) === Matrixvariate
        @test variate_form(vague(SampleList, (3, 4))) === Matrixvariate

        @test ndims(vague(SampleList))    === 1
        @test ndims(vague(SampleList, 2)) === 2
        @test ndims(vague(SampleList, 2, 2)) === (2, 2)
        @test ndims(vague(SampleList, (3, 4))) === (3, 4)

        for nsamples in [ 10, 100, 1000 ]
            @test variate_form(vague(SampleList; nsamples = nsamples))    === Univariate
            @test variate_form(vague(SampleList, 2; nsamples = nsamples)) === Multivariate
            @test variate_form(vague(SampleList, 2, 2; nsamples = nsamples)) === Matrixvariate
            @test variate_form(vague(SampleList, (3, 4); nsamples = nsamples)) === Matrixvariate

            @test length(vague(SampleList; nsamples = nsamples))    === nsamples
            @test length(vague(SampleList, 2; nsamples = nsamples)) === nsamples
            @test length(vague(SampleList, 2, 2; nsamples = nsamples)) === nsamples
            @test length(vague(SampleList, (3, 4); nsamples = nsamples)) === nsamples

            @test ndims(vague(SampleList; nsamples = nsamples))    === 1
            @test ndims(vague(SampleList, 2; nsamples = nsamples)) === 2
            @test ndims(vague(SampleList, 2, 2; nsamples = nsamples)) === (2, 2)
            @test ndims(vague(SampleList, (3, 4); nsamples = nsamples)) === (3, 4)

            @test size(vague(SampleList; nsamples = nsamples))    === (nsamples, )
            @test size(vague(SampleList, 2; nsamples = nsamples)) === (nsamples, )
            @test size(vague(SampleList, 2, 2; nsamples = nsamples)) === (nsamples, )
            @test size(vague(SampleList, (3, 4); nsamples = nsamples)) === (nsamples, )
        end

    end

    @testset "SampleListMeta" begin
        @test_throws ErrorException get_meta(SampleList([ 0.1 ]))

        rng = StableRNG(1234)

        for uweights in [ rand(rng, 2), rand(rng, 2) ],
            entropy in [ -75.6, 234.2 ],
            logproposal in [ NormalMeanVariance(-3.0, 1.0), Gamma(2.0, 4.0), (x) -> x + 1.0 ],
            logintegrand in [ NormalMeanPrecision(-1.0, 2.0), Beta(3.0, 4.0), (x) -> log(abs(x)) ]

            meta = SampleListMeta(uweights, entropy, logproposal, logintegrand)
            sl   = SampleList([ 0.1, 0.1 ], [ 0.5, 0.5 ], meta)

            @test get_meta(sl) === meta
            @test get_unnormalised_weights(sl) == uweights
            @test get_entropy(sl) == entropy
            @test get_logproposal(sl) == logproposal
            @test get_logintegrand(sl) == logintegrand
            @test is_meta_present(sl) === true

            some_random_numbers = rand(rng, 100)

            __call(dist::Distribution, x) = logpdf(dist, x)
            __call(dist::Function, x)     = dist(x)

            @test map(e -> call_logproposal(sl, e), some_random_numbers) == map(e -> __call(logproposal, e), some_random_numbers)
            @test map(e -> call_logintegrand(sl, e), some_random_numbers) == map(e -> __call(logintegrand, e), some_random_numbers)
        end

    end

    @testset "Iteration utilities" begin
        rng = StableRNG(42)

        uni_distribution = Uniform(-10rand(rng), 10rand(rng))

        μ  = rand(rng, 3)
        L1 = rand(rng, 3, 3)
        Σ = L1' * L1
        mv_distribution = MvNormal(μ, Σ)

        L2 = rand(rng, 3, 3)
        W  = L2' * L2
        mvx_distribution = Wishart(3, W)

        # Entity to entity
        f1(e) = e .+ 1
        f2(e) = exp.(e)

        # Entity to Number
        f3(e::Number) = e + 1
        f3(e::AbstractVector) = norm(e .+ 1)
        f3(e::AbstractMatrix) = det(e .+ 1)

        # Entity to Vector
        f4(e::Number)         = @SVector [ e, e ]
        f4(e::AbstractVector) = reverse(e)
        f4(e::AbstractMatrix) = diag(e)

        # Entity to Matrix
        f5(e::Number)         = @SMatrix [ e + 1 e; e e + 1 ]
        f5(e::AbstractVector) = SMatrix{length(e), length(e)}(Diagonal(ones(length(e))))
        f5(e::AbstractMatrix) = [ e[1, 1] + 1 e[1, 2]; e[2, 1] e[2, 2] + 1 ]

        for N in (500, 1000, 5_000)
            for distribution in (uni_distribution, mv_distribution, mvx_distribution)
                samples = [ rand(rng, distribution) for _ in 1:N ]
                weights = ones(N) ./ N
                samplelist = SampleList(samples, weights)

                @test collect(samplelist) == collect(zip(samples, weights))
                @test map(i -> samplelist[i], 1:N) == collect(zip(samples, weights))

                for f in (f1, f2, f3, f4, f5)
                    @test all(map(e -> all(e[1] .≈ e[2]), zip(collect(transform_samples(f, samplelist)), collect(zip(map(f, samples), weights)))))
                    @test all(map(e -> all(e[1] .≈ e[2]), zip(map(i -> (f(samplelist[i][1]), samplelist[i][2]), 1:N), collect(zip(f.(samples), weights)))))
                end

                iter = N:-1:1
                index = 0

                old_weights = copy(weights)

                transform_weights!(w -> w * iter[ index += 1 ], samplelist)

                newweights = map(prod, zip(old_weights, iter))
                newweights ./= sum(newweights)

                @test get_weights(samplelist) ≈ newweights

            end
        end

    end

    @testset "prod approximations" begin

        rng = StableRNG(1234)

        posdefm(rng, s) = begin L = rand(rng, s, s); L'*L end

        sizes  = [ 2_500, 5_000, 10_000 ]
        inputs = [
            (x = NormalMeanPrecision(3.0, 7.0), y = NormalMeanVariance(-4.0, 6.0), mean_tol = [ 1e-1, 1e-1, 1e-1 ], cov_tol = [ 1e-1, 1e-1, 1e-1 ], entropy_tol = [ 1e-1, 1e-1, 1e-1 ]),
            (x = NormalMeanVariance(3.0, 7.0), y = NormalWeightedMeanPrecision(4.0, 6.0), mean_tol = [ 1e-1, 1e-1, 1e-1 ], cov_tol = [ 1e-1, 1e-1, 1e-1 ], entropy_tol = [ 1e-1, 1e-1, 1e-1 ]),
            (x = GammaShapeRate(3.0, 7.0), y = GammaShapeScale(4.0, 6.0), mean_tol = [ 1e-1, 1e-1, 1e-1 ], cov_tol = [ 1e-1, 1e-1, 1e-1 ], entropy_tol = [ 3e-1, 3e-1, 3e-1 ]),
            (x = MvNormalMeanCovariance(10rand(rng, 4), posdefm(rng, 4)), y = MvNormalMeanPrecision(10rand(rng, 4), posdefm(rng, 4)), mean_tol = [ 2e-1, 2e-1, 2e-1 ], cov_tol = [ 6e-1, 6e-1, 6e-1 ], entropy_tol = [ 4e-1, 4e-1, 4e-1 ]),
            (x = Wishart(10.0, posdefm(rng, 3)), y = Wishart(5.0, posdefm(rng, 3)), mean_tol = [ 7e-1, 7e-1, 7e-1 ], cov_tol = [ 5e-1, 5e-1, 5e-1 ], entropy_tol = [ 2e-1, 2e-1, 2e-1 ]),
        ]

        for (i, N) in enumerate(sizes)
            for input in inputs
                analytical = prod(ProdAnalytical(), input[:x], input[:y])
                approximation = approximate_prod_with_sample_list(rng, input[:x], input[:y], N)

                @test is_meta_present(approximation) === true
                @test length(approximation) === N

                μᵣ, Σᵣ = mean_cov(analytical)
                μₐ, Σₐ = mean_cov(approximation)

                @test norm(μᵣ .- μₐ) < input[:mean_tol][i]
                @test norm(Σᵣ .- Σₐ) < input[:cov_tol][i]
                @test abs(entropy(analytical) - entropy(approximation)) < input[:entropy_tol][i]

                # Second order approximation here
                if (variate_form(input[:x]) === Univariate)

                    analytical2 = prod(ProdAnalytical(), analytical, input[:x])
                    approximation2 = approximate_prod_with_sample_list(rng, input[:x], approximation, N)

                    @test is_meta_present(approximation2) === true
                    @test length(approximation2) === N

                    μᵣ, Σᵣ = mean_cov(analytical2)
                    μₐ, Σₐ = mean_cov(approximation2)

                    if !(abs(entropy(analytical2) - entropy(approximation2)) < input[:entropy_tol][i])
                        @show analytical2
                    end

                    @test norm(μᵣ .- μₐ) < input[:mean_tol][i]
                    @test norm(Σᵣ .- Σₐ) < input[:cov_tol][i]
                    @test abs(entropy(analytical2) - entropy(approximation2)) < input[:entropy_tol][i]
                end
            end
        end

    end

end

end

ext: jl


text:
export area

"""
    area(polygon::Polygon, radius::Float64=1.0)

Return the area of the spherical polygon. Without given radius a unit sphere is assumed.

Source: en.wikipedia.org/wiki/Spherical_trigonometry # Area and spherical excess

The method has been altered to allow for the usage of an internal point given with the Polygon.
This method only works with polygons that fit within half a sphere as it uses great circle distances between
subsequent points.
"""
function area(polygon::Polygon, radius::Float64=1.0)
    point₁ = polygon.points[1]
    area_polygon_unit_sphere = 0.0
    for point₂ in polygon.points[2:end]
        n_intersectionsᵢ = length(intersection_points(Arc(polygon.inside_point,
        midpoint(point₁, point₂)), polygon))
        isodd(n_intersectionsᵢ) ? area_polygon_unit_sphere +=
        area(polygon.inside_point, point₁, point₂) : area_polygon_unit_sphere -=
        area(polygon.inside_point, point₁, point₂)

        point₁ = point₂
    end
    return area_polygon_unit_sphere * radius^2
end

"""
    area(angular_distance₁₂::Float64, angular_distance₂₃::Float64, angular_distance₁₃::Float64,
    radius::Float64=1.0)

Return the area of the spherical triangle 123 based on the angular distances [deg] 12-23-13. Without given
radius a unit sphere is assumed.

Source: mathworld.wolfram.com/SphericalTriangle.html
"""
function area(angular_distance₁₂::Float64, angular_distance₂₃::Float64, angular_distance₁₃::Float64,
    radius::Float64=1.0)
    return deg2rad(spherical_excess(angular_distance₁₂, angular_distance₁₃, angular_distance₂₃))*radius^2
end

"""
    area(point₁::Point, point₂::Point, point₃::Point, radius::Float64=1.0)

Return the area of the spherical triangle 123 based on the points point₁ - point₂ - point₃. Without given radius
a unit sphere is assumed.

Source: mathworld.wolfram.com/SphericalTriangle.html
"""
area(point₁::Point, point₂::Point, point₃::Point, radius::Float64=1.0) =
area(angular_distance(point₁, point₂), angular_distance(point₂, point₃), angular_distance(point₃, point₁),
radius)


text:
using LAP_julia: filt_onebyone!

img_sizes = [200, 400, 700]
point_counts = [30, 50, 100, 200]
fhs = [5, 10, 30]

filt_onebyone!()


text:
function line_score(delta2::Int, tl::Int, br::Int)::Float64
    (tl += 1) == br && return 0.0
    delta = sqrt(delta2)
    abs(delta - br) - abs(delta - tl)
end

function grow_corner_diags(n::Int, leq=false)::BitMatrix
    ds = Matrix{Int}(undef, n+1, n+1)
    ds[:, 1] = ds[1, :] = 0:n
    diags = BitMatrix(undef, n, n)
    pred = leq ? >=(0) : >(0)
    for i in CartesianIndices(diags)
        tl = ds[i]
        l = ds[i+onex]
        t = ds[i+oney]
        br = min(l, t) + 1
        if pred(line_score(sum(i.I.^2), tl, br))
            ds[i + onexy] = tl + 1
            diags[i] = true
        else
            ds[i + onexy] = br
            diags[i] = false
        end
    end
    diags
end

gamma_weight(opp, adj, x::Int, y::Int, ::Val{:unweighted})::Float64 = 1.0
gamma_weight(opp, adj, x::Int, y::Int, ::Val{:angular})::Float64 = atan(opp, adj)
gamma_weight(opp, adj, x::Int, y::Int, ::Val{:minmax})::Float64 = /(minmax(x, y)...)^-1
gamma_weight(opp, adj, x::Int, y::Int, ::Nothing=nothing)::Float64 =
    gamma_weight(opp, adj, x, y, Val(:angular)) * gamma_weight(opp, adj, x, y, Val(:minmax))

function gamma_score(v::Vertex, tl, br, W=nothing)::Float64
    vsum = sum(v.I)
    v2 = v.I.^2
    summands = Float64[]

    for i in 1:v[1]-1
        delta2 = i^2 + v2[2]
        s = line_score(delta2, tl[i], br[i])
        s != 0 && push!(summands, gamma_weight(v[2], delta2 - 0.25, i, v[2], W) * s)
    end

    delta2 = sum(v2)
    s = line_score(delta2, tl[v[1]], br[v[1]])
    s != 0 && push!(summands, gamma_weight(vsum - 0.5, 2 * delta2 - vsum, v.I..., W) * s)

    for i in v[1]+1:vsum-1
        j = vsum - i
        delta2 = v2[1] + j^2
        s = line_score(delta2, tl[i], br[i])
        s != 0 && push!(summands, gamma_weight(v[1], delta2 - 0.25, v[1], j, W) * s)
    end

    sum_kbn(summands)
end

const Buffer = Array{Int, 3}

mutable struct State
    diags::BitMatrix
    blocked::BitMatrix
    position::Int
    vertices::Vector{Vertex}
    buffer::Buffer
    children::SubArray{Int, 2, Buffer}
    parents::SubArray{Int, 2, Buffer}
    grandparents::SubArray{Int, 2, Buffer}

    function State(n::Int; blocked=nothing)
        isnothing(blocked) && (blocked = falses(n, n))
        checksquare(blocked) == n ||
            throw(DimensionMismatch("size(blocked) $(size(blocked)) != ($n, $n)"))
        buffer = Buffer(undef, 2n+1, n+2, 3)
        children = view(buffer, 1:2, 1:2, 2) .= [0 1; 1 0]
        parents = view(buffer, 1:1, 1:1, 1) .= 0
        new(falses(n, n), blocked, 0, sizehint!(Vertex[], n), buffer, children, parents)
    end
end

Base.isempty(s::State)::Bool = s.position == size(s.buffer, 1) - 2

function Base.popfirst!(s::State)
    @assert !isempty(s)

    n = size(s.buffer, 2) - 2
    i = s.position += 1

    s.grandparents = i > n ? view(s.parents, :, 2:2n+1-i) : s.parents
    s.parents = s.children
    s.children = children = view(s.buffer, 1:i+2, 1:(i < n ? i+2 : 2n-i), mod1(i+2, 3))

    width = size(s.grandparents, 2)
    if i < n
        children[:, 1] .= 0:size(children, 1)-1
        children[:, width+2] .= size(children, 1)-1:-1:0
        children = view(children, :, 2:width+1)
        v = Vertex(i, 1)
    else
        v = Vertex(n, i-n+1)
    end
    resize!(s.vertices, width) .= (v + Vertex(-1, 1) * j for j in 0:width-1)

    (s.grandparents, s.parents, children)
end

function propagate_distances!(i::Int, parents, children)
    depth = size(parents, 1) - 1
    l = view(parents, :, i)
    t = view(parents, :, i+1)
    br = view(children, :, i)
    br[1] = l[1] + 1
    br[depth+2] = t[depth+1] + 1
    view(br, 2:depth+1) .= min.(view(l, 2:depth+1), view(t, 1:depth)) .+ 1
end

function score!(s::State, grandparents, parents, children; W=nothing)::Vector{Float64}
    scores = Vector{Float64}(undef, length(s.vertices))
    Threads.@threads for i in eachindex(s.vertices)
        tl = view(grandparents, :, i)
        br = propagate_distances!(i, parents, children)
        scores[i] = s.blocked[s.vertices[i]] ? -1 : gamma_score(s.vertices[i], tl, br, W)
    end
    scores
end

function setdiags!(s::State, diag_indices, grandparents, children)
    s.diags[s.vertices[diag_indices]] .= true
    Threads.@threads for i in diag_indices
        children[2:s.position+1, i] .= view(grandparents, 1:s.position, i) .+ 1
    end
end

function step!(s::State, (grandparents, parents, children); W=nothing)
    scores = score!(s, grandparents, parents, children, W=W)
    setdiags!(s, findall(>=(0.0), scores), grandparents, children)
end

function grow_gamma_diags(n::Int; W=nothing, margin=0, blocked=nothing)::BitMatrix
    m = n + margin
    s = State(m, blocked=blocked)
    for _ in 1:2m-1
        step!(s, popfirst!(s), W=W)
    end
    s.diags[(margin+1)*onexy:m*onexy]
end

max_stars_repo_name: moshelooks/eugrid

max_issues_repo_name: moshelooks/eugrid

max_forks_repo_name: moshelooks/eugrid


text:
# This file is a part of SimilaritySearch.jl

export closestpair

"""
    closestpair(idx::AbstractSearchContext; parallel=false, pools=getpools(idx))

Finds the closest pair among all elements in `idx`. If the index `idx` is approximate then pair of points could be also an approximation.

# Arguments:
- `idx`: the search structure that indexes the set of points

# Keyword Arguments:
- `parallel`: If true then the algorithm uses all available threads to compute the closest pair
- `pools`: The pools needed for the index. Only used for special cases, default values should work in most cases. See [`getpools`](@ref) for more information.
"""
function closestpair(idx::AbstractSearchContext; parallel=false, pools=getpools(idx))
    parallel ? parallel_closestpair(idx, Threads.nthreads(), pools) : sequential_closestpair(idx, pools)
end

function search_hint(idx::AbstractSearchContext, i::Integer, pools)
    res = getknnresult(2)
    search(idx, idx[i], res; pools)
    argmin(res) == i ? (argmax(res), maximum(res)) : (argmin(res), minimum(res))
end

function search_hint(G::SearchGraph, i::Integer, pools)
    res = getknnresult(2)
    search(G.search_algo, G, G[i], res, first(G.links[i]), pools)
    argmin(res) == i ? (argmax(res), maximum(res)) : (argmin(res), minimum(res))
end

function parallel_closestpair(idx::AbstractSearchContext, parallel_block, pools)
    n = length(idx)
    parallel_block = min(n, parallel_block)
    B = [(zero(Int32), zero(Int32), typemax(Float32)) for _ in 1:parallel_block]

    for block in Iterators.partition(1:n, parallel_block)
        Threads.@threads for i in 1:length(block)
            @inbounds objID = block[i]
            id_, d_ = search_hint(idx, objID, pools)
            @inbounds if d_ < last(B[i])
                B[i] = (objID, id_, d_)
            end
        end
    end

    _, i = findmin(last, B)
    B[i]
end

function sequential_closestpair(idx::AbstractSearchContext, pools)
    mindist = typemax(Float32)
    I = J = zero(Int32)

    for i in eachindex(idx)
        id_, d_ = search_hint(idx, i, pools)
        if d_ < mindist
            I, J, mindist = i, id_, d_
        end
    end

    (I, J, mindist)
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
using Test
using MusicManipulations: scale_identification

tester = getnotes(readMIDIFile(joinpath(@__DIR__, "recording_uwe_2.mid")),1)
@test scale_identification(tester) == "A♯ Major/G minor"


text:
module Accounting

using Compat
import Compat: String

import Currencies: currency
using Currencies
using DataStructures
import Base.push!

include("debitcredit.jl")
include("accounts.jl")
include("entries.jl")
include("ledger.jl")
include("reports.jl")

export Split, Entry, Ledger
export Asset, Liability, Equity, Revenue, Expense, Trading
export asset!, liability!, equity!, revenue!, expense!, trading!
export debit, credit, balances, currency
export register!, transfer!

end  # module


text:
module WiagDataSetup

using MySQL
using Infiltrator
using DataFrames

dbwiag = nothing

function setDBWIAG(;pwd = missing, host = "127.0.0.1", user = "wiag", db = "wiag")
    global dbwiag
    if ismissing(pwd)
        println("Passwort für User ", user)
        pwd = readline()
    end
    dbwiag = DBInterface.connect(MySQL.Connection, host, user, pwd, db = db)
end

"""
    updatenamevariant()

obsolete: see namelookup
"""
function updatenamevariant(fieldsrc::AbstractString, tablename::AbstractString)::Int
    dbwiag = DBInterface.connect(MySQL.Connection, "localhost", "wiag", "Wogen&Wellen", db="wiag");

    DBInterface.execute(dbwiag, "DELETE FROM " * tablename);

    # Do not use the same DB connection with an open cursor to insert data
    dfsrc = DBInterface.execute(dbwiag,
                                "SELECT wiagid, " * fieldsrc
                                * " FROM person") |> DataFrame;

    tblid = 1
    for row in eachrow(dfsrc)
        id, fns = row
        # println("id: ", id)
        # println(fns)

        if ismissing(fns) || fns == "" continue end
        for nv in split(fns, r",|;")
            insertstmt = ("INSERT INTO " * tablename * " VALUES ("
                          * string(tblid) * ","
                          * string(id) * ","
                          * "'" * strip(nv) * "')")
            # println(insertstmt)
            DBInterface.execute(dbwiag, insertstmt)
            tblid += 1;
        end
    end
    return tblid;
end

"""
    fillera(tblera::AbstractString, tblperson::AbstractString, tbloffice::AbstractString, colnameid::AbstractString, checkisonline = true)::Int

Compute earliest and latest date for each person, identified by `colnameid` and `colnameidinoffice`.

"""
function fillera(tblera::AbstractString,
                 tblperson::AbstractString,
                 tbloffice::AbstractString,
                 colid = "id",
                 colidinoffice = "id_person";
                 checkisonline = false)::Int
    global dbwiag
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    msg = 1000

    DBInterface.execute(dbwiag, "DELETE FROM " * tblera);

    sqlselect = "SELECT " * colid * " as idperson, date_birth, date_death " *
        " FROM " * tblperson

    if checkisonline
        sqlselect *= " WHERE status = 'online'"
    end

    dfperson = DBInterface.execute(dbwiag, sqlselect) |> DataFrame;

    tblid = 0;

    # get office data
    sqlselect = "SELECT " * colidinoffice * ", date_start, date_end " * " FROM " * tbloffice
    dfoffice = DBInterface.execute(dbwiag, sqlselect) |> DataFrame

    insstmt = DBInterface.prepare(dbwiag, "INSERT INTO " * tblera * " VALUES (?, ?, ?)")

    for row in eachrow(dfperson)
        idperson, datebirth, datedeath = row[[:idperson, :date_birth, :date_death]]

        dfioffice = filter([Symbol(colidinoffice)] => isequal(idperson), dfoffice)

        erastart, eraend = extremaera(datebirth, datedeath, dfioffice)

        # it is a bit slower to call the database in each step, but needs less code
        DBInterface.execute(insstmt, [idperson, erastart, eraend]);
        tblid += 1
        if tblid % msg == 0
            @info tblid
        end
    end

    DBInterface.close!(insstmt)

    return tblid
end

"""
    updatecnera(tblera::AbstractString, tblperson::AbstractString, tblonline::AbstractString, tbloffice::AbstractString, colnameid::AbstractString)::Int

Compute earliest and latest date for each person, identified by `colnameid` and `colnameidinoffice`.

"""
function updatecnera(tblera::AbstractString,
                     tblonline::AbstractString,
                     tblperson::AbstractString,
                     tbloffice::AbstractString,
                     colid::AbstractString,
                     colidinoffice = "id_person")::Int
    # TODO take into account offices in bishop's database
    global dbwiag
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    msg = 1000

    @info "update '" * tblera * "'"

    sqlperson = "SELECT po.id as idonline, po." * colid * " as idperson " * ", date_birth, date_death " *
        " FROM " * tblonline * " AS po " *
        " JOIN " * tblperson * " AS dh ON dh.id = po." * colid

    dfperson = DBInterface.execute(dbwiag, sqlperson) |> DataFrame;

    tblid = 0;

    # get office data
    sqloffice = "SELECT " * colidinoffice * ", date_start, date_end " * " FROM " * tbloffice
    dfoffice = DBInterface.execute(dbwiag, sqloffice) |> DataFrame
    # @infiltrate

    sqlstart = "UPDATE " * tblera * " SET era_start=?" *
        " WHERE id_online=? AND (era_start is NULL OR era_start > ?)"
    updstmtstart = DBInterface.prepare(dbwiag, sqlstart)

    sqlend = "UPDATE " * tblera * " SET era_end=?" *
        " WHERE id_online=? AND (era_end is NULL OR era_end < ?)"
    updstmtend = DBInterface.prepare(dbwiag, sqlend)

    for row in eachrow(dfperson)

        idonline, idperson, datebirth, datedeath = row[[:idonline, :idperson, :date_birth, :date_death]]

        dfioffice = filter([Symbol(colidinoffice)] => isequal(idperson), dfoffice)

        erastart, eraend = eraextrema(datebirth, datedeath, dfioffice)

        # it is a bit slower to call the database in each step, but needs less code
        if !ismissing(erastart)
            DBInterface.execute(updstmtstart, [erastart, idonline, erastart])
        end
        if !ismissing(eraend)
            DBInterface.execute(updstmtend, [eraend, idonline, eraend])
        end

        tblid += 1
        if tblid % msg == 0
            @info tblid
        end
    end

    DBInterface.close!(updstmtstart)
    DBInterface.close!(updstmtend)

    return tblid
end

"""
    eraextrema(datebirth, datedeath, dfoffice)

find extrema in [`datebirth`, `datedeath`] and `dfoffice`
"""
function eraextrema(datebirth, datedeath, dfoffice)

    erastart = Inf
    eraend = -Inf

    vcand = parsemaybe(datebirth, :lower)
    if !ismissing(vcand) erastart = vcand end

    vcand = parsemaybe(datedeath, :upper)
    if ismissing(vcand)
        # we may have a date like "13. Jhd" only in datebirth
        vcand = parsemaybe(datebirth, :upper)
    end

    if !ismissing(vcand)
        eraend = vcand
    end

    for oc in eachrow(dfoffice)
        datestart = oc[:date_start]
        dateend = oc[:date_end]

        vcand = parsemaybe(datestart, :lower)
        if !ismissing(vcand) && vcand < erastart
            erastart = vcand
        end

        vcand = parsemaybe(dateend, :upper)
        if ismissing(vcand)
            # we may have a date like "13. Jhd"
            vcand = parsemaybe(datestart, :upper)
        end

        if !ismissing(vcand) && vcand > eraend
            eraend = vcand
        end

    end

    if erastart == Inf && eraend != -Inf
        erastart = eraend
    elseif erastart != Inf && eraend == -Inf
        eraend = erastart
    end

    if erastart == Inf
        erastart = missing
    end
    if eraend == -Inf
        eraend = missing
    end

    return erastart, eraend
end

"""
    fillpersondate(tblperson::AbstractString)::Int

parse `date_birth` and `date_death`.

"""
function fillpersondate(tblperson::AbstractString)::Int

    global dbwiag
    if isnothing(dbwiag)
        setDBWIAG()
    end

    msg = 1000

    sqlupdstart = "UPDATE " * tblperson *
        " SET numdate_birth = ?, date_hist_first = ?" *
        " WHERE id = ?";

    updstartstmt = DBInterface.prepare(dbwiag, sqlupdstart)

    sqlupdend = "UPDATE " * tblperson *
        " SET numdate_death = ?, date_hist_last = ?" *
        " WHERE id = ?";

    updendstmt = DBInterface.prepare(dbwiag, sqlupdend)

    sqldf = "SELECT id, date_birth, date_death FROM " * tblperson

    dfperson = DBInterface.execute(dbwiag, sqldf) |> DataFrame;

    tblid = 0
    for row in eachrow(dfperson)
        id, date_birth, date_death = row

        numdate_start = parsemaybe(date_birth, :lower)
        numdate_end = parsemaybe(date_death, :upper)
        if ismissing(numdate_end)
            numdate_end = parsemaybe(date_birth, :upper)
        end

        if !ismissing(numdate_start)
            DBInterface.execute(updstartstmt, [numdate_start, numdate_start, id])
        end

        if !ismissing(numdate_end)
            DBInterface.execute(updendstmt, [numdate_end, numdate_end, id])
        end

        tblid += 1
        if tblid % msg == 0
            @info tblid
        end
    end

    DBInterface.close!(updstartstmt)
    DBInterface.close!(updendstmt)

    return tblid
end

"""
    fillofficedateowntable(tblofficedate::AbstractString,
                       tbloffice::AbstractString,
                       colid::AbstractString,
                       colidperson = "wiagid",
                       colidinoffice = "wiagid_person",
                       tblperson = nothing)::Int

Extract dates as integer values. (2021-05-06 obsolete)
"""
function fillofficedateowntable(tblofficedate::AbstractString,
                                tbloffice::AbstractString;
                                colid::AbstractString = "id",
                                colidperson = "wiagid",
                                colidinoffice = "wiagid_person",
                                tblperson = nothing)::Int

    global dbwiag
    if isnothing(dbwiag)
        setDBWIAG()
    end

    msg = 1000

    sql = "SELECT o." * colid * ", date_start, date_end, id_monastery FROM " * tbloffice * " as o"

    if !isnothing(tblperson)
        sql *= " JOIN " * tblperson * " as p ON p." * colidperson * " = " * "o." * colidinoffice
    end

    dfoffice = DBInterface.execute(dbwiag, sql) |> DataFrame;

    tblid = 0;
    # do not clear the table
    @info "Append to table " tblofficedate
    sqli = "INSERT INTO " * tblofficedate * "(id_office, date_start, date_end) VALUES (?, ?, ?)"

    insstmt = DBInterface.prepare(dbwiag, sqli)

    for row in eachrow(dfoffice)
        id, date_start, date_end, id_monastery = row

        numdate_start = parsemaybe(date_start, :lower)
        numdate_end = parsemaybe(date_end, :upper)
        if ismissing(numdate_end)
            numdate_end = parsemaybe(date_start, :upper)
        end

        # push!(csqlvalues, "(" * id * ", " * numdate_start * ", " * numdate_end * ")")
        da = [id, numdate_start, numdate_end]
        DBInterface.execute(insstmt, da)
        tblid += 1
        if tblid % msg == 0
            @info tblid
        end
        # if tblid > 25 break end
    end

    DBInterface.close!(insstmt)
    #sqlvalues = join(csqlvalues, ", ")
    #DBInterface.execute(dbwiag, "INSERT INTO " * tblofficedate * " VALUES " * sqlvalues)

    return tblid
end

function fillofficedate(tbloffice::AbstractString;
                        colid::AbstractString = "id",
                        colidperson = "id_dh",
                        colidinoffice = "id_canon",
                        tblperson = nothing)::Int

    global dbwiag
    if isnothing(dbwiag)
        setDBWIAG()
    end

    msg = 1000

    sql = "SELECT o." * colid * ", date_start, date_end, id_monastery FROM " * tbloffice * " as o"

    # e.g. take only care of persons that are 'online'
    if !isnothing(tblperson)
        sql *= " JOIN " * tblperson * " as p ON p." * colidperson * " = " * "o." * colidinoffice
    end

    dfoffice = DBInterface.execute(dbwiag, sql) |> DataFrame;

    tblid = 0;
    # do not clear the table
    @info "Update table " tbloffice
    sqlstart = "UPDATE " * tbloffice * " SET numdate_start = ? WHERE id = ?"
    updstartstmt = DBInterface.prepare(dbwiag, sqlstart)

    sqlend = "UPDATE " * tbloffice * " SET numdate_end = ? WHERE id = ?"
    updendstmt = DBInterface.prepare(dbwiag, sqlend)

    for row in eachrow(dfoffice)
        id, date_start, date_end, id_monastery = row

        numdate_start = parsemaybe(date_start, :lower)
        numdate_end = parsemaybe(date_end, :upper)
        if ismissing(numdate_end)
            numdate_end = parsemaybe(date_start, :upper)
        end

        # push!(csqlvalues, "(" * id * ", " * numdate_start * ", " * numdate_end * ")")
        if !ismissing(numdate_start)
            DBInterface.execute(updstartstmt, [numdate_start, id])
        end

        if !ismissing(numdate_end)
            DBInterface.execute(updendstmt, [numdate_end, id])
        end

        tblid += 1
        if tblid % msg == 0
            @info tblid
        end
        # if tblid > 25 break end
    end

    DBInterface.close!(updstartstmt)
    DBInterface.close!(updendstmt)
    #sqlvalues = join(csqlvalues, ", ")
    #DBInterface.execute(dbwiag, "INSERT INTO " * tblofficedate * " VALUES " * sqlvalues)

    return tblid
end

"""
    fillofficelocation(tbloffice::AbstractString, tblmonasterylocation::AbstractString, tblplace::AbstractString; tblonline::AbstractString = nothing, colid::AbstractString)::Int

Find locations for offices that are related to a monastery or write the value of field diocese to location
"""
function fillofficelocation(tbloffice::AbstractString,
                            tblmonasterylocation::AbstractString,
                            tblplace::AbstractString,
                            tblonline::AbstractString = nothing,
                            colid::AbstractString = nothing)

    colnameofficeid = "id_office";

    global dbwiag
    if isnothing(dbwiag)
        setDBWIAG()
    end

    msg = 1000

    if isnothing(tblonline)
        sqlo = "SELECT o.id, id_monastery, location, diocese," *
            " o.numdate_start, o.numdate_end" *
            " FROM " * tbloffice * " as o"
    else
        sqlo = "SELECT o.id, id_monastery, location, diocese," *
            " o.numdate_start, o.numdate_end" *
            " FROM " * tbloffice * " as o" *
            " JOIN " * tblonline * " as l ON o.id_canon = l." * colid;
    end

    dfoffice = DBInterface.execute(dbwiag, sqlo) |> DataFrame;

    sqlml = "SELECT place_id, location_name, location_begin_tpq as loc_start, location_end_tpq as loc_end" *
        " FROM " * tblmonasterylocation * " WHERE location_name IS NOT NULL AND wiagid_monastery = ?"
    mlstmt = DBInterface.prepare(dbwiag, sqlml)

    sqlmlnn = "SELECT place_id, location_name, location_begin_tpq as loc_start, location_end_tpq as loc_end" *
        " FROM " * tblmonasterylocation * " WHERE location_name IS NULL AND wiagid_monastery = ?"
    mlnonamestmt = DBInterface.prepare(dbwiag, sqlmlnn)

    sqlp = "SELECT place_name FROM " * tblplace * " WHERE id_places IN (?)"
    plstmt = DBInterface.prepare(dbwiag, sqlp)

    sqlupd = "UPDATE " * tbloffice * " SET location_show = ? WHERE id = ?"
    updstmt = DBInterface.prepare(dbwiag, sqlupd)

    ntest = 100000
    ir = 0
    for row in eachrow(dfoffice)
        places = String[];
        id, id_monastery, location, diocese, date_start, date_end =
            row[[:id, :id_monastery, :location, :diocese, :numdate_start, :numdate_end]]

        # use values of `location` or `diocese` if present
        if !ismissing(location) && location != ""
            push!(places, location)
        elseif ismissing(id_monastery) || id_monastery == ""
            if !ismissing(diocese)
                push!(places, diocese)
            end
        else
            ffilter(loc_start, loc_end) = filterlocbydate(loc_start, loc_end, date_start, date_end)

            # monasteries where a location name is given
            ml = DBInterface.execute(mlstmt, [id_monastery]) |> DataFrame;
            nloc = size(ml, 1)
            if nloc == 1
                push!(places, ml[1, :location_name])
            elseif nloc > 1
                mlfilter = filter([:loc_start, :loc_end] => ffilter, ml)
                places = mlfilter[:, :location_name]
            else
                # monasteries where no location name is given
                mlnn = DBInterface.execute(mlnonamestmt, [id_monastery]) |> DataFrame;
                # filter only if there is more than one alternative
                if size(mlnn, 1) > 1
                    mlfilter = filter([:loc_start, :loc_end] => ffilter, mlnn)
                ids_place = mlfilter[:, :place_id]
                else
                    ids_place = mlnn[:, :place_id]
                end
                if length(ids_place) > 0
                    dfp = DBInterface.execute(plstmt, ids_place) |> DataFrame;
                    places = dfp[:, :place_name]
                else
                    @warn "No place found for office", id
                end
            end
        end
        if length(places) > 0
            locationshow = String(strip(places[1]))
            DBInterface.execute(updstmt, [locationshow, id])
            ir += 1
            if ir % msg == 0
                @info ir
            end
            if ir > ntest
                break
            end
        end

    end

    DBInterface.close!(updstmt)
    DBInterface.close!(plstmt)
    DBInterface.close!(mlnonamestmt)
    DBInterface.close!(mlstmt)

    return ir;
end

"""
    filterlocbydate(loc_start, loc_end, date_start, date_end)

Check if dates for the location and the office are compatibel with each other
"""
function filterlocbydate(loc_start, loc_end, date_start, date_end)
    loc_startint = parsemaybe(loc_start, :lower)
    if !ismissing(loc_startint) && !ismissing(date_end) && loc_startint > date_end
        return false
    end
    loc_endint = parsemaybe(loc_end, :upper)
    if !ismissing(loc_endint) && !ismissing(date_start) && loc_endint < date_start
        return false
    end
    return true
end

"""
    fillnamelookup(tbllookup::AbstractString,
                   tblperson::AbstractString,
                   colnameid::AbstractString = "id";
                   checkisonline = false)::Int

Fill `tablename` with combinations of givenname and familyname and their variants.
"""
function fillnamelookup(tbllookup::AbstractString,
                        tblperson::AbstractString,
                        colnameid::AbstractString = "id";
                        checkisonline = false)::Int
    msg = 200
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    DBInterface.execute(dbwiag, "DELETE FROM " * tbllookup)

    sql = "SELECT " * colnameid * " as id_person, " *
        "givenname, prefix_name, familyname, givenname_variant, familyname_variant " *
        "FROM " * tblperson * " person " *
        (checkisonline ? "WHERE status = 'online' " : "")

    dfperson = DBInterface.execute(dbwiag, sql) |> DataFrame

    # SQL
    # INSERT INTO dsttable VALUES (NULL, 'id_person1', 'givenname1', 'prefix_name1', 'familyname1'),
    # ('NULL', 'id_person2', 'givenname2', 'prefix_name2', 'familyname2');
    #
    # structure
    # gn[:] prefix fn|fnv
    # gn[1] prefix fn|fnv
    # gnv[:] prefix fn|fnv
    # gnv[1] prefix fn|fnv

    # In the web application choose a version with or without prefix.

    imsg = 0
    csqlvalues = String[]
    appendtosqlrow(row) = append!(csqlvalues, row)
    for row in eachrow(dfperson)
        idperson = row[:id_person]
        gn = row[:givenname]
        prefix = row[:prefix_name]
        fn = row[:familyname]
        gnv = row[:givenname_variant]
        fnv = row[:familyname_variant]

        fillnamelookupgn(idperson, gn, prefix, fn, fnv) |> appendtosqlrow

        if !ismissing(gnv) && gnv != ""
            # sets of givennames
            cgnv = split(gnv, r", *")
            for gnve in cgnv
                fillnamelookupgn(idperson, gnve, prefix, fn, fnv) |> appendtosqlrow
            end
        end
        imsg += 1

        if imsg % msg == 0
            @info imsg
        end
    end
    sqlvalues = join(csqlvalues, ", ")

    irowout = length(csqlvalues)

    DBInterface.execute(dbwiag, "INSERT INTO " * tbllookup * " VALUES " * sqlvalues)

    return irowout

end

"""
    fillcnnamelookup(tbllookup::AbstractString,
                     tblonline::AbstractString,
                     tblcanon::AbstractString,
                     tblcanongs::AbstractString)::Int

Fill `tablename` (for canons) with combinations of givenname and familyname and their variants.
"""
function fillcnnamelookup(tbllookup::AbstractString,
                          tblonline::AbstractString,
                          tblcanon::AbstractString,
                          tblcanongs::AbstractString)::Int
    msg = 400
    if isnothing(dbwiag)
        error("There is no valid database connection. Use `setDBWIAG'.")
    end

    DBInterface.execute(dbwiag, "DELETE FROM " * tbllookup)
    @info "clear " tbllookup

    sql = "SELECT id, id_dh, id_gs FROM " * tblonline

    dfonline = DBInterface.execute(dbwiag, sql) |> DataFrame

    sqldh = "SELECT id as id_person, " *
        "givenname, prefix_name, familyname, givenname_variant, familyname_variant " *
        "FROM " * tblcanon *
        " WHERE ID = ?"

    stmtdh = DBInterface.prepare(dbwiag, sqldh)
    sqlgs = "SELECT id as id_person, " *
        "givenname, prefix_name, familyname, givenname_variant, familyname_variant " *
        "FROM " * tblcanongs *
        " WHERE ID = ?"

    stmtgs = DBInterface.prepare(dbwiag, sqlgs)

    sqllookup = "INSERT INTO " * tbllookup * " VALUES (NULL, ?, ?, ?, ?, ?, ?)"

    stmtlookup = DBInterface.prepare(dbwiag, sqllookup)

    # SQL
    # INSERT INTO dsttable VALUES (NULL, 'id_person1', 'givenname1', 'prefix_name1', 'familyname1'),
    # ('NULL', 'id_person2', 'givenname2', 'prefix_name2', 'familyname2');
    #
    # structure
    # gn[:] prefix fn|fnv
    # gn[1] prefix fn|fnv
    # gnv[:] prefix fn|fnv
    # gnv[1] prefix fn|fnv

    # In the web application choose a version with or without prefix.

    imsg = 0

    for online in eachrow(dfonline)
        id_online, id_person = online[[:id, :id_dh]]
        nvars = 0
        if !ismissing(id_person)
            dfperson = DBInterface.execute(stmtdh, [id_person]) |> DataFrame
            if size(dfperson, 1) != 1
                @warn "No exact match for " id_online
            else
                nvars = insertlookuprows(stmtlookup, id_online, dfperson[1, :])
            end
        end
        id_online, id_person = online[[:id, :id_gs]]
        if !ismissing(id_person)
            dfperson = DBInterface.execute(stmtgs, [id_person]) |> DataFrame
            if size(dfperson, 1) != 1
                @warn "No exact match for " id_online
            else
                nvars += insertlookuprows(stmtlookup, id_online, dfperson[1, :])
            end
        end
        # @infiltrate id_online == "WIAG-Pers-CANON-80886-001"

        imsg += 1
        if imsg % msg == 0
            @info imsg
        end
    end

    return imsg
end

function insertlookuprows(stmt, id_online, row)
    gn = row[:givenname]
    prefix = row[:prefix_name]
    fn = row[:familyname]
    gnv = row[:givenname_variant]
    fnv = row[:familyname_variant]

    csqlvalues = String[]

    function insertvars(variants)
        for variant in variants
            DBInterface.execute(stmt, vcat([id_online], variant))
        end
    end

    nvars = 0
    vars = makevariantsgn(gn, prefix, fn, fnv)
    nvars += length(vars)
    insertvars(vars)

    if !ismissing(gnv) && gnv != ""
        # set of givennames
        cgnv = split(gnv, r", *")
        for gnve in cgnv
            vars = makevariantsgn(gnve, prefix, fn, fnv)
            nvars += length(vars)
            insertvars(vars)
        end
    end
    # @infiltrate id_online == "WIAG-Pers-CANON-80886-001"

    return nvars

end

"""
    sqlstring(s::AbstractString)::AbstractString

remove labels in data fields ("Taufname: Karl") and escape apostrophes
"""
function sqlstring(s::AbstractString)::AbstractString
    poslabel = findfirst(':', s)
    if !isnothing(poslabel)
        s = s[poslabel + 1:end]
    end
    s = replace(strip(s), "'" => "''")
    return s
end

"""
    makevariantsgn(gn, prefix, fn, fnv)

return an array of variants
"""
function makevariantsgn(gn, prefix, fn, fnv)
    csql = Vector{Vector{Union{String, Missing}}}()

    function getvalue(s)
        if ismissing(s) || s == "NULL" || s == ""
            return missing
        else
            return sqlstring(s)
        end
    end

    function pushcsql(gni, fni)
        sgni = getvalue(gni)
        sfni = getvalue(fni)
        sprefix = getvalue(prefix)
        # skip the prefix if it is contained in the variant
        if !ismissing(sprefix) && !ismissing(sfni) &&
            (occursin(sprefix * " ", sfni) || occursin(" " * sprefix * " ", sfni))
            sprefix = missing
        end
        values = [
            sgni,
            sprefix,
            sfni,
            ismissing(sfni) ? missing : sgni * " " * sfni,
            ismissing(sfni) || ismissing(sprefix) ? missing : sgni * " " * sprefix * " " * sfni
        ]
        push!(csql, values)
    end

    # pushcsql(gn, fn)
    cgn = split(gn);
    for gnsingle in cgn
        pushcsql(gnsingle, fn)
    end

    # more than one givenname; write complete name
    if length(cgn) > 1
        pushcsql(gn, fn)
    end

    # familyname variants
    if !ismissing(fnv) && strip(fnv) != ""
        cfnv = split(fnv, r", *")
        for fnve in cfnv
            pushcsql(gn, fnve)
            # more than one givenname
            for gnsingle in cgn
                pushcsql(gnsingle, fnve)
            end
        end
    end

    return csql
end

# parse time data
const rgpcentury = "([1-9][0-9]?)\\. (Jahrh|Jh)"
const rgpyear = "([1-9][0-9][0-9]+)"
const rgpyearfc = "([1-9][0-9]+)"

# turn of the century
const rgxtcentury = Regex("([1-9][0-9]?)\\./" * rgpcentury, "i")

# quarter
const rgx1qcentury = Regex("(1\\.|erstes) Viertel +(des )?" * rgpcentury, "i")
const rgx2qcentury = Regex("(2\\.|zweites) Viertel +(des )?" * rgpcentury, "i")
const rgx3qcentury = Regex("(3\\.|drittes) Viertel +(des )?" * rgpcentury, "i")
const rgx4qcentury = Regex("(4\\.|viertes) Viertel +(des )?" * rgpcentury, "i")

# begin, middle end
const rgx1tcentury = Regex("Anfang (des )?" * rgpcentury, "i")
const rgx2tcentury = Regex("Mitte (des )?" * rgpcentury, "i")
const rgx3tcentury = Regex("Ende (des )?" * rgpcentury, "i")

# half
const rgx1hcentury = Regex("(1\\.|erste) Hälfte +(des )?" * rgpcentury, "i")
const rgx2hcentury = Regex("(2\\.|zweite) Hälfte +(des )?" * rgpcentury, "i")

# between
const rgxbetween = Regex("zwischen " * rgpyear * " und " * rgpyear)

# early, late
const rgxearlycentury = Regex("frühes " * rgpcentury, "i")
const rgxlatecentury = Regex("spätes " * rgpcentury, "i")

# around, ...
const rgpmonth = "(Januar|Februar|März|April|Mai|Juni|Juli|August|September|Oktober|November|Dezember|Jan\\.|Feb\\.|Mrz\\.|Apr\\.|Jun\\.|Jul\\.|Aug\\.|Sep\\.|Okt\\.|Nov\\.|Dez\\.)"
const rgxbefore = Regex("(vor|bis|spätestens|spät\\.|v\\.)( [1-9][0-9]?\\.)? " * rgpmonth * "? ?" * rgpyear, "i")
# add 'circa'
const rgxaround = Regex("(um|circa|ca\\.|wahrscheinlich|wohl|etwa|evtl\\.) " * rgpyear, "i")
const rgxafter = Regex("(nach|frühestens|seit|ab) " * rgpyear, "i")

const rgxcentury = Regex("^ *" * rgpcentury)
const rgxyear = Regex("^( *|erwählt *)" * rgpyear)
const rgxyearfc = Regex("^( *|erwählt *)" * rgpyearfc)
const stripchars = ['†', '[', ']', ' ', '(', ')']
const rgxbelegt = r"belegt(.*)"

"""
    parsemaybe(s, Symbol::dir)

Parse `s` for an earliest or latest date.
"""
function parsemaybe(s, dir::Symbol)::Union{Missing, Int}
    if !(dir in [:lower, :upper])
        error("parameter dir must be :lower or :upper got ", dir)
    end

    year = missing
    if ismissing(s) || s == ""
        return year
    end

    # strip 'belegt'
    # prog: use replace instead
    rgm = match(rgxbelegt, s)
    if !isnothing(rgm)
        s = rgm[1]
    end

    # handle special cases
    s = strip(s, stripchars)

    if strip(s) == "?" return year end

    # turn of the century
    rgm = match(rgxtcentury, s)
    if !isnothing(rgm) && !isnothing(rgm[2])
        year = parse(Int, rgm[2])
        return year * 100
    end

    # quarter
    rgxq = [rgx1qcentury, rgx2qcentury, rgx3qcentury, rgx4qcentury]
    for (q, rgx) in enumerate(rgxq)
        rgm = match(rgx, s)
        if !isnothing(rgm) && !isnothing(rgm[3])
            century = parse(Int, rgm[3])
            if dir == :lower
                year = (century - 1) * 100 + (q - 1) * 25 + 1
                return year
            elseif dir == :upper
                year = (century - 1) * 100 + q * 25
                return year
            end
        end
    end

    # begin, middle, end
    rgxq = [rgx1tcentury, rgx2tcentury, rgx3tcentury]
    for (q, rgx) in enumerate(rgxq)
        rgm = match(rgx, s)
        if !isnothing(rgm) && !isnothing(rgm[3])
            century = parse(Int, rgm[2])
            if dir == :lower
                year = (century - 1) * 100 + (q - 1) * 33 + 1
                return year
            elseif dir == :upper
                year = (century - 1) * 100 + q * 33 + (q == 3 ? 1 : 0)
                return year
            end
        end
    end

    # half
    rgxq = [rgx1hcentury, rgx2hcentury]
    for (q, rgx) in enumerate(rgxq)
        rgm = match(rgx, s)
        if !isnothing(rgm) && !isnothing(rgm[3])
            century = parse(Int, rgm[3])
            if dir == :lower
                year = (century - 1) * 100 + (q - 1) * 50 + 1
                return year
            elseif dir == :upper
                year = (century - 1) * 100 + q * 50
                return year
            end
        end
    end

    # between
    rgm = match(rgxbetween, s)
    if !isnothing(rgm) && !isnothing(rgm[1]) && !isnothing(rgm[2])
        if dir == :lower
            year = parse(Int, rgm[1])
            return year
        elseif dir == :upper
            year = parse(Int, rgm[2])
            return year
        end
    end

    # early, late
    rgm = match(rgxearlycentury, s)
    if !isnothing(rgm) && !isnothing(rgm[1])
        century = parse(Int, rgm[1])
        if dir == :lower
            year = (century - 1) * 100 + 1
            return year
        elseif dir == :upper
            year = (century - 1) * 100 + 20
            return year
        end
    end

    rgm = match(rgxlatecentury, s)
    if !isnothing(rgm) && !isnothing(rgm[1])
        century = parse(Int, rgm[1])
        if dir == :lower
            year = century * 100 - 19
            return year
        elseif dir == :upper
            year = century * 100
            return year
        end
    end

    # before, around, after
    rgm = match(rgxbefore, s)
    if !isnothing(rgm)
        year = parse(Int, rgm[4])
        if dir == :lower
            year -= 50
        end
        return year
    end

    rgm = match(rgxafter, s)
    if !isnothing(rgm)
        year = parse(Int, rgm[2])
        if dir == :upper
            year += 50
        end
        return year
    end

    rgm = match(rgxaround, s)
    if !isnothing(rgm)
        year = parse(Int, rgm[2])
        if dir == :lower
            year -= 5
        elseif dir == :upper
            year += 5
        end
        return year
    end

    # century
    rgm = match(rgxcentury, s)
    if !isnothing(rgm) && !isnothing(rgm[1])
        century = parse(Int, rgm[1])
        if dir == :lower
            year = (century - 1) * 100 + 1
        elseif dir == :upper
            year = century * 100
        end
        return year
    end

    # plain year
    rgm = match(rgxyear, s)
    if !isnothing(rgm) && !isnothing(rgm[2])
        year = parse(Int, rgm[2])
        return year
    end

    # first century
    rgm = match(rgxyearfc, s)
    if !isnothing(rgm) && !isnothing(rgm[2])
        @info "First century date " s
        year = parse(Int, rgm[2])
        return year
    end

    # handle other special cases
    if strip(s) == "?" return year end

    ssb = strip(s, ['(', ')'])
    if ssb != s
        return parsemaybe(ssb, dir)
    end

    @warn "Could not parse " s
    return year

end

"""
    sqlvalue(data)

return quoted string or "NULL"
"""
function sqlvalue(data)::String
    value = ismissing(data) ? "NULL" : "'" * string(data) * "'"
    return value
end

end

max_stars_repo_name: ghertko/WiagDataSetup

max_issues_repo_name: ghertko/WiagDataSetup

max_forks_repo_name: ghertko/WiagDataSetup


text:
@testset "295.find-median-from-data-stream.jl" begin
    @testset "295-case 1" begin
        heap = MedianFinder()
        add_num!(heap, 1)
        add_num!(heap, 2)
        @test find_median(heap) == 1.5
        add_num!(heap, 3)
        @test find_median(heap) == 2
    end
    @testset "295-case 2" begin
        heap = MedianFinder()
        add_num!(heap, 5)
        add_num!(heap, 4)
        @test find_median(heap) == 4.5
        add_num!(heap, 2)
        @test find_median(heap) == 4
        add_num!(heap, 6)
        @test find_median(heap) == 4.5
    end
end