text:
using CSV
using DataFrames
using DataStructures
using JSON
const Path = Vector{String}

function parserange(str::AbstractString)
  parts = split(str, ":")
  if 2 == length(parts)
    b = parse(Float64, parts[1])
    e = parse(Float64, parts[2])
    b:e
  elseif 3 == length(parts)
    b = parse(Float64, parts[1])
    s = parse(Float64, parts[2])
    e = parse(Float64, parts[3])
    b:s:e
  else
    nothing
  end
end

function findranges(dict::AbstractDict{String,T} where T)
  paths = Vector{Path}()
  for key in keys(dict)
    val = dict[key]
    if isa(val, AbstractString)
      if nothing != parserange(val)
        push!(paths, Path([key]))
      end
    elseif isa(val, AbstractDict)
      append!(paths, map(path->insert!(copy(path), 1, key), findranges(dict[key])))
    end
  end
  paths
end

getbypath(dict::AbstractDict{T} where T<:AbstractString, path::Path) = getbypath(dict[path[1]], path[2:end])
getbypath(val, ::Path) = val

function setbypath!(dict::AbstractDict{T} where T<:AbstractString, path::Path, value::Real)
  if length(path) > 1
    setbypath!(dict[path[1]], path[2:end], value)
  else
    dict[only(path)] = value
  end
end

make_job_script(;
  cmd_dir::AbstractString,
  image_path::AbstractString,
  julia_path::AbstractString="julia",
  num_threads::Integer=1,
  )="""
#!/bin/bash
set -Eeuxo pipefail

cd $cmd_dir

JOB_IDX=`expr \$PBS_ARRAY_INDEX + 1`
JOB_DIR=`head -n "\${JOB_IDX}" jobdirs.txt | tail -n1`
cd "\${JOB_DIR}"

test -f _SUCCESS && exit 0

mkdir -p output

$julia_path -O3 --threads $num_threads -J $image_path \\
  -e 'MocosSimLauncher.launch(["params_experiment.json", "--output-summary", "output/summary.jld2"])' \\
  1> stdout.log \\
  2> stderr.log \\
  &

JULIA_PID=\$!
echo JULIA_PID=\$JULIA_PID
pidstat -r -p \$JULIA_PID 1 > memory.log &
pidstat -u -p \$JULIA_PID 1 > cpu.log &

wait \$JULIA_PID && touch _SUCCESS && exit 0

"""

function main()
  @assert length(ARGS) > 0 "JSON file needed"

  json = JSON.parsefile(ARGS[1], dicttype=OrderedDict)
  memory_gb = length(ARGS) > 1 ? ARGS[2] : 16
  num_threads = length(ARGS) > 2 ? ARGS[3] : 1
  workdir = length(ARGS) > 3 ? ARGS[4] : splitext(ARGS[1])[1]
  jobname = basename(workdir)

  rangepaths = findranges(json) |> sort
  ranges = map(x->getbypath(json, x) |> parserange, rangepaths)

  @info "ranges are" rangepaths ranges

  @assert length(ranges) == 2 "we support only 2D grids now"
  xpath, ypath = rangepaths
  xrange, yrange = ranges

  df = DataFrame([String[], Float64[], Float64[]], [:path, Symbol(join(xpath, "_")), Symbol(join(ypath, "_"))])

  for (x,xval) in enumerate(xrange)
    for (y,yval) in enumerate(yrange)
      subdir = joinpath(workdir, "grid_$(x-1)_$(y-1)")
      mkpath(subdir)
      specific_json = deepcopy(json)
      setbypath!(specific_json, xpath, xval)
      setbypath!(specific_json, ypath, yval)

      open(joinpath(subdir, "params_experiment.json"),"w") do f
        JSON.print(f, specific_json)
      end
      push!(df, ("grid_$(x-1)_$(y-1)", xval, yval))
    end
  end
  CSV.write(joinpath(workdir, "parameters_map.csv"), df)
  CSV.write(joinpath(workdir, "jobdirs.txt"), select(df, :path), writeheader=false)

  @info "generated $(nrow(df)) jobs"

  open(joinpath(workdir, "template.json"),"w") do f
    JSON.print(f, json, 2)
  end

  launcher_path = "/home/tomoz/MocosSimLauncher/"
  script = make_job_script(cmd_dir=abspath(workdir), image_path=joinpath(launcher_path, "sysimage.img"), num_threads=num_threads)
  write(joinpath(workdir, "script.sh"), script)
  num_jobs = nrow(df)

  cd(workdir)
  mkpath("task-logs")
  cd("task-logs")

  command = `qsub
    -J 0-$num_jobs
    -N $jobname
    -l walltime=48:00:00
    -l select=1:ncpus=$(num_threads):mem=$(memory_gb)gb
    -q "covid-19"
    ../script.sh
  `

  @info "executing command" command
  run(command)

end

main()


text:
module MathPlay

# Write your package code here.

export colatz

function colatz(x,n)
    res = BitVector(undef,n)
    for i = 1:n
        res[i] = isodd(x)
        x = iseven(x) ? x ÷ 2 : (3x + 1) ÷ 2
    end
    return res
end

end

ext: jl


text:
#= """
Produces a fresh temporary dir, which is deleted on teardown

Usage:

