text:
using FileIO
using LasIO
using Test
using FixedPointNumbers
using Dates

include("read_write_laspoints.jl")
include("stream.jl")

workdir = dirname(@__FILE__)
# source: http://www.liblas.org/samples/
filename = "libLAS_1.2.las" # point format 0
testfile = joinpath(workdir, filename)
writefile = joinpath(workdir, "libLAS_1.2-out.las")

"Find the centroid of all points in a LAS file"
function centroid(io, header)
    x_sum = 0.0
    y_sum = 0.0
    z_sum = 0.0
    n = Int(header.records_count)

    for i = 1:n
        p = read(io, LasPoint0)
        x = xcoord(p, header)
        y = ycoord(p, header)
        z = zcoord(p, header)

        x_sum += x
        y_sum += y
        z_sum += z
    end

    x_avg = x_sum / n
    y_avg = y_sum / n
    z_avg = z_sum / n

    x_avg, y_avg, z_avg
end

# reading point by point
open(testfile) do io
    # magic bytes
    @test String(read(io, 4)) == "LASF"
    header = read(io, LasHeader)

    seek(io, header.data_offset)
    x_avg, y_avg, z_avg = centroid(io, header)

    @test x_avg ≈ 1442694.2739025319
    @test y_avg ≈ 377449.24373880465
    @test z_avg ≈ 861.60254888088491

    seek(io, header.data_offset)
    p = read(io, LasPoint0)

    @test xcoord(p, header) ≈ 1.44013394e6
    @test xcoord(1.44013394e6, header) ≈ p.x
    @test ycoord(p, header) ≈ 375000.23
    @test ycoord(375000.23, header) ≈ p.y
    @test zcoord(p, header) ≈ 846.66
    @test zcoord(846.66, header) ≈ p.z
    @test intensity(p) === 0x00fa
    @test scan_angle(p) === Int8(0)
    @test user_data(p) === 0x00
    @test pt_src_id(p) === 0x001d
    @test return_number(p) === 0x00
    @test number_of_returns(p) === 0x00
    @test scan_direction(p) === false
    @test edge_of_flight_line(p) === false
    @test classification(p) === 0x02
    @test synthetic(p) === false
    @test key_point(p) === false
    @test withheld(p) === false

    # raw bytes composed of bit fields
    @test flag_byte(p) === 0x00
    @test raw_classification(p) === 0x02

    # recompose bytes with bit fields
    @test flag_byte(return_number(p),number_of_returns(p),scan_direction(p),edge_of_flight_line(p)) === p.flag_byte
    @test raw_classification(classification(p),synthetic(p),key_point(p),withheld(p)) === p.raw_classification

    # TODO GPS time, colors (not in this test file, is point data format 0)
end

# reading complete file into memory
# test if output file matches input file
header, pointdata = load(testfile)
n = length(pointdata)
save(writefile, header, pointdata)
@test hash(read(testfile)) == hash(read(writefile))
rm(writefile)

# testing a las file version 1.0 point format 1 file with VLRs
srsfile = joinpath(workdir, "srs.las")
srsfile_out = joinpath(workdir, "srs-out.las")
srsheader, srspoints = load(srsfile)
for record in srsheader.variable_length_records
    @test record.reserved === 0xaabb
    @test record.user_id == "LASF_Projection"
    @test typeof(record.description) == String
    if record.record_id == 34735
        @test record.data.key_directory_version === UInt16(1)
        @test record.data.key_reversion === UInt16(1)
        @test record.data.minor_revision === UInt16(0)
        @test record.data.number_of_keys === UInt16(length((record.data.keys)))
        @test typeof(record.data.keys) == Vector{LasIO.KeyEntry}
    end
end
@test srsheader.version_major == 1
@test srsheader.version_minor == 0
@test srsheader.data_format_id == 1
@test srsheader.n_vlr == 3
@test isa(srsheader.variable_length_records, Vector{LasVariableLengthRecord})
for vlr in srsheader.variable_length_records
    @test vlr.reserved === 0xaabb
    @test vlr.user_id == "LASF_Projection"
    @test vlr.description == ""
end

@test srsheader.variable_length_records[1].record_id == LasIO.id_geokeydirectorytag
@test srsheader.variable_length_records[2].record_id == LasIO.id_geodoubleparamstag
@test srsheader.variable_length_records[3].record_id == LasIO.id_geoasciiparamstag
@test typeof(srsheader.variable_length_records[1].data) == LasIO.GeoKeys
@test typeof(srsheader.variable_length_records[2].data) == LasIO.GeoDoubleParamsTag
@test typeof(srsheader.variable_length_records[3].data) == LasIO.GeoAsciiParamsTag

@test LasIO.epsg_code(header) === nothing
@test LasIO.epsg_code(srsheader) === UInt16(32617)
# set the SRS. Note: this will not change points, but merely set SRS-metadata.
epsgheader = deepcopy(srsheader)
LasIO.epsg_code!(epsgheader, 32633)  # set to WGS 84 / UTM zone 33N, not the actual SRS
@test epsgheader.variable_length_records[1].record_id == LasIO.id_geokeydirectorytag
@test count(LasIO.is_srs, srsheader.variable_length_records) == 3
@test count(LasIO.is_srs, epsgheader.variable_length_records) == 1

save(srsfile_out, srsheader, srspoints)
@test hash(read(srsfile)) == hash(read(srsfile_out))
rm(srsfile_out)

@test boundingbox(header) === (xmin = 1.44e6, ymin = 375000.03, zmin = 832.1800000000001,
    xmax = 1.44499996e6, ymax = 379999.99, zmax = 972.6700000000001)
@test LasIO.boundingbox_unscaled(header) === (xmin = Int32(144000000), ymin = Int32(37500003), zmin = Int32(83218), xmax = Int32(144499996), ymax = Int32(37999999), zmax = Int32(97267))


text:
mutable struct NodeData
  exerciseValue::Float64
  cumulatedCashFlows::Float64
  values::Vector{Float64}
  controlValue::Float64
  isValid::Bool
end

NodeData() = NodeData(0.0, 0.0, Vector{Float64}(), 0.0, false)

ext: jl


text:
module AnyBoost

using Formatting
using MicroLogging # CoreLogging
using DataFrames
using Distributions

include("Tbasis.jl")
include("Tlosses.jl")
include("Tconstraintsets.jl")
include("Tproblem.jl")
include("Tensemble.jl")
include("Tlrensemble.jl")

include("Tnextbasissolver.jl")
include("SubNextBasis.jl")
include("SubLassoOnA.jl")
include("SubLineSearch.jl")
include("SubBackProp.jl")

include("AlgADFW.jl")
include("AlgBoosting.jl")
include("AlgStagewise.jl")
include("AlgNN.jl")

end


text:
get_layout(gl::GridLayout) = gl
get_layout(f::Union{Figure, GridPosition}) = f.layout
get_layout(l::Union{MakieLayout.Layoutable, GridSubposition}) = get_layout(l.parent)

# Wrap layout updates in an update block to avoid triggering multiple updates
function update(f, fig)
    layout = get_layout(fig)
    block_updates = layout.block_updates
    layout.block_updates = true
    output = f(fig)
    layout.block_updates = block_updates
    block_updates || update!(layout)
    return output
end

function Makie.plot!(fig, s::OneOrMoreLayers;
                     axis=NamedTuple(), palettes=NamedTuple())
    if isa(fig, Union{Axis, Axis3}) && !isempty(axis)
        @warn("Axis got passed, but also axis attributes. Ignoring axis attributes $axis.")
    end
    grid = update(f -> compute_axes_grid(f, s; axis, palettes), fig)
    foreach(plot!, grid)
    return grid
end

function Makie.plot(s::OneOrMoreLayers;
                    axis=NamedTuple(), figure=NamedTuple(), palettes=NamedTuple())
    fig = Figure(; figure...)
    grid = plot!(fig, s; axis, palettes)
    return FigureGrid(fig, grid)
end

"""
    draw(s; axis=NamedTuple(), figure=NamedTuple, palettes=NamedTuple())

Draw a [`AlgebraOfGraphics.Layer`](@ref) or [`AlgebraOfGraphics.Layers`](@ref) object `s`.
The output can be customized by giving axis attributes to `axis`, figure attributes
to `figure`, or custom palettes to `palettes`.
Legend and colorbar are drawn automatically. For finer control, use [`draw!`](@ref),
[`legend!`](@ref), and [`colorbar!`](@ref) independently.
"""
function draw(s::OneOrMoreLayers;
              axis=NamedTuple(), figure=NamedTuple(), palettes=NamedTuple(),
              facet=NamedTuple(), legend=NamedTuple(), colorbar=NamedTuple())
    return update(Figure(; figure...)) do f
        grid = plot!(f, s; axis, palettes)
        fg = FigureGrid(f, grid)
        facet!(fg; facet)
        colorbar!(fg; colorbar...)
        legend!(fg; legend...)
        resize_to_layout!(fg)
        return fg
    end
end

"""
    draw!(fig, s; axis=NamedTuple(), palettes=NamedTuple())

Draw a [`AlgebraOfGraphics.Layer`](@ref) or [`AlgebraOfGraphics.Layers`](@ref) object `s` on `fig`.
`fig` can be a figure, a position in a layout, or an axis if `s` has no facet specification.
The output can be customized by giving axis attributes to `axis` or custom palettes
to `palettes`.
"""
function draw!(fig, s::OneOrMoreLayers;
               axis=NamedTuple(), palettes=NamedTuple(), facet=NamedTuple())
    return update(fig) do f
        ag = plot!(f, s; axis, palettes)
        facet!(f, ag; facet)
        return ag
    end
end


text:
module CcallMacros
export @ccall, @cdef, @disable_sigint, @check_syserr

struct NoType end
struct CcallError{T <: AbstractString} <: Exception
    msg::T
end
Base.showerror(io::IO, e::CcallError) = print(io, "CcallError: ", e.msg)

getcfunc(f) = QuoteNode(f)
getcfunc(f::Expr) = :(($(f.args[2]), $(f.args[1])))
hashead(symbol, _) = false
hashead(expr::Expr, head) = expr.head === head

"""
Determine if there are varargs and return a vector of all arguments
based on the call signature.

returns a tuple of `(hasvarargs, arguments)`
"""
function getargs(call)
    # no arguments
    length(call.args) > 2 &&
        return [], []

    firstarg = call.args[2]
    hashead(firstarg, :parameters) &&
        return call.args[3:end], firstarg.args

    return call.args[2:end], []
end

"""
       mkarg(arg)

`mkarg` takes an argument and returns a tuple of symbol and type. The
argument must have a type annotation or an error will be thrown.
"""
function mkarg(arg)
    !hashead(arg, :(::)) &&
        throw(CcallError("args in @ccall need type annotations. '$arg' doesn't have one."))
    return (arg=arg.args[1], type=arg.args[2])
end

"""
get the type specified on the varargs. throw if there is more than one type.
"""
function getvarargtype(varargs)
    vararg_types = Set(a.type for a in varargs)
    length(vararg_types) > 1 &&
        throw(CcallError("varargs @ccall with different argument types not yet supported"))
    return pop!(vararg_types)
end

"""
    parsecall(expression)

`parsecall` is an implementation detail of `@ccall

it takes and expression like `:(printf("%d"::Cstring, value::Cuint)::Cvoid)`
returns: a tuple of `(function_name, return_type, arg_types, args)`

The above input outputs this:

    (:printf, :Cvoid, :((Cstring, Cuint)), ["%d", :value])

Note that the args are in an array, not a quote block and have to be
appended to the ccall in a separate step.
"""
function parsecall(expr::Expr)
    # setup and check for errors
    !hashead(expr, :(::)) &&
        throw(CcallError("@ccall needs a function signature with a return type"))
    rettype = expr.args[2]

    call = expr.args[1]
    !hashead(call, :call) &&
        throw(CcallError("@ccall has to take a function call"))

    # get the function symbols
    func = getcfunc(call.args[1])
    normalargs, varargs = getargs(call)

    # separate annotations from names
    normalargs = mkarg.(normalargs)
    args = Any[a.arg for a in normalargs]
    types = Any[a.type for a in normalargs]
    argtypes = :(())
    argtypes.args = types

    isempty(varargs) && return func, rettype, argtypes, args

    # vararg handling. to be changed if rebased on foreigncall.
    varargs = mkarg.(varargs)
    append!(args, (a.arg for a in varargs))
    vararg_type = getvarargtype(varargs)
    push!(types, :($vararg_type...))
    return func, rettype, argtypes, args
end