```
@pytest function(tempdir_fixture)
  cp(stuff, tempdir_fixture)
  # test stuff in the tempdir_fixture, don't worry about cleanup or clashes
end
```
""" =#
# FIXME this is not documentable, hmm...
@fixture function tempdir_fixture()
  the_directory = Base.mktempdir()
  @yield the_directory  # next lines are fixture teardown
  rm(the_directory, recursive=true)
end


text:
### mixture GP types.

# Structure for storing debug parameters from a MixtureGPType object when it is being queried.
mutable struct MixtureGPDebugType{T}

    w_tilde_set::Vector{Vector{T}}
    u_set::Vector{Vector{T}}
    v_set::Vector{Vector{T}}

    region_inds_set::Vector{Vector{Int}}
    p_region_ind_set::Vector{Int}

    hps_keep_flags_set::Vector{BitVector}
    zs_set::Vector{Vector{Vector{T}}}
    ts_set::Vector{Vector{T}}
end

function MixtureGPDebugType(dummy_val::T) where T

    w_tilde_set = Vector{Vector{T}}(undef, 0)
    u_set = Vector{Vector{T}}(undef, 0)
    v_set = Vector{Vector{T}}(undef, 0)

    region_inds_set = Vector{Vector{Int}}(undef, 0)
    p_region_ind_set = Vector{Int}(undef, 0)

    hps_keep_flags_set = Vector{BitVector}(undef, 0)
    zs_set = Vector{Vector{Vector{T}}}(undef, 0)
    ts_set = Vector{Vector{T}}(undef, 0)

    return MixtureGPDebugType{T}(w_tilde_set, u_set, v_set,
    region_inds_set, p_region_ind_set, hps_keep_flags_set,
    zs_set, ts_set)
end

# single kernel for all partitions.
mutable struct MixtureGPType{T}
    #
    X_parts::Vector{Vector{Vector{T}}}
    c_set::Vector{Vector{T}}

    # can probably delte the following.
    σ²_set::Vector{T}
    U_set::Vector{Matrix{T}}
    L_set::Vector{LowerTriangular{T, Matrix{T}}}

    #boundary_labels::Vector{Tuple{Int,Int}}
    hps::Vector{HyperplaneType{T}}

    #MixtureGPType{KT,T}(X_parts::Vector{Vector{Vector{T}}}, θ::KT) where {KT,T} = new{KT,T}(X_parts, θ)
end

function MixtureGPType(X_parts::Vector{Vector{Vector{T}}},
    hps::Vector{HyperplaneType{T}}) where T

    N = length(X_parts)

    c_set = Vector{Vector{T}}(undef, N)
    σ²_set = Vector{T}(undef, N)
    K_set = Vector{Matrix{T}}(undef, N)
    L_set = Vector{LowerTriangular{T, Matrix{T}}}(undef, N)

    return MixtureGPType(X_parts, c_set, σ²_set, K_set, L_set, hps)
end

function fitmixtureGP!(η::MixtureGPType{T},
    y_parts::Vector{Vector{T}},
    θ::KT, σ²) where {T,KT}

    N_parts = length(η.X_parts)

    ## get boundary labels.
    #M = round(Int, N_parts*(N_parts-1)/2)
    # boundary_labels = Vector{Tuple{Int,Int}}(undef, M)
    # a = 0

    # for j = 2:N_parts
    #     for i = 1:j-1
    #         a += 1
    #         boundary_labels[a] = (i,j)
    #     end
    # end
    # N_boundaries = length(boundary_labels)

    ### Set up kernel matrices.
    #K_XX_set = Vector{Matrix{T}}(undef, N_parts)

    for n = 1:N_parts

        X = η.X_parts[n]
        y = y_parts[n]

        # train GPs.
        U = constructkernelmatrix(X, θ)
        η.U_set[n] = copy(U) # store.

        # add observation model's noise.
        for i = 1:size(U,1)
            U[i,i] += σ²
        end

        c = U\y

        # cholesky.
        chol_U = cholesky(U)

        # store.
        η.L_set[n] = chol_U.L
        η.c_set[n] = c
        η.σ²_set[n] = σ²
    end

    return η
end

function querymixtureGP(xq::Vector{T},
    η::MixtureGPType{T},
    root,
    levels,
    radius::T,
    δ::T, # tolerance for partition search.
    θ::KT, σ²,
    weight_θ;
    debug_flag = false) where {KT, T}

    Xq = Vector{Vector{T}}(undef, 1)
    Xq[1] = xq

    return querymixtureGP(Xq, η::MixtureGPType{T}, root, levels,
    radius, δ, θ, σ², weight_θ; debug_flag = debug_flag)
end

function querymixtureGP(Xq::Vector{Vector{T}},
    η::MixtureGPType{T},
    root,
    levels,
    radius::T,
    δ::T, # tolerance for partition search.
    θ::KT, σ²,
    weight_θ;
    debug_flag = false) where {KT, T}

    Yq = Vector{T}(undef, 0)
    Vq = Vector{T}(undef, 0)

    debug_vars = MixtureGPDebugType(one(T))

    querymixtureGP!(Yq, Vq,
    Xq, η, root, levels, radius, δ, θ, σ², weight_θ,
    debug_vars; debug_flag = debug_flag)

    return Yq, Vq, debug_vars
end

function querymixtureGP!(Yq::Vector{T},
    Vq::Vector{T},
    Xq::Vector{Vector{T}},
    η::MixtureGPType{T},
    root,
    levels,
    radius::T,
    δ::T, # tolerance for partition search.
    θ::KT, σ²,
    weight_θ,
    debug_vars::MixtureGPDebugType{T};
    debug_flag = false)::Nothing where {KT, T}

    c_set = η.c_set
    L_set = η.L_set
    X_parts = η.X_parts
    hps = η.hps

    # set up output buffers.
    Nq = length(Xq)
    resize!(Yq, Nq)
    resize!(Vq, Nq)

    # debug.
    fill!(Yq, -20.0)
    fill!(Yq, -20.0)
    if debug_flag
        resize!(debug_vars.w_tilde_set, Nq)
        resize!(debug_vars.u_set, Nq)
        resize!(debug_vars.v_set, Nq)
        resize!(debug_vars.region_inds_set, Nq) # Vector{Vector{Int}}
        resize!(debug_vars.p_region_ind_set, Nq) # Vector{Int}

        resize!(debug_vars.hps_keep_flags_set, Nq)
        resize!(debug_vars.zs_set, Nq)
        resize!(debug_vars.ts_set, Nq)
    end

    # intermediates.
    kq = Vector{T}(undef, 0)
    w = Vector{T}(undef, 0) # mixture weights for GPs.
    u = Vector{T}(undef, 0) # means of GPs.
    v = Vector{T}(undef, 0) # variances of GPs..

    for j = 1:Nq

        xq = Xq[j]

        # find the region for p.
        p_region_ind = PatchMixtureKriging.findpartition(xq, root, levels)

        # get active partition list, then use it to get u, v, w.
        region_inds, ts, zs, hps_keep_flags = findneighbourpartitions(xq, radius, root,
            levels, hps, p_region_ind; δ = δ)
        t_kept = ts[hps_keep_flags]

        N_regions = length(region_inds)
        resize!(w, N_regions+1)
        resize!(u, N_regions+1)
        resize!(v, N_regions+1)

        # debug.
        fill!(w, -8.70)
        fill!(u, -8.70)

        for i = 1:N_regions

            r = region_inds[i]

            # i is index for objects returned by findneighbourpartitions()
            # r is the global region index, which is used by objects from η.

            w[i] = PatchMixtureKriging.evalkernel(abs(t_kept[i]), weight_θ)

            u[i], v[i] = queryinner!(kq, xq, X_parts[r], θ, c_set[r], L_set[r])
        end

        # u, v, w for the region that contains xq.
        w[end] = one(T)

        u[end], v[end] = queryinner!(kq, xq, X_parts[p_region_ind], θ, c_set[p_region_ind], L_set[p_region_ind])

        #debug.
        if debug_flag
            # w_set[j] = zeros(T, N_total_regions)

            # for r = 1:N_regions
            #     w_set[j][r] = w[r]
            # end
            # w_set[j][p_region_ind] = w[end]

            debug_vars.w_tilde_set[j] = copy(w)
            debug_vars.u_set[j] = copy(u)
            debug_vars.v_set[j] = copy(v)

            debug_vars.region_inds_set[j] = region_inds
            debug_vars.p_region_ind_set[j] = p_region_ind

            debug_vars.hps_keep_flags_set[j] = hps_keep_flags
            debug_vars.zs_set[j] = zs
            debug_vars.ts_set[j] = ts
        end

        # normalize to get a convex combination.
        sum_w = sum(w)
        for i = 1:length(w)
            w[i] = w[i]/sum_w
        end

        ## apply mixture.
        Yq[j] = dot(w, u)

        #Vq[j] = dot(w,diagm(v)*w)
        Vq[j] = dot(w, v .* w) # faster.

        # # hack.
        # Yq[j] = u[end]
        # Vq[j] = v[end]

        # ## debug.
        # if debug_flag
        #     println("p_region_ind = ", p_region_ind)
        #     println("region_inds = ", region_inds)
        #     println("ts = ", ts)
        #     println("zs = ", zs)
        #     println("hps_keep_flags = ", hps_keep_flags)
        #     println("w = ", w)
        #     println("u = ", u)
        #     println("v = ", v)
        #     #println("kq = ", kq)
        #     println()
        # end
    end

    return nothing
end

function queryinner!(kq::Vector{T}, xq, X, θ, c, L; min_v = 1e-12) where T

    @assert length(c) == length(X)

    resize!(kq, length(X))

    for i = 1:length(X)
        #kq[i] = evalkernel(norm(xq-X[i]), θ) # figure this out later.
        kq[i] = evalkernel(xq, X[i], θ) # figure this out later.
    end

    # mean.
    μq = dot(kq, c)

    ## variance.
    v = L\kq
    vq = clamp(evalkernel(xq, xq, θ) - dot(v,v), min_v, Inf)
    #vq = -1.23 # debug: placeholder.

    return μq, vq
end

function queryinner(xq::Vector{T}, X, θ, c, L) where T
    return queryinner!(Vector{T}(undef, 0), xq, X, θ, c, L)
end

function fetchhyperplanes(root::BinaryNode{PartitionDataType{T}}) where T

    hps = Vector{HyperplaneType{T}}(undef, 0)
    for node in AbstractTrees.PreOrderDFS(root)

        if isdefined(node.data.hp, :v)
            # non-leaf node.
            push!(hps, node.data.hp)
        end
    end

    return hps
end

"""
Given a point p and a list of hyperplanes hps, find the regions within adjacent to p within radius.
"""
function findneighbourpartitions(p::Vector{T},
    radius::T,
    root,
    levels,
    hps,
    p_region_ind::Int;
    δ::T = 1e-10) where T

    region_inds = Vector{Int}(undef, length(hps))
    j = 0

    keep_flags = falses(length(hps)) # for debug.
    zs = Vector{Vector{T}}(undef, length(hps)) # for debug.
    ts = Vector{T}(undef, length(hps)) # for debug.

    for i = 1:length(hps)

        # parse.
        u = hps[i].v
        c = hps[i].c

        # intersection with hyperplane.
        t = -dot(u, p) + c
        z = p + t .* u

        zs[i] = z # debug.
        ts[i] = t

        if norm(z-p) < radius

            # get two points along ray normal to the plane, on either side of plane.
            z1 = p + (t+δ) .* u
            z2 = p + (t-δ) .* u

            # get z1's, z2's region.
            z1_region_ind = PatchMixtureKriging.findpartition(z1, root, levels)
            z2_region_ind = PatchMixtureKriging.findpartition(z2, root, levels)

            # println("i = ", i)
            # println("z1 = ", z1)
            # println("z2 = ", z2)
            # println("z1_region_ind = ", z1_region_ind)
            # println("z2_region_ind = ", z2_region_ind)
            # println()

            # do not keep if both are in p's region.
            # keep if either z1 or z2 is in p's region.
            #if !(z2_region_ind == p_region_ind && z1_region_ind == p_region_ind) &&
            #    (z2_region_ind == p_region_ind || z1_region_ind == p_region_ind)
            if xor(z2_region_ind == p_region_ind, z1_region_ind == p_region_ind)

                keep_flags[i] = true

                j += 1
                region_inds[j] = z1_region_ind
                if z1_region_ind == p_region_ind
                    region_inds[j] = z2_region_ind
                end

            end
        end
    end

    resize!(region_inds, j)

    return region_inds, ts, zs, keep_flags
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelimitedFiles
using PyPlot
pygui(true)

times_tde = readdlm("./scripts/performance/results/times_tde.csv")
times_GA = readdlm("./scripts/performance/results/times_GA.csv")
times_mdop = readdlm("./scripts/performance/results/times_mdop.csv")
times_pec = readdlm("./scripts/performance/results/times_pec.csv")
Nss = readdlm("./scripts/performance/results/times.csv")

## Plot results

lw = 2.5      # linewidth
ms = 7        # markersize
fsa = 14      # fontsize
axislabelsize = 10
lwa = 1.5                 # axis line width
ticklabelsize = 10
fsl = 10        # legendfontsize

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
font0 = Dict(
        "font.size" => fsa,
        "font.weight" => "bold",
        "axes.labelweight" => "bold",
        "axes.titleweight" => "bold",
        "axes.labelsize" => axislabelsize,
        "axes.linewidth" => lwa,
        "xtick.labelsize" => ticklabelsize,
        "ytick.labelsize" => ticklabelsize,
        "legend.fontsize" => fsl)
merge!(rcParams, font0)

figure()
plot(Nss, times_tde./100000, label="TDE", marker="o", linestyle="solid", linewidth=lw, markersize=ms)
plot(Nss, times_GA/100000, label="G&A", marker="D", linestyle="solid", linewidth=lw, markersize=ms)
plot(Nss, times_mdop/100000, label="MDOP", marker="s", linestyle="solid", linewidth=lw, markersize=ms)
plot(Nss, times_pec/100000, label="PECUZAL", marker="x", linestyle="solid", linewidth=lw, markersize=ms)
plot(Nss, Nss.*log.(Nss), label="N log(N)", linestyle="dashed", linewidth=lw)
legend()
xlabel("time series length N")
ylabel("time [μs]")
title("  Median computation time \n of reconstruction methods")
yscale("log")
grid()


text:
using VertexSafeGraphs

using Test
import Random

import LightGraphs
const LG = LightGraphs

@testset "Graph construction and basic interface" begin
    nv = 20
    g1 = VSafeGraph(nv)
    @test LG.nv(g1) == nv
    @test LG.nv(g1.g) == nv
    @test LG.ne(g1) == 0
    @test LG.add_edge!(g1, 1, 2)
    g1_copy = copy(g1)
    @test LG.ne(g1) == 1 == LG.ne(g1_copy)
    @test LG.add_edge!(g1, 1, 3)
    @test LG.ne(g1) == 2 == LG.ne(g1_copy) + 1

    g2_inner = LG.CompleteGraph(nv)
    g2 = VSafeGraph(g2_inner)
    @test LG.nv(g2) == LG.nv(g2_inner)
    @test LG.ne(g2) == LG.ne(g2_inner)

    @test all(sort(collect(LG.vertices(g2))) .== sort(collect(LG.vertices(g2_inner))))

    @test LG.has_edge(g2, 1, 2)
    @test LG.has_edge(g2, LG.SimpleEdge(1, 2))
    @test LG.edgetype(g2) == LG.edgetype(g2_inner)

    g3 = VSafeGraph(LG.CompleteDiGraph(30))
    @test !LG.add_edge!(g3, 1, 2) # no possible addition on a complete graph
    @test !LG.add_edge!(g3, LG.SimpleEdge(1, 2))
    @test LG.is_directed(g3)
    @test !LG.is_directed(g2)
    @test LG.is_directed(typeof(g3))
    @test !LG.is_directed(typeof(g2))
    @test !LG.is_directed(VSafeGraph)
end

@testset "Vertex deletion" begin
    Random.seed!(33)
    nv = 45
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    @test LG.ne(inner) == LG.ne(g)
    @test LG.nv(inner) == LG.nv(g)
    nrm = 0
    for _ in 1:15
        removed_ok = LG.rem_vertex!(g, rand(1:nv))
        if !removed_ok
            continue
        end
        nrm += 1
        @test LG.nv(inner) == nv
        @test LG.nv(g) == nv - nrm
        @test length(g.deleted_vertices) == nrm

        @test LG.ne(inner) == LG.ne(g)
    end
    @test LG.add_vertex!(g)
    @test LG.nv(g) == nv - nrm + 1
    @test LG.nv(g.g) == nv + 1
end

@testset "Edge addition" begin
    Random.seed!(45)
    nv = 10
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    nea = 0         #number of edges added
    nrv = 0         #number of removed vertices
    neo = LG.ne(g)  #original number of edges
    for _ in 1:5
        removed_ok = LG.rem_vertex!(g, rand(1:nv))
        if !removed_ok
            continue
        end
        nrv += 1
    end

    ne = length(LG.edges(g))   #new number of edges
    for _ in 1:10
        v1 = rand(1:nv)
        v2 = rand(1:nv)
        while v2 == v1
            v2 = rand(1:nv)
        end
        added_ok = LG.add_edge!(g, v1, v2)
        if added_ok
            nea += 1
        end
        @test LG.ne(g) == ne + nea
	    @test LG.nv(g) == nv - nrv
    end
end

@testset "Neighbor sets" begin
    Random.seed!(99)
    nv = 10
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    nrv = 0
    for _ in 1:5
        rv = rand(1:nv)
        tv = rand(1:nv)
        while !LG.has_vertex(g, tv)
            tv = rand(1:nv)
        end
	removed_ok = LG.rem_vertex!(g, rv)
        if !removed_ok
            continue
        end
        nrv += 1
        @test_throws ArgumentError LG.inneighbors(g, rv)
        @test_throws ArgumentError LG.outneighbors(g, rv)
        @test !(rv in LG.outneighbors(g, tv))
        @test !(rv in LG.inneighbors(g, tv))
        @test LG.nv(g) == nv - nrv
    end
end

@testset "Generic LG algorithms work" begin
    nv = 45
    inner = LG.CompleteGraph(nv)
    g = VSafeGraph(inner)
    removed_ok = LG.rem_vertex!(g, rand(1:nv))
    @test removed_ok
    # LG broken here
    @test_throws BoundsError LG.pagerank(g)
    @test_throws BoundsError LG.kruskal_mst(g)
end


text:
using BenchmarkTools

function load_input(input_path)
  map(line -> parse(Int32, line), readlines(input_path))
end

function solve_1(data)
  increases = 0
  for i = 2:length(data)
    if data[i] > data[i-1]
      increases += 1
    end
  end

  return increases
end

function solve_2(data)
  increases = 0
  for i = 4:length(data)
    if data[i] > data[i-3]
      increases += 1
    end
  end

  return increases
end

function solve(data)
  println(@btime solve_1($data))
  println(@btime solve_2($data))
end

test_data = load_input("test_input.txt")
println("Test input:")
solve(test_data)

data = load_input("input.txt")
println("\nReal input:")
solve(data)

max_stars_repo_name: wbthomason/advent

max_issues_repo_name: wbthomason/advent

max_forks_repo_name: wbthomason/advent


text:
using Pkg
using Dates
const HOME = ENV["HOME"]
Pkg.activate("$HOME/HelioseismicKernels")

using MPI
MPI.Init()
using HelioseismicKernels
const comm = MPI.COMM_WORLD
const nworkers = MPI.Comm_size(comm)

try
    HelioseismicKernels.Cω(comm, Point2D(π/2, 0), Point2D(π/2, deg2rad(45)), los_earth(),
    ν_ind_range = 1:nworkers)

    tstart = time()
    HelioseismicKernels.Cω(comm, Point2D(π/2, 0), Point2D(π/2, deg2rad(45)), los_earth())
    tend = time()
    Δt = Time(0) + Second(round(Int, tend - tstart))
    if MPI.Comm_rank(comm) == 0
        println("Time taken = $Δt")
    end
finally
    MPI.Finalize()
end


text:
__precompile__(true)

module CQLdriver
using DataFrames, Dates, StructArrays, JuliaDB
export DataFrames, cqlinit, cqlclose, cqlwrite, cqlread, cqlexec, cqlreadtodf, CassSession, CassCluster

include("cqlwrapper.jl")
const CQL_OK = 0x0000

function Base.size(result::Ptr{CassResult})
    rows = cql_result_row_count(result)
    cols = cql_result_column_count(result)
    return (Int(rows)::Int, Int(cols)::Int)
end

Base.size(data::IndexedTable) = length(data), length(colnames(data))

# Table helpers
cass_combine(tbl::AbstractDataFrame, other_tbl::AbstractDataFrame) = hcat(tbl, other_tbl)
cass_combine(tbl::IndexedTable, other_tbl::IndexedTable) = merge(tbl, other_tbl)

cass_tbl_select(tbl::IndexedTable, colindex::Int, rowindex::Int) = columns(tbl)[colindex][rowindex]
cass_tbl_select(tbl::AbstractDataFrame, colindex::Int, rowindex::Int) = tbl[rowindex, colindex]

cass_tbl_slice(tbl::IndexedTable, row_start::Int, row_end::Int) = view(tbl, row_start:row_end)
cass_tbl_slice(tbl::AbstractDataFrame, row_start::Int, row_end::Int) = view(tbl, row_start:row_end, :)

cass_col_names(tbl::IndexedTable) = colnames(tbl)
cass_col_names(tbl::AbstractDataFrame) = names(tbl)

cass_get_lastindex(tbl::IndexedTable) = lastindex(tbl)
cass_get_lastindex(tbl::AbstractDataFrame) = lastindex(tbl, 1)

"""
    function cqlinit(hosts; username, password, threads, connections, queuesize, bytelimit, requestlimit)
Change the performance characteristics of your CQL driver
# Arguments
- `hosts::String`: a string with ipv4 addreses of hosts
- `username::String`: provide username for authenticated connections
- `password::String`: provide password for authenticated connections
- `threads::Int64`: set number of IO threads that handle query requests (default 1)
- `connections::Int64`: set number of connections per thread (default 2)
- `queuesize::Int64`: set queuesize that stores pending requests (default 4096)
- `bytelimit::Int64`: set max number of bytes pending on connection (default 65536 - 64KB)
- `requestlimit::Int64`: set max number of requests pending on connection (default 128 * connections)
- `whitelist::String`: set whiteslist of hosts - will only connect to these hosts and all other connections will be ignored
## Return
- `session::Ptr{CassSession}`: a pointer to the active session
- `cluster::Ptr{CassCluster}`: a pointer to the cluster
- `err::UInt`: a 16 bit integer with an error code. No error returns 0
"""
function cqlinit(hosts::String; username = "", password = "", whitelist = "", blacklist="", threads = 0, connections = 0, queuesize = 0, bytelimit = 0, requestlimit = 0, heartbeatinterval = 0, connecttimeo = 0, idletimeo = 0, reconnecttimeo = 0)
    cluster = cql_cluster_new()
    session = cql_session_new()

    err = CQL_OK

    if username != ""
	    cql_cluster_set_credentials(cluster, username, password)
    end
    if threads != 0
        err = cql_cluster_set_concurrency(cluster, threads) | err
    end
    if connections != 0
        err = cql_cluster_set_connections_per_host(cluster, connections) | err
    end
    if queuesize != 0
        err = cql_cluster_set_queue_size(cluster, queuesize) | err
    end
    if bytelimit != 0
        err = cql_cluster_set_write_bytes_high_water_mark(cluster, bytelimit) | err
    end
    if requestlimit != 0
        err = cql_cluster_set_pending_requests_high_water_mark(cluster, requestlimit) | err
    end
    if whitelist != ""
        cql_cluster_set_whitelist_filtering(cluster, whitelist)
    end
    if blacklist != ""
        cql_cluster_set_blacklist_filtering(cluster, blacklist)
    end
    if heartbeatinterval != 0
        cql_set_cluster_heartbeat_interval(cluster, heartbeatinterval)
    end
    if connecttimeo != 0
        cql_set_cluster_timeout(cluster, connecttimeo)
    end
    if idletimeo !=0
        cql_set_cluster_idle_timeout(cluster, idletimeo)
    end
    if reconnecttimeo !=0
        cql_set_cluster_reconnect_timeout(cluster, reconnecttimeo)
    end

    cql_cluster_set_contact_points(cluster, hosts)
    future = cql_session_connect(session, cluster)
    err = cqlfuturecheck(future, "Session Connect") | err
    @show err
    cql_future_free(future)
    return session::Ptr{CassSession}, cluster::Ptr{CassCluster}, err::UInt16
end

"""
    function cql_future_check(future, caller)
Check if a future contains any errors
# Arguments
- `future::Ptr{CassFuture}`: a pointer to a future
- `caller::String`: a string to help identify where this function is called from
# Return
- `err::UInt`: a 16 bit integer with an error code. No error returns 0
"""
function cqlfuturecheck(future::Ptr{CassFuture}, caller::String = "")
    #println("incqlfuturecheck")
    err = cql_future_error_code(future)
    #@show err
    # only prints valid messages for client errors
    if err != CQL_OK
        println("Error in CQL operation: ", caller)
        str = Vector{UInt8}(undef, 256)
        strref = Ref{Ptr{UInt8}}(pointer(str))
        siz = Ref{Csize_t}(sizeof(str))
        cql_future_error_message(future, strref, siz)
        println(unsafe_string(strref[], siz[]))
    end
    return err::UInt16
end

"""
    cql_val_type(result, idx)
Takes a CassResult and returns the type in a given column
# Arguments
- `result::Ptr{CassResult}`: a valid result from a query
- `idx::Int64`: the column to check
# Return
- `typ`: the type of the value in the specified column
"""
function cqlvaltype(result::Ptr{CassResult}, idx::Int64)
# http://datastax.github.io/cpp-driver/api/cassandra.h/#enum-CassValueType
    val = cql_result_column_type(result, idx)
    val == 0x0009 ? typ = Int32    : # INTEGER
    val == 0x0002 ? typ = Int64    : # BIGINT
    val == 0x0005 ? typ = Int64    : # COUNTER
    val == 0x0007 ? typ = Float64  : # DOUBLE
    val == 0x0008 ? typ = Float32  : # FLOAT
    val == 0x000A ? typ = String   : # TEXT
    val == 0x000D ? typ = String      : # VARCHAR
    val == 0x000B ? typ = DateTime : # TIMESTAMP
    val == 0x0014 ? typ = Int8     : # TINYINT
    val == 0x0013 ? typ = Int16    : # SMALLINT
    val == 0x0011 ? typ = Date     : # DATE
    val == 0x0004 ? typ = Bool     : # BOOLEAN
    val == 0x000C ? typ = UInt128  : # UUID
    val == 0x000F ? typ = UInt128  : # TIMEUUID
    val == 0x000E ? typ = BigInt   : # VARINT
    val == 0x0010 ? typ = IPAddr   : # INET
    val == 0x0006 ? typ = BigFloat : # DECIMAL
    val == 0x0012 ? typ = Missing      : # TIME
    val == 0x0001 ? typ = Missing      : # ASCII
    val == 0x0003 ? typ = Missing      : # BLOB
    val == 0xFFFF ? typ = Missing      : # UNKNOWN
    val == 0x0000 ? typ = Missing      : # CUSTOM
    val == 0x0015 ? typ = Missing      : # DURATION
    val == 0x0020 ? typ = Missing      : # LIST
    val == 0x0021 ? typ = Missing      : # MAP
    val == 0x0022 ? typ = Missing      : # SET
    val == 0x0030 ? typ = Missing      : # UDT
    val == 0x0031 ? typ = Missing      : # TUPLE
    typ = Missing
    if typ == Missing
        typ = UInt8
        println("Warning, unsupported datatype: $(num2hex(val))
        https://docs.datastax.com/en/developer/cpp-driver/2.8/api/cassandra.h/#enum-CassValueType")
    end
    return Union{typ, Missing}
end

"""
    function cqlgetvalue(val, t, strlen)
retrieve value using the correct type
# Arguments
- `val::Ptr{CassValue}`: a returned value from a query
- `t::DataType`: the type of the value being extracted
- `strlen::Int`: for string values specify max-length of output
# Return
- `out`: the return value, can by of any type
"""
function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Int64, Missing}}, strlen::Int)
    num = Ref{Clonglong}(0)
    err = cql_value_get_int64(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Bool, Missing}}, strlen::Int)
    num = Ref{Cint}(0)
    err = cql_value_get_bool(val, num)
    return ifelse(err == CQL_OK, Bool(num[]), missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Int32, Missing}}, strlen::Int)
    num = Ref{Cint}(0)
    err = cql_value_get_int32(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Int16, Missing}}, strlen::Int)
    num = Ref{Cshort}(0)
    err = cql_value_get_int16(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Int8, Missing}}, strlen::Int)
    num = Ref{Cshort}(0)
    err = cql_value_get_int8(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{String, Missing}}, strlen::Int)
    str = Vector{UInt8}(undef, strlen)
    strref = Ref{Ptr{UInt8}}(pointer(str))
    siz = Ref{Csize_t}(sizeof(str))
    err = cql_value_get_string(val, strref, siz)
    return ifelse(err == CQL_OK, unsafe_string(strref[], siz[]), missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Float64, Missing}}, strlen::Int)
    num = Ref{Cdouble}(0)
    err = cql_value_get_double(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Float32, Missing}}, strlen::Int)
    num = Ref{Cfloat}(0)
    err = cql_value_get_float(val, num)
    return ifelse(err == CQL_OK, num[], missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{Date, Missing}}, strlen::Int)
    num = Ref{Cuint}(0)
    err = cql_value_get_uint32(val, num)
    s = string(num[])
    l = length(s)
    o = ifelse(l == 8, s[1:4]*"-"*s[5:6]*"-"*s[7:8], "")
    return ifelse(err == CQL_OK, Date(o), missing)
end

function cqlgetvalue(val::Ptr{CassValue}, T::Type{Union{DateTime, Missing}}, strlen::Int)
    unixtime = Ref{Clonglong}(0)
    err = cql_value_get_int64(val, unixtime)
    return ifelse(err == CQL_OK, Dates.unix2datetime(unixtime[]/1000), missing)
end

cqlgetvalue(val::Ptr{CassValue}, T::Type{Union}, strlen::Int) = missing

"""
    function cqlstrprep(table, data)
create a prepared query string for use with batch inserts
# Arguments
- `table::String`: name of the table on the server
- `columns::Array{String}`: name of the columns on the server
- `data::Array{Any,1}`: an array of data to be inserted
# Return
- `out::String`: a valid INSERT or UPDATE query
"""
function cqlstrprep(table::String, data::Union{IndexedTable, AbstractDataFrame}; update::Union{IndexedTable, AbstractDataFrame, Nothing}=nothing, counter::Bool=false)
    out = ""
    if update == nothing
        datacolnames = string.(cass_col_names(data))
        cols, vals = "", ""

        for c in datacolnames
            cols = cols * c * ","
            vals = vals * "?,"
        end
        out = "INSERT INTO " * table * " (" * cols[1:end-1] * ") VALUES (" * vals[1:end-1] * ")"
    else write == :update
        datacolnames = string.(cass_col_names(data))
        updtcolnames = string.(cass_col_names(update))
        cols, vals = "", ""
        for c in datacolnames

            cols = cols * c * "=" * ifelse(counter, c*"+?, ", "?, ")
        end
        for u in updtcolnames
            vals = vals * u * "=? AND "
        end
        out = "UPDATE " * table * " SET " * cols[1:end-2] * " WHERE " * vals[1:end-5]
    end
    return out::String
end

"""
    function cqlstatementbind(statement, pos, typ, data)