"""
    @ccall(call expression)

convert a julia-style function definition to a ccall:

    @ccall printf("%d"::Cstring, 10::Cint)::Cint

same as:

    ccall(:printf, Cint, (Cstring, Cint), "%d", 10)

All arguments must have type annotations and the return type must also
be annotated.

varargs are supported with the following convention:

    @ccall printf("%d, %d, %d"::Cstring ; 1::Cint, 2::Cint, 3::Cint)::Cint

Mind the semicolon. Note that, as with the current ccall API, all
varargs must be of the same type.

Using functions from other libraries is supported by prefixing
the function name with the name of the C library, like this:

    const glib = "libglib-2.0"
    @ccall glib.g_uri_escape_string(
        uri::Cstring, ":/"::Cstring, true::Cint
    )::Cstring

The string literal could also be used directly before the symbol of
the function name, if desired `"libglib-2.0".g_uri_escape_string(...`
"""
macro ccall(expr)
    func, rettype, argtypes, args = parsecall(expr)
    output = :(ccall($func, $rettype, $argtypes))
    append!(output.args, args)
    esc(output)
end

"""
    nolinenum(expr)

remove `LineNumberNodes`
"""
nolinenum(s) = s
nolinenum(e::Expr) =
    Expr(e.head, (nolinenum(a) for a in e.args if !isa(a, LineNumberNode))...)

getsym(arg) = hashead(arg, :(::)) ? arg.args[1] : arg
getmacrocall(expr) = begin
    hashead(expr, :macrocall) ? Tuple(expr.args) : (nothing, nothing, expr)
end

function cdef(funcname, expr)
    macrocall, lnnode, expr = getmacrocall(expr)
    func, rettype, argtypes, args = parsecall(expr)
    realargs = getsym.(args)
    call = :(ccall($func, $realret, $argtypes))
    append!(call.args, realargs)
    if macrocall != nothing
        call = Expr(:macrocall, macrocall, lnnode, call)
    end
    definition = :($funcname())
    append!(definition.args, args)
    esc(:($definition = $call))
end

"""
define a _very_ thin wrapper function on a ccall. Mostly for wrapping
libraries quickly as a foundation for a higher-level interface.

    @cdef mkfifo(path::Cstring, mode::Cuint)::Cint

becomes:

   mkfifo(path, mode) = ccall(:mkfifo, Cint, (Cstring, Cuint), path, mode)
"""
macro cdef(funcname, expr)
    cdef(funcname, expr)
end

macro cdef(expr)
    _, _, inner = getmacrocall(expr)
    func, _, _, _ = parsecall(inner)
    name = func isa QuoteNode ? func.value : func.args[1].value
    cdef(name, expr)
end

"""
disable SIGINT while expr is being executed. Mostly useful for calling
C functions that call back into Julia in a concurrent context because
memory corruption can occur and crash the whole program.
"""
macro disable_sigint(expr)
    out = quote
        disable_sigint() do
            $expr
        end
    end
    esc(out)
end

"""
throw a system error if the expression returns a non-zero exit status.
"""
macro check_syserr(expr, message=nothing)
    if message == nothing
        message = nolinenum(expr) |> string
    end
    return quote
        err = $(esc(expr))
        systemerror($message, err != 0)
    end
end

end # module


text:
using BenchmarkTools
using DataFrames
using PooledArrays
using Random

@show Threads.nthreads()

Random.seed!(1234)
ref_dfi = DataFrame(rand(1:10^4, 10^7, 4), :auto)
ref_dfs = string.(ref_dfi)
ref_dfp = mapcols(PooledArray, ref_dfs)

res = DataFrame(rows=Int[],cols=Int[], type=String[], op=String[], time=Float64[])

for x in (10, 10^6-1, 10^6, 10^7), y in 1:4
    dfi = ref_dfi[1:x, 1:y]
    dfs = ref_dfs[1:x, 1:y]
    dfp = ref_dfp[1:x, 1:y]

    @show (x, y) # ping that the process is alive
    push!(res, (x, y, "integer", "copy", @belapsed DataFrame($dfi)))
    push!(res, (x, y, "string", "copy", @belapsed DataFrame($dfs)))
    push!(res, (x, y, "pooled", "copy", @belapsed DataFrame($dfp)))
    push!(res, (x, y, "integer", ":", @belapsed $dfi[:, :]))
    push!(res, (x, y, "string", ":", @belapsed $dfs[:, :]))
    push!(res, (x, y, "pooled", ":", @belapsed $dfp[:, :]))
    push!(res, (x, y, "integer", "1:end-5", @belapsed $dfi[1:end-5, :]))
    push!(res, (x, y, "string", "1:end-5", @belapsed $dfs[1:end-5, :]))
    push!(res, (x, y, "pooled", "1:end-5", @belapsed $dfp[1:end-5, :]))
    push!(res, (x, y, "integer", "1:5", @belapsed $dfi[1:5, :]))
    push!(res, (x, y, "string", "1:5", @belapsed $dfs[1:1:5, :]))
    push!(res, (x, y, "pooled", "1:5", @belapsed $dfp[1:1:5, :]))
end

res.time *= 1_000

@show Threads.nthreads()
@show unstack(res, [:cols, :type, :op], :rows, :time)


text:
"""
	disjp : example from Bazinga.jl

	Mathematical program with disjunctive constraints, from [Meh19].

	Original formulations:

	minimize      f(x)
	subject to    x ∈ X1 ∪ X2
    with          f(x) := (x1 - 1)^2 + (x2 - 2)^2 + (x3 + 2)^2
                  X1 := { (x1,x2,x3) | x1 ≥ 4, x1 + (x2 - 2)^2 + (x3 + 2)^2 ≥ 5 }
                  X2 := { (x1,x2,x3) | x1^2 + x2^2 ≤ x3, (x1 - 1)^2 + x2^2 + x3 ≥ 1, x2 ≤ 0 }

	Reformulation as a constrained structured problem in the form

	minimize     f(x)
	subject to   c(x) in S

	where
	    c(x) := [ x1 - 4                           ]
       	        [ x1 + (x2 - 2)^2 + (x3 + 2)^2 - 5 ]
   	            [ x3 - x1^2 - x2^2                 ]
       	        [ (x1 - 1)^2 + x2^2 + x3 - 1       ]
                [ - x2                             ]
	    S    := { (c12,c345) | c12 ≥ 0  or  c345 ≥ 0 }

    References:
    [Meh19]     Mehlitz, "A comparison of first-order methods for the numerical
                solution of or-constrained optimization problems" (2019).
                ....
                arXiv:1905.01893v1 [math.OC] 6 May 2019
"""

using Bazinga, OptiMo
using Random, LinearAlgebra
using Printf

###################################################################################
# problem definition
###################################################################################
mutable struct DISJP <: AbstractOptiModel
  meta::OptiModelMeta
end

function DISJP()
    name = "DisjProgExample"
    nvar = 3
    ncon = 5
    meta = OptiModelMeta( nvar, ncon, x0=ones(Float64,nvar), name=name)
    return DISJP( meta )
end

# necessary methods:
# obj, grad!: cons!, jprod!, jtprod!, proj!, prox!, objprox!
function OptiMo.obj( prob::DISJP, x::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    return (x[1] - 1)^2 + (x[2] - 2)^2 + (x[3] + 2)^2
end

function OptiMo.grad!( prob::DISJP, x::AbstractVector, dfx::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    dfx .= [2*(x[1] - 1); 2*(x[2] - 2); 2*(x[3] + 2)]
    return nothing
end

function OptiMo.cons!( prob::DISJP, x::AbstractVector, cx::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x
    OptiMo.@lencheck prob.meta.ncon cx
    cx .= [ x[1] - 4;
            x[1] + (x[2] - 2)^2 + (x[3] + 2)^2 - 5;
            x[3] - x[1]^2 - x[2]^2;
            (x[1] - 1)^2 + x[2]^2 + x[3] - 1;
            - x[2] ]
    return nothing
end

function OptiMo.jprod!( prob::DISJP, x::AbstractVector, v::AbstractVector, Jv::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x v
    OptiMo.@lencheck prob.meta.ncon Jv
    Jv .= [ v[1];
            v[1] + 2*(x[2]-2)*v[2] + 2*(x[3]+2)*v[3];
            -2*x[1]*v[1] - 2*x[2]*v[2] + v[3];
            2*(x[1]-1)*v[1] + 2*x[2]*v[2] + v[3];
            -v[2] ]
    return nothing
end

function OptiMo.jtprod!( prob::DISJP, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x Jtv
    OptiMo.@lencheck prob.meta.ncon v
    Jtv .= [ v[1] + v[2] - 2*x[1]*v[3];
             2*(x[2]-2)*v[2] - 2*x[2]*v[3] + 2*x[2]*v[4] - v[5];
             2*(x[3]+2)*v[2] + v[3] + v[4] ]
    return nothing
end

function OptiMo.prox!( prob::DISJP, x::AbstractVector, a::Real, z::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x
    return nothing
end

function OptiMo.objprox!( prob::DISJP, x::AbstractVector, a::Real, z::AbstractVector )
    OptiMo.@lencheck prob.meta.nvar x z
    z .= x
    return 0.0
end

function OptiMo.proj!( prob::DISJP, cx::AbstractVector, px::AbstractVector )
    OptiMo.@lencheck prob.meta.ncon cx px
    px .= cx
    if any(cx[1:2] .< 0) && any(cx[3:5] .< 0)
        px .= max.(cx, 0)
        if norm(px[1:2]-cx[1:2], 2) < norm(px[3:5]-cx[3:5], 2)
            px[3:5] .= cx[3:5]
        else
            px[1:2] .= cx[1:2]
        end
    end
    return nothing
end

###################################################################################
# problem build
###################################################################################
problem = DISJP()

###################################################################################
# solver build
###################################################################################
solver = Bazinga.ALPX( max_iter=50,
                       max_sub_iter=1000,
                       tol_optim=1e-6,
                       tol_cviol=1e-8,
                       verbose=true,
                       subsolver=:zerofpr,
                       subsolver_verbose=false )

out = solver( problem )
print( out )
x = out.x;

#=
R = eltype( problem.meta.x0 )
nvar = problem.meta.nvar

global c22 = 0
global c44 = 0
global cun = 0

xmin = -5.0
xmax = +5.0

ntests = 1e+3
atol = 1e-4
ndots = round(Int,cbrt(ntests^2))

for i=1:ntests

    x0 = xmin .+ (xmax - xmin) .* rand(R,nvar);

    out = solver( problem, x0=x0 )

    x = out.x;

    if norm(x - [2; -2], Inf) <= atol
        global c22 += 1
    elseif norm(x - [4; 4], Inf) <= atol
        global c44 += 1
    else
        global cun += 1
        @printf "(%6.4f,%6.4f) from (%6.4f,%6.4f)\n" x[1] x[2] x0[1] x0[2]
    end

    @printf "."
    if mod(i,ndots) == 0
        @printf "\n"
    end

end
=#


text:
#Boundary types
struct Periodic<:FlatBoundary
    s::Symbol
    min::Real
    max::Real
    addSymbols::Dict{String,Array{Symbol,1}}
end

function Periodic(s::Symbol,min::Real,max::Real;additional::Array{Symbol,1}=Array{Symbol,1}())
    return Periodic(s,min,max,Dict("add"=>additional))
end

function returnBound_(b::Periodic,p::Program_)

    posm = p.update["Local"][b.s]

    up = quote localVCopy[ic1_,$posm] += $(b.max-b.min) end
    for i in b.addSymbols["add"]
        pos = p.update["Local"][i]
        push!(up.args,:(localVCopy[ic1_,$pos] += $(b.max-b.min)))
    end

    down = quote localVCopy[ic1_,$posm] -= $(b.max-b.min) end
    for i in b.addSymbols["add"]
        pos = p.update["Local"][i]
        push!(down.args,:(localVCopy[ic1_,$pos] -= $(b.max-b.min)))
    end

    code = quote
                if localVCopy[ic1_,$posm] < $(b.min)
                    $up
                elseif localVCopy[ic1_,$posm] >= $(b.max)
                    $down
                end
           end

    return code
end

struct Bound<:NonPeriodicFlat
    s::Symbol
    min::Real
    max::Real
    addSymbols::Dict{String,Array{Symbol,1}}
end

function Bound(s::Symbol,min::Real,max::Real;
                stop::Array{Symbol,1}=Array{Symbol,1}(),
                stopMin::Array{Symbol,1}=Array{Symbol,1}(),
                stopMax::Array{Symbol,1}=Array{Symbol,1}(),
                bounce::Array{Symbol,1}=Array{Symbol,1}(),
                bounceMin::Array{Symbol,1}=Array{Symbol,1}(),
                bounceMax::Array{Symbol,1}=Array{Symbol,1}(),
                reflect::Array{Symbol,1}=Array{Symbol,1}(),
                reflectMin::Array{Symbol,1}=Array{Symbol,1}(),
                reflectMax::Array{Symbol,1}=Array{Symbol,1}()
                )

        d = Dict{String,Array{Symbol,1}}()

        d["stop"]=stop
        d["stopMin"]=stopMin
        d["stopMax"]=stopMax
        d["bounce"]=bounce
        d["bounceMin"]=bounceMin
        d["bounceMax"]=bounceMax
        d["reflect"]=reflect
        d["reflectMin"]=reflectMin
        d["reflectMax"]=reflectMax

    return Bound(s,min,max,d)
end

function returnBound_(b::Bound,p::Program_)

    min = quote end
    for i in b.addSymbols["stop"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundStopMin(localVCopy[ic1_,$pos],$(b.min))))
    end
    for i in b.addSymbols["stopMin"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundStopMin(localVCopy[ic1_,$pos],$(b.min))))
    end
    for i in b.addSymbols["bounce"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundBounceMin(localVCopy[ic1_,$pos],$(b.min))))
    end
    for i in b.addSymbols["bounceMin"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundBounceMin(localVCopy[ic1_,$pos],$(b.min))))
    end
    for i in b.addSymbols["reflect"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundReflect(localVCopy[ic1_,$pos])))
    end
    for i in b.addSymbols["reflectMin"]
        pos = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$pos] = boundReflect(localVCopy[ic1_,$pos])))
    end

    max = quote end
    for i in b.addSymbols["stop"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundStopMax(localVCopy[ic1_,$pos],$(b.max))))
    end
    for i in b.addSymbols["stopMax"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundStopMax(localVCopy[ic1_,$pos],$(b.max))))
    end
    for i in b.addSymbols["bounce"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundBounceMax(localVCopy[ic1_,$pos],$(b.max))))
    end
    for i in b.addSymbols["bounceMax"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundBounceMax(localVCopy[ic1_,$pos],$(b.max))))
    end
    for i in b.addSymbols["reflect"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundReflect(localVCopy[ic1_,$pos])))
    end
    for i in b.addSymbols["reflectMax"]
        pos = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$pos] = boundReflect(localVCopy[ic1_,$pos])))
    end

    pos = p.update["Local"][b.s]

    if !emptyquote_(min) && !emptyquote_(max)

        code = quote
            if localVCopy[ic1_,$pos] < $(b.min)
                $min
            elseif localVCopy[ic1_,$pos] >= $(b.max)
                $max
            end
       end

       return code
    elseif !emptyquote_(min)
        code = quote
            if localVCopy[ic1_,$pos] < $(b.min)
                $min
            end
       end

        return code
    elseif !emptyquote_(max)
        code = quote
            if localVCopy[ic1_,$pos] >= $(b.max)
                $max
            end
       end

        return code
    else

        return quote end
    end

end

export boundStopMin, boundStopMax, boundBounceMin, boundBounceMax, boundReflect

boundStopMin(x,min) = max(x,min)
boundStopMax(x,max) = min(x,max)

boundBounceMin(x,min) = min + (min-x)
boundBounceMax(x,max) = max - (x-max)

boundReflect(x) = -x

function returnBound_(b::Array{<:FlatBoundary,1},p::Program_)

    code = quote end
    for i in b
        if i.s in keys(p.update["Local"])
            codeb = returnBound_(i,p)

            if !emptyquote_(codeb)
                push!(code.args,codeb)
            end
        end
    end

    return code

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
const Rmath = "libRmath-julia"

function rbeta(alpha::Float64, beta::Float64)
  ccall( (:rbeta, Rmath), Float64, (Float64, Float64), alpha, beta )
end

function rbeta(n::Int, alpha::Float64, beta::Float64)
  res = zeros(Float64, n)
  for i in 1:n
    res[i] = rbeta(alpha, beta)
  end
  res
end

#function pbeta(x::Float64, alpha::Float64, beta::Float64)
#  ccall( (:pbeta, Rmath), Float64, (Float64, Float64, Float64, Int64, Int64), x, alpha, beta, 0, 0)
#end

# Incomplete beta function.
#function betainc(x, a, b)
#  pbeta(float(x), float(a), float(b))
#end


text:
using RecipesBase

@recipe function f(tree::Phylo.AbstractTree; treetype = :dendrogram, showtips = true, tipfont = (7,))

    #linecolor --> :black
    grid --> false
    framestyle --> :none
    legend --> false
    colorbar --> true
    size --> (1000, 1000)

    d, h, n = _findxy(tree)
    adj = 0.03maximum(values(d))
    tipannotations = map(x->(d[x] + adj, h[x], x), getleafnames(tree))

    x, y = Float64[], Float64[]
    for node ∈ n
        if hasinbound(tree, node)
            push!(x, d[getparent(tree, node)], d[getparent(tree, node)], d[node], NaN)
            push!(y, h[getparent(tree, node)], h[node], h[node], NaN)
        end
    end

    marker_x, marker_y = Float64[], Float64[]
    if any(x->occursin(r"marker", String(x)), keys(plotattributes))
        f = nodenamefilter(!isleaf, tree)
        append!(marker_x, getindex.(Ref(d), f))
        append!(marker_y, getindex.(Ref(h), f))
    end

    if treetype == :dendrogram
        Dendrogram(x, y, tipannotations, marker_x, marker_y, showtips, tipfont)
    elseif treetype == :fan
        Fan(x, y, tipannotations, marker_x, marker_y, showtips, tipfont)
    else
        throw(ArgumentError("Unsupported `treetype`; valid values are `:dendrogram` or `:fan`"))
    end
end

struct Dendrogram; x; y; tipannotations; marker_x; marker_y; showtips; tipfont; end
struct Fan; x; y; tipannotations; marker_x; marker_y; showtips; tipfont; end

@recipe function f(dend::Dendrogram)

    sa = get(plotattributes, :series_annotations, nothing)
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing

        lc = _extend(get(plotattributes, :linecolor, nothing), dend.x)
        lc !== nothing && (linecolor := lc)
        la = _extend(get(plotattributes, :linealpha, nothing), dend.x)
        la !== nothing && (linealpha := la)
        lz = _extend(get(plotattributes, :line_z, nothing), dend.x)
        lz !== nothing && (line_z := lz)

        dend.x, dend.y
    end
    if !isempty(dend.marker_x) || sa !== nothing
        @series begin
            seriestype := :scatter
            sa !== nothing && (series_annotations := sa)
            dend.marker_x, dend.marker_y
        end
    end
    dend.showtips && (annotations := map(x -> (x[1], x[2], (x[3], :left, dend.tipfont...)), dend.tipannotations))
    [],[]
end

@recipe function f(fan::Fan)
    adjust(y) = 2pi*y / (length(fan.tipannotations) + 1)

    sa = get(plotattributes, :series_annotations, nothing)
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing

        x, y = _circle_transform_segments(fan.x, adjust(fan.y))
        lc = _extend(get(plotattributes, :linecolor, nothing), x)
        lc !== nothing && (linecolor := lc)
        la = _extend(get(plotattributes, :linealpha, nothing), x)
        la !== nothing && (linealpha := la)
        lz = _extend(get(plotattributes, :line_z, nothing), x)
        lz !== nothing && (line_z := lz)
        x, y
    end
    if !isempty(fan.marker_x) || sa !== nothing
        @series begin
            seriestype := :scatter
            a !== nothing && (series_annotations := sa)
            _xcirc.(adjust(fan.marker_y), fan.marker_x), _ycirc.(adjust(fan.marker_y), fan.marker_x)
        end
    end
    aspect_ratio := 1
    mx = maximum(filter(isfinite, fan.x))
    if fan.showtips
        xlim --> (1.5 .* (-mx, mx))
        ylim --> (1.5 .* (-mx, mx))
        annotations := map(x -> (_tocirc(x[1], adjust(x[2]))..., (x[3], :left,
            rad2deg(adjust(x[2])), fan.tipfont...)), fan.tipannotations)
    end
    [],[]
end

function _extend(tmp, x)
    tmp isa AbstractVector && abs(length(tmp) - count(isnan, x)) < 2 || return nothing
    ret = similar(x, eltype(tmp))
    j = 1 + length(tmp) - count(isnan, x)
    for i in eachindex(x)
        ret[i] = tmp[j]
        isnan(x[i]) && (j += 1)
    end
    return ret
end

function Base.sort!(tree::AbstractTree)
    function loc!(clade::String)
        if isleaf(tree, clade)
            return 1
        end

        sizes = map(loc!, getchildren(tree, clade))
        node = getnode(tree, clade)
        node.outbounds .= node.outbounds[sortperm(sizes)]
        sum(sizes) + 1
    end

    loc!(first(nodenamefilter(isroot, tree)))
    tree
end

function _findxy(tree::Phylo.AbstractTree)

    # two convenience recursive functions using captured variables
    function findheights!(clade::String)
        if !in(clade, keys(height))
            for subclade in getchildren(tree, clade)
                findheights!(subclade)
            end
        end
        if !isleaf(tree, clade)
            ch_heights = [height[child] for child in getchildren(tree, clade)]
            height[clade] = (maximum(ch_heights) + minimum(ch_heights)) / 2.
        end
    end

    function finddepths!(clade::String, parentdepth::Float64 = 0.0)
        mydepth = parentdepth
        push!(names, clade)
        if hasinbound(tree, clade)
             mydepth += getlength(tree, getinbound(tree, clade))
        end
        depth[clade] = mydepth
        for ch in getchildren(tree, clade)
            finddepths!(ch, mydepth)
        end
    end

    root = getnodename(tree, getroot(tree))
    height = Dict(tip => float(i) for (i, tip) in enumerate(nodefuture(tree, root)) if isleaf(tree, tip))
    sizehint!(height, nnodes(tree))
    findheights!(root)

    depth = Dict{String, Float64}()
    names = String[]
    sizehint!(depth, nnodes(tree))
    sizehint!(names, nnodes(tree))
    finddepths!(root)

    depth, height, names
end

function _find_tips(depth, height, tree)
    x, y, l = Float64[], Float64[], String[]
    for k in keys(depth)
        if isleaf(tree, k)
            push!(x, depth[k])
            push!(y, height[k])
            push!(l, k)
        end
    end
    x, y, l
end

function _p_circ(start_θ, end_θ, r=1)
    steps = range(start_θ, stop=end_θ, length = 1+ceil(Int, 60abs(end_θ - start_θ)))
    retx = Array{Float64}(undef, length(steps))
    rety = similar(retx)
    for u in eachindex(steps)
        retx[u] = _xcirc(steps[u], r)
        rety[u] = _ycirc(steps[u], r)
    end
    retx, rety
end

_xcirc(x, r) = r*cos(x)
_ycirc(y, r) = r*sin(y)
_tocirc(x, y) = _xcirc(y, x), _ycirc(y, x)

function _circle_transform_segments(xs, ys)
    retx, rety = Float64[], Float64[]
    function _transform_seg(_x, _y)
        tmpx, tmpy = _p_circ(_y[1], _y[2], _x[1])
        append!(retx, tmpx)
        append!(rety, tmpy)
        push!(retx, _xcirc(_y[3], _x[3]), NaN)
        push!(rety, _ycirc(_y[3], _x[3]), NaN)
    end
    i = 1
    while !(i === nothing) && i < length(xs)
        j = findnext(isnan, xs, i) - 1
        _transform_seg(view(xs,i:j), view(ys, i:j))
        i = j + 2
    end
    retx, rety
end

# a function to update a value successively from the root to the tips
function map_depthfirst(FUN, start, tree, eltype = nothing)
    root = first(nodenamefilter(isroot, tree))
    eltype === nothing && (eltype = typeof(FUN(start, root)))
    ret = Vector{eltype}()
    function local!(val, node)
        push!(ret, val)
        for ch in getchildren(tree, node)
            local!(FUN(val, node), ch)
        end
    end
    local!(start, root)
    ret
end


text:
const kMinimalTargetExponent = -60
const kMaximalTargetExponent = -32

function roundweed(buffer,len,rest,tk,unit,kappa,too_high::Uint64,unsafe_interval::Uint64)
    small = too_high - unit
    big = too_high + unit
    while rest < small &&
            unsafe_interval - rest >= tk &&
            (rest + tk < small ||
            small - rest >= rest + tk - small)
        buffer[len-1] -= 1
        rest += tk
    end
    if rest < big &&
        unsafe_interval - rest >= tk &&
        (rest + tk < big ||
        big - rest > rest + tk - big)
        return false, kappa
    end
    return (2 * unit <= rest) && (rest <= unsafe_interval - 4 * unit), kappa
end