Bind data to a column in a statement for use with batch inserts
# Arguments
- `statement::Ptr{CassStatement}`: pointer to a statement
- `pos::Int`: what column to put data into
- `typ::DataType, data)`: the datatype of the data
# Return
- `Void`:
"""
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Missing) = nothing # default to unset_value
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::String) = cql_statement_bind_string(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Bool) = cql_statement_bind_bool(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Int8) = cql_statement_bind_int8(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Int16) = cql_statement_bind_int16(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Int32) = cql_statement_bind_int32(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Int64) = cql_statement_bind_int64(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Float32) = cql_statement_bind_float(statement, pos, data)
cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Float64) = cql_statement_bind_double(statement, pos, data)

function cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::Date)
    d = parse(UInt32, replace(string(data),"-" => ""))
    cql_statement_bind_uint32(statement, pos, d)
end

function cqlstatementbind(statement::Ptr{CassStatement}, pos::Int, data::DateTime)
    d = convert(Int64, Dates.datetime2unix(data)*1000)
    cql_statement_bind_int64(statement, pos, d)
end

"""
function cqlclose(session, cluster)
Decommission a connection and free its resources
# Arguments
- `session::Ptr{CassSession}`: the current active session
- `cluster::Ptr{CassCluster}`: the cluster associated with the active session
# Return
- `Void`:
"""
function cqlclose(session::Ptr{CassSession}, cluster::Ptr{CassCluster})
    cql_session_free(session)
    cql_cluster_free(cluster)
end

function _cqlresultscheck(session::Ptr{CassSession}, statement::Ptr{CassStatement}, retries::Int)
    future = nothing
    while(true)
        future = cql_session_execute(session, statement)
        err = cqlfuturecheck(future, "Session Execute")
        err == CQL_OK && break
        if (err != CQL_OK) & (retries == 0)
            cql_statement_free(statement)
            cql_future_free(future)
            @show err
            println(typeof(err))
            return err
        end
        sleep(1)
        retries -= 1
        cql_future_free(future)
    end
    return CQL_OK, future
end

function _cqlprocessresult!(result::Ptr{CassResult}, output_arr::StructArray{NT}, statement::Ptr{CassStatement}, types::Vector, cols::Int, strlen::Int) where {NT}
    iterator = cql_iterator_from_result(result)
    @inbounds for r = eachindex(output_arr)
        cql_iterator_next(iterator)
        row = cql_iterator_get_row(iterator)
        output_arr[r] = NT(Tuple([cqlgetvalue(cql_row_get_column(row, c-1), types[c], strlen) for c = 1:cols]))
    end
    morepages = cql_result_has_more_pages(result)
    cql_statement_set_paging_state(statement, result)
    cql_iterator_free(iterator)
    cql_result_free(result)

    return morepages
end

function _cqlprocessresultdf!(result::Ptr{CassResult}, output::DataFrame, statement::Ptr{CassStatement}, types::Vector, cols::Int, strlen::Int, rows::Int)
    iterator = cql_iterator_from_result(result)
    @inbounds for r = 1:rows
        cql_iterator_next(iterator)
        row = cql_iterator_get_row(iterator)
        push!(output, [cqlgetvalue(cql_row_get_column(row, c-1), types[c], strlen) for c = 1:cols])
        #output_arr[r] = NT(Tuple([cqlgetvalue(cql_row_get_column(row, c-1), types[c], strlen) for c = 1:cols]))
    end
    morepages = cql_result_has_more_pages(result)
    cql_statement_set_paging_state(statement, result)
    cql_iterator_free(iterator)
    cql_result_free(result)

    return morepages
end

"""
    function cqlread(session, query; pgsize, retries, strlen)
    function cqlread(session, queries, concurrency; strlen)
Query the server for the contents of a table
- `session::Ptr{CassSession}`: pointer to the active session
- `query::String`: a valid SELECT query
- `queries::Array{String}`: an array of valid queries
- `concurrency::Int=500`: how many queries to execute
- `pgsize::Int=10000`: how many lines to pull at a time
- `retries::Int=5`: number of times to retry pulling a page of data
- `timeout::Int=10000`: time to wait for response in milliseconds
- `strlen::Int=128`: the maximum number of characters in a string
# Return
- `err::UInt16`: status of the query
- `output::DataFrame`: a dataframe with named columns
- `outputs::Array{DataFrame}`: an array of dataframe results
"""

function cqlreadtodf(session::Ptr{CassSession}, query::String; pgsize::Int=10000, retries::Int=8, timeout::Int=20000, strlen::Int=128)
    statement = cql_statement_new(query, 0)
    cql_statement_set_request_timeout(statement, timeout)
    cql_statement_set_paging_size(statement, pgsize)

    morepages = true
    err = CQL_OK

    # process first page
    err, future = _cqlresultscheck(session, statement, retries)
    if err != CQL_OK
        @show err
        return err::UInt16, DataFrame()
    end

    # get result
    result = cql_future_get_result(future)
    cql_future_free(future)
    rows, cols = size(result)

    # define all the types we will need
    types = [cqlvaltype(result, c-1) for c = 1:cols]
    names = Array{Symbol}(UndefInitializer(), cols)

    @inbounds for c = eachindex(names)
        str = zeros(UInt8, strlen)
        strref = Ref{Ptr{UInt8}}(pointer(str))
        siz = Ref{Csize_t}(sizeof(str))
        errcol = cql_result_column_name(result, c-1, strref, siz)
        names[c] = Symbol(ifelse(errcol == CQL_OK, unsafe_string(strref[], siz[]), string("C",c)))
    end

    outDF = DataFrame(types, names)

    morepages = _cqlprocessresultdf!(result, outDF, statement, types, cols, strlen, rows)

    while(morepages)
        err, future = _cqlresultscheck(session, statement, retries)
        if err != CQL_OK
            println("**")
            @show err
            return err::UInt16, outDF::DataFrame
        end

        # get result
        result = cql_future_get_result(future)
        cql_future_free(future)
        rows, cols = size(result)
        morepages = _cqlprocessresultdf!(result, outDF, statement, types, cols, strlen,rows)

    end
    cql_statement_free(statement)
    return err::UInt16, outDF::DataFrame
end

function cqlread(session::Ptr{CassSession}, query::String; pgsize::Int=10000, retries::Int=5, timeout::Int=10000, strlen::Int=128)
    statement = cql_statement_new(query, 0)
    #cql_statement_set_request_timeout(statement, timeout)
    cql_statement_set_paging_size(statement, pgsize)

    morepages = true
    err = CQL_OK

    # process first page
    err, future = _cqlresultscheck(session, statement, retries)
    if err != CQL_OK
        return err::UInt16, StructArray()
    end

    # get result
    result = cql_future_get_result(future)
    cql_future_free(future)
    rows, cols = size(result)

    # define all the types we will need
    types = [cqlvaltype(result, c-1) for c = 1:cols]
    names = Array{Symbol}(UndefInitializer(), cols)

    @inbounds for c = eachindex(names)
        str = zeros(UInt8, strlen)
        strref = Ref{Ptr{UInt8}}(pointer(str))
        siz = Ref{Csize_t}(sizeof(str))
        errcol = cql_result_column_name(result, c-1, strref, siz)
        names[c] = Symbol(ifelse(errcol == CQL_OK, unsafe_string(strref[], siz[]), string("C",c)))
    end

    NT = NamedTuple{Tuple(names), Tuple{types...}}
    SA_Type = StructArray{NT}

    output = SA_Type(undef, rows)

    morepages = _cqlprocessresult!(result, output, statement, types, cols, strlen)

    while(morepages)
        err, future = _cqlresultscheck(session, statement, retries)
        if err != CQL_OK
            return err::UInt16, output::SA_Type
        end

        # get result
        result = cql_future_get_result(future)
        cql_future_free(future)
        rows, cols = size(result)
        output_arr = SA_Type(undef, rows)

        morepages = _cqlprocessresult!(result, output_arr, statement, types, cols, strlen)

        output = vcat(output, output_arr)
    end
    cql_statement_free(statement)
    return err::UInt16, output::SA_Type
end

function cqlread(session::Ptr{CassSession}, queries::Array{String}; concurrency::Int=500, retries::Int=5, timeout::Int=10000, strlen::Int=128)
    out = Array{StructArray}(UndefInitializer(), 0)
    err = CQL_OK

    for query in 1:concurrency:length(queries)
        concurrency = ifelse(length(queries)-query < concurrency, length(queries)-query+1, concurrency)

        futures = Array{Ptr{CassFuture}}(UndefInitializer(), 0)
        for c in 1:concurrency
            statement = cql_statement_new(queries[query+c-1], 0)
          #  cql_statement_set_request_timeout(statement, timeout)
            push!(futures, cql_session_execute(session, statement))
            cql_statement_free(statement)
        end

        results = Array{Ptr{CassResult}}(UndefInitializer(), 0)
        for f in 1:length(futures)
            retry = retries
            future = futures[f]
            while(true)
                futerr = cqlfuturecheck(future, "Async Read")
                if futerr == CQL_OK
                    push!(results, cql_future_get_result(future))
                    cql_future_free(future)
                    break
                end
                if (futerr != CQL_OK) & (retry == 0)
                    err = futerr
                    cql_future_free(future)
                    break
                end
                sleep(1)
                retry -= 1
                statement = cql_statement_new(queries[query+f-1], 0)
                cql_statement_set_request_timeout(statement, timeout)
                future = cql_session_execute(session, statement)
                cql_statement_free(statement)
            end
        end

        for result in results
            rows, cols = size(result)
            iterator = cql_iterator_from_result(result)
            sa, types = cqlbuildstructarray(result, strlen)
            NT = eltype(sa)
            println(NT)
            @inbounds for r = eachindex(sa)
                cql_iterator_next(iterator)
                row = cql_iterator_get_row(iterator)
                sa[r] = NT(Tuple([cqlgetvalue(cql_row_get_column(row, c-1), types[c], strlen) for c = 1:cols]))
            end
            push!(out, sa)
            cql_iterator_free(iterator)
            cql_result_free(result)
        end
    end
    return err::UInt16, out
end

function cqlbuilddf(result::Ptr{CassResult}, strlen::Int)
    rows, cols = size(result)
    types = Array{Union}(UndefInitializer(), cols)
    for c in 1:cols
        types[c] = cqlvaltype(result, c-1)
    end
    names = Array{Symbol}(UndefInitializer(), cols)
    for c in 1:cols
        str = zeros(UInt8, strlen)
        strref = Ref{Ptr{UInt8}}(pointer(str))
        siz = Ref{Csize_t}(sizeof(str))
        errcol = cql_result_column_name(result, c-1, strref, siz)
        names[c] = Symbol(ifelse(errcol == CQL_OK, unsafe_string(strref[], siz[]), string("C",c)))
    end
    output = DataFrame(types, names, 0)
    return output::DataFrame, types::Array{Union}
end

function cqlbuildstructarray(result::Ptr{CassResult}, strlen::Int)
    rows, cols = size(result)
    types = [cqlvaltype(result, c-1) for c = 1:cols]
    names = Array{Symbol}(UndefInitializer(), cols)
    for c in 1:cols
        str = zeros(UInt8, strlen)
        strref = Ref{Ptr{UInt8}}(pointer(str))
        siz = Ref{Csize_t}(sizeof(str))
        errcol = cql_result_column_name(result, c-1, strref, siz)
        names[c] = Symbol(ifelse(errcol == CQL_OK, unsafe_string(strref[], siz[]), string("C",c)))
    end
    NT = NamedTuple{Tuple(names), Tuple{types...}}
    SA_Type = StructArray{NT}
    output = SA_Type(undef, rows)
    return output::SA_Type, types
end

"""
    function cqlbatchwrite(session, table, data; retries, update, counter)