function roundweed(buffer,len,rest,tk,unit,kappa,null1::Int=0,null2::Int=0)
    unit >= tk && return false, kappa
    tk - unit <= unit && return false, kappa
    tk - rest > rest && (tk - 2 * rest >= 2 * unit) && return true, kappa
    if rest > unit && (tk - (rest - unit) <= (rest - unit))
        buffer[len-1] += 1
        for i = (len-1):-1:2
            buffer[i] != 0x30 + 10 && break
            buffer[i] = 0x30
            buffer[i-1] += 1
        end
        if buffer[1] == 0x30 + 10
            buffer[1] = 0x31
            kappa += 1
        end
        return true, kappa
    end
    return false, kappa
end

const SmallPowersOfTen = [
        0, 1, 10, 100, 1000, 10000, 100000,
        1000000, 10000000, 100000000, 1000000000]

function BiggestPowerTen(n,n_bits)
    guess = (n_bits + 1) * 1233 >> 12
    guess += 1
    i = SmallPowersOfTen[guess+1]
    return n < i ? (SmallPowersOfTen[guess], guess-1) : (i,guess)
end

function DigitGen(low::Float,w::Float,high::Float,
                  buffer,requested_digits=1000,shortest=true)
    unit::Uint64 = 1
    one = Float(unit << -w.e, w.e)
    if shortest
        too_low = Float(low.s-unit,low.e)
        too_high = Float(high.s+unit,high.e)
        unsafe_interval = too_high - too_low
        ref = too_high
    else
        ref = w
        unsafe_interval = w
    end
    integrals = ref.s >> -one.e
    fractionals = ref.s & (one.s-1)
    divisor, divisor_exp_plus_one = BiggestPowerTen(integrals, 64 + one.e)
    kappa = divisor_exp_plus_one
    len = 1
    rest = 0
    check = false
    while kappa > 0
        digit = div(integrals,divisor)
        buffer[len] = 0x30 + digit
        len += 1
        requested_digits -= 1
        integrals %= divisor
        kappa -= 1
        rest = (uint64(integrals) << -one.e) + fractionals
        check = shortest ? (rest < unsafe_interval.s) : (requested_digits == 0)
        check && break
        divisor = div(divisor,10)
    end
    if check
        a,b = shortest ? ((too_high - w).s, unsafe_interval.s) : (0,0)
        r, kappa = roundweed(buffer, len, rest, uint64(divisor) << -one.e,
                    unit,kappa,a,b)
        return r, kappa, len
    end
    while true
        !shortest && (requested_digits <= 0 || fractionals < unit) && break
        fractionals *= 10
        unit *= 10
        unsafe_interval = Float(unsafe_interval.s*10,unsafe_interval.e)
        digit = fractionals >> -one.e
        buffer[len] = 0x30 + digit
        len += 1
        requested_digits -= 1
        fractionals &= one.s - 1
        kappa -= 1
        shortest && fractionals < unsafe_interval.s && break
    end
    if shortest
        a,b = (too_high - w).s*unit,unsafe_interval.s
    else
        requested_digits != 0 && return false, kappa, len
        a,b = 0, 0
    end
    r, kappa = roundweed(buffer,len,fractionals,one.s,
                         unit,kappa,a,b)
    return r, kappa, len
end

function fastdtoa(v,mode,requested_digits=1000,buffer=Array(Uint8,100))
    f = normalize(float64(v))
    shortest = mode != PRECISION
    bound_minus, bound_plus = normalizedbound(v) #
    ten_mk_minimal_binary_exponent =
        kMinimalTargetExponent - (f.e + FloatSignificandSize) #e
    ten_mk_maximal_binary_exponent =
        kMaximalTargetExponent - (f.e + FloatSignificandSize) #e
    cp = GetCachedPowerForBinaryExponentRange(
           ten_mk_minimal_binary_exponent,
           ten_mk_maximal_binary_exponent)
    scaled_w = f * cp #scale
    scaled_bound_minus = bound_minus * cp
    scaled_bound_plus = bound_plus * cp
    r, kappa, len = DigitGen(scaled_bound_minus,scaled_w,
                             scaled_bound_plus,buffer,
                             requested_digits,shortest)
    decimal_exponent = -cp.de + kappa
    return r, len, decimal_exponent+len-1, buffer
end


text:
# functions: return_boundary_nodes,
# Matern_3d_Grid, Laplace_3D_grid,
# Parallel_Matern_3DGrid

"""
  return_boundary_nodes(xpoints, ypoints, zpoints)

...
# Arguments

  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
...

...
# Outputs
  - `BoundaryNodes3D::Vector{Int64}`: vector containing the indices of coordinates
  on the boundary of the rectangular 3D volume
...

"""
function return_boundary_nodes(xpoints, ypoints, zpoints)
    BoundaryNodes3D =[]
    xneighbors = []
    yneighbors = []
    zneighbors = []
    counter = 0
    for k = 1:zpoints
        for j = 1:ypoints
            for i = 1:xpoints
                counter=counter+1
                if(k == 1 || k == zpoints || j == 1|| j == ypoints || i == 1 || i == xpoints)
                    BoundaryNodes3D = push!(BoundaryNodes3D, counter)
                    if(k == 1 || k == zpoints)
                        push!(zneighbors, 1)
                    else
                        push!(zneighbors, 2)
                    end
                    if(j == 1 || j == ypoints)
                        push!(yneighbors, 1)
                    else
                        push!(yneighbors, 2)
                    end
                    if(i == 1 || i == xpoints)
                        push!(xneighbors, 1)
                    else
                        push!(xneighbors, 2)
                    end
                end
            end
        end
    end
    return BoundaryNodes3D, xneighbors, yneighbors, zneighbors
end

"""

  Matern3D_Grid(xpoints, ypoints, zpoints, imgg, epsilon, radius, h, k, l, m)

...
# Arguments
  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg`: the matrix containing the image
  - `epsilon`: Matern parameter epsilon
  - `radius::Vector`: the tuple containing the punch radii
  - `h::Float`: grid spacing along the x-axis
  - `k::Float`: grid spacing along the y-axis
  - `l::Float`: grid spacing along the z-axis
  - `m::Int` : Matern parameter

# Outputs
  - tuple containing the restored image and the punched image.
...

"""
function Matern3D_Grid(xpoints, ypoints, zpoints, imgg, epsilon, radius, h, k, l, m)
    A3D = nablasq_3d_grid(length(xpoints), length(ypoints), length(zpoints), h, k, l)
    sizeA = size(A3D, 1)
    for i = 1:sizeA
        A3D[i, i] = A3D[i, i] + epsilon^2
    end
    A3DMatern = A3D
    for i = 1:m - 1
        A3DMatern = A3DMatern * A3D
    end
    discard = punch_holes_nexus(xpoints, ypoints, zpoints, radius)
    punched_image = copy(imgg)
    punched_image[discard] .= 1
    totalsize = prod(size(imgg))
    C = sparse(I, totalsize, totalsize)
    rhs_a = punched_image[:]
    for i in discard
        C[i, i] = 0
        rhs_a[i] = 0
    end
    Id = sparse(I, totalsize, totalsize)
    u = ((C - (Id - C) * A3DMatern)) \ rhs_a
    return (u, punched_image[:])
end

"""

  Laplace3D_Grid(xpoints, ypoints, zpoints, imgg, radius, h, k, l)

...
# Arguments
  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg`: the matrix containing the image
  - `radius::Vector`: the tuple containing the punch radii
  - `h::Float`: grid spacing along the x-axis
  - `k::Float`: grid spacing along the y-axis
  - `l::Float`: grid spacing along the z-axis

# Outputs
  - tuple containing the restored image and the punched image.
...

"""
function Laplace3D_Grid(xpoints, ypoints, zpoints, imgg, radius, h, k, l)
    A3D = nablasq_3d_grid(length(xpoints), length(ypoints), length(zpoints), h, k, l)
    discard = punch_holes_nexus(xpoints, ypoints, zpoints, radius)
    punched_image = copy(imgg)
    punched_image[discard] .= 1
    totalsize = prod(size(imgg))
    C = sparse(I, totalsize, totalsize)
    rhs_a = punched_image[:]
    for i in discard
        C[i,i] = 0
        rhs_a[i] = 0
    end
    Id = sparse(I, totalsize, totalsize)
    u = ((C - (Id - C) * A3D)) \ rhs_a
    return u, punched_image[:]
end

"""

  Parallel_Laplace3D_Grid(xpoints, ypoints, zpoints, imgg, radius, h, k, l,
          xmin, xmax, ymin, ymax, zmin, zmax)

Compute the spherically-punched, Laplace-interpolated 3D data

...
# Arguments
  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg`: the matrix containing the image
  - `radius::Vector`: the tuple containing the punch radii
  - `h::Float`: grid spacing along the x-axis
  - `k::Float`: grid spacing along the y-axis
  - `l::Float`: grid spacing along the z-axis
  - `xmin::Int64`: Vishwas should fill in the next six fields.
  - `xmax::Int64`:
  - `ymin::Int64`:
  - `ymax::Int64`:
  - `zmin::Int64`:
  - `zmax::Int64`:

# Outputs
  - array containing the restored image

# Example

```julia-repl
(xmin, xmax) = (ymin, ymax) = (zmin,zmax) = (0.0, 1.0)
xpoints = ypoints = zpoints = -0.2:0.2:1.2
h = k = l = 0.2
imgg = randn(8,8,8)
radius = 0.2
restored = Parallel_Laplace3D_Grid(xpoints, ypoints, zpoints, imgg, radius,
                                 h, k, l, xmin, xmax, ymin, ymax, zmin, zmax)
```

...

"""
function Parallel_Laplace3D_Grid(xpoints::Union{StepRangeLen{T},Vector{T}},
                                ypoints::Union{StepRangeLen{T},Vector{T}},
                                zpoints::Union{StepRangeLen{T},Vector{T}},
                                imgg::Array{P,3},
                                radius::Union{Q,Tuple{Q,Q,Q}},
                                h::Float64, k::Float64, l::Float64,
                                xmin::R, xmax::R, ymin::R, ymax::R, zmin::R, zmax::R
                                ) where{T<:Number,P<:Number,Q<:Number,R<:Number}
  #
  fun(x,y,z,w) = Int(round((x -y)/z) - w )
  ran(x,y,z,w) = fun(x, y, z, w) .+ (0, 2*w + 1)
  #
  rad = (typeof(radius) <: Tuple) ? radius : (radius, radius, radius)
  radius_x, radius_y, radius_z = rad
  stride_h = Int64(round(radius_x / h))
  stride_k = Int64(round(radius_y / k))
  stride_l = Int64(round(radius_z / l))
  cartesian_product_boxes = [(ran(i, xpoints[1], h, stride_h)...,
                              ran(j, ypoints[1], k, stride_k)...,
                              ran(kk, zpoints[1], l, stride_l)...)
                         for i in xmin:xmax for j in ymin:ymax for kk in zmin:zmax]
  #
  z3d_restored = copy(imgg)
  #
  Threads.@threads for i = 1:length(cartesian_product_boxes)
    i1, i2, j1, j2, k1, k2 = cartesian_product_boxes[i]
    restored_img = Laplace3D_Grid(xpoints[i1 + 1:i2], ypoints[j1 + 1:j2],
                           zpoints[k1 + 1:k2],
                           imgg[i1 + 1:i2, j1 + 1:j2, k1 + 1:k2],
                           radius, h, k, l)[1]
    z3d_restored[i1 + 1:i2, j1 + 1:j2, k1 + 1:k2] = reshape(restored_img,
                           (2 * stride_h + 1, 2 * stride_k + 1, 2 * stride_l + 1))
  end
  #
  return z3d_restored[:]
  #
end