Write a set of rows to a table as a prepared batch
# Arguments
- `session::Ptr{CassSession}`: pointer to the active session
- `table::String`: the name of the table you want to write to
- `data::DataFrame`: a DataFrame with named columns
- `retries::Int=5`: number of retries per batch insert
- `update::DataFrame`: the arguments for WHERE during an UPDATE
- `counter::Bool`: for updating the counter datatype
# Return
- `err::UInt16`: status of the batch insert
"""
function cqlbatchwrite(session::Ptr{CassSession}, cass_table::String, data::AbstractDataFrame; retries::Int=5, update::Union{AbstractDataFrame, Nothing}=nothing, counter::Bool=false)
    query = cqlstrprep(cass_table, data, update=update, counter=counter)
    future = cql_session_prepare(session, query)
    cql_future_wait(future)
    err = cqlfuturecheck(future, "Session Prepare")
    if err != CQL_OK
        cql_future_free(future)
        return err::UInt16
    end

    prep = cql_future_get_prepared(future)
    cql_future_free(future)
    batchtype = ifelse(!counter, 0x00, 0x02)
    batch = cql_batch_new(batchtype)
    rows, cols = size(data)
    frame = data
    if update != nothing
        urows, ucols = size(update)
        cols += ucols
        frame = hcat(data, update)
    end
    for r in 1:rows
        statement = cql_prepared_bind(prep)
        for c in 1:cols
            cqlstatementbind(statement, c-1, frame[r,c])
        end
        cql_batch_add_statement(batch, statement)
        cql_statement_free(statement)
    end
    while(true)
        future = cql_session_execute_batch(session, batch)
        cql_future_wait(future)
        err = cqlfuturecheck(future, "Execute Batch")
        cql_future_free(future)
        err == CQL_OK && break
        retries == 0 && break
        retries -= 1
        sleep(1)
    end
    cql_prepared_free(prep)
    cql_batch_free(batch)
    return err::UInt16
end

function cqlbatchwrite(session::Ptr{CassSession}, table::String, data::IndexedTable; retries::Int=5, update::Union{IndexedTable, Nothing}=nothing, counter::Bool=false)
    query = cqlstrprep(table, data, update=update, counter=counter)
    future = cql_session_prepare(session, query)
    cql_future_wait(future)
    err = cqlfuturecheck(future, "Session Prepare")
    if err != CQL_OK
        cql_future_free(future)
        return err::UInt16
    end

    prep = cql_future_get_prepared(future)
    cql_future_free(future)
    batchtype = ifelse(!counter, 0x00, 0x02)
    batch = cql_batch_new(batchtype)
    rows, cols = size(data)
    frame = data
    if update != nothing
        urows, ucols = size(update)
        cols += ucols
        frame = merge(data, update)
    end
    for r in 1:rows
        statement = cql_prepared_bind(prep)
        for c in 1:cols
            cqlstatementbind(statement, c-1, columns(frame)[c][r])
        end
        cql_batch_add_statement(batch, statement)
        cql_statement_free(statement)
    end
    while(true)
        future = cql_session_execute_batch(session, batch)
        cql_future_wait(future)
        err = cqlfuturecheck(future, "Execute Batch")
        cql_future_free(future)
        err == CQL_OK && break
        retries == 0 && break
        retries -= 1
        sleep(1)
    end
    cql_prepared_free(prep)
    cql_batch_free(batch)
    return err::UInt16
end

"""
    function cqlrowwrite(session, table, data; retries, update, counter)
Write one row of data to a table
# Arguments
- `session::Ptr{CassSession}`: pointer to the active session
- `cass_table::String`: the name of the table you want to write to
- `data::Union{IndexedTable, DataFrame}`: a DataFrame or IndexedTable with named columns
- `retries::Int=5`: number of retries per batch insert
- `update::DataFrame`: the arguments for WHERE during an UPDATE
- `counter::Bool`: for updating the counter datatype
# Return
- `err::UInt16`: status of the insert
"""

function cqlrowwrite(session::Ptr{CassSession}, cass_table::String, data::Union{DataFrame, IndexedTable}; retries::Int=5, update::Union{DataFrame, IndexedTable, Nothing} = nothing, counter::Bool = false)
    err = CQL_OK
    query = cqlstrprep(cass_table, data, update=update, counter=counter)
    rows, cols = size(data)
    frame = data
    if update != nothing
        urows, ucols = size(update)
        cols += ucols
        frame = cass_combine(data, update)
    end
    statement = cql_statement_new(query, cols)
    for c in 1:cols
        cqlstatementbind(statement, c-1, cass_tbl_select(frame, c, 1))
    end

    while(true)
        future = cql_session_execute(session, statement)
        cql_future_wait(future)
        err = cqlfuturecheck(future, "Execute Statement")
        cql_future_free(future)
        err == CQL_OK && break
        retries == 0 && break
        retries -= 1
        sleep(1)
    end
    cql_statement_free(statement)
    return err::UInt16
end

"""
    function cqlwrite(session, table, data; batchsize, retries, update, counter)
Write to a table
# Arguments
- `session::Ptr{CassSession}`: pointer to the active session
- `table::String`: the name of the table you want to write to
- `data::DataFrame`: a DataFrame with named columns
- `retries::Int=5`: number of retries per batch insert
- `batchsize::Int=1000`: number of rows to write per batch
- `update::DataFrame`: the arguments for WHERE during an UPDATE
- `counter::Bool`: for updating the counter datatype
# Return
- `err::UInt16`: status of the insert
"""
function cqlwrite(s::Ptr{CassSession}, cass_table::String, data::Union{DataFrame, IndexedTable}; update::Union{DataFrame, IndexedTable, Nothing}=nothing, batchsize::Int=500, retries::Int=5, counter::Bool=false, returnanyerror=false)
    rows, cols = size(data)
    rows == 0 && return 0x9999
    err = CQL_OK
    if rows == 1
        err = cqlrowwrite(s, cass_table, data, retries=retries, update=update, counter=counter)
    elseif rows <= batchsize
        err = cqlbatchwrite(s, cass_table, data, retries=retries, update=update, counter=counter)
    else
        pages = (rows ÷ batchsize)
        errs = zeros(UInt16, pages)
        @sync for p in 1:pages
            to = p * batchsize
            fr = to - batchsize + 1
            if p < pages
                @async errs[p] = cqlbatchwrite(s, cass_table, cass_tbl_slice(data, fr, to), retries=retries, update=update == nothing ? nothing : cass_tbl_slice(update, fr, to), counter=counter)
            else
                @async errs[p] = cqlbatchwrite(s, cass_table, cass_tbl_slice(data, fr, cass_get_lastindex(data)), retries=retries, update=update==nothing ? nothing : cass_tbl_slice(update, fr, cass_get_lastindex(update)), counter=counter)
            end
        end
        reterr = union(errs)
        if returnanyerror
            err = reterr[end]
        else
            err = reterr[1]
        end
    end
    return err::UInt16
end

function cqlwrite(s::Ptr{CassSession}, cass_table::String, data::JuliaDB.DIndexedTable; paritionsize::Int=100000, batchsize::Int=500, retries::Int=5, counter::Bool=false)
    errs = Vector{UInt16}()
    for tbl = Iterators.partition(data, paritionsize)
        push!(errs, cqlwrite(s, cass_table, tbl; batchsize=batchsize, retries=retries, counter=counter))
    end
    return union(errs)[1]
end

"""
    function cqlexec(session, statement)
Execute arbitrary command to the CQL database
# Arguments
- `session::Ptr{CassSession}`: pointer to the active session
- `statement::String`: a valid CQL command
# Return
- `err::UInt16`: status of the command
"""
function cqlexec(session::Ptr{CassSession}, cmd::String)
    err = CQL_OK
    statement = cql_statement_new(cmd, 0)
    future = cql_session_execute(session, statement)
    cql_future_wait(future)
    err = cqlfuturecheck(future, "Execute Statement")
    cql_future_free(future)
    cql_statement_free(statement)
    return err::UInt16
end

end


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

df = DataFrame(pb = rand(1:50, 50),
    t = rand(1:50, 50),
    h =[0.0 for i in 1:50])

plot(df, x = :pb, y = :t, color = :h, Geom.rectbin)


text:
@system Weather begin
    vp(context): vapor_pressure ~ ::VaporPressure

    PFD: photon_flux_density ~ preserve(u"μmol/m^2/s", parameter) #Quanta
    CO2: carbon_dioxide ~ preserve(u"μmol/mol", parameter)
    RH: relative_humidity ~ preserve(u"percent", parameter)
    T_air: air_temperature ~ preserve(u"°C", parameter)
    Tk_air(T_air): absolute_air_temperature ~ track(u"K")
    wind: wind_speed ~ preserve(u"m/s", parameter)
    P_air: air_pressure => 100 ~ preserve(u"kPa", parameter)

    VPD(T_air, RH, D=vp.D): vapor_pressure_deficit => D(T_air, RH) ~ track(u"kPa")
    VPD_Δ(T_air, Δ=vp.Δ): vapor_pressure_saturation_slope_delta => Δ(T_air) ~ track(u"kPa/K")
    VPD_s(T_air, P_air, s=vp.s): vapor_pressure_saturation_slope => s(T_air, P_air) ~ track(u"K^-1")
end


text:
facts("E calculus") do
    context("Get symbols in expression") do
        @fact get_all_syms(5) --> Set{Symbol}([])
        @fact get_all_syms(:x) --> Set([:x])
        @fact get_all_syms(:(x + 1)) --> Set([:x])
        @fact get_all_syms(:(x + y)) --> Set([:x, :y])
        @fact get_all_syms(:(x + (y * x))) --> Set([:x, :y])
        @fact get_all_syms(:(2x + (y * x))) --> Set([:x, :y])
    end

    context("Get indices in expression") do
        @fact get_all_inds(5) --> Set{Symbol}([])
        @fact get_all_inds(:x) --> Set{Symbol}([])
        @fact get_all_inds(:(x[i] + 1)) --> Set([:i])
        @fact get_all_inds(:(x[i] + y[j, k])) --> Set([:i, :j, :k])
        @fact get_all_inds(:(log(x[i]) + y[j, k])) --> Set([:i, :j, :k])
    end

    context("Strip indexing from expression") do
        @fact strip_inds(5) --> 5
        @fact strip_inds(:x) --> :x
        @fact strip_inds(:(x[i] + 1)) --> :(x + 1)
        @fact strip_inds(:(x[i] + y[j, k])) --> :(x + y)
        @fact strip_inds(:(log(x[i]) + y[j, k])) --> :(log(x) + y)
    end

    context("Basic identities") do
        @fact _simplify(1) --> 1
        @fact _simplify(ones(5)) --> ones(5)
        @fact _simplify(:x) --> :(x)
        @fact _simplify(:(E(x))) --> :(E(x))
        @fact _simplify(:(x + y)) --> :(x + y)

        @fact _simplify_call(Val{:f}, [:(x + y + z)]) --> :(f(x + y + z))
        @fact _simplify_call(Val{:f}, [:(log(x))]) --> :(f(log(x)))
        @fact _simplify_call(Val{:f}, [:(log(x)), :y, 5]) --> :(f(log(x), y, 5))
        @fact _simplify_compose(Val{:f}, Val{:g}, [:(x + y + z)]) --> :(f(g(x + y + z)))
        @fact _simplify_compose(Val{:E}, Val{:E}, [:x, :y, :z]) --> :(E(x, y, z))
        @fact _simplify_compose(Val{:E}, Val{:+}, [:x, :y, :z]) --> :(E(x) + E(y) + E(z))
        @fact _simplify_compose(Val{:E}, Val{:*}, [:x, :y, :z]) --> :(E(x) * E(y) * E(z))

        @fact _simplify_inside(Val{Symbol("'")}, [:x]) --> :(x')
        @fact _simplify_inside(Val{Symbol("'")}, [:(E(x + y))]) --> :((E(x) + E(y))')
    end

    context("Commutative operations") do
        @fact _simplify(:(E(x'))) --> :(E(x)')
        @fact _simplify(:(E(sum(x)))) --> :(sum(E(x)))
        @fact _simplify(:(E(dot(x, y)))) --> :(dot(E(x), E(y)))
        @fact _simplify(:(E(prod(x)))) --> :(prod(E(x)))
        @fact _simplify(:(C(sum(x)))) --> :(sum(C(x)))
        @fact _simplify(:(V(sum(x)))) --> :(sum(V(x)))
    end

    context("+ and -") do
        @fact _simplify(:(E(x + y))) --> :(E(x) + E(y))
        @fact _simplify(:(E(x + y + z))) --> :(E(x) + E(y) + E(z))
        @fact _simplify(:(E(x - y))) --> :(E(x) - E(y))
        @fact _simplify(:(E(x - y + z))) --> :(E(x) - E(y) + E(z))
        @fact _simplify(:(E(x .+ y))) --> :(E(x) .+ E(y))
        @fact _simplify(:(C(x + y))) --> _simplify(:(C(x) + C(y)))
        @fact _simplify(:(V(x + y))) --> _simplify(:(V(x) + V(y)))
    end

    context("*") do
        @fact _simplify(:(E(2x))) --> :(2 * E(x))
        @fact _simplify(:(E(2x * 3))) --> :((2 * E(x)) * 3)
        @fact _simplify(:(E(2x * y))) --> :((2 * E(x)) * E(y))
        @fact _simplify(:(E(2 * x * y))) --> :(2 * E(x) * E(y))
        @fact _simplify(:(E(2 * x * y * x))) --> :(2 * E(x * y * x))
        @fact _simplify(:(E(2 * x * y * x * z))) --> :(2 * E(x * y * x) * E(z))
        @fact _simplify(:(E((x * y) * x))) --> :(E((x * y) * x))
        @fact _simplify(:(E((x * y) * x * (w * z)))) --> :(E((x * y) * x) * (E(w) * E(z)))

        # @fact _simplify(:(C(x * y))) --> :((C(x) + E(x) * E(x)') * (C(y) + E(y) * E(y)') - (E(x) * E(x)') * (E(y) * E(y)'))
        @fact _simplify(:(V(dot(x, y)))) --> :(sum(C(x) .* C(y)) + E(y)' * C(x) * E(y) + E(x)' * C(y) * E(x))
        @fact _simplify(:(V(x * y))) --> :((V(x) + E(x)^2) * (V(y) + E(y)^2) - E(x)^2 * E(y)^2)
    end

    context("^") do
        @fact _simplify(:(E(x^2))) --> :(V(x) + E(x)^2)
        @fact _simplify(:(E((x + y)^2))) --> :((V(x) + V(y)) + (E(x) + E(y))^2)
        @fact _simplify(:((E(x + y)^2))) --> :((E(x) + E(y))^2)
    end

    context("Other special cases") do
        @fact _simplify(:(E(log(x)))) --> :(Elog(x))
    end

    context("Matrix linear combinations") do
        @fact _simplify(:(E(B[k]))) --> :(E(B[k]))
        @fact _simplify(:(E(B[k] * C[k]))) --> :(E(B[k]) * E(C[k]))
        @fact _simplify(:(E(sum(B[k] * C[k])))) --> :(sum(E(B[k]) * E(C[k])))
        @fact _simplify(:(E(A + sum(B[k] * C[k])))) --> :(E(A) + sum(E(B[k]) * E(C[k])))
        @fact _simplify(:(E(prod(B[k] * C[k])))) --> :(prod(E(B[k]) * E(C[k])))
    end

    context("macro expansion") do
        x ~ Normal(rand(), rand())
        y ~ Normal(rand(), rand())

        xy = @simplify E(x.data[1] + y.data[1])
        @fact xy --> E(x.data[1]) + E(y.data[1])

        xy = @simplify E(x.data[1] * y.data[1] + 5)
        @fact xy --> E(x.data[1]) * E(y.data[1]) + 5
    end
end


text:
struct BoolDomainViewNot <: BoolDomainView
    orig            ::AbstractBoolDomain
end

struct BoolVarViewNot <: BoolVarView
    x               ::AbstractBoolVar
    domain          ::BoolDomainViewNot
    id              ::String

    """
        BoolVarViewNot(x::AbstractBoolVar, id::String)

    Create a *fake* variable `y`, such that `y = ¬x`. This variable behaves like an usual one.
    """
    function BoolVarViewNot(x::AbstractBoolVar, id::String)
        dom = BoolDomainViewNot(x.domain)
        return new(x, dom, id)
    end
end

assignedValue(x::BoolVarViewNot) = !assignedValue(x.x)

"""
    isempty(dom::BoolDomainView)