"""

  Parallel_Matern3D_Grid(xpoints, ypoints, zpoints, imgg, epsilon, radius, h, k, l,
          xmin, xmax, ymin, ymax, zmin, zmax, m)

Compute the spherically-punched, Matern-interpolated 3D data

...
# Arguments
  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg`: the matrix containing the image
  - `epsilon`: one of the matern parameters
  - `radius::Vector`: the tuple containing the punch radii
  - `h::Float`: grid spacing along the x-axis
  - `k::Float`: grid spacing along the y-axis
  - `l::Float`: grid spacing along the z-axis
  - `xmin::Int64`: Vishwas should fill in the next six fields.
  - `xmax::Int64`:
  - `ymin::Int64`:
  - `ymax::Int64`:
  - `zmin::Int64`:
  - `zmax::Int64`:
  - `m::Int64`: The matern order parameter

# Outputs
  - array containing the restored image

# Example

```julia-repl
(xmin, xmax) = (ymin, ymax) = (zmin,zmax) = (0.0, 1.0)
xpoints = ypoints = zpoints = -0.2:0.2:1.2
h = k = l = 0.2
imgg = randn(8,8,8)
m = 2
epsilon = 0.0
radius = 0.2
restored = Parallel_Matern3D_Grid(xpoints, ypoints, zpoints, imgg, epsilon, radius,
                                  h, k, l, xmin, xmax, ymin, ymax, zmin, zmax, m)
```

...

"""
function Parallel_Matern3D_Grid(xpoints::Union{StepRangeLen{T},Vector{T}},
                                ypoints::Union{StepRangeLen{T},Vector{T}},
                                zpoints::Union{StepRangeLen{T},Vector{T}},
                                imgg::Array{P,3}, epsilon::Q,
                                radius::Union{Q,Tuple{Q,Q,Q}},
                                h::Float64, k::Float64, l::Float64,
                                xmin::R, xmax::R, ymin::R, ymax::R, zmin::R, zmax::R,
                                m::Int) where{T<:Number,P<:Number,Q<:Number,R<:Number}
  #
  fun(x,y,z,w) = Int(round((x -y)/z) - w )
  ran(x,y,z,w) = fun(x, y, z, w) .+ (0, 2*w + 1)
  #
  rad = (typeof(radius) <: Tuple) ? radius : (radius, radius, radius)
  radius_x, radius_y, radius_z = rad
  stride_h = Int64(round(radius_x / h))
  stride_k = Int64(round(radius_y / k))
  stride_l = Int64(round(radius_z / l))
  cartesian_product_boxes = [(ran(i, xpoints[1], h, stride_h)...,
                              ran(j, ypoints[1], k, stride_k)...,
                              ran(kk, zpoints[1], l, stride_l)...)
                         for i in xmin:xmax for j in ymin:ymax for kk in zmin:zmax]
  #
  z3d_restored = copy(imgg)
  #
  Threads.@threads for i = 1:length(cartesian_product_boxes)
    i1, i2, j1, j2, k1, k2 = cartesian_product_boxes[i]
    restored_img = Matern3D_Grid(xpoints[i1 + 1:i2], ypoints[j1 + 1:j2],
                           zpoints[k1 + 1:k2],
                           imgg[i1 + 1:i2, j1 + 1:j2, k1 + 1:k2],
                           epsilon, radius, h, k, l, m)[1]
    z3d_restored[i1 + 1:i2, j1 + 1:j2, k1 + 1:k2] = reshape(restored_img,
                           (2 * stride_h + 1, 2 * stride_k + 1, 2 * stride_l + 1))
  end
  #
  return z3d_restored[:]
  #
end

# function Parallel_Matern3D_Grid(xpoints, ypoints, zpoints, imgg, epsilon, radius, h, k, l, xmin, xmax, ymin, ymax, zmin, zmax, m)
#   xbegin = xpoints[1];
#   ybegin = ypoints[1];
#   zbegin = zpoints[1];
#   cartesian_product_boxes = [];
#   stride = Int(round(radius/h));
#   z3d_restored = copy(imgg);
#   for i = xmin:xmax+1
#       i1 = Int(round((i-xbegin) /h))-stride;
#       i2 = i1+2*stride+1;
#       for j = ymin:ymax+1
#           j1 = Int(round((j-ybegin)/h))-stride;
#           j2 = j1+2*stride+1;
#           for k = zmin:zmax+1
#               k1 = Int(round((k-ybegin)/h)) - stride;
#               k2 = k1+2*stride+1;
#               append!(cartesian_product_boxes,[(i1,i2,j1,j2,k1,k2)]);
#           end
#       end
#   end

#   Threads.@threads for i = 1:length(cartesian_product_boxes)
#     i1 = cartesian_product_boxes[i][1];
#     i2 = cartesian_product_boxes[i][2];
#     j1 = cartesian_product_boxes[i][3];
#     j2 = cartesian_product_boxes[i][4];
#     k1 = cartesian_product_boxes[i][5];
#     k2 = cartesian_product_boxes[i][6];
#     z3temp = imgg[i1+1:i2,j1+1:j2,k1+1:k2];
#     restored_img, punched_image = Matern3D_Grid(xpoints[i1+1:i2], ypoints[j1+1:j2], zpoints[k1+1:k2], z3temp, epsilon, radius, h, k, l, m);
#     restored_img_reshape = reshape(restored_img, (2*stride+1,2*stride+1,2*stride+1));
#     z3d_restored[i1+1:i2, j1+1:j2, k1+1:k2] = restored_img_reshape;
#   end
#   return z3d_restored[:];
# end

max_stars_repo_name: rayosborn/nxpeaks

max_issues_repo_name: rayosborn/nxpeaks

max_forks_repo_name: axmas-anl/nxrefine


text:
export estimateCoilSensitivities, mergeChannels, espirit, estimateCoilSensitivitiesFixedPoint

"""
    `s = estimateCoilSensitivities(I::AbstractArray{T,5})`

Estimates the coil sensitivity based on a reconstruction where the data
from each coil has been reconstructed individually.
Returns a 5D array.
"""
function estimateCoilSensitivities(I::AbstractArray{T,5}, thresh=1.e-2) where T
  nx,ny,nz,ne,numChan = size(I)

  I_sum = sqrt.( sum(abs.(I).^2, dims=5) ).+eps()
  I_max = maximum(abs.(I_sum))
  msk = zeros(size(I_sum))
  msk[findall(x->x>thresh*I_max,I_sum)] .= 1

  s = zeros(eltype(I),size(I))
  for i=1:numChan
      s[:,:,:,:,i] = msk .* I[:,:,:,:,i] ./ I_sum
  end

  return s
end

"""
    `I4 = mergeChannels(I::AbstractArray{T,5})`

Merge the channels of a multi-coil reconstruction.
Returns a 4D array.
"""
mergeChannels(I::AbstractArray{T,5}) where T = sqrt.(sum(abs.(I).^2,dims=5))

"""
    `maps = espirit(acqData::AcquisitionData, ksize::NTuple{2,Int64}, ncalib::Int64
               ; eigThresh_1::Float64=0.02, eigThresh_2::Float64=0.95)`

obtains coil sensitivities from a calibration area using ESPIRiT
adapted from the MATLAB code by Uecker et al. for the paper'
M. Uecker, P. Lai, MJ Murphy, P. Virtue, M Elad, JM Pauly, SS Vasanawala and M Lustig, "ESPIRiT- an
eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA", Magn Reson Med, 2013

the matlab code can be found at: [http://people.eecs.berkeley.edu/~mlustig/Software.html]

# Arguments
* `acqData::AcquisitionData`  - AcquisitionData
* `ksize::NTuple{2,Int64}`    - size of the k-space kernel
* `ncalib::Int64`             - number of calibration points in each dimension
* `eigThresh_1::Float64=0.02` - threshold for the singular values of the calibration matrix (relative to the largest value)
* `eigThresh_2::Float64=0.95` - threshold of the image space kernels (if no singular value > `eigThresh_2` exists)
                                , the corresponding pixel has a sensitivity of 0.

Returns a 4D array.
"""
function espirit(acqData::AcquisitionData, ksize::NTuple{2,Int64}, ncalib::Int64
                ; eigThresh_1::Float64=0.02, eigThresh_2::Float64=0.95)

  if !isCartesian(trajectory(acqData,1))
    @error "espirit does not yet support non-cartesian sampling"
  end

  nx,ny = acqData.encodingSize[1:2]
  numChan, numSl = numChannels(acqData), numSlices(acqData)
  maps = zeros(ComplexF64,acqData.encodingSize[1],acqData.encodingSize[2],numSl,numChan)

  for slice = 1:numSl
    # form zeropadded array with kspace data
    kdata = zeros(ComplexF64,nx*ny,numChan)
    for coil = 1:numChan
      kdata[acqData.subsampleIndices[1],coil] .= kData(acqData,1,coil,slice)
    end
    kdata = reshape(kdata,nx,ny,numChan)

    calib = crop(kdata,(ncalib,ncalib,numChan))

    maps[:,:,slice,:] .= espirit(calib,(nx,ny),ksize,eigThresh_1=eigThresh_1,eigThresh_2=eigThresh_2)
  end
  return maps
end

"""
`maps = espirit(calibData::Array{ComplexF64,3}, imsize::NTuple{2,Int64}, ksize::NTuple{2,Int64}
                ; eigThresh_1::Float64=0.02, eigThresh_2::Float64=0.95)`
Returns a 3D array.
"""
function espirit(calibData::Array{ComplexF64,3}, imsize::NTuple{2,Int64}, ksize::NTuple{2,Int64}
                ; eigThresh_1::Float64=0.02, eigThresh_2::Float64=0.95)
  sx,sy = imsize
  nc = size(calibData,3)

  # compute calibration matrix, perform SVD and convert right singular vectors
  # into k-space kernels
  k, S = dat2Kernel(calibData,ksize)
  idx = findlast(x->x>=S[1]*eigThresh_1,S) # last singular vector to keep

  # crop kernels and compute eigenvalue decomposition in images space
  # to get sensitivity maps
  M, W = kernelEig(k[:,:,:,1:idx],(sx,sy))

  # sensitivity maps correspond to eigen vectors with a singular value of 1
  msk = zeros(size(W,1),size(W,2),1)
  msk[findall(x->x>eigThresh_2, abs.(W[:,:,end]))] .= 1
  maps = M[:,:,:,end] .* repeat(msk,1,1,nc)

  return maps
end

#
# computes ESPIRiT step I:
#   form calibration matrix from data, perform SVD and convert right singular vectors
#   into k-space kernels
#
function dat2Kernel(data::Array{T,3}, ksize::NTuple{2,Int64}) where T
  sx,sy,nc = size(data)
  imSize=(sx,sy)

  tmp = im2row(data,ksize)
  tsx,tsy,tsz = size(tmp)
  A = reshape(tmp, tsx, tsy*tsz)
  usv = svd(A)
  kernel = reshape(usv.V,ksize[1],ksize[2],nc,size(usv.V,2))
  return kernel, usv.S
end

# computes ESPIRiT step II:
#   eigenvalue decomposition of a k-space kernel in image space.
#
# inputs:
#         kernel - kspace kernels (4D)
#         imSize - size of the image to compute maps
#
# outputs :
#         eigenvecs: images representing the eigenvectors (sx,sy,numChan,numChan)
#         eigenvals: images representing th eigenvalues (sx,sy,numChan)
function kernelEig(kernel::Array{T,4},imsize::NTuple{2,Int64}) where T
  nx,ny,nc,nv = size(kernel)
  ksize = (nx,ny)

  # rotate kernel to order by maximum variance
  k = permutedims(kernel,[1,2,4,3])
  k = reshape(k,nx*ny*nv,nc)

  if size(k,1) < size(k,2)
    usv = svd(k,full=true)
  else
    usv = svd(k)
  end

  k = k*usv.V
  kernel = reshape(k, nx,ny,nv,nc)
  kernel = permutedims(kernel,[1,2,4,3])
  kern2 = zeros(ComplexF64, imsize[1],imsize[2],size(kernel,3),size(kernel,4))
  for n = 1:size(kernel,4)
    kern2[:,:,:,n] .= fft2c( zeropad( conj.(kernel[end:-1:1,end:-1:1,:,n])*sqrt(imsize[1]*imsize[2]), (imsize[1],imsize[2],size(kernel,3)) ) )
  end
  kern2 = kern2/sqrt(prod(ksize))

  eigenVecs = zeros(ComplexF64,imsize[1],imsize[2],nc,min(nc,nv))
  eigenVals = zeros(ComplexF64,imsize[1],imsize[2],min(nc,nv))

  for n = 1:prod(imsize)
    x,y = Tuple( CartesianIndices((imsize[1],imsize[2]))[n] )
    mtx = kern2[x,y,:,:]

    cdv = svd(mtx)
    ph = transpose( repeat( exp.(-1im*angle.(cdv.U[1,:])), 1,size(cdv.U,1) ) )
    C = usv.V*(cdv.U .* ph)
    D = real.(cdv.S)
    eigenVals[x,y,:] .= D[end:-1:1]
    eigenVecs[x,y,:,:] .= C[:,end:-1:1]
  end

  return eigenVecs, eigenVals
end

#
# rearrange data from the image data to form kernels for dat2Kernel
#
function im2row(img::Array{T,3}, winSize::NTuple{2,Int64}) where T
  sx,sy,sz = size(img)
  res = zeros(T,(sx-winSize[1]+1)*(sy-winSize[2]+1),prod(winSize),sz)
  cnt = 0
  for y = 1:winSize[2]
    for x = 1:winSize[1]
      cnt += 1
      res[:,cnt,:] .= reshape(img[x:sx-winSize[1]+x,y:sy-winSize[2]+y,:]
                              ,(sx-winSize[1]+1)*(sy-winSize[2]+1),sz)
    end
  end
  return res
end

#
# crop the central area (of size s) of an array A
#
function crop(A::Array{T,3},s::NTuple{3,Int64}) where T
  nx,ny,nz = size(A)
  idx_x = div(nx,2)-div(s[1],2)+1:div(nx,2)-div(s[1],2)+s[1]
  idx_y = div(ny,2)-div(s[2],2)+1:div(ny,2)-div(s[2],2)+s[2]
  idx_z = div(nz,2)-div(s[3],2)+1:div(nz,2)-div(s[3],2)+s[3]
  return A[idx_x,idx_y,idx_z]
end

#
# zeropoad an array A to a target size s
#
function zeropad(A::Array{T,3},s::NTuple{3,Int64}) where T
  nx,ny,nz = size(A)
  idx_x = div(s[1],2)-div(nx,2)+1:div(s[1],2)-div(nx,2)+nx
  idx_y = div(s[2],2)-div(ny,2)+1:div(s[2],2)-div(ny,2)+ny
  idx_z = div(s[3],2)-div(nz,2)+1:div(s[3],2)-div(nz,2)+nz

  res = zeros(T,s[1],s[2],s[3])
  res[idx_x,idx_y,idx_z] .= A

  return res
end

#
#  apply 2d fft to all slices of a higher dimensional array
#
function fft2c(x::Array{T}) where T
  s = size(x)
  x = reshape(x,s[1],s[2],:)
  res = zeros(ComplexF64,size(x))
  for n=1:size(x,3)
    res[:,:,n] = 1/sqrt(s[1]*s[2])*fftshift(fft(ifftshift(x[:,:,n])))
  end
  return reshape(res,s)
end

#
#  apply 2d ifft to all slices of a higher dimensional array
#
function ifft2c(x::Array{T}) where T
  s = size(x)
  x = reshape(x,s[1],s[2],:)
  res = zeros(ComplexF64,size(x))
  for n=1:size(x,3)
    res[:,:,n] = sqrt(s[1]*s[2])*fftshift(fft(ifftshift(x[:,:,n])))
  end
  return reshape(res,s)
end

function estimateCoilSensitivitiesFixedPoint(acqData::AcquisitionData;
                                             iterations=1, outerIterations=3, cmap = nothing)

  N = acqData.encodingSize

  params = Dict{Symbol, Any}()
  params[:reco] = "standard" # "multiCoil"
  #params[:solver] = "cgnr" #solver
  #params[:iterations] = iterations
  params[:reconSize] = (N,N)
  params[:alpha] = 1.25
  #params[:iterationsInner] = 5
  params[:normalizeReg] = true

  if cmap != nothing
    params[:correctionMap] = cmap
    params[:K] = 15
    params[:m] = 3.0
  end

  Ireco = reconstruction(acqData, params)
  s = estimateCoilSensitivities(IrecoCorr).data

  #if !isempty(coilsens)
  #  params[:senseMaps] = coilsens
  #end

end


text:
# using Pkg
# Pkg.activate(@__DIR__)
# Pkg.instantiate()

using CompScienceMeshes, BEAST
using LinearAlgebra

# Γ = meshcuboid(1.0, 1.0, 1.0, 0.25)
Γ = meshsphere(1.0, 0.3)
X = raviartthomas(Γ)
Y = buffachristiansen(Γ, ibscaled=true)

sol = 1.0
T = TDMaxwell3D.singlelayer(speedoflight=sol, numdiffs=1)
K = TDMaxwell3D.doublelayer(speedoflight=sol, numdiffs=1)

Δt, Nt = 0.6, 200
δ = timebasisdelta(Δt, Nt)
T0 = timebasiscxd0(Δt, Nt)
T1 = timebasisshiftedlagrange(Δt,Nt,1)
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
T3 = timebasisshiftedlagrange(Δt, Nt, 3)

duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)
H = direction × E

@hilbertspace j m
@hilbertspace k l

pmchwt = @discretise(
    2.0T[k,j] + 2.0K[k,m] -
    2.0K[l,j] + 2.0T[l,m] == H[k] + E[l],
    j∈X⊗T2, m∈Y⊗T2, k∈X⊗δ, l∈Y⊗δ)

u = solve(pmchwt)

# nX = numfunctions(X)
# uj = u[1:nX]
# um = u[nX+1:end]

# Θ, Φ = range(0.0,stop=2π,length=100), 0.0
# ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
#
# # Don't forgt the far field comprises two contributions
# κ, η = 1.0, 1.0
# ffm = potential(MWFarField3D(κ*im), ffpoints, um, X)
# ffj = potential(MWFarField3D(κ*im), ffpoints, uj, X)
# ff = η*im*κ*ffj + im*κ*cross.(ffpoints, ffm)
#
# using Plots
# plot(xlabel="theta")
# plot!(Θ,norm.(ffm),label="far field")
#
# import PlotlyJS
# using LinearAlgebra
# fcrj, _ = facecurrents(uj,X)
# PlotlyJS.plot(patch(Γ, norm.(fcrj)))


text:
using Test
using LyapunovExponents: is_semi_unitary

@testset "Utils: is_semi_unitary" begin
    @test is_semi_unitary(eye(3))
    @test ! is_semi_unitary(zeros(Int, 2, 2))
    @test ! is_semi_unitary(ones(Int, 2, 2))
    @test is_semi_unitary(qr([[1, 2] [2, 1]])[1])
end


text:
fsort(x::AbstractVector{Union{T, Missing}}; rev = false) where T = begin
    nmissing = sum(ismissing, x)
    cx = similar(x)
    cx[1:length(x) - nmissing] .= fsort(collect(skipmissing(x)); rev = rev)
    cx[length(x)+1:end] .= missing
    cx
end


text:
using .MPI
import Libdl

# Check whether the HDF5 libraries were compiled with parallel support.
HAS_PARALLEL[] = Libdl.dlopen(libhdf5) do lib
    Libdl.dlsym(lib, :H5Pget_fapl_mpio, throw_error=false) !== nothing
end

# Low-level MPI handles.
const MPIHandle = Union{MPI.MPI_Comm, MPI.MPI_Info}

# MPI.jl wrapper types.
const MPIHandleWrapper = Union{MPI.Comm, MPI.Info}

const H5MPIHandle = let csize = sizeof(MPI.MPI_Comm)
    @assert csize in (4, 8)
    csize == 4 ? Hmpih32 : Hmpih64
end

h5_to_mpi(handle::H5MPIHandle) = reinterpret(MPI.MPI_Comm, handle)

mpi_to_h5(handle::MPIHandle) = reinterpret(H5MPIHandle, handle)
mpi_to_h5(mpiobj::MPIHandleWrapper) = mpi_to_h5(mpiobj.val)

# Set MPIO properties in HDF5.
# Note: HDF5 creates a COPY of the comm and info objects.
function h5p_set_fapl_mpio(fapl_id, comm::MPI.Comm, info::MPI.Info)
    h5comm = mpi_to_h5(comm)
    h5info = mpi_to_h5(info)
    h5p_set_fapl_mpio(fapl_id, h5comm, h5info)
end

h5p_set_fapl_mpio(fapl_id, comm::Hmpih32, info::Hmpih32) =
    h5p_set_fapl_mpio32(fapl_id, comm, info)
h5p_set_fapl_mpio(fapl_id, comm::Hmpih64, info::Hmpih64) =
    h5p_set_fapl_mpio64(fapl_id, comm, info)

# Retrieves the copies of the comm and info MPIO objects from the HDF5 property list.
function h5p_get_fapl_mpio(fapl_id)
    h5comm, h5info = h5p_get_fapl_mpio(fapl_id, H5MPIHandle)
    comm = MPI.Comm(h5_to_mpi(h5comm))
    info = MPI.Info(h5_to_mpi(h5info))
    comm, info
end

function h5p_get_fapl_mpio(fapl_id, ::Type{h}) where {h}
    comm, info = Ref{h}(), Ref{h}()
    h5p_get_fapl_mpio(fapl_id, comm, info)
    comm[], info[]
end

h5p_get_fapl_mpio(fapl_id, comm::Ref{Hmpih32}, info::Ref{Hmpih32}) =
    h5p_get_fapl_mpio32(fapl_id, comm, info)
h5p_get_fapl_mpio(fapl_id, comm::Ref{Hmpih64}, info::Ref{Hmpih64}) =
    h5p_get_fapl_mpio64(fapl_id, comm, info)

function check_hdf5_parallel()
    has_parallel() && return
    error(
        "HDF5.jl has no parallel support." *
        " Make sure that you're using MPI-enabled HDF5 libraries, and that" *
        " MPI was loaded before HDF5." *
        " See HDF5.jl docs for details."
    )
end

"""
    h5open(filename, [mode="r"], comm::MPI.Comm, [info::MPI.Info]; pv...)

Open or create a parallel HDF5 file using the MPI-IO driver.

Equivalent to `h5open(filename, mode; fapl_mpio=(comm, info), pv...)`.
Throws an informative error if the loaded HDF5 libraries do not include parallel
support.

See the [HDF5 docs](https://portal.hdfgroup.org/display/HDF5/H5P_SET_FAPL_MPIO)
for details on the `comm` and `info` arguments.
"""
function h5open(
        filename::AbstractString, mode::AbstractString,
        comm::MPI.Comm, info::MPI.Info = MPI.Info(); pv...
    )
    check_hdf5_parallel()
    h5open(filename, mode; fapl_mpio=(comm, info), pv...)
end

h5open(filename::AbstractString, comm::MPI.Comm, args...; pv...) =
    h5open(filename, "r", comm, args...; pv...)


text:
"""
LDA correlation according to Vosko Wilk,and Nusair, (DOI 10.1139/p80-159)
"""
function lda_c_vwn!(ρ; E=nothing, Vρ=nothing)
    x0  = sqrt(5913592351)
    x1  = 2^(2/3)*3^(1/6)
    x2  = pi^(1/6)
    x3  = ρ.^(1/6)
    x4  = x2.*x3
    x5  = 468052052839577158149*x0/928058053502134909615000000
    x6  = 6^(1/3)
    x7  = 12500*x4.*x6
    x8  = pi^(1/3)
    x9  = ρ.^(1/3)
    x10 = x8.*x9
    x11 = 1 ./(46593*x1.*x10 + x7 + 323380*sqrt(pi)*sqrt.(ρ))
    x12 = x1./x2
    x13 = x12./x3
    x14 = x13 .+ 46593/12500
    x15 = x6./x8
    x16 = x15./(2*x9)
    x17 = 46593*x13/25000 .+ x16 .+ 16169/1250
    x18 = 1 ./x17
    x19 = x13/2 .+ 5249/50000
    x20 = x19.^2
    x21 = x14.^(-2)
    x22 = x12./ρ.^(7/6)
    x23 = x15./(6*ρ.^(4/3))
    x24 = (15531*x22/50000 + x23)./x17.^2

    if E !== nothing
        @. E .= x5*atan(x0*x4/(12500*x1 + 46593*x4)) + 310907*log(x11*x7)/10000000 + 76037485627899*log(x11*x4*(25000*x1 + 5249*x4)^2/100000)/78468213432500000
    end
    if Vρ !== nothing
        @. Vρ .= x5*atan(x0/(12500*x14)) + ρ*(310907*6^(2/3)*x10*x17*(x16*x24 - x18*x23)/30000000 + 76037485627899*x17*(-x18*x19*x22/6 + x20*x24)/(78468213432500000*x20) + 156017350946525719383*x21*x22/(3923410671625000000000*(5913592351*x21/156250000 + 1))) + 310907*log(x16*x18)/10000000 + 76037485627899*log(x18*x20)/78468213432500000
    end

    (E=E, Vρ=Vρ)
end


text:
# Implementation follows Rasmussen & Williams, Gaussian Processes for Machine
# Learning, the MIT Press, 2006. In the following referred to as 'RW'.
# Online text:
# - http://www.gaussianprocess.org/gpml/chapters/RW3.pdf
# - http://www.gaussianprocess.org/gpml/chapters/RW5.pdf

struct LaplaceApproximation{Tkw}
    newton_kwargs::Tkw
end

LaplaceApproximation(; newton_kwargs...) = LaplaceApproximation((; newton_kwargs...))

"""
    posterior(la::LaplaceApproximation, lfx::LatentFiniteGP, ys)

Construct a Gaussian approximation `q(f)` to the posterior `p(f | y)` using the
Laplace approximation. Solves for a mode of the posterior using Newton's
method.
"""
function AbstractGPs.posterior(la::LaplaceApproximation, lfx::LatentFiniteGP, ys)
    dist_y_given_f, K, newton_kwargs = _check_laplace_inputs(lfx, ys; la.newton_kwargs...)
    f_opt = newton_inner_loop(dist_y_given_f, ys, K; newton_kwargs...)
    # We could call _newton_inner_loop directly and immediately use the
    # returned cache, but this would make the posterior call
    # non-differentiable, at only a marginal computational saving.
    cache = _laplace_train_intermediates(dist_y_given_f, ys, K, f_opt)
    f_post = ApproxPosteriorGP(la, lfx.fx, cache)
    return f_post  # TODO return LatentGP(f_post, lfx.lik, lfx.fx.Σy)