Return `true` iff `dom` is an empty set.
"""
Base.isempty(dom::BoolDomainView) = isempty(dom.orig)

reset_domain!(dom::BoolDomainView) = reset_domain!(dom.orig)

"""
    rootVariable(x::BoolVarView)

Return the "true" variable behind `x`.
"""
rootVariable(x::BoolVarView) = rootVariable(x.x)

"""
    length(dom::BoolDomainView)

Return the size of `dom`.
"""
Base.length(dom::BoolDomainView) = length(dom.orig)

"""
    Base.in(value::Bool, dom::BoolDomainViewNot)

Check if a boolean is in the domain.
"""
Base.in(value::Bool, dom::BoolDomainViewNot) = (!value) in dom.orig

"""
    remove!(dom::BoolDomainViewNot, value::Bool)

Remove `value` from `dom`.
"""
remove!(dom::BoolDomainViewNot, value::Bool) = map(x -> !x, remove!(dom.orig, !value))

"""
    removeAll!(dom::BoolDomainViewNot)

Remove every value from `dom`. Return the removed values.
"""
removeAll!(dom::BoolDomainViewNot) = map(x -> !x, removeAll!(dom.orig))

"""
    assign!(dom::BoolDomainViewNot, value::Bool)

Remove everything from the domain but `value`. Return the removed values. Return the pruned values.
"""
assign!(dom::BoolDomainViewNot, value::Bool) = map(x -> !x, assign!(dom.orig, !value))

"""
    Base.iterate(dom::BoolDomainViewNot, state=1)

Iterate over the domain in an efficient way. The order may not be consistent.
WARNING: Do **NOT** update the domain you are iterating on.
"""
function Base.iterate(dom::BoolDomainViewNot, state=1)
    returned = iterate(dom.orig, state)
    if isnothing(returned)
        return nothing
    end

    value, newState = returned
    return !value, newState
end

parentValue(::BoolVarViewNot, v::Bool) = ! v


text:
using Interpolations, Distributions, ForwardDiff

const TINY = 1e-30
const _FWHM = 2 * sqrt(2 * log(2.0))
const PI4 = 4e-6 * pi

Base.float(d::ForwardDiff.Dual{T}) where T = ForwardDiff.Dual{T}(float(d.value), d.partials)
Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(prevfloat(float(d.value)), d.partials)
Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(nextfloat(float(d.value)), d.partials)
function Base.ldexp(x::T, e::Integer) where T<:ForwardDiff.Dual
    if e >=0
        x * (1<<e)
    else
        x / (1<<-e)
    end
end

"""
    abeles(q::Array{Float64, 1}, layers::Array{Float64, 2})

Performs the Abeles optical matrix calculation.

### Parameters
- `q::Array{Float64, 1}` : q-vector values.
- `layers::Array{Float64, 2}` : an Nx4 array, where N is the number of layers in the system, 1st item in a given row is the thickness, the 2nd the SLD, the 3rd the imaginary SLD, and the 4th the roughness with the layer above.

### Returns
- `::Array{Float64, 1}` : unsmeared reflectometry values for the given q-vectors.
"""
function abeles(q::Array{Float64, 1}, layers)
    nlayers = size(layers, 1) - 2
    npnts = length(q)
    reflectivity = Array{Any}(undef, npnts)

    sld = Array{Any}(undef, nlayers + 2)
    thick = Array{Any}(undef, size(layers, 1))
    for i = 1:size(layers, 1)
        thick[i] = layers[i, 1] * 1im
    end
    rough = Array{Any}(undef, nlayers + 1)

    for i = 2:nlayers+2
        sld[i] = PI4 * (layers[i, 2] - layers[1, 2] + 1im * (
            abs(layers[i, 3]) + TINY))
        rough[i - 1] = -2 * layers[i, 4] ^ 2
    end

    Threads.@threads for j = 1:npnts
        mr = Array{Any}(undef, (2, 2))
        mi = Array{Any}(undef, (2, 2))

        q2 = q[j] ^ 2. / 4. + 0im
        kn = q[j] / 2. + 0im

        for i = 1:nlayers + 1
            knn = sqrt(q2 - sld[i + 1])
            rj = (kn - knn) / (kn + knn) * exp(kn * knn * rough[i])

            if i == 1
                mr[1, 1] = 1. + 0im
                mr[1, 2] = rj
                mr[2, 2] = 1. + 0im
                mr[2, 1] = rj
            else
                beta = exp(kn * thick[i])
                mi[1, 1] = beta
                mi[2, 2] = (1. + 0im) / beta
                mi[2, 1] = rj * mi[1, 1]
                mi[1, 2] = rj * mi[2, 2]

                p0 = mr[1, 1] * mi[1, 1] + mr[2, 1] * mi[1, 2]
                p1 = mr[1, 1] * mi[2, 1] + mr[2, 1] * mi[2, 2]
                mr[1, 1] = p0
                mr[2, 1] = p1

                p0 = mr[1, 2] * mi[1, 1] + mr[2, 2] * mi[1, 2]
                p1 = mr[1, 2] * mi[2, 1] + mr[2, 2] * mi[2, 2]

                mr[1, 2] = p0
                mr[2, 2] = p1
            end
            kn = knn
        end

        r = mr[1, 2] / mr[1, 1]
        reflectivity[j] = real(r * conj(r))
    end
    return reflectivity
end

"""
    same_convolution(a::Array{Float64, 1}, b::Array{Float64, 1})