end

"""
    approx_lml(la::LaplaceApproximation, lfx::LatentFiniteGP, ys)

Compute an approximation to the log of the marginal likelihood (also known as
"evidence"), which can be used to optimise the hyperparameters of `lfx`.

This should become part of the AbstractGPs API (see JuliaGaussianProcesses/AbstractGPs.jl#221).
"""
function approx_lml(la::LaplaceApproximation, lfx::LatentFiniteGP, ys)
    return laplace_lml(lfx, ys; la.newton_kwargs...)
end

"""
    build_laplace_objective(build_latent_gp, xs, ys; kwargs...)

Construct a closure that computes the minimisation objective for optimising
hyperparameters of the latent GP in the Laplace approximation. The returned
closure passes its arguments to `build_latent_gp`, which must return the
`LatentGP` prior.

# Keyword arguments

- `newton_warmstart=true`: (default) begin Newton optimisation at the mode of
  the previous call to the objective
- `newton_callback`: called as `newton_callback(fnew, cache)` after each Newton step
- `newton_maxiter=100`: maximum number of Newton steps.
"""
function build_laplace_objective(build_latent_gp, xs, ys; kwargs...)
    # TODO assumes type of `xs` will be same as `mean(lfx.fx)`
    f = similar(xs, length(xs))  # will be mutated in-place to "warm-start" the Newton steps
    return build_laplace_objective!(f, build_latent_gp, xs, ys; kwargs...)
end

function build_laplace_objective!(
    f,
    build_latent_gp,
    xs,
    ys;
    newton_warmstart=true,
    newton_callback=nothing,
    newton_maxiter=100,
)
    initialize_f = true

    function objective(args...)
        lf = build_latent_gp(args...)
        lfx = lf(xs)
        ignore_derivatives() do
            # Zygote does not like the try/catch within @info etc.
            @debug "Objective arguments: $args"
            # Zygote does not like in-place assignments either
            if initialize_f
                f .= mean(lfx.fx)
            end
        end
        f_opt, lml = laplace_f_and_lml(
            lfx, ys; f_init=f, maxiter=newton_maxiter, callback=newton_callback
        )
        ignore_derivatives() do
            if newton_warmstart
                f .= f_opt
                initialize_f = false
            end
        end
        return -lml
    end

    return objective
end

"""
    laplace_f_and_lml(lfx::LatentFiniteGP, ys; newton_kwargs...)

Compute a mode of the posterior and the Laplace approximation to the log
marginal likelihood.
"""
function laplace_f_and_lml(lfx::LatentFiniteGP, ys; newton_kwargs...)
    dist_y_given_f, K, newton_kwargs = _check_laplace_inputs(lfx, ys; newton_kwargs...)
    f_opt = newton_inner_loop(dist_y_given_f, ys, K; newton_kwargs...)
    lml = laplace_lml(dist_y_given_f, ys, K, f_opt)
    return f_opt, lml
end

"""
    laplace_lml(lfx::LatentFiniteGP, ys; newton_kwargs...)

Compute the Laplace approximation to the log marginal likelihood.
"""
function laplace_lml(lfx::LatentFiniteGP, ys; newton_kwargs...)
    dist_y_given_f, K, newton_kwargs = _check_laplace_inputs(lfx, ys; newton_kwargs...)
    return laplace_lml(dist_y_given_f, ys, K; newton_kwargs...)
end

function laplace_lml(dist_y_given_f, ys, K; f_init, maxiter, newton_kwargs...)
    f_opt = newton_inner_loop(dist_y_given_f, ys, K; f_init, maxiter, newton_kwargs...)
    return laplace_lml(dist_y_given_f, ys, K, f_opt)
end

function laplace_lml(dist_y_given_f, ys, K, f_opt)
    cache = _laplace_train_intermediates(dist_y_given_f, ys, K, f_opt)
    return _laplace_lml(f_opt, cache)
end

function _check_laplace_inputs(
    lfx::LatentFiniteGP, ys; f_init=nothing, maxiter=100, newton_kwargs...
)
    fx = lfx.fx
    @assert mean(fx) == zero(mean(fx))  # might work with non-zero prior mean but not checked
    @assert length(ys) == length(fx)  # LaplaceApproximation currently does not support multi-latent likelihoods
    dist_y_given_f = lfx.lik
    K = cov(fx)
    if isnothing(f_init)
        f_init = mean(fx)
    end
    return dist_y_given_f, K, (; f_init, maxiter, newton_kwargs...)
end

struct LaplaceCache{
    Tm<:AbstractMatrix,Tv<:AbstractVector,Td<:Diagonal,Tf<:Real,Tc<:Cholesky
}
    K::Tm  # kernel matrix
    f::Tv  # mode of posterior p(f | y)
    W::Td  # diagonal matrix of ∂²/∂fᵢ² loglik
    Wsqrt::Td  # sqrt(W)
    loglik::Tf  # ∑ᵢlog p(yᵢ|fᵢ)
    d_loglik::Tv  # ∂/∂fᵢloglik
    B_ch::Tc  # cholesky(I + Wsqrt * K * Wsqrt)
    a::Tv  # K⁻¹ f
end

function _laplace_train_intermediates(dist_y_given_f, ys, K, f)
    # Ψ = log p(y|f) + log p(f)
    #   = loglik + log p(f)
    # dΨ/df = d_loglik - K⁻¹ f
    # at fhat: d_loglik = K⁻¹ f

    # d²Ψ/df² = d2_loglik - K⁻¹
    #         = -W - K⁻¹

    ll, d_ll, d2_ll = loglik_and_derivs(dist_y_given_f, ys, f)

    # inner loop iteration of RW Algorithm 3.1, lines 4-7
    W = -Diagonal(d2_ll)
    Wsqrt = sqrt(W)
    B = I + Wsqrt * K * Wsqrt
    B_ch = cholesky(Symmetric(B))
    b = W * f + d_ll
    a = b - Wsqrt * (B_ch \ (Wsqrt * K * b))

    #return (; K, f, W, Wsqrt, loglik=ll, d_loglik=d_ll, B_ch, a)
    return LaplaceCache(K, f, W, Wsqrt, ll, d_ll, B_ch, a)
end

"""
    loglik_and_derivs(dist_y_given_f, ys, f)

`dist_y_given_f` must be a scalar function from a Real to a Distribution object.
`ys` and `f` are vectors of observations and latent function values, respectively.
"""
function loglik_and_derivs(dist_y_given_f, ys::AbstractVector, f::AbstractVector{<:Real})
    l(_f, _y) = logpdf(dist_y_given_f(_f), _y)
    dl(_f, _y) = ForwardDiff.derivative(__f -> l(__f, _y), _f)
    d2l(_f, _y) = ForwardDiff.derivative(__f -> dl(__f, _y), _f)

    ls = map(l, f, ys)
    d_ll = map(dl, f, ys)
    d2_ll = map(d2l, f, ys)

    ll = sum(ls)
    return ll, d_ll, d2_ll
end

function _newton_step(dist_y_given_f, ys, K, f)
    cache = _laplace_train_intermediates(dist_y_given_f, ys, K, f)
    # inner loop iteration of RW Algorithm 3.1, line 8
    fnew = K * cache.a
    return fnew, cache
end

function _laplace_lml(f, cache)
    # inner loop iteration of RW Algorithm 3.1, line 10
    # -a' * f / 2 + loglik - sum(log.(diag(B_ch.L)))
    return -cache.a' * f / 2 + cache.loglik - sum(log.(diag(cache.B_ch.L)))
end

function _newton_inner_loop(dist_y_given_f, ys, K; f_init, maxiter, callback=nothing)
    @assert maxiter >= 1
    f = f_init
    cache = nothing
    for i in 1:maxiter
        @debug "  - Newton iteration $i: f[1:3]=$(f[1:3])"
        fnew, cache = _newton_step(dist_y_given_f, ys, K, f)
        if !isnothing(callback)
            callback(fnew, cache)
        end

        if isapprox(f, fnew)
            @debug "  + converged"
            #f = fnew
            break  # converged
        else
            f = fnew
        end
    end
    return f, something(cache)
end

function ChainRulesCore.frule(Δargs, ::typeof(_newton_inner_loop), args...; kwargs...)
    return _newton_inner_loop_no_derivatives()
end
function ChainRulesCore.rrule(::typeof(_newton_inner_loop), args...; kwargs...)
    return _newton_inner_loop_no_derivatives()
end
function _newton_inner_loop_no_derivatives()
    # It's important that we don't simply call _newton_inner_loop and pass the
    # resulting cache directly to _laplace_lml. This would result in the wrong
    # gradients. Instead, in newton_inner_loop we return only f_opt, which will
    # have the correct gradients thanks to the custom frule/rrule for
    # newton_inner_loop. The f_opt form of laplace_lml will then explicitly
    # call _laplace_train_intermediates, which allows the AD system to
    # correctly compute the gradients.
    return error(
        "Do not try to compute the derivatives of _newton_inner_loop directly. " *
        "Instead, call newton_inner_loop, which has the correct frule/rrule.",
    )
end

# Currently, we have a separate function that returns only f_opt to simplify frule/rrule
"""
    newton_inner_loop(dist_y_given_f, ys, K; f_init, maxiter, callback=nothing)

Find a mode of `p(f | y)` using Newton's method.
"""
function newton_inner_loop(dist_y_given_f, ys, K; f_init, maxiter, callback=nothing)
    f_opt, _ = _newton_inner_loop(dist_y_given_f, ys, K; f_init, maxiter, callback)
    return f_opt
end

function ChainRulesCore.frule(
    (Δself, Δdist_y_given_f, Δys, ΔK),
    ::typeof(newton_inner_loop),
    dist_y_given_f,
    ys,
    K;
    kwargs...,
)
    f_opt, cache = _newton_inner_loop(dist_y_given_f, ys, K; kwargs...)

    # f = K grad_log_p_y_given_f(f)
    # fdot = Kdot grad_log_p_y_given_f(f) + K grad2_log_p_y_given_f(f) fdot
    # fdot (I - K grad2_log_p_y_given_f(f)) = Kdot grad_log_p_y_given_f(f)
    # fdot = (I - K grad2_log_p_y_given_f(f))⁻¹ Kdot grad_log_p_y_given_f(f)
    # (I - K grad2_log_p_y_given_f(f)) = (I + K W) = (√W)⁻¹ (I + √W K √W) √W = (√W)⁻¹ B √W
    # fdot = (√W)⁻¹ B⁻¹ √W Kdot grad_log_p_y_given_f(f)
    ∂f_opt = cache.Wsqrt \ (cache.B_ch \ (cache.Wsqrt * (ΔK * cache.d_loglik)))

    return f_opt, ∂f_opt
end