Performs a convolution of one-dimensional arrays equivalent to the [`np.convolve`](https://numpy.org/doc/stable/reference/generated/numpy.convolve.html), where the `mode` is `'same'`.

### Parameters
- `a::Array{Float64, 1}` : first array to convolve.
- `b::Array{Float64, 1}` : second array to convolve.

### Returns
- `::Array{Float64, 1}` : discrete, linear convolution of `a` and `b`.
"""
function same_convolution(a::Array{Any, 1}, b::Array{Float64, 1})

    m = length(a)
    n = length(b)

    padleft = ceil(Int32, n/2) - 1
    padright = floor(Int32, n/2)

    zero_padleft = Array{Any}(undef, padleft)
    for i in 1:padleft
        zero_padleft[i] = 0.
    end
    zero_padright = Array{Any}(undef, padleft)
    for i in 1:padleft
        zero_padright[i] = 0.
    end
    a = append!(zero_padleft, a)
    a = append!(a, zero_padright)

    output = Array{Any}(undef, m)
    for i in 1:m
        output[i] = 0.
    end
    for i in 1:m
        for j in 1:n
            output[i] += a[i+j-1] * b[n-j+1]
        end
    end
    return output
end

"""
    gauss(x::Array{Float64, 1}, s::Float64)

A Gaussian kernel for resolution smearing.

### Parameters
- `x::Array{Float64, 1}` : the kernal positions.
- `s::Float64` : the width of the kernel.

### Returns
- `::Array{Float64, 1}`: probabilities for `x`.
"""
function gauss(x::Array{Float64, 1}, s::Float64)
    g = zeros(size(x, 1))
    for i = 1:size(x, 1)
        g[i] = 1. / s / sqrt(2 * pi) * exp(-0.5 * x[i] ^ 2 / s / s)
    end
    return g
end

"""
    pointwise_smearing(q::Array{Float64, 1}, w::Array{Any, 2}, dq::Array{Float64, 1}; quad_order::Int32=17)

Perform the reflectometry calculation with a pointwise smearing.

### Parmaeters
- `q::Array{Float64, 1}` : q-vector values.
- `w::Array{Any, 2}` : an Nx4 array, where N is the number of layers in the system, 1st item in a given row is the thickness, the 2nd the SLD, the 3rd the imaginary SLD, and the 4th the roughness with the layer above.
- `dq::Array{Float64, 1}` : `dq` values for each `q` value.
- `quad_order{Int32, 1}` : The order of the Gaussian quadrature integration for the resolution integration.

### Returns
- `::Array{Float64, 1}` : smeared reflectometry values for the given q-vectors, with a pointise smearing.
"""
function pointwise_smearing(q::Array{Float64, 1}, w::Array{Any, 2}, dq::Array{Float64, 1}; quad_order::Int32=17)
end

"""
    constant_smearing(q::Array{Float64, 1}, w::Array{Any, 2}, resolution::Any, scale::Any, bkg::Any)

Perform the reflectometry calculation with a constant convolutional smearing.

### Parameters
- `q::Array{Float64, 1}` : q-vector values.
- `w::Array{Any, 2}` : an Nx4 array, where N is the number of layers in the system, 1st item in a given row is the thickness, the 2nd the SLD, the 3rd the imaginary SLD, and the 4th the roughness with the layer above.
- `resolution::Any` : the percentage resolution (dq/q) to be used.
- `scale::Any` : the multiplicative scale factor assocaited with the reflectometry profile.
- `bkg::Any` : the uniform background to add to the profile.

### Returns
- `::Array{Float64, 1}` : smeared reflectometry values for the given q-vectors.
"""
function constant_smearing(q::Array{Float64, 1}, w::Array{Any, 2}, resolution, scale, bkg)
    if resolution < 0.5
        return abeles(q, w)
    end

    resolution /= 100
    gaussnum = 51
    gaussgpoint = (gaussnum - 1) / 2

    lowq = minimum(q[:])
    highq = maximum(q[:])
    if lowq <= 0.
        lowq = 1e-6
    end

    start = log10(lowq) - 6 * resolution / _FWHM
    finish = log10(highq * (1 + 6 * resolution / _FWHM))
    interpnum = round(abs(1 * (abs(start - finish)) / (1.7 * resolution / _FWHM / gaussgpoint)))
    xtemp = range(start, finish, length=Int(interpnum))
    xlin = zeros(size(xtemp, 1))
    for i = 1:size(xtemp, 1)
        xlin[i] = 10. ^ xtemp[i]
    end

    gauss_x = collect(range(-1.7 * resolution, 1.7 * resolution, length=gaussnum))
    gauss_y = gauss(gauss_x, resolution / _FWHM)

    rvals = abeles(xlin, w)
    smeared_rvals = same_convolution(rvals, gauss_y)
    smeared_rvals *= gauss_x[2] - gauss_x[1]

    itp = LinearInterpolation(xlin, smeared_rvals)
    smeared_output = itp(q)
    for i = 1:size(smeared_output, 1)
        smeared_output[i] *= scale
        smeared_output[i] += bkg
    end
    return smeared_output
end


text:
const ASet{T} = Base.AbstractSet{T}
const ODict = IdDict

struct ObjectIdSet{T} <: ASet{T}
  dict::IdDict{Any, Any}
  ObjectIdSet{T}() where T = new(IdDict{Any, Any}())
end

Base.eltype(::ObjectIdSet{T}) where T= T

ObjectIdSet() = ObjectIdSet{Any}()

Base.push!(s::ObjectIdSet{T}, x::T) where T= (s.dict[x] = nothing; s)
Base.delete!(s::ObjectIdSet{T}, x::T) where T= (delete!(s.dict, x); s)
Base.in(x, s::ObjectIdSet) = haskey(s.dict, x)

(::Type{ObjectIdSet{T}})(xs) where T= push!(ObjectIdSet{T}(), xs...)

ObjectIdSet(xs) = ObjectIdSet{eltype(xs)}(xs)

Base.collect(s::ObjectIdSet) = collect(keys(s.dict))
Base.similar(s::ObjectIdSet, T::Type) = ObjectIdSet{T}()

@forward ObjectIdSet.dict Base.length

@iter xs::ObjectIdSet -> keys(xs.dict)

const OSet = ObjectIdSet

struct ObjectArraySet{T} <: ASet{T}
  xs::Vector{T}
  ObjectArraySet{T}() where T = new(T[])
end

Base.in(x::T, s::ObjectArraySet{T}) where T= any(y -> x ≡ y, s.xs)
Base.push!(s::ObjectArraySet, x) = (x ∉ s && push!(s.xs, x); s)

function Base.delete!(s::ObjectArraySet, x)
  i = findfirst(s.xs, x)
  i ≠ 0 && deleteat!(s.xs, i)
  return s
end

(::Type{ObjectArraySet{T}})(xs) where T= push!(ObjectArraySet{T}(), xs...)

ObjectArraySet(xs) = ObjectArraySet{eltype(xs)}(xs)

Base.collect(xs::ObjectArraySet) = xs.xs
Base.similar(s::ObjectArraySet, T::Type) = ObjectArraySet{T}()

@forward ObjectArraySet.xs Base.length

@iter xs::ObjectArraySet -> xs.xs

const OASet{T} = ObjectArraySet{T}


text:
using LRUCache
using Test
using Random
using Base.Threads

# Test insertion and reinsertion ordering
@testset "Single-threaded insertion and reinsertion" begin
    for cache in [LRU(; maxsize = 100), LRU{Int, Int}(; maxsize = 100)]
        r = 1:100
        for i in reverse(r)
            cache[i] = i
        end
        @test collect(cache) == collect(i=>i for i in r)

        for i = 1:10:100
            @test haskey(cache, i)
            @test !haskey(cache, 100+i)
        end

        # reinsert in random order
        p = randperm(100)
        for i in reverse(p)
            cache[i] = i
        end
        @test collect(cache) == collect(i=>i for i in p)
    end
end

@testset "Multi-threaded insertion and reinsertion" begin
    cache = LRU{Int, Int}(; maxsize = 100)

    r = 1:100
    @threads for i in reverse(r)
        cache[i] = i
    end

    @threads for i = 1:10:100
        @test haskey(cache, i)
        @test !haskey(cache, 100+i)
    end

    @test Set(cache) == Set(i=>i for i=1:100)

    @threads for i in 100:-1:1
        cache[i] = 2*i
    end

    @test Set(cache) == Set(i=>2*i for i=1:100)
end

@testset "Single-threaded getting and setting" begin
    cache = LRU{Int, Int}(; maxsize = 50*sizeof(Int), by = sizeof)

    for i in 100:-1:1
        @test 2*i == get(cache, i, 2*i)
        @test 2*i == get(()->2*i, cache, i)
        @test i == (iseven(i) ? get!(cache, i, i) : get!(()->i, cache, i))
        @test i == get!(cache, i, 2*i)
        @test i == get!(()->error("this should not happen"), cache, i)
        @test i == get(cache, i, 2*i)
        @test i == get(()->2*i, cache, i)
    end
    for i in 1:50
        @test haskey(cache, i)
        @test !haskey(cache, i+50)
        @test_throws KeyError getindex(cache, i+50)
    end

    @test collect(cache) == collect(i=>i for i in 1:50)

    p = randperm(50)
    for i in reverse(p)
        cache[i] = i
    end
    @test collect(cache) == collect(i=>i for i in p)

    p10 = p[1:10]
    resize!(cache; maxsize = 10*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    resize!(cache; maxsize = 50*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    for i in 50:-1:1
        @test get!(cache, i, 2*i) == (i in p10 ? i : 2*i)
    end

    for i in reverse(p10)
        @test cache[i] == i
    end
    resize!(cache; maxsize = 10*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    delete!(cache, p10[1])
    @test !haskey(cache, p10[1])

    @test p10[10] == pop!(cache, p10[10])
    @test !haskey(cache, p10[10])

    @test length(empty!(cache)) == 0
    @test isempty(cache)
    for i in p10
        @test !haskey(cache, p10)
    end
end

@testset "Multi-threaded getting and setting" begin
    cache = LRU{Int, Int}(; maxsize = 50*sizeof(Int), by = sizeof)

    @threads for i in 50:-1:1
        @test 2*i == get(cache, i, 2*i)
        @test 2*i == get(()->2*i, cache, i)
        @test i == (iseven(i) ? get!(cache, i, i) : get!(()->i, cache, i))
        @test i == get!(cache, i, 2*i)
        @test i == get!(()->error("this should not happen"), cache, i)
        @test i == get(cache, i, 2*i)
        @test i == get(()->2*i, cache, i)
    end

    p = randperm(50)
    p10 = p[1:10]
    @threads for i in p10
        @test cache[i] == i
    end
    resize!(cache; maxsize = 10*sizeof(Int))
    @test Set(cache) == Set(i=>i for i in p10)

    resize!(cache; maxsize = 50*sizeof(Int))
    @test Set(cache) == Set(i=>i for i in p10)

    delete!(cache, p10[1])
    @test !haskey(cache, p10[1])
    @test_throws KeyError getindex(cache, p10[1])

    @test p10[10] == pop!(cache, p10[10])
    @test !haskey(cache, p10[10])
    @test_throws KeyError getindex(cache, p10[1])
end

@testset "Recursive lock in get(!)" begin
    cache = LRU{Int,Int}(; maxsize = 100)
    p = randperm(100)
    cache[1] = 1

    f!(cache, i) = get!(()->(f!(cache, i-1) + 1), cache, i)
    @threads for i = 1:100
        f!(cache, p[i])
    end

    @threads for i = 1:100
        @test haskey(cache, i)
        @test cache[i] == i
    end
end

@testset "Eviction callback" begin
    # Julia 1.0 and 1.1 crash with multiple threads on this
    # combination of @threads and Channel.
    if VERSION >= v"1.2" || Threads.nthreads() == 1
        resources = Channel{Matrix{Float64}}(11)
        for _ = 1:11
            put!(resources, zeros(5, 5))
        end
        callback = (key, value) -> put!(resources, value)
        cache = LRU{Int,Matrix{Float64}}(; maxsize = 10, finalizer = callback)

        @threads for i = 1:100
            cache[i ÷ gcd(i, 60)] = take!(resources)
        end
        # Note: It's not ideal to rely on the Channel internals but there
        # doesn't seem to be a public way to check how much is occupied.
        @test length(resources.data) == 1
        cache[101] = take!(resources)
        @test length(resources.data) == 1
        pop!(cache, 101)
        @test length(resources.data) == 2
        cache[101] = take!(resources)
        @test length(resources.data) == 1
        delete!(cache, 101)
        @test length(resources.data) == 2
        get!(cache, 101, take!(resources))
        @test length(resources.data) == 1
        get!(() -> take!(resources), cache, 102)
        @test length(resources.data) == 1
        get!(() -> take!(resources), cache, 102)
        @test length(resources.data) == 1
        resize!(cache, maxsize = 5)
        @test length(resources.data) == 6
        empty!(cache)
        @test length(resources.data) == 11
    end
end

include("originaltests.jl")


text:
import BasisFunctions: approximate

export ErrorStyle
abstract type ErrorStyle end
export ResidualStyle, OversampledResidual, RandomPoints, FNAStyle
struct ResidualStyle <: ErrorStyle end
struct OversampledResidual <: ErrorStyle end
struct RandomPoints <: ErrorStyle end
struct FNAStyle{CREL,EREL} <: ErrorStyle
    FNAStyle() = new{true,true}()
end

export AdaptiveStrategy
abstract type AdaptiveStrategy end
export GreedyStyle, OptimalStyle, SimpleStyle
struct GreedyStyle <: AdaptiveStrategy end
struct OptimalStyle <: AdaptiveStrategy end
struct SimpleStyle <: AdaptiveStrategy end
struct OptimalStyleFirstFase <: AdaptiveStrategy end

emptylogbook() = Array{Any,1}(undef,0)
addlogentry!(log, entry) = push!(log, entry)

# Any error criterion returns a result of the form (converged,error):
# - converged: indicates whether the approximation satisfies all convergence criteria
# - error: a criterion-specific measure of the error

# Verify the average pointwise error in a limited number of random points.
function errormeasure(::RandomPoints, platform, tolerance, f, F, args...; numrandompts = 50, verbose=false, options...)
    g = randomgrid(support(F), numrandompts)
    z = sample(g, f)
    max_error = LinearAlgebra.generic_normInf(z-F(g))
    converged = max_error < tolerance
    verbose && println("Errormeasure: Maximum error in $numrandompts random points did $(converged ? "" : "not ")converge "*@sprintf("%1.3e (%1.3e)",max_error,tolerance))
    converged, max_error
end

# Measure the norm of the residual of the system Ax=B.
function errormeasure(::ResidualStyle, platform, tolerance, f, F, n, A, B, C, S, L; options...)
    residual = norm(A*C-B)
    residual <= tolerance, residual
end

# Measure the norm of the residual using oversampling (as specified by the residualoversampling factor)
function errormeasure(::OversampledResidual, platform, tolerance, f, F, args...; residualoversampling=2, M = nothing, L = nothing, options...)
    # Note that we leave out any M or L parameter that the user may have passed,
    # so that we can specify an oversampling factor to the discretization routine.
    A, B = full_discretization(f, platform, length(F), samplingstyle=OversamplingStyle(), oversamplingfactor=residualoversampling)
    residual = norm(A*coefficients(F)-B)
    residual <= error, residual
end

# Measure the error FNA style:
# - the residual has to meet a threshold
# - the norm of the coefficients has to be smaller than a value times the estimated norm of the right hand side
function errormeasure(::FNAStyle{CREL,EREL}, platform, δ, f, F, n, A, B, C, S, L; optimizefase=false,
        numrandompts = 3, FNAη = Inf, FNAδ = δ, verbose=false, Sbest=nothing, Bbest=nothing,options...) where {CREL, EREL}
    residual = norm(A*C-B)

    # Note: we pass on options, because it may contain a measure
    # Q = discrete_normalization(platform, n, L; S=S, options...)
    if optimizefase
        w = BasisFunctions.quadweights(grid(dest(Sbest)), measure(platform))
        normF = abs(sqrt(sum(w .* Bbest.^2)))
    else
        w = BasisFunctions.quadweights(grid(dest(S)), measure(platform))
        normF = abs(sqrt(sum(w .* B.^2)))
    end
    normErr = residual
    normCoef = norm(C)
    verbose && println( @sprintf("Errormeasure: approximate ||f|| %1.3e",normF))

    errConverged = EREL ? (normErr < FNAδ*normF) : (normErr < FNAδ)
    coefConverged = CREL ? (normCoef < FNAη*normF) : (normErr < FNAη)

    verbose && println( "Errormeasure: error did $(errConverged ? "" : "not ")converge: " * @sprintf("error=%1.3e, ||f||=%1.3e, δ=%1.3e",normErr,normF,FNAδ))
    verbose && println( "Errormeasure: coefficients did $(coefConverged ? "" : "not ")converge: " * @sprintf("||c||=%1.3e, ||f||=%1.3e, η=%1.3e",normCoef,normF,FNAη))

    converged = errConverged && coefConverged
    if converged
        if optimizefase
            converged, residual
        else
            # If N is too small the norm of `f` might not be good enough. User random points for robustness.
            verbose && println( "Errormeasure: Check random points for FNA robustness. ")
            converged, _ = errormeasure(RandomPoints(), platform, normF*FNAδ, f, F; numrandompts = numrandompts, verbose=verbose, options...)
            converged, residual
        end
    else
        converged, residual
    end
end

ErrorStyle(platform::Platform) = ErrorStyle(platform, DictionaryStyle(platform))
ErrorStyle(platform, ::BasisStyle) = RandomPoints()
ErrorStyle(platform, ::FrameStyle) = ResidualStyle()
ErrorStyle(platform, ::UnknownDictionaryStyle) = RandomPoints()

function adaptive_approximate_step!(logbook, f, platform, n, δ; criterion = ErrorStyle(platform), verbose=false, options...)
    fullverbose = (length(logbook) == 0) && verbose
    F, A, B, C, S, L = approximate(f, platform, n; verbose=fullverbose, options...)
    converged, error = errormeasure(criterion, platform, δ, f, F, n, A, B, C, S, L; verbose=verbose,options...)
    addlogentry!(logbook, (length(F), n, error,))
    F, A, B, C, S, L, converged, error
end

function initial_weight(platform, nprev, n, error)
    dict = dictionary(platform, n)
    ScalingOperator(dict, error)
end

function restrict_weight(platform, W, nbest, n)
    dictbest = dictionary(platform, nbest)
    dict = dictionary(platform, n)
    DiagonalOperator(dict, restriction(dictbest, dict)*diag(W))
end

function next_weight(platform, W, nprev, n, error)
    dict_small = dictionary(platform, nprev)
    dict_large = dictionary(platform, n)
    E = extension(dict_small, dict_large)
    D = E*diag(W)
    for i in eachindex(D)
        # might not be robust
        if D[i] == 0
            D[i] = error
        end
    end
    DiagonalOperator(dict_large, D)
end

# Dispatch on the style of the adaptive algorithm
function approximate(ap::AdaptiveApproximation;
            adaptivestyle = OptimalStyle(), criterion = ErrorStyle(platform(ap)),
            verbose = false, options...)

    # Inform the user that an adaptive computation is starting
    verbose && println("\nAdaptive style: $(adaptivestyle)")
    verbose && println("Error Criterium: $criterion")

    # Do the actual computation
    F, logbook, n, tol, error, iterations, converged = adaptive_approximation(adaptivestyle, ap_fun(ap), platform(ap); criterion=criterion, verbose=verbose, options...)

    # Report on the final convergence
    if !converged
        @warn "Adaptive: convergence to desired tolerance not reached after $(iterations) iterations."
        if verbose
            @printf "     Tolerance: %1.3e\n" tol
            @printf "     Minimal error: %1.3e\n\n" error
        end
    elseif verbose
        @printf "     Tolerance: %1.3e\n" tol
        println("\nAdaptive: Tolerance met using $(length(F)) degrees of freedom (n=$n) in $(iterations) iterations.\n")
    end

    return F
end

function adaptive_approximation(algorithm::Symbol, args...; options...)
    if algorithm == :greedy
        style = GreedyStyle()
    elseif algorithm == :simple
        style = SimpleStyle()
    elseif algorithm == :optimal
        style = OptimalStyle()
    else
        error("Unknown adaptive algorithm.")
    end
    adaptive_approximation(style, args...; options...)
end

# "SimpleStyle: the number of degrees of freedom is doubled until the tolerance is achieved."
# "GreedyStyle: Greedy scheme to get decreasing coefficients."
for (STYLE,nextparam,ret) in zip((:GreedyStyle,:SimpleStyle,:OptimalStyleFirstFase),(:param_increment,:param_double,:param_double),
        (Meta.parse("F, logbook, n, δ, error, iterations, converged"),
        Meta.parse("F, logbook, n, δ, error, iterations, converged"),
        Meta.parse("F, A, B, C, S, L, logbook, n, nprev, δ, error, iterations, converged, (weightedAZ ? (W) : ())...")) )
    @eval function adaptive_approximation(::$STYLE, f, platform;
                p0 = param_first(platform), maxlength = 2^12, maxiterations = 10,
                threshold = 1e-12, δ=100*threshold, verbose=false, weightedAZ=false, options...)

        iterations = 0
        logbook = emptylogbook()
        n = p0
        nprev = p0
        # First approximation
        verbose && println("Adaptive: initial value of n is $n")
        iterations += 1
        F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ; threshold=threshold, verbose=verbose, options...)
        verbose && @printf "Adaptive: N = %d, err = %1.3e, ||x|| = %1.3e\n\n" length(F) error norm(coefficients(F))
        converged && (return $ret)

        # Preparation for next step
        nprev = n
        n = $nextparam(platform, n)
        if weightedAZ
            W = initial_weight(platform, nprev, n, error)
        end

        # Next Steps
        while (length(F) <= maxlength) && (iterations <= maxiterations)
            iterations += 1

            # If weightedAZ add the weight
            if weightedAZ
                F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                    threshold=threshold, verbose=verbose,
                    weightedAZ = true, AZ_Cweight = W,
                    options...)
            else
                F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                    threshold=threshold, verbose=verbose,
                    options...)
            end
            verbose && @printf "Adaptive: N = %d, err = %1.3e, ||x|| = %1.3e\n\n" length(F) error norm(coefficients(F))
            converged && (return $ret)

            # Preparation for next step
            nprev = n
            n = $nextparam(platform, n)
            if weightedAZ
                W = next_weight(platform, W, nprev, n, error)
            end
        end
        @warn "Adaptive: No convergence with $n dofs and $iterations iterations"
        return $ret
    end
end

# "OptimalStyle: the number of degrees of freedom is chosen adaptively and optimally."
function adaptive_approximation(::OptimalStyle, f, platform;
        maxlength = 2^12, maxiterations = 100,
        threshold = 1e-12, δ=100*threshold, verbose=false, weightedAZ=false, stoptolerance=δ, options...)

    # First let the size grow until the tolerance is reached
    R = adaptive_approximation(OptimalStyleFirstFase(), f, platform;
            weightedAZ = weightedAZ, maxlength = maxlength, maxiterations = maxiterations,
            threshold = threshold, δ=δ, verbose=verbose, options...)

    if weightedAZ
        Fbest, Abest, Bbest, Cbest, Sbest, Lbest, logbook, upper_n, lower_n, δ, error, iterations, converged, W = R
    else
        Fbest, Abest, Bbest, Cbest, Sbest, Lbest, logbook, upper_n, lower_n, δ, error, iterations, converged = R
    end
    n = nbest = upper_n
    F, A, B, C, S, L = Fbest, Abest, Bbest, Cbest, Sbest, Lbest

    # The first phase was not successful
    if !converged
        verbose && println("Adaptive: The first fase was not successful. Optimal n lies beyond $nbest.\n")
        return Fbest, logbook, upper_n, δ, error, iterations, converged
    end

    # The first phase was successful
    verbose && println("Adaptive: Tolerance first met using $(length(Fbest)) degrees of freedom (n=$nbest). Optimal n lies in [$lower_n,$upper_n].\n")

    # Next, bisect until minimal n is found that achieves the tolerance
    while hasparam_inbetween(platform, lower_n, upper_n, stoptolerance) && (iterations <= maxiterations)
        n = param_inbetween(platform, lower_n, upper_n)

        if weightedAZ
            F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                threshold=threshold, verbose=verbose,
                weightedAZ = true, AZ_Cweight = restrict_weight(platform, W, nbest, n),
                optimizefase=true,Sbest=Sbest,Bbest=Bbest,
                options...)
        else
            F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                threshold=threshold, verbose=verbose,
                optimizefase=true,Sbest=Sbest,Bbest=Bbest,
                options...)
        end

        if !converged
            # At this stage, we know that lower_n does not meet the criterium, but the
            # new value of n does not either. It could be that lower_n==n, for example
            # when lower_n = 3 and upper_n = 4, and the param_inbetween method above
            # just returns 3 again. In order to avoid this, we increment n before
            # assigning to lower_n
            lower_n = param_increment(platform, n)
        else
            upper_n = n
        end
        verbose && @printf "Adaptive (phase 2): N = %d, err = %1.3e, ||x|| = %1.3e.\n" length(F) error norm(coefficients(F))
        verbose && println("Optimal param in [$(lower_n),$(upper_n)].\n")
        iterations += 1
    end
    if n != lower_n
        # We should take the upper_n. For simplicity, we recompute the approximation.
        n = upper_n
        if weightedAZ
            F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                threshold=threshold, verbose=verbose,
                weightedAZ = true, AZ_Cweight = restrict_weight(platform, W, nbest, n),
                optimizefase=true,Sbest=Sbest,Bbest=Bbest,
                options...)
        else
            F, A, B, C, S, L, converged, error = adaptive_approximate_step!(logbook, f, platform, n, δ;
                threshold=threshold, verbose=verbose,
                optimizefase=true,Sbest=Sbest,Bbest=Bbest,
                options...)
        end

    end
    if !converged
        @warn "Adaptive: No convergence with $n dofs and $iterations iterations"
    else
        verbose && println("Adaptive: Optimal n is $(n)!")
    end

    return F, logbook, n, δ, error, iterations, converged
end


text:
using StanDump, Test

using StanDump: dump, _dump, StanDumpIO

function stanrepr(x, internal = false; options...)
    io = IOBuffer()
    sd = StanDumpIO(io ; options...)
    (internal ? _dump : dump)(sd, x)
    String(take!(sd.io))
end

@testset "integer" begin
    @test stanrepr(1) == "1"
    @test stanrepr(typemax(Int32)+1) == string(typemax(Int32)+1) * "L"
    @test stanrepr(-99) == "-99"
    @test_throws ArgumentError stanrepr(BigInt(typemax(Int64))+1, true)
end

@testset "float" begin
    @test stanrepr(1/7) == string(1/7)
    @test stanrepr(-pi) == string(Float64(-pi))
    @test stanrepr(Inf) == "Inf"
    @test stanrepr(-Inf) == "-Inf"
    @test stanrepr(Inf32) == "Inf"
    @test stanrepr(-Inf32) == "-Inf"
    @test stanrepr(NaN) == "NaN"
    @test stanrepr(1//2) == "0.5"
end

@testset "unhandled" begin
    @test_throws ArgumentError stanrepr(:s)      # standalone symbol
    @test_throws ArgumentError stanrepr(nothing) # unknown type
    @test_throws ArgumentError stanrepr(:s__ => 1) # invalid name
    @test_throws ArgumentError stanrepr(Symbol("1s") => 1) # invalid name
end

@testset "definition and formatting" begin
    _repr(x; kwargs...) = stanrepr(x, true; kwargs...)
    @test _repr(:A => 99) == "A <- 99\n" # defaults
    # non-compact
    @test _repr(:A => 99; def_arrow = false, def_newline = false, compact = false) ==
        "A = 99\n" # same as above
    @test _repr(:A => 99; def_arrow = true, def_newline = false, compact = false) ==
        "A <- 99\n"
    @test _repr(:A => 99; def_arrow = false, def_newline = true, compact = false) ==
        "A =\n99\n"
    @test _repr(:A => 99; def_arrow = true, def_newline = true, compact = false) ==
        "A <-\n99\n"
    # compact
    @test _repr(:A => 99; def_arrow = false, def_newline = false, compact = true) ==
        "A=99\n"
    @test _repr(:A => 99; def_arrow = true, def_newline = false, compact = true) ==
        "A<-99\n"
    @test _repr(:A => 99; def_arrow = false, def_newline = true, compact = true) ==
        "A=\n99\n" # same as above
    @test _repr(:A => 99; def_arrow = true, def_newline = true, compact = true) ==
        "A<-\n99\n"
end

@testset "vector" begin
    @test stanrepr([1, 2]) == "c(1, 2)"
    @test stanrepr(2:10) == "2:10"
    @test stanrepr(10:2) == "integer(0)"
    @test stanrepr(Int[]) == "integer(0)"
    @test stanrepr([3.0, 7.0]) == "c(3.0, 7.0)"
    @test stanrepr(range(0, 1; length = 3)) == "c(0.0, 0.5, 1.0)"
end

@testset "array" begin
    @test stanrepr([1 2 3; 4 5 6], compact = true) ==
        "structure(c(1,4,2,5,3,6),.Dim=c(2,3))"
end

@testset "dict" begin
    d = Dict(:a => [1, 2], :b => 9.0)
    s = "a <- c(1, 2)\nb <- 9.0\n"
    let io = IOBuffer()
        stan_dump(io, d)
        @test String(take!(io)) == s
    end
end

@testset "general" begin
    io = IOBuffer()
    sd = StanDumpIO(io; compact = true)
    stan_dump(sd, (a = 1, b = 2)) # multiple arguments
    @test String(take!(io)) == "a<-1\nb<-2\n"
end

@testset "accidental varnames (from pairs)" begin
    @test_throws ArgumentError stan_dump(stdout, 1:3)
end