function ChainRulesCore.rrule(::typeof(newton_inner_loop), dist_y_given_f, ys, K; kwargs...)
    f_opt, cache = _newton_inner_loop(dist_y_given_f, ys, K; kwargs...)

    # f = K (∇log p(y|f))                               (RW 3.17)
    # δf = δK (∇log p(y|f)) + K δ(∇log p(y|f))
    #    = δK (∇log p(y|f)) + K ∇(∇log p(y|f)) δf
    # δf (I - K ∇∇log p(y|f)) = δK (∇log p(y|f))
    # δf (I + K W) = δK (∇log p(y|f))
    # δf = (I + K W)⁻¹ δK (∇log p(y|f))                 (RW 5.24)
    # (I + K W) = (√W)⁻¹ (I + √W K √W) √W = (√W)⁻¹ B √W
    # δf = (√W)⁻¹ B⁻¹ √W δK (∇log p(y|f))

    # ∂f_opt = cache.Wsqrt \ (cache.B_ch \ (cache.Wsqrt * (ΔK * cache.d_loglik)))

    # Re⟨Δf, δf⟩ = Re⟨Δf, Wsqrt⁻¹ B⁻¹ Wsqrt δK d_loglik⟩
    #            = Re⟨Wsqrt' B⁻¹' Wsqrt⁻¹' Δf d_loglik', δK⟩
    #
    # ΔK = Wsqrt' * cache.B_ch' \ Wsqrt' \ Δf_opt * cache.d_loglik'

    function newton_pullback(Δf_opt)
        ∂self = NoTangent()

        ∂dist_y_given_f = @not_implemented(
            "gradient of Newton's method w.r.t. likelihood parameters"
        )

        ∂ys = @not_implemented("gradient of Newton's method w.r.t. observations")

        # ∂K = df/dK Δf
        ∂K = @thunk(cache.Wsqrt * (cache.B_ch \ (cache.Wsqrt \ Δf_opt)) * cache.d_loglik')

        return (∂self, ∂dist_y_given_f, ∂ys, ∂K)
    end

    return f_opt, newton_pullback
end

"""
    laplace_f_cov(cache)

Compute the covariance of `q(f)` from the results of the training computation
that are stored in a `LaplaceCache`.
"""
function laplace_f_cov(cache)
    # (K⁻¹ + W)⁻¹
    # = (√W⁻¹) (√W⁻¹ (K⁻¹ + W) √W⁻¹)⁻¹ (√W⁻¹)
    # = (√W⁻¹) (√W⁻¹ K⁻¹ √W⁻¹ + I)⁻¹ (√W⁻¹)
    # ; (I + C⁻¹)⁻¹ = I - (I + C)⁻¹
    # = (√W⁻¹) (I - (I + √W K √W)⁻¹) (√W⁻¹)
    # = (√W⁻¹) (I - B⁻¹) (√W⁻¹)
    B_ch = cache.B_ch
    Wsqrt_inv = inv(cache.Wsqrt)
    return Wsqrt_inv * (I - inv(B_ch)) * Wsqrt_inv
end

function LaplaceResult(fnew, cache)
    f = cache.f
    f_cov = laplace_f_cov(cache)
    q = MvNormal(f, AbstractGPs._symmetric(f_cov))
    lml_approx = _laplace_lml(f, cache)

    return (; fnew, f_cov, q, lml_approx, cache...)
end

"""
    laplace_steps(lfx::LatentFiniteGP, ys; newton_kwargs...)

For demonstration purposes: returns an array of all the intermediate
approximations of each Newton step.

If you are only interested in the actual posterior, use
`posterior(::LaplaceApproximation, ...`.

TODO figure out how to get the `@ref` to work to point to the LaplaceApproximation-specific `posterior` docstring...
"""
function laplace_steps(lfx::LatentFiniteGP, ys; newton_kwargs...)
    dist_y_given_f, K, newton_kwargs = _check_laplace_inputs(lfx, ys; newton_kwargs...)

    res_array = []

    function store_result!(fnew, cache)
        return push!(res_array, LaplaceResult(fnew, cache))
    end

    _ = newton_inner_loop(dist_y_given_f, ys, K; newton_kwargs..., callback=store_result!)

    return res_array
end

const LaplacePosteriorGP = ApproxPosteriorGP{<:LaplaceApproximation}

function _laplace_predict_intermediates(cache, prior_at_x, xnew)
    k_x_xnew = cov(prior_at_x.f, prior_at_x.x, xnew)
    f_mean = mean(prior_at_x.f, xnew) + k_x_xnew' * cache.d_loglik  # RW (3.21)
    L = cache.B_ch.L
    v = L \ (cache.Wsqrt * k_x_xnew)  # RW (3.29)
    return f_mean, v
end

function StatsBase.mean_and_var(f::LaplacePosteriorGP, x::AbstractVector)
    f_mean, v = _laplace_predict_intermediates(f.data, f.prior, x)
    f_var = var(f.prior.f, x) - vec(sum(v .^ 2; dims=1))  # RW (3.29) |> diag
    return f_mean, f_var
end

function StatsBase.mean_and_cov(f::LaplacePosteriorGP, x::AbstractVector)
    f_mean, v = _laplace_predict_intermediates(f.data, f.prior, x)
    f_cov = cov(f.prior.f, x) - v' * v  # RW (3.29)
    return f_mean, f_cov
end

function Statistics.mean(f::LaplacePosteriorGP, x::AbstractVector)
    d_loglik = f.data.d_loglik
    return mean(f.prior.f, x) + cov(f.prior.f, f.prior.x, x)' * d_loglik
end

function Statistics.var(f::LaplacePosteriorGP, x::AbstractVector)
    return last(mean_and_var(f, x))
end

function Statistics.cov(f::LaplacePosteriorGP, x::AbstractVector)
    return last(mean_and_cov(f, x))
end

function Statistics.cov(f::LaplacePosteriorGP, x::AbstractVector, y::AbstractVector)
    L = f.data.B_ch.L
    vx = L \ (f.data.Wsqrt * cov(f.prior.f, f.prior.x, x))
    vy = L \ (f.data.Wsqrt * cov(f.prior.f, f.prior.x, y))
    return cov(f.prior.f, x, y) - vx' * vy
end


text:
export climscatter, climscatter!, climsurface!, climsurface, climplot

"""
    climplot(A::ClimArray; kwargs...) → fig, ax, el, cb
Main plotting function that dispatches to [`climscatter!`](@ref) if `A` has
an [`CoordinateSpace`](@ref) dimension, or to [`climsurface!`](@ref) for [`OrthogonalSpace`](@ref).

Return the figure, axis, plotted element, and colorbar.

Optionally you can provide keyword `scatter = true` to force using the scatterplot.

Plotting from ClimateBase.jl also works with `Observable`s that enclose a `ClimArray`.
You can update the values of the observable with another `ClimArray` with the same spatial
dimension, and the plot will be updated. See documentation online for examples.
"""
function climplot(A, args...; scatter = spacestructure(A) == CoordinateSpace(),
    source = "+proj=longlat +datum=WGS84", dest = "+proj=eqearth",
    colorbar = true, name = string(DimensionalData.name(A)), kwargs...)

    # TODO: Perhaps setting custom colorrange leads to better plots?
    # vmin = haskey(kwargs, :vmin) ? kwargs[:vmin] : quantile(data, 0.025)
    # vmax = haskey(kwargs, :vmax) ? kwargs[:vmax] : quantile(data, 0.975)

    fig = GeoMakie.Figure()
    ax = GeoMakie.GeoAxis(fig[1,1]; source, dest, title = GeoMakie.Observable(""))

    coastplot = lines!(ax, GeoMakie.coastlines(); color = :black, overdraw = true, coastkwargs...)
    translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements

    # TODO: @assert A has only space dimension

    if scatter
        el = climscatter!(ax, A; kwargs...)
    else
        el = climsurface!(ax, A; kwargs...)
    end

    if colorbar
        cb = GeoMakie.Colorbar(fig[1, 2], el; label = name)
    else
        cb = nothing
    end
    return fig, ax, el, cb
end

##########################################################################################
# # Scatter
##########################################################################################
# Notice that `A` is not declared as `ClimArray`, but assumed to be.
# Duck-typing for Observables.
function climscatter!(ax, A; colormap = :dense, kwargs...)
    if hasdim(A, Coord)
        lonlat = dims(A, Coord).val
    elseif dims(A)[1] isa Lon
        londim = dims(A, Lon)
        latdim = dims(A, Lat)
        lonlat = [GeoMakie.Point2f0(l,lat) for lat in latdim for l in londim]
    else
        error("Unknown spatial dimensions for input.")
    end
    data = GeoMakie.lift(A -> vec(A.data), A)
    GeoMakie.scatter!(ax, lonlat; color = data, colormap, kwargs...)
end

##########################################################################################
# # Surface
##########################################################################################
# TODO: Change this to `contourf`

# Notice that `A` is not declared as `ClimArray`, but assumed to be.
# Duck-typing for Observables.
function climsurface!(ax, A; colormap = :dense, kwargs...)
    # TODO: @assert A has only space
    if hasdim(A, Coord)
        # TODO: @pkeil this is for you
        error("Surface plots for `Coord` arrays are not supported yet!")
    end
    # TODO: This needs to depend on lon_0 or use the "meridian cut" function
    B = GeoMakie.lift(A -> longitude_circshift(A), A)
    lon = dims(B, Lon).val
    lat = dims(B, Lat).val
    data = GeoMakie.lift(B -> B.data, B)
    GeoMakie.surface!(ax, lon, lat, data; shading = false, colormap, kwargs...)
end

##########################################################################################
# # Observables overloads
##########################################################################################
# These overloads allows me to write generic code that does not care if input
# is observable or not. This leads to simple, clean, small, elegant code for
# animating fields by simply replacing them with their observables.

GeoMakie.lift(f, x::ClimArray) = f(x)
DimensionalData.dims(o::GeoMakie.Observable, args...) = dims(o.val, args...)
DimensionalData.name(o::GeoMakie.Observable, args...) = name(o.val, args...)
spacestructure(A::GeoMakie.Observable) = spacestructure(A.val)


text:
using Documenter, NeuralPDE

makedocs(
    sitename="NeuralPDE.jl",
    authors="#",
    clean=true,
    doctest=false,
    modules=[NeuralPDE],

    format=Documenter.HTML(# analytics = "",
                             assets=["assets/favicon.ico"],
                             canonical="#"),
    pages=[
        "NeuralPDE.jl: Scientific Machine Learning (SciML) for Partial Differential Equations" => "index.md",
        "Tutorials" => Any[
            "examples/ode.md",
            "examples/pinns_example.md",
            "examples/100_HJB.md",
            "examples/blackscholes.md",
            "examples/kolmogorovbackwards.md",
            "examples/optimal_stopping_american.md",
            "examples/nnrode_example.md",
        ],
        "Neural-Enhanced PDE Solvers" => Any[
            "solvers/ode.md",
            "solvers/pinn.md",
            "solvers/deep_fbsde.md",
            "solvers/kolmogorovbackwards_solver.md",
            "solvers/optimal_stopping.md"
        ]
    ]
)

deploydocs(
   repo="github.com/SciML/NeuralPDE.jl.git";
   push_preview=true
)


text:
############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenStrict,
    x::AbstractVector{R}
) where {T<:AbstractFloat,F<:FlattenTypes,R<:Real}
    len = length(x)
    buffer_flat = zeros(T, size(x))
    function flatten_to_Real_Vec(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(buffer_flat, v)
    end
    buffer_unflat = zeros(R, size(x))
    function unflatten_to_Real_Vec(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(buffer_unflat, v)
    end
    return flatten_to_Real_Vec, unflatten_to_Real_Vec
end

function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenFlexible,
    x::AbstractVector{R}
) where {T<:AbstractFloat,F<:FlattenTypes,R<:Real}
    len = length(x)
    function flatten_to_Real_Vec_AD(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return identity(v)
    end
    function unflatten_to_Real_Vec_AD(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return identity(v)
    end
    return flatten_to_Real_Vec_AD, unflatten_to_Real_Vec_AD
end

############################################################################################
#Export
export
    construct_flatten


text:
output_clause(v) = join(v, " ") * " 0"

"""
Create string for DIMACS output.
Each clause is terminated by 0.

TODO: If clauses are too long they should be split up into separate lines.
"""
function dimacs_output(p::SATProblem)

    header = """
    c DIMACS created by sat.jl
    p cnf $(p.num_variables) $(length(p.clauses))
    """

    body = join(output_clause.(p.clauses), "\n")

    return header * body * "\n"
end

parse_line(s) = parse.(Int, split(s)[2:end])

# exit codes from CryptoMiniSat:
const SAT_status = Dict(10 => :sat,
                         15 => :unknown,
                         20 => :unsat)

function solve(p::SATProblem; solver=cryptominisat)
    s = dimacs_output(p)

    input = "_input.cnf"
    write(input, s)

    return @time call_solver(solver, input=input)
end

abstract type AbstractSATSolver end

"""
Object representing an external SAT solver binary (executable).

- `name`: Name of the executable
- `options`: Options to pass in
- `status`: Dictionary containing exit codes and their interpretation
"""
struct ExternalSATSolver <: AbstractSATSolver
    name::String
    options::String
    status::Dict{Int, Symbol}
end

const cryptominisat =
    ExternalSATSolver("cryptominisat5", "",
                Dict(10 => :sat,
                     15 => :unknown,
                     20 => :unsat) )

"""
Call a solver executable (binary) that takes `.cnf` input.
The solver's standard output is captured to a file.
"""
function call_solver(solver; input="_input.cnf", output="_out.txt")

    # run the solver:
    cmd = `$(solver.name) $input`

    pipe = pipeline(cmd, stdout=output)
    out = run(pipe, wait=false);
    wait(out)

    status = solver.status[out.exitcode]

    status ∈ (:unsat, :unknown) && return status, Int[]

    output = filter(line -> !startswith(line, "c"), readlines(output))

    # TODO: check for correct initial letter

    results = reduce(vcat, parse_line.(output[2:end]))
    pop!(results)  # remove final 0

    return status, results
end

"""
Check that the clause is satisfied by the assignments
in `results`
"""
function satisfies(clause::Vector{Int}, results)
    return any(x -> results[abs(x)] == x, clause)
end

"""
Check that all clases in `p` are satisfied by the
assignments in `results`
"""
function satisfies(p, results)
    return all(satisfies.(p.clauses, Ref(results)))
end