text:
function test_diagnosticsteps(dev::Device=CPU(); nsteps=100, freq=1, ndata=ceil(Int, (nsteps+1)/freq))
  prob = Problem(nx=6, Lx=2π, dev=dev)
  getone(prob) = 1
  diagnostic = Diagnostic(getone, prob, nsteps=nsteps, freq=freq, ndata=ndata)
  stepforward!(prob, diagnostic, nsteps)
  expectedsteps = cat([0], freq:freq:nsteps, dims=1)

  return diagnostic.steps[1:diagnostic.i] == expectedsteps
end

function test_scalardiagnostics(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2, nsteps=100, freq=1, ndata=ceil(Int, (nsteps+1)/freq))
  k1 = 2π/Lx
  dt = 1e-2 / (kappa*k1^2) # time-scale for diffusive decay dynamics are resolved
  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="RK4", dev=dev)
  c0 = @. sin(k1 * prob.grid.x)

  ct(t) = exp(-kappa*k1^2*t)

  t1 = dt*nsteps
  c0 = @. sin(k1*prob.grid.x)
  c1 = @. ct(t1) * c0 # analytical solution

  set_c!(prob, c0)

  variance(prob) = parsevalsum2(prob.sol, prob.grid) / Lx
  diagnostic = Diagnostic(variance, prob, nsteps=nsteps, freq=freq, ndata=ndata)

  stepforward!(prob, diagnostic, nsteps)

  da = @. 0.5*ct(diagnostic[:t])^2
  return isapprox(diagnostic[:data], da)
end

function test_basicdiagnostics(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2)
  k1 = 2π/Lx
  dt = 1e-2 / (kappa*k1^2) # time-scale for diffusive decay dynamics are resolved
  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="RK4", dev=dev)
  c0 = @. cos(k1 * prob.grid.x)
  set_c!(prob, c0)

  getsol(prob) = prob.sol
  soldiag = Diagnostic(getsol, prob, nsteps=1)

  stepforward!(prob, soldiag, 1)

  return isapprox(soldiag.data[soldiag.i], prob.sol)
end

function test_extenddiagnostic(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2)
  k1 = 2π/Lx
  dt = 1e-2 / (kappa*k1^2) # time-scale for diffusive decay dynamics are resolved
  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="RK4", dev=dev)
  c0 = @. sin(k1 * prob.grid.x)
  set_c!(prob, c0)

  getsol(prob) = prob.sol

  soldiagnostic1 = Diagnostic(getsol, prob, nsteps=1)
  soldiagnostic2 = Diagnostic(getsol, prob, nsteps=4)

  nsteps_initially1 = length(soldiagnostic1.t)
  nsteps_extend1 = 12
  FourierFlows.extend!(soldiagnostic1, nsteps_extend1)

  nsteps_initially2 = length(soldiagnostic2.t)
  FourierFlows.extend!(soldiagnostic2)

  return length(soldiagnostic1.t) == nsteps_initially1 + nsteps_extend1 && length(soldiagnostic2.t) == 2*nsteps_initially2
end

function test_incrementdiagnostic(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2)
  k1 = 2π/Lx
  dt = 1e-2 / (kappa*k1^2) # time-scale for diffusive decay dynamics are resolved
  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="ETDRK4", dev=dev)
  c0 = @. sin(k1 * prob.grid.x)
  set_c!(prob, c0)

  dummydiagnostic1(prob) = 2.0
  dummydiagnostic2(prob) = 2.0im

  diagnostic1 = Diagnostic(dummydiagnostic1, prob, nsteps=10)
  diagnostic2 = Diagnostic(dummydiagnostic2, prob, nsteps=10)

  diagnostics = [diagnostic1, diagnostic2]

  FourierFlows.increment!(diagnostics)

  return diagnostic1.data[2]==2.0 && diagnostic2.data[2]==2.0im
end

function test_getindex(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2)
  k1 = 2π/Lx
  dt = 1e-2 / (kappa*k1^2) # time-scale for diffusive decay dynamics are resolved
  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="ETDRK4", dev=dev)
  c0 = @. sin(k1 * prob.grid.x)
  set_c!(prob, c0)

  dummydiagnostic1(prob) = 2.0
  dummydiagnostic2(prob) = 2.0im

  nsteps = 10

  diagnostic1 = Diagnostic(dummydiagnostic1, prob, nsteps=nsteps)
  diagnostic2 = Diagnostic(dummydiagnostic2, prob, nsteps=nsteps)

  diagnostics = [diagnostic1, diagnostic2]

  ntimes=7
  for j=1:ntimes
    FourierFlows.increment!(diagnostics)
  end

  return getindex(diagnostic1, 2:5) == 2.0*ones(4) && getindex(diagnostic2, 2:5) == 2.0im*ones(4)
end


text:
using Paillier

rngprime(rng::AbstractRNG,N) = Paillier.nbit_prime_of_size(rng,N)
rngint(rng::AbstractRNG,N) = Paillier.n_bit_random_number(rng,N)

function hexstrtoint(str)
    hexstr = lowercase(join(split(replace(str,"\n"=>"")," ")))
    return parse(BigInt,str,base=16)
end

macro hex_str(str)
    esc(hexstrtoint(str))
end

ext: jl


text:
export conv
import NNlib
using AutoGrad

"""
    conv(w, x; kwargs...)

Return the convolution of filter `w` with tensor `x`, optionally with bias/residual
addition, activation and/or scaling:

    y = activation(alpha * conv(w,x) + beta * z .+ bias)

All tensors should have the same number of dimensions. 3-5 dimensional tensors are
supported. The sizes should be:

    tensor  default       channelmajor
    ------  -------       ------------
    x       (X...,Cx,N)   (Cx,X...,N)
    y       (Y...,Cy,N)   (Cy,Y...,N)
    w       (W...,Cw,Cy)  (Cw,W...,Cy) # Cw=Cx÷group
    bias    (1...,Cy,1)   (Cy,1...)
    z       (Y...,Cy,N)   (Cy,Y...,N)  # same as y

where `(X...)`, `(Y...)`, `(W...)` are spatial dimensions, `Cx`,`Cy`,`Cw` are the number of
input/output/filter channels, and `N` is the number of instances. Both `Cx` and `Cy` have to
be an exact multiple of `group` and `Cw` must be `Cx÷group`. The `channelmajor` option
corresponds to `CUDNN_TENSOR_NHWC` format in cuDNN.

The arguments `padding`, `stride` and `dilation` can be specified as `n-2` dimensional
integer vectors, tuples or a single integer which is assumed to be repeated `n-2` times
where `n` is the number of `x` dimensions. If any of the entries is larger than the
corresponding `x` dimension, the `x` dimension is used instead. For a description of
different types of convolution see:
https://towardsdatascience.com/a-comprehensive-introduction-to-different-types-of-convolutions-in-deep-learning-669281e58215

Keyword arguments:
* `activation = nothing`: apply activation function if provided
* `normalization = nothing`: apply normalization function if provided
* `alpha = 1, beta = 0`: scaling parameters
* `bias = nothing`: add bias if provided
* `channelmajor = false`: assume channel-major format tensors if specified
* `crosscorrelation = true`: if true apply cross-correlation rather than convolution (flipped kernel)
* `dilation = 1`: dilation factor
* `groups = 1`: number of groups to be used
* `padding = 0`: padding assumed around `x`
* `stride = 1`: how far to shift the convolution window at each step
* `z = nothing`: add `beta*z` to the result if specified
"""
function conv(
    w, x;
    activation = nothing,
    normalization = nothing,
    alpha = 1,
    beta = 0,
    bias = nothing,
    channelmajor = false,
    crosscorrelation = true,
    dilation = 1,
    groups = 1,
    padding = 0,
    stride = 1,
    z = nothing,
    o...
)
    if channelmajor
        @warn "channelmajor is not implemented on the CPU yet, using slow permutedims, see NNlib#267" maxlog=1
        w,x,z,bias = (x->(x===nothing ? nothing : permutedims(x,(2,3,1,4)))).((w,x,z,bias))
    end
    N = ndims(w)
    stride = NNlib.expand(Val(N-2), stride)
    padding = NNlib.expand(Val(N-2), padding)
    dilation = NNlib.expand(Val(N-2), dilation)
    if groups == 1
        cdims = NNlib.DenseConvDims(size(x), size(w); stride, padding, dilation, flipkernel=crosscorrelation)
        y = NNlib.conv(x, w, cdims)
    else
        @warn "groups > 1 is not implemented on the CPU yet, using slow manual implementation, see NNlib#267" maxlog=1
        @assert size(x,3) == size(w,3)*groups
        @assert size(w,4) % groups == 0
        ys,dx,dy = Array{Any}(undef,groups), size(w,3), size(w,4)÷groups
        for i = 1:groups
            xi = view(x, :, :, (1+(i-1)*dx):(i*dx), :)
            wi = view(w, :, :, :, (1+(i-1)*dy):(i*dy))
            cdims = NNlib.DenseConvDims(size(xi), size(wi); stride, padding, dilation, flipkernel=crosscorrelation)
            ys[i] = NNlib.conv(xi, wi, cdims)
        end
        y = cat(ys...; dims=3)
    end
    if alpha != 1; y = alpha * y; end
    if beta != 0 && z !== nothing; y = y + beta * z; end
    if bias !== nothing; y = y .+ bias; end
    if normalization ∉ (nothing, identity); y = normalization(y); end
    if activation ∉ (nothing, identity); y = activation(y); end
    if channelmajor; y = permutedims(y, (3,1,2,4)); end
    return y
end

@primitive1 NNlib.conv(x, w, cdims; o...),dy,y  NNlib.∇conv_data(dy, w, cdims; o...)  NNlib.∇conv_filter(x, dy, cdims; o...)
@primitive1 NNlib.∇conv_data(dy, w, cdims; o...),ddx,dx  NNlib.conv(ddx, w, cdims; o...)  NNlib.∇conv_filter(ddx, dy, cdims; o...)
@primitive1 NNlib.∇conv_filter(x, dy, cdims; o...),ddw,dw  NNlib.∇conv_data(dy, ddw, cdims; o...)  NNlib.conv(x, ddw, cdims; o...)


text:
struct Torrent
    file::String
    content::Dict{String,Any}
end

function Torrent(file::AbstractString)
    content = open(bdecode, file)
    return Torrent(file, content)
end

function create_torrent(path; block_size=256*1024, meta::Dict=Dict{String,Any}())
    path = abspath(path)
    # Compute file information
    info = Dict{String,Any}()
    if isfile(path)
        files = [path]
        ## Populate single-file torrent keys
        info["name"] = basename(path)
        info["length"] = filesize(path)
    elseif isdir(path)
        ## Populate multi-file torrent keys
        info["name"] = splitpath(path)[end]
        files = String[]
        for (root, _, _files) in walkdir(path)
            append!(files, joinpath.(root, _files))
        end
        sort!(files; by=filesize, rev=true)
        info["files"] = [Dict{String,Any}(
            "length"=>filesize(f),
            "path" => splitpath(relpath(f, path)),
            ) for f in files]
    else
        throw(ArgumentError("path `$path` is not an existing file or folder."))
    end

    info["piece length"] = block_size # TODO: Better default block_size to make e.g. 1000 pieces or something.

    # Compute pieces hash
    total_size = sum(filesize(x) for x in files)
    npieces = ceil(Int, total_size/block_size)
    pieces = sizehint!(UInt8[], npieces * 20)
    buffer = Vector{UInt8}(undef, block_size)
    buffer_ptr = 1
    for file in files
        open(file) do io
            while !eof(io)
                nb = readbytes!(io, @view(buffer[buffer_ptr:end]))
                buffer_ptr += nb
                if buffer_ptr === block_size + 1
                    append!(pieces, SHA.sha1(buffer))
                    buffer_ptr = 1
                end
            end
        end
    end
    ## Compute the hash of the remaining bytes (if there are any in the buffer)
    if buffer_ptr > 1
        append!(pieces, SHA.sha1(@view(buffer[1:buffer_ptr-1])))
    end
    info["pieces"] = String(pieces)

    # Create content dictionary
    content = Dict{String, Any}("info" => info)
    ## Add some default values
    content["creation date"] = round(Int, time())
    content["created by"] = "BitTorrent.jl"
    ## Add user meta
    merge!(content, meta)

    # Write to disk
    torrent_file = joinpath(splitpath(path)...) * ".torrent"
    open(io -> bencode(io, content), torrent_file, "w")

    return Torrent(torrent_file, content)
end

# npieces(t::Torrent) = div(sizeof(t.content["info"]["pieces"]), 20)
# piece_length(t::Torrent) = t.content["info"]["piece length"]
# creation_date(t::Torrent) = t.content["info"]["creation date"]

function info_hash(t::Torrent)
    io = IOBuffer()
    bencode(io, t.content["info"])
    sha1 = bytes2hex(SHA.sha1(seekstart(io)))
    return sha1
end

function magnet_uri(t::Torrent)
    btih = info_hash(t)
    return "magnet:?xt=urn:btih:$(btih)"
end


text:
# Designed to be used with PkgBenchmark
# i.e.
# using PkgBenchmark
# import PNGFiles
# benchmarkpkg(PNGFiles)

using BenchmarkTools
using PNGFiles
using PNGFiles.ImageCore

include("../test/test_images/synth_images.jl")

SUITE = BenchmarkGroup()

SUITE["synth_imgs"] = BenchmarkGroup()
for (name, img) in synth_imgs
    SUITE["synth_imgs"][name] = BenchmarkGroup()
    SUITE["synth_imgs"][name]["save"] = @benchmarkable PNGFiles.save("tst.png", $img) teardown=(rm("tst.png"))
    SUITE["synth_imgs"][name]["load"] = @benchmarkable PNGFiles.load("tst.png") setup=(PNGFiles.save("tst.png", $img)) teardown=(rm("tst.png"))
end


text:
for F in (:(+), :(-), :(*), :(/), :(\), :(^))
    @eval begin
        $F(x::T, y::T) where {T<:MarkableFloat64} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::T) where {T<:MarkableFloat32} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::T) where {T<:MarkableFloat16} = T($F(Float16(x), Float16(y)))

        $F(x::T1, y::T2) where {T1<:MarkableFloat64, T2<:MarkableFloat64} =
            promote_type(T1,T2)($F(Float64(x), Float64(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat32, T2<:MarkableFloat32} =
            promote_type(T1,T2)($F(Float32(x), Float32(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat16, T2<:MarkableFloat16} =
            promote_type(T1,T2)($F(Float16(x), Float16(y)))

        $F(x::T, y::R) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::R) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::R) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))

        $F(x::R, y::T) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::R, y::T) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::R, y::T) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))

        $F(x::T, y::R) where {T<:MarkableFloat64, R<:Integer} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::R) where {T<:MarkableFloat32, R<:Integer} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::R) where {T<:MarkableFloat16, R<:Integer} = T($F(Float16(x), Float16(y)))

        $F(x::R, y::T) where {T<:MarkableFloat64, R<:Integer} = T($F(Float64(x), Float64(y)))
        $F(x::R, y::T) where {T<:MarkableFloat32, R<:Integer} = T($F(Float32(x), Float32(y)))
        $F(x::R, y::T) where {T<:MarkableFloat16, R<:Integer} = T($F(Float16(x), Float16(y)))
    end
end

for F in (:ceil, :floor, :trunc, :round)
    @eval begin
        $F(x::T) where {T<:MarkableFloat64} = T($F(Float64(x)))
        $F(x::T) where {T<:MarkableFloat32} = T($F(Float32(x)))
        $F(x::T) where {T<:MarkableFloat16} = T($F(Float16(x)))

        $F(::Type{I}, x::T) where {T<:MarkableFloat64, I<:Integer} = $F(I, Float64(x))
        $F(::Type{I}, x::T) where {T<:MarkableFloat32, I<:Integer} = $F(I, Float32(x))
        $F(::Type{I}, x::T) where {T<:MarkableFloat16, I<:Integer} = $F(I, Float16(x))

        $F(::Type{Integer}, x::T) where {T<:MarkableFloat64} = $F(Int, Float64(x))
        $F(::Type{Integer}, x::T) where {T<:MarkableFloat32} = $F(Int, Float32(x))
        $F(::Type{Integer}, x::T) where {T<:MarkableFloat16} = $F(Int, Float16(x))
    end
end

for F in (:div, :fld, :cld, :rem, :mod)
    @eval begin
        $F(x::T, y::T) where {T<:MarkableFloat64} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::T) where {T<:MarkableFloat32} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::T) where {T<:MarkableFloat16} = T($F(Float16(x), Float16(y)))

        $F(x::T1, y::T2) where {T1<:MarkableFloat64, T2<:MarkableFloat64} =
            promote_type(T1,T2)($F(Float64(x), Float64(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat32, T2<:MarkableFloat32} =
            promote_type(T1,T2)($F(Float32(x), Float32(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat16, T2<:MarkableFloat16} =
            promote_type(T1,T2)($F(Float16(x), Float16(y)))

        $F(x::T, y::R) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::T, y::R) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::T, y::R) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))

        $F(x::R, y::T) where {T<:MarkableFloat64, R<:AbstractFloat} = T($F(Float64(x), Float64(y)))
        $F(x::R, y::T) where {T<:MarkableFloat32, R<:AbstractFloat} = T($F(Float32(x), Float32(y)))
        $F(x::R, y::T) where {T<:MarkableFloat16, R<:AbstractFloat} = T($F(Float16(x), Float16(y)))
    end
end

for F in (:divrem, :fldmod)
    @eval begin
        $F(x::T, y::T) where {T<:MarkableFloat64} = T.($F(Float64(x), Float64(y)))
        $F(x::T, y::T) where {T<:MarkableFloat32} = T.($F(Float32(x), Float32(y)))
        $F(x::T, y::T) where {T<:MarkableFloat16} = T.($F(Float16(x), Float16(y)))

        $F(x::T1, y::T2) where {T1<:MarkableFloat64, T2<:MarkableFloat64} =
            promote_type(T1,T2).($F(Float64(x), Float64(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat32, T2<:MarkableFloat32} =
            promote_type(T1,T2).($F(Float32(x), Float32(y)))
        $F(x::T1, y::T2) where {T1<:MarkableFloat16, T2<:MarkableFloat16} =
            promote_type(T1,T2).($F(Float16(x), Float16(y)))

        $F(x::T, y::R) where {T<:MarkableFloat64, R<:AbstractFloat} = T.($F(Float64(x), Float64(y)))
        $F(x::T, y::R) where {T<:MarkableFloat32, R<:AbstractFloat} = T.($F(Float32(x), Float32(y)))
        $F(x::T, y::R) where {T<:MarkableFloat16, R<:AbstractFloat} = T.($F(Float16(x), Float16(y)))

        $F(x::R, y::T) where {T<:MarkableFloat64, R<:AbstractFloat} = T.($F(Float64(x), Float64(y)))
        $F(x::R, y::T) where {T<:MarkableFloat32, R<:AbstractFloat} = T.($F(Float32(x), Float32(y)))
        $F(x::R, y::T) where {T<:MarkableFloat16, R<:AbstractFloat} = T.($F(Float16(x), Float16(y)))
    end
end


text:
immutable Levy <: ContinuousUnivariateDistribution
    μ::Float64
    c::Float64

    function Levy(μ::Real, c::Real)
        c >= zero(c) || error("scale must be non-negative")
        @compat new(Float64(μ), Float64(c))
    end

    Levy(μ::Real) = new(μ, 1.0)
    Levy() = new(0.0, 1.0)
end

@distr_support Levy d.location Inf

#### Parameters

location(d::Levy) = d.μ
params(d::Levy) = (d.μ, d.c)

#### Statistics

mean(d::Levy) = Inf
var(d::Levy) = Inf
skewness(d::Levy) = NaN
kurtosis(d::Levy) = NaN

mode(d::Levy) = d.c / 3.0 + d.μ

function entropy(d::Levy)
    c = scale(d)
    (1.0 - 3.0 * digamma(1.0) + log(16.0 * pi * c * c)) / 2.0
end

function median(d::Levy)
    μ, c = params(d)
    μ + c / (2.0 * erfcinv(0.5)^2)
end

#### Evaluation

function pdf(d::Levy, x::Float64)
    μ, c = params(d)
    z = x - μ
    (sqrt(c) / sqrt2π) * exp((-c) / (2.0 * z)) / z^1.5
end

function logpdf(d::Levy, x::Float64)
    μ, c = params(d)
    z = x - μ
    0.5 * (log(c) - log2π - c / z - 3.0 * log(z))
end

cdf(d::Levy, x::Float64) = erfc(sqrt(d.c / (2.0 * (x - d.μ))))
ccdf(d::Levy, x::Float64) = erf(sqrt(d.c / (2.0 * (x - d.μ))))

quantile(d::Levy, p::Float64) = d.μ + d.c / (2.0 * erfcinv(p)^2)
cquantile(d::Levy, p::Float64) = d.μ + d.c / (2.0 * erfinv(p)^2)

mgf(d::Levy, t::Real) = t == zero(t) ? 1.0 : NaN

function cf(d::Levy, t::Real)
    μ, c = params(d)
    exp(im * μ * t - sqrt(-2.0 * im * c * t))
end

#### Sampling

rand(d::Levy) = d.μ + d.c / randn()^2


text:
using PyCall

function __init__()
    py"""
    import basis_set_exchange as bse

    def dump_json(basis, path=None):
        if path is None:
            path = "."
        bse.fileio.write_json_basis(f"{path}/{basis}.json", bse.get_basis(basis))
    """
end

dump_json(basis::String, path::String) = py"dump_json"(basis, path)
dump_json(basis::String) = py"dump_json"(basis)


text:
"""
# LocalRegistry

Create and maintain local registries for Julia packages. This package
is intended to provide a simple but manual workflow for maintaining
small local registries (private or public) without making any
assumptions about how the registry or the packages are hosted.

Registry creation is done by the function `create_registry`.

Registration of new and updated packages is done by the function `register`.
"""
module LocalRegistry

using RegistryTools: RegistryTools, gitcmd, Compress
using UUIDs: uuid4
using LibGit2
using Pkg: Pkg, TOML

export create_registry, register

"""
    create_registry(name, repo; description = nothing)
    create_registry(path, repo; description = nothing)

Create a registry with the given `name` or at the local directory
`path`, and with repository URL `repo`. Optionally add a description
of the purpose of the registry with the keyword argument
`description`. The first argument is interpreted as a path if it has
more than one path component and otherwise as a name. If a path is
given, the last path component is used as the name of the registry. If
a name is given, it is created in the standard registry location. In
both cases the registry path must not previously exist.

Note: This will only prepare the registry locally. Review the result
and `git push` it manually.
"""
function create_registry(name_or_path, repo; description = nothing,
                         gitconfig::Dict = Dict(), uuid = nothing)
    if length(splitpath(name_or_path)) > 1
        path = abspath(expanduser(name_or_path))
    else
        path = joinpath(first(DEPOT_PATH), "registries", name_or_path)
    end
    name = basename(path)
    if isempty(name)
        # If `path` ends with a slash, `basename` becomes empty and
        # the slash needs to be peeled off with `dirname`.
        name = basename(dirname(path))
    end
    if isdir(path) || isfile(path)
        throw("$path already exists.")
    end
    mkpath(path)

    # The only reason for the `uuid` function argument is to allow
    # deterministic testing of the package.
    if isnothing(uuid)
        uuid = string(uuid4())
    end

    registry_data = RegistryTools.RegistryData(name, uuid, repo = repo,
                                               description = description)
    RegistryTools.write_registry(joinpath(path, "Registry.toml"), registry_data)

    git = gitcmd(path, gitconfig)
    run(`$git init -q`)
    run(`$git add Registry.toml`)
    run(`$git commit -qm 'Create registry.'`)
    run(`$git remote add origin $repo`)
    @info "Created registry in directory $(path)"

    return path
end

"""
    register(package, registry)

Register the new `package` in the registry `registry`. The version
number and other information is obtained from the package's
`Project.toml`.

    register(package)

Register a new version of `package`.  The version number is obtained
from the package's `Project.toml`.

Note: In both cases this will only update the registry locally. Review
the result and `git push` it manually. The package must live in a git
working copy, e.g. having been cloned by `Pkg.develop`.

`package` can be specified in the following ways:
* By package name. The package must be available in the active `Pkg`
  environment.
* By path. This is distinguished from package name by having more than
  one path component. A path in the current working directory can be
  specified by starting with `"./"`.
* By module. It needs to first be loaded by `using` or `import`.

`registry` can be specified by name or by path in the same way as
`package`. If omitted or `nothing`, it will be automatically located
by `package`.

*Keyword arguments*

    register(package, registry; commit = true, repo = nothing, gitconfig = Dict())

* `commit`: If `false`, only make the changes to the registry but do not commit.
* `repo`: Specify the package repository explicitly. Otherwise looked up as the `git remote` of the package.
* `gitconfig`: Optional configuration parameters for the `git` command.
"""
function register(package::Union{Module, AbstractString},
                  registry::Union{Nothing, AbstractString} = nothing;
                  repo = nothing, commit = true, gitconfig::Dict = Dict())
    # Find and read the `Project.toml` for the package.
    package_path = find_package_path(package)
    pkg = Pkg.Types.read_project(joinpath(package_path, "Project.toml"))
    if isnothing(pkg.name)
        throw("$(package) does not have a Project.toml file")
    end

    registry_path = find_registry_path(registry, pkg)
    if LibGit2.isdirty(LibGit2.GitRepo(registry_path))
        throw("Registry repo is dirty. Stash or commit files.")
    end

    # Compute the tree hash for the package.
    tree_hash = get_tree_hash(package_path, gitconfig)

    # Use the `repo` argument or check for the git remote if not provided.
    if isnothing(repo)
        package_repo = get_remote_repo(package_path, gitconfig)
    else
        package_repo = repo
    end

    @info "Registering package" package_path registry_path package_repo uuid=pkg.uuid version=pkg.version tree_hash
    clean_registry = true

    git = gitcmd(registry_path, gitconfig)
    try
        package_path = get_package_path(registry_path, pkg)
        update_package_data(pkg, registry_path, package_repo, package_path, tree_hash)
        if commit
            commit_registry(pkg, package_path, package_repo, tree_hash, git)
        end
    finally
        if clean_registry
            run(`$git reset --hard`)
        end
    end

    return
end

# If the package is provided as a module, directly find the package
# path from the loaded code. This works both if the module is loaded
# from the current package environment or found in LOAD_PATH.
function find_package_path(package::Module)
    return dirname(dirname(pathof(package)))
end

# A string argument is either interpreted as a path or as a package
# name, decided by the number of components returned by `splitpath`.
#
# If the package is given by name, it must be available in the current
# package environment as a developed package.
function find_package_path(package_name::AbstractString)
    if length(splitpath(package_name)) > 1
        if !isdir(package_name)
            error("Package path $(package_name) does not exist.")
        end
        return abspath(expanduser(package_name))
    end

    ctx = Pkg.Types.Context()
    if !haskey(ctx.env.project.deps, package_name)
        error("Unknown package $package_name.")
    end
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    if isnothing(pkg_path)
        error("Package must be developed to be registered.")
    end
    return pkg_path
end

function find_registry_path(registry::AbstractString, pkg)
    if length(splitpath(registry)) > 1
        return abspath(expanduser(registry))
    end

    all_registries = Pkg.Types.collect_registries()
    matching_registries = filter(r -> r.name == registry, all_registries)
    if isempty(matching_registries)
        error("Registry $(registry) not found.")
    end

    return first(matching_registries).path
end

function find_registry_path(registry::Nothing, pkg)
    all_registries = Pkg.Types.collect_registries()

    matching_registries = filter(all_registries) do reg_spec
        reg_data = Pkg.Types.read_registry(joinpath(reg_spec.path,
                                                    "Registry.toml"))
        haskey(reg_data["packages"], string(pkg.uuid))
    end

    if isempty(matching_registries)
        error("Package $(pkg.name) not found in any registry. Please specify in which registry you want to register it.")
    elseif length(matching_registries) > 1
        error("Package $(pkg.name) is registered in more than one registry, please specify in which you want to register the package.")
    end

    return first(matching_registries).path
end

function get_tree_hash(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    return read(`$git log --pretty=format:%T -1`, String)
end

function get_remote_repo(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    remote_name = split(readchomp(`$git remote`), '\n')
    length(remote_name) > 1 && throw("Repository has multiple remotes.")
    remote_name[1] == "" && throw("Repository does not have a remote.")
    return readchomp(`$git remote get-url $(remote_name[1])`)
end

# Find the package in the repository and return its path. If not
# previously registered, first create its directory.
function get_package_path(registry_path, pkg::Pkg.Types.Project)
    @debug("find package in registry")
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = RegistryTools.parse_registry(registry_file)

    uuid = string(pkg.uuid)
    if haskey(registry_data.packages, uuid)
        package_data = registry_data.packages[uuid]
        if package_data["name"] != pkg.name
            err = "Changing package names not supported yet"
            @debug(err)
            throw(err)
        end
        package_path = joinpath(registry_path, package_data["path"])
    else
        @debug("Package with UUID: $uuid not found in registry, checking if UUID was changed")
        if pkg.name in (v["name"] for (k, v) in registry_data.packages)
            err = "Changing UUIDs is not allowed"
            @debug(err)
            throw(err)
        end

        @debug("Creating directory for new package $(pkg.name)")
        first_letter = uppercase(pkg.name[1:1])
        relative_package_path = joinpath(first_letter, pkg.name)
        package_path = joinpath(registry_path, relative_package_path)
        mkpath(package_path)

        # Do not use joinpath to construct the stored path. Should be
        # a forward slash independent of platform.
        registry_data.packages[uuid] = Dict("name" => pkg.name,
                                            "path" => string(first_letter, "/",
                                                             pkg.name))
        RegistryTools.write_registry(registry_file, registry_data)
    end

    return package_path
end

function update_package_data(pkg::Pkg.Types.Project, registry_path,
                             package_repo, package_path, tree_hash)
    # Package file.
    RegistryTools.update_package_file(pkg, package_repo, package_path)

    # Versions file.
    update_versions_file(pkg, package_path, tree_hash)

    # Dependencies file.
    check_dependencies(pkg, registry_path)
    update_deps_file(pkg, package_path)

    # Compatibilities file.
    update_compat_file(pkg, package_path)

    return
end

# Same code as in RegistryTools but without the RegBranch version checks.
function update_versions_file(pkg::Pkg.Types.Project,
                              package_path::AbstractString,
                              tree_hash::AbstractString)
    versions_file = joinpath(package_path, "Versions.toml")
    versions_data = isfile(versions_file) ? TOML.parsefile(versions_file) : Dict()
    versions = sort!([VersionNumber(v) for v in keys(versions_data)])

    version_info = Dict{String,Any}("git-tree-sha1" => string(tree_hash))
    versions_data[string(pkg.version)] = version_info

    open(versions_file, "w") do io
        TOML.print(io, versions_data; sorted=true, by=x->VersionNumber(x))
    end
    return
end

function check_dependencies(pkg::Pkg.Types.Project, registry_path)
    @debug("Verifying package name and uuid in deps")
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = TOML.parsefile(registry_file)
    for (k, v) in pkg.deps
        u = string(v)
        if haskey(registry_data["packages"], u)
            name_in_reg = registry_data["packages"][u]["name"]
            if name_in_reg != k
                err = "Error in `[deps]`: UUID $u refers to package '$name_in_reg' in registry but deps file has '$k'"
                throw(err)
            end
        elseif haskey(BUILTIN_PKGS, k)
            if BUILTIN_PKGS[k] != u
                err = "Error in `[deps]`: UUID $u for package $k should be $(BUILTIN_PKGS[k])"
                throw(err)
            end
        else
#             err = "Error in `[deps]`: Package '$k' with UUID: $u not found in registry or stdlib"
#             throw(err)
        end
    end
end

function update_deps_file(pkg::Pkg.Types.Project, package_path)
    deps_file = joinpath(package_path, "Deps.toml")
    if isfile(deps_file)
        deps_data = Compress.load(deps_file)
    else
        deps_data = Dict()
    end

    deps_data[pkg.version] = pkg.deps
    Compress.save(deps_file, deps_data)
end

function update_compat_file(pkg::Pkg.Types.Project, package_path)
    compat_file = joinpath(package_path, "Compat.toml")
    if isfile(compat_file)
        compat_data = Compress.load(compat_file)
    else
        compat_data = Dict()
    end

    d = Dict()
    for (n, v) in pkg.compat
        spec = Pkg.Types.semver_spec(v)
        # The call to map(versionrange, ) can be removed once support
        # for older Julia versions than 1.3 are dropped.
        ranges = map(r->versionrange(r.lower, r.upper), spec.ranges)
        ranges = Pkg.Types.VersionSpec(ranges).ranges # this combines joinable ranges
        d[n] = length(ranges) == 1 ? string(ranges[1]) : map(string, ranges)
    end
    compat_data[pkg.version] = d
    Compress.save(compat_file, compat_data)
end

function versionrange(lo::Pkg.Types.VersionBound, hi::Pkg.Types.VersionBound)
    lo.t == hi.t && (lo = hi)
    return Pkg.Types.VersionRange(lo, hi)
end

function commit_registry(pkg::Pkg.Types.Project, package_path, package_repo, tree_hash, git)
    @debug("commit changes")
    message = """
    New version: $(pkg.name) v$(pkg.version)

    UUID: $(pkg.uuid)
    Repo: $(package_repo)
    Tree: $(string(tree_hash))
    """
    run(`$git add --all`)
    if occursin("nothing to commit, working tree clean",
                read(`$git status`, String))
        @info("nothing to commit")
    else
        run(`$git commit -qm $message`)
    end
end

const BUILTIN_PKGS = Dict("Pkg"=>"44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
                          "Statistics"=>"10745b16-79ce-11e8-11f9-7d13ad32a3b2",
                          "Test"=>"8dfed614-e22c-5e08-85e1-65c5234f0b40",
                          "CRC32c"=>"8bf52ea8-c179-5cab-976a-9e18b702a9bc",
                          "Random"=>"9a3f8284-a2c9-5f02-9a11-845980a1fd5c",
                          "Libdl"=>"8f399da3-3557-5675-b5ff-fb832c97cbdb",
                          "UUIDs"=>"cf7118a7-6976-5b1a-9a39-7adc72f591a4",
                          "Distributed"=>"8ba89e20-285c-5b6f-9357-94700520ee1b",
                          "Serialization"=>"9e88b42a-f829-5b0c-bbe9-9e923198166b",
                          "DelimitedFiles"=>"8bb1440f-4735-579b-a4ab-409b98df4dab",
                          "LinearAlgebra"=>"37e2e46d-f89d-539d-b4ee-838fcccc9c8e",
                          "FileWatching"=>"7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee",
                          "SharedArrays"=>"1a1011a3-84de-559e-8e89-a11a2f7dc383",
                          "Base64"=>"2a0f44e3-6c83-55bd-87e4-b1978d98bd5f",
                          "SparseArrays"=>"2f01184e-e22b-5df5-ae63-d93ebab69eaf",
                          "Profile"=>"9abbd945-dff8-562f-b5e8-e1ebf5ef1b79",
                          "Mmap"=>"a63ad114-7e13-5084-954f-fe012c677804",
                          "Unicode"=>"4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5",
                          "InteractiveUtils"=>"b77e0a4c-d291-57a0-90e8-8db25a27a240",
                          "Future"=>"9fa8497b-333b-5362-9e8d-4d0656e87820",
                          "Sockets"=>"6462fe0b-24de-5631-8697-dd941f90decc",
                          "Printf"=>"de0858da-6303-5e67-8744-51eddeeeb8d7")

end


text:
"""
    compare_summary(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1}; kwargs...)

# Arguments
- `df::Array{DataFrame,1}`: List of DataFrames to compare.
    These must all share the same column names.
- `d::Array{Dict,1}`: Array of dictonaries of DataFrames to compare.
- `id::Array{Symbol,1}`: List of IDs that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and IDs.

# Keywords
- `tol::Float64 = 1E-6`: Tolerance used when determining whether values are equal.
    Default values is `1E-6`.
- `complete_summary::Bool = false`: Should we include the full summary DataFrame or only the
    problematic rows (i.e., ones where either `equal_keys = false` or `equal_values = false`)?

# Returns
- `df::DataFrame`: Summary DataFrame including the index columns and values from the original
    `df`, marked to indicate which input DataFrame they came from, the `reldiff` between
    the values, and booleans `equal_key` (are values for the index present in both DataFrames?)
    and `equal_value` (are the values the same within the given `tol`?)
- `d::Dict`: Dictionary of summary DataFrames
"""
function compare_summary(
    df::Array{DataFrame,1},
    id::Array{Symbol,1};
    tol::Float64=DEFAULT_TOL,
    complete_summary::Bool=false
)
    df = copy.(df)
    rel = length(df) == 2 ? :reldiff : :maxreldiff

    # Do some checks on the indices before comparing.
    idx = findindex.(df)
    # !!!! If columns are different, print a warning message, but make them consistent and
    # press on!
    length(unique(sort.(idx))) > 1 && @error("Cannot compare DataFrames with different indices.")
    length(unique(idx)) > 1 && @warn("Comparing DataFrames with different index orders.")

    df = indexjoin(df; id=id, indicator=true, fillmissing=true)
    idx = findindex(df)
    val = setdiff(findvalue(df), id)

    # Are all of the keys here?
    df[!,:equal_keys] .= prod.(eachrow(df[:,id]))

    if !isempty(val)
        # Are there discrepancies between PRESENT values (within the specified tolerance)?
        # All values in a row x will be considered "equal" if (max(x) - x_i) / mean(x) < tol
        df_comp = abs.(df[:,val])
        df_comp = (maximum.(skipmissing.(eachrow(df_comp))) .- df_comp) ./
            Statistics.mean.(skipmissing.(eachrow(df_comp)))

        df[!,rel] .= maximum.(skipmissing.(eachrow(df_comp)))
        df[!,:equal_values] .= df[:,rel] .<= tol

        # What if some zeros were include but not others?
        # df[all.(eachrow(.|(ismissing.(df[:,val]), df[:,val].==0))), :equal_values] .= true
        select!(df, [idx; val; id; rel; :equal_keys; :equal_values])
        ii = df[:,:equal_keys] .* df[:,:equal_values]
    else
        select!(df, [idx; id; :equal_keys])
        ii = df[:,:equal_keys]
    end

    !complete_summary && (df = df[.!ii,:])
    return df
end

function compare_summary(
    d::Array{Dict{Symbol,DataFrame},1},
    id::Array{Symbol,1};
    tol::Float64=DEFAULT_TOL,
    complete_summary::Bool=false
)
    keys_comp = intersect(collect.(keys.(d))...)

    length(keys_comp) == 0 && @warn("Cannot compare Dictionaries that share no common keys.")

    d = Dict(k => compare_summary([d[m][k] for m in keys(d)], id;
        tol=tol, complete_summary=complete_summary) for k in keys_comp)
    return d
end

"""
    compare_values(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})

# Arguments
- df_lst::Array{DataFrame,1}: List of DataFrames to compare.
    These must all share the same column names.
- `inds::Array{Symbol,1}`: List of IDs that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and IDs.

# Keyword Argument
- `tol::Float64 = 1E-6`: Tolerance used when determining whether values are equal.
    Default values is `1E-6`.
"""
function compare_values(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1}; tol=DEFAULT_TOL)
    df_lst = copy.(df_lst)
    df = compare_summary(copy.(df_lst), inds; tol=tol)
    df = df[.!df[:,:equal_values],:]

    # size(df, 1) > 0 && @warn("Inconsistent values:", df)
    return df
end

"""
    compare_keys(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})

# Arguments
- df_lst::Array{DataFrame,1}: List of DataFrames to compare.
    These must all share the same column names.
- `inds::Array{Symbol,1}`: List of IDs that describe each DataFrame and track which
    values/keys are present in each DataFrame. There must be an equal number of input
    DataFrames and IDs.
"""
function compare_keys(df_lst::Array{DataFrame,1}, inds::Array{Symbol,1})
    df_lst = copy.(df_lst)
    N = length(inds)
    cols = intersect(find_oftype.(df_lst, Not(AbstractFloat))...)

    # Determine whether to consider case when comparing keys. Only consider case if there
    # instances of the same keys with differing cases in the same DataFrame.
    d_unique = Dict(col => Dict(inds[ii] => sort(unique(df_lst[ii][:,col]))
        for ii in 1:N) for col in cols)
    d_lower = Dict(col => Dict(inds[ii] => lowercase.(d_unique[col][inds[ii]])
        for ii in 1:N) for col in cols);
    CHECKCASE = Dict(col => any(length.(unique.(values(d_lower[col]))) .!==
        length.(values(d_unique[col]))) for col in cols)

    d_all = Dict(col => CHECKCASE[col] ? sort(unique([values(d_unique[col])...;])) :
        sort(unique([values(d_lower[col])...;])) for col in cols)

    df = DataFrame()

    for col in cols
        df_temp = DataFrame(key=fill(col, size(d_all[col])))
        d_check = CHECKCASE[col] ? d_unique[col] : d_lower[col]

        [df_temp[!,ind] = [v in d_check[ind] ? d_unique[col][ind][v .== d_check[ind]][1] :
            missing for v in d_all[col]] for ind in inds]
        df_temp = unique(df_temp[length.(unique.(eachrow(df_temp[:,inds]))) .> 1, :])

        df = [df; df_temp]
    end

    # size(df, 1) > 0 && @warn("Inconsistent keys:", df)
    return df
end

"""
    benchmark_against(df_calc, df_bench; kwargs...)
"""
function benchmark_against(
    df_calc::DataFrame,
    df_bench::DataFrame;
    key=missing,
    tol=DEFAULT_TOL,
    small=missing
)
    # # Make sure that all input DataFrames have the same columns.
    # # !!!! Todo: keep units where they appear!
    # col = intersect(propertynames(df_calc), propertynames(df_bench))
    # select!(df_calc, col)
    # select!(df_bench, col)

    # Remove very small numbers. These might be zero or missing in the other DataFrame,
    # and we're not splitting hairs here.
    if small !== missing
        df_calc = df_calc[abs.(df_calc[:,:value] .> small), :]
        df_bench = df_bench[abs.(df_bench[:,:value] .> small), :]
    end

    !ismissing(key) && println("  Comparing keys and values for ", key)

    df_comp = compare_summary([df_calc, df_bench], [:calc,:bench]; tol=tol)
    # key == :utd && (df_comp = edit_with(df_comp, Drop(:yr,2002,"<")))

    # If the dataframes are in agreement, store this value as "true".
    # Otherwise, store the comparison dataframe rows that are not in agreement.
    return size(df_comp, 1) == 0 ? true : df_comp
end

function benchmark_against(calc::Dict, bench::Dict; tol=DEFAULT_TOL, small=DEFAULT_SMALL)
    keys_comp = intersect(keys(calc), keys(bench))

    length(keys_comp) == 0 && @warn("Cannot compare Dictionaries that share no common keys.")

    return Dict(k => benchmark_against(calc[k], bench[k]; key=k, tol=tol) for k in keys_comp)
end

"""
    verify_over(df::DataFrame, col::Any; tol = 1E-6)
"""
function verify_over(df::DataFrame, col::Any; tol=DEFAULT_TOL)
    df = combine_over(df, col)
    df = df[(df[:,:value] .- 1.0) .> tol, :]
    return size(df, 1) == 0 ? true : df
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
#!/usr/bin/julia

a = time()

s = 1
for i=1:100000000
    global s*=i
end

b = time();

print(b-a)
print("\n")


text:
### Function to analyze the geometry
###
function voronoi_perimeter(ver, region)
    # compute the perimetet  in the voronoi region and vertices

    perimeter = 0.
    sizeRegion = length(region)

    if -1 in region
        perimeter = 0.
    else
        for index in 1:sizeRegion
            index2 = (index%sizeRegion) + 1

            if region[index]!= -1 || region[index2] != -1

                x0 = ver[region[index]+1,1]
                y0 = ver[region[index]+1,2]

                x1 = ver[region[index2]+1,1]
                y1 = ver[region[index2]+1,2]
                perimeter += sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
            end
        end
    end
    return(perimeter)
end
##
function voronoi_area(ver, region)
    # compute the area   in the voronoi region and vertices

    area = 0.
    sizeRegion = length(region)

    if -1 in region
        area = 0.
    else
        for index in 1:sizeRegion
            x0 = ver[region[index]+1,1]
            y0 = ver[region[index]+1,2]

            index2 = (index%sizeRegion) + 1
            x1 = ver[region[index2]+1,1]
            y1 = ver[region[index2]+1,2]

            area += 0.5 * (x0*y1 - x1*y0)
        end
    end
    return(abs(area))
end
## Voronoi tesselation using the scipy function
##
function voronoi_python(pts, verbose = true)
    let
        ndat = length(pts)
        peri = zeros(ndat)
        area = zeros(ndat)

        vor = 0
        spatial = pyimport("scipy.spatial")

        try
            # vor = spatial[:Voronoi](pts)
            vor = spatial.Voronoi(pts)
        catch
            println("## Voronoi error...")
            return([1000,1000],[1000,1000])   ## arbitrary values for peri and area
        end

        if verbose println("## Voronoi tesselation done.") end

        ver = vor.vertices
        reg = vor.regions
        pt  = vor.point_region

        for i in 1:ndat
            region  =  reg[pt[i]+1]
            peri[i] =  voronoi_perimeter(ver,region)
            area[i] =  voronoi_area(ver,region)
        end
        println("ending voronoi ...")
        return(peri , area)
    end
end
##
## Julia voronoi implementation
function voronoi(pts, verbose = false)

        ndat = length(pts)
        peri = zeros(ndat)
        area = zeros(ndat)

        xx= [] ; yy=[]
        for i in 1:ndat
            push!(xx,pts[i][1])
            push!(yy,pts[i][2])
        end

        xmin= minimum(xx)
        xmax= maximum(xx)
        ymin= minimum(yy)
        ymax= maximum(yy)

        rect = Rectangle(gb.Point2(xmin,ymin),gb.Point2(xmax,ymax))
        ptarr= [gb.Point2(xx[i], yy[i]) for i in 1:ndat]

        tess = voronoicells(ptarr, rect)
        area= voronoiarea(tess)

        return(area , area)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
@generated function check_size_parameters(::Type{S}) where {S}
    if !all(x -> isa(x, Int) || x <: Symmetry, S.parameters)
        return :(throw(ArgumentError("Tensor's size parameter must be a tuple of `Int`s or `Symmetry` (e.g. `Tensor{Tuple{3,3}}` or `Tensor{Tuple{3, @Symmetry{3,3}}}`).")))
    end
end

@generated function flatten_tuple(x::Tuple{Vararg{Tuple, N}}) where {N}
    exps = [Expr(:..., :(x[$i])) for i in 1:N]
    :(tuple($(exps...)))
end

@inline fill_tuple(x, ::Val{N}) where {N} = ntuple(i -> x, Val(N))
@inline fill_tuple(f::Function, ::Val{N}) where {N} = ntuple(i -> f(), Val(N))

@inline convert_ntuple(::Type{T}, x::NTuple{N, T}) where {N, T} = x
@generated function convert_ntuple(::Type{T}, x::NTuple{N, Any}) where {N, T}
    exps = [:(convert(T, x[$i])) for i in 1:N]
    quote
        @_inline_meta
        tuple($(exps...))
    end
end

@inline promote_ntuple_eltype(x::NTuple{N, T}) where {N, T} = T
@generated function promote_ntuple_eltype(x::NTuple{N, Any}) where {N}
    T = promote_type(x.parameters...)
    quote
        @_inline_meta
        $T
    end
end

tuple_sort(x::NTuple{N, Any}; kwargs...) where {N} = Tuple(sort(SVector(x); kwargs...))


text:
using Serialization
## 1. Julia-only serialization. Use it for technical data (e.g. functions, codes, etc.)

"Return serialized object as an array of bytes"
function serialized(x)
    buf = IOBuffer()
    serialize(buf, x)
    return buf.data
end

"Return object deserialized from array of bytes"
function deserialized(x::Vector{UInt8})
    return deserialize(IOBuffer(x))
end


text:
module day9

readInput() = parse.(Int, eachline(joinpath(@__DIR__, "input.txt")))

function isValid(target, list)
    for i in 1:length(list) - 1, j in i + 1:length(list)
        list[i] + list[j] == target && return true
    end
    false
end

function part1(input, width=25)
    i = width + 1
    while isValid(input[i], view(input, (i - width):(i - 1)))
        i += 1
    end
    input[i]
end

function findRange(target, list)
    i, j = 1, 2
    sum = list[i] + list[j]
    while sum != target
        if sum > target && j - 1 > 1
            sum -= list[i]
            i += 1
        elseif sum < target
            j += 1
            sum += list[j]
        end
    end
    return i, j
end

const resultP1 = part1(readInput()) # To leave part1 out of the benchmark

function part2(input, target=resultP1)
    i, j = findRange(target, input)
    +(extrema(view(input, i:j))...)
end

end


text:
@testset "constructors" begin

    # We can't use test_rrule here (as it's currently implemented) because the elements of
    # the array have arbitrary values. The only thing we can do is ensure that we're getting
    # `ZeroTangent`s back, and that the forwards pass produces the correct thing still.
    # Issue: https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/202
    @testset "undef" begin
        val, pullback = rrule(Array{Float64}, undef, 5)
        @test size(val) == (5, )
        @test val isa Array{Float64, 1}
        @test pullback(randn(5)) == (NoTangent(), NoTangent(), NoTangent())
    end
    @testset "from existing array" begin
        test_rrule(Array, randn(2, 5))
        test_rrule(Array, Diagonal(randn(5)))
        test_rrule(Matrix, Diagonal(randn(5)))
        test_rrule(Matrix, transpose(randn(4)))
        test_rrule(Array{ComplexF64}, randn(3))
    end
end

@testset "vect" begin
    test_rrule(Base.vect)
    @testset "homogeneous type" begin
        test_rrule(Base.vect, (5.0, ), (4.0, ))
        test_rrule(Base.vect, 5.0, 4.0, 3.0)
        test_rrule(Base.vect, randn(2, 2), randn(3, 3))
    end
    @testset "inhomogeneous type" begin
        test_rrule(
            Base.vect, 5.0, 3f0;
            atol=1e-6, rtol=1e-6, check_inferred=VERSION>=v"1.6",
        ) # tolerance due to Float32.
        test_rrule(Base.vect, 5.0, randn(3, 3); check_inferred=false)
        test_rrule(Base.vect, (5.0, 4.0), (y=randn(3),); check_inferred=false)
    end
end

@testset "reshape" begin
    test_rrule(reshape, rand(4, 5), (2, 10))
    test_rrule(reshape, rand(4, 5), 2, 10)
    test_rrule(reshape, rand(4, 5), 2, :)
end

@testset "repeat" begin

    test_rrule(repeat, rand(4, ))
    test_rrule(repeat, rand(4, 5))
    test_rrule(repeat, rand(4, 5); fkwargs = (outer=(1,2),))
    test_rrule(repeat, rand(4, 5); fkwargs = (inner=(1,2), outer=(1,3)))

    test_rrule(repeat, rand(4, ), 2; check_inferred=VERSION>=v"1.6")
    test_rrule(repeat, rand(4, 5), 2; check_inferred=VERSION>=v"1.6")
    test_rrule(repeat, rand(4, 5), 2, 3; check_inferred=VERSION>=v"1.6")
    test_rrule(repeat, rand(1,2,3), 2,3,4; check_inferred=VERSION>v"1.6")
    test_rrule(repeat, rand(0,2,3), 2,0,4; check_inferred=VERSION>v"1.6")
    test_rrule(repeat, rand(1,1,1,1), 2,3,4,5; check_inferred=VERSION>v"1.6")

    if VERSION>=v"1.6"
        # These are cases where repeat itself fails in earlier versions
        test_rrule(repeat, rand(4, 5); fkwargs = (inner=(2,4), outer=(1,1,1,3)))
        test_rrule(repeat, rand(1,2,3), 2,3)
        test_rrule(repeat, rand(1,2,3), 2,3,4,2)
        test_rrule(repeat, fill(1.0), 2)
        test_rrule(repeat, fill(1.0), 2, 3)

        # These fail for other v1.0 related issues (add!!)
        # v"1.0": fill(1.0) + fill(1.0) != fill(2.0)
        # v"1.6: fill(1.0) + fill(1.0) == fill(2.0) # Expected
        test_rrule(repeat, fill(1.0); fkwargs = (inner=2,))
        test_rrule(repeat, fill(1.0); fkwargs = (inner=2, outer=3,))

    end

    @test rrule(repeat, [1,2,3], 4)[2](ones(12))[2] == [4,4,4]
    @test rrule(repeat, [1,2,3], outer=4)[2](ones(12))[2] == [4,4,4]

end

@testset "hcat" begin
    test_rrule(hcat, randn(3, 2), randn(3), randn(3, 3); check_inferred=VERSION>v"1.1")
    test_rrule(hcat, rand(), rand(1,2), rand(1,2,1); check_inferred=VERSION>v"1.1")
    test_rrule(hcat, rand(3,1,1,2), rand(3,3,1,2); check_inferred=VERSION>v"1.1")

    # mix types
    test_rrule(hcat, rand(2, 2), rand(2, 2)'; check_inferred=VERSION>v"1.1")
end

@testset "reduce hcat" begin
    mats = [randn(3, 2), randn(3, 1), randn(3, 3)]
    test_rrule(reduce, hcat, mats)

    vecs = [rand(3) for _ in 1:4]
    test_rrule(reduce, hcat, vecs)

    mix = AbstractVecOrMat[rand(4,2), rand(4)]  # this is weird, but does hit the fast path
    test_rrule(reduce, hcat, mix)

    adjs = vec([randn(2, 4), randn(1, 4), randn(3, 4)]')  # not a Vector
    # test_rrule(reduce, hcat, adjs ⊢ map(m -> rand(size(m)), adjs))
    dy = 1 ./ reduce(hcat, adjs)
    @test rrule(reduce, hcat, adjs)[2](dy)[3] ≈ rrule(reduce, hcat, collect.(adjs))[2](dy)[3]

    # mix types
    mats = [randn(2, 2), rand(2, 2)']
    test_rrule(reduce, hcat, mats; check_inferred=VERSION>v"1.1")
end

@testset "vcat" begin
    test_rrule(vcat, randn(2, 4), randn(1, 4), randn(3, 4); check_inferred=VERSION>v"1.1")
    test_rrule(vcat, rand(), rand(); check_inferred=VERSION>v"1.1")
    test_rrule(vcat, rand(), rand(3), rand(3,1,1); check_inferred=VERSION>v"1.1")
    test_rrule(vcat, rand(3,1,2), rand(4,1,2); check_inferred=VERSION>v"1.1")

    # mix types
    test_rrule(vcat, rand(2, 2), rand(2, 2)'; check_inferred=VERSION>v"1.1")
end

@testset "reduce vcat" begin
    mats = [randn(2, 4), randn(1, 4), randn(3, 4)]
    test_rrule(reduce, vcat, mats)

    vecs = [rand(2), rand(3), rand(4)]
    test_rrule(reduce, vcat, vecs)

    mix = AbstractVecOrMat[rand(4,1), rand(4)]
    test_rrule(reduce, vcat, mix)
end

@testset "cat" begin
    test_rrule(cat, rand(2, 4), rand(1, 4); fkwargs=(dims=1,), check_inferred=VERSION>v"1.1")
    test_rrule(cat, rand(2, 4), rand(2); fkwargs=(dims=Val(2),), check_inferred=VERSION>v"1.1")
    test_rrule(cat, rand(), rand(2, 3); fkwargs=(dims=[1,2],), check_inferred=VERSION>v"1.1")
    test_rrule(cat, rand(1), rand(3, 2, 1); fkwargs=(dims=(1,2),), check_inferred=false) # infers Tuple{Zero, Vector{Float64}, Any}

    test_rrule(cat, rand(2, 2), rand(2, 2)'; fkwargs=(dims=1,), check_inferred=VERSION>v"1.1")
end

@testset "hvcat" begin
    test_rrule(hvcat, 2, rand(ComplexF64, 6)...; check_inferred=VERSION>v"1.1")
    test_rrule(hvcat, (2, 1), rand(), rand(1,1), rand(2,2); check_inferred=VERSION>v"1.1")
    test_rrule(hvcat, 1, rand(3)' ⊢ rand(1,3), transpose(rand(3)) ⊢ rand(1,3); check_inferred=VERSION>v"1.1")
    test_rrule(hvcat, 1, rand(0,3), rand(2,3), rand(1,3,1); check_inferred=VERSION>v"1.1")

    # mix types (adjoint and transpose)
    test_rrule(hvcat, 1, rand(3)', transpose(rand(3)) ⊢ rand(1,3); check_inferred=VERSION>v"1.1")
end

@testset "reverse" begin
    # Forward
    test_frule(reverse, rand(5))
    test_frule(reverse, rand(5), 2, 4)
    test_frule(reverse, rand(5), fkwargs=(dims=1,))

    test_frule(reverse, rand(3,4), fkwargs=(dims=2,))
    if VERSION >= v"1.6"
        test_frule(reverse, rand(3,4))
        test_frule(reverse, rand(3,4,5), fkwargs=(dims=(1,3),))
    end

    # Reverse
    test_rrule(reverse, rand(5))
    test_rrule(reverse, rand(5), 2, 4)
    test_rrule(reverse, rand(5), fkwargs=(dims=1,))

    test_rrule(reverse, rand(3,4), fkwargs=(dims=2,))
    if VERSION >= v"1.6"
        test_rrule(reverse, rand(3,4))
        test_rrule(reverse, rand(3,4,5), fkwargs=(dims=(1,3),))

        # Structured
        y, pb = rrule(reverse, Diagonal([1,2,3]))
        @test unthunk(pb(rand(3,3))[2]) isa Diagonal
    end
end

@testset "circshift" begin
    # Forward
    test_frule(circshift, rand(10), 1)
    test_frule(circshift, rand(10), (1,))
    test_frule(circshift, rand(3,4), (-7,2))

    # Reverse
    test_rrule(circshift, rand(10), 1)
    test_rrule(circshift, rand(10) .+ im, -2)
    test_rrule(circshift, rand(10), (1,))
    test_rrule(circshift, rand(3,4), (-7,2))
end

@testset "fill" begin
    test_frule(fill, 12.3, 4)
    test_frule(fill, 5.0, (6, 7))

    test_rrule(fill, 44.4, 4)
    test_rrule(fill, 55 + 0.5im, 5)
    test_rrule(fill, 3.3, (3, 3, 3))
end


text:
include("utils/load_parameters.jl")
include("components/RCPSSPScenario.jl")
include("components/CO2emissions.jl")
include("components/CO2cycle.jl")
include("components/CO2forcing.jl")
include("components/CH4emissions.jl")
include("components/CH4cycle.jl")
include("components/CH4forcing.jl")
include("components/N2Oemissions.jl")
include("components/N2Ocycle.jl")
include("components/N2Oforcing.jl")
include("components/LGemissions.jl")
include("components/LGcycle.jl")
include("components/LGforcing.jl")
include("components/SulphateForcing.jl")
include("components/TotalForcing.jl")
include("components/ClimateTemperature.jl")
include("components/PermafrostSiBCASA.jl")
include("components/PermafrostJULES.jl")
include("components/PermafrostTotal.jl")

scenario = "NDCs"
use_permafrost = true

function climatemodel(scenario::String, use_permafrost::Bool=true, use_seaice::Bool=false)
    m = Model()

    set_dimension!(m, :year, collect(2015:2300)) # NEW - for annual version
    set_dimension!(m, :time, [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300])
    set_dimension!(m, :region, ["EU", "USA", "OECD","USSR","China","SEAsia","Africa","LatAmerica"])

    #add all the components
    scenario = addrcpsspscenario(m, scenario)
    climtemp = addclimatetemperature(m, use_seaice)
    if use_permafrost
        permafrost_sibcasa = add_comp!(m, PermafrostSiBCASA)
        permafrost_jules = add_comp!(m, PermafrostJULES)
        permafrost = add_comp!(m, PermafrostTotal)
    end
    co2emit = add_comp!(m,co2emissions)
    co2cycle = addco2cycle(m, use_permafrost)
    add_comp!(m, co2forcing)
    ch4emit = add_comp!(m, ch4emissions)
    ch4cycle = addch4cycle(m, use_permafrost)
    add_comp!(m, ch4forcing)
    n2oemit = add_comp!(m, n2oemissions)
    add_comp!(m, n2ocycle)
    add_comp!(m, n2oforcing)
    lgemit = add_comp!(m, LGemissions)
    add_comp!(m, LGcycle)
    add_comp!(m, LGforcing)
    sulfemit = add_comp!(m, SulphateForcing)
    totalforcing = add_comp!(m, TotalForcing)

    #connect parameters together
    set_param!(m, :ClimateTemperature, :y_year_ann, collect(2015:2300))
    set_param!(m, :ClimateTemperature, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.])
    set_param!(m, :ClimateTemperature, :y_year_0, 2015.)
    connect_param!(m, :ClimateTemperature => :fant_anthroforcing, :TotalForcing => :fant_anthroforcing)

    if use_permafrost
        permafrost_sibcasa[:rt_g] = climtemp[:rt_g_globaltemperature]
        permafrost_jules[:rt_g] = climtemp[:rt_g_globaltemperature]
        permafrost[:perm_sib_ce_co2] = permafrost_sibcasa[:perm_sib_ce_co2]
        permafrost[:perm_sib_e_co2] = permafrost_sibcasa[:perm_sib_e_co2]
        permafrost[:perm_sib_ce_ch4] = permafrost_sibcasa[:perm_sib_ce_ch4]
        permafrost[:perm_jul_ce_co2] = permafrost_jules[:perm_jul_ce_co2]
        permafrost[:perm_jul_e_co2] = permafrost_jules[:perm_jul_e_co2]
        permafrost[:perm_jul_ce_ch4] = permafrost_jules[:perm_jul_ce_ch4]
    end

    co2emit[:er_CO2emissionsgrowth] = scenario[:er_CO2emissionsgrowth]
    co2cycle[:y_year] = [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.]
    co2cycle[:y_year_0] = 2015.
    co2cycle[:e_globalCO2emissions] = co2emit[:e_globalCO2emissions]
    co2cycle[:rt_g_globaltemperature] = climtemp[:rt_g_globaltemperature]
    if use_permafrost
        co2cycle[:permte_permafrostemissions] = permafrost[:perm_tot_e_co2]
    end

    connect_param!(m, :co2forcing => :c_CO2concentration, :CO2Cycle => :c_CO2concentration)

    ch4emit[:er_CH4emissionsgrowth] = scenario[:er_CH4emissionsgrowth]
    ch4cycle[:y_year] = [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.]
    ch4cycle[:y_year_0] = 2015.
    ch4cycle[:e_globalCH4emissions] = ch4emit[:e_globalCH4emissions]
    ch4cycle[:rtl_g0_baselandtemp] = climtemp[:rtl_g0_baselandtemp]
    ch4cycle[:rtl_g_landtemperature] = climtemp[:rtl_g_landtemperature]
    if use_permafrost
        ch4cycle[:permtce_permafrostemissions] = permafrost[:perm_tot_ce_ch4]
    end

    connect_param!(m, :ch4forcing => :c_CH4concentration, :CH4Cycle => :c_CH4concentration)
    connect_param!(m, :ch4forcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    n2oemit[:er_N2Oemissionsgrowth] = scenario[:er_N2Oemissionsgrowth]
    set_param!(m, :n2ocycle, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.])
    set_param!(m, :n2ocycle, :y_year_0, 2015.)
    connect_param!(m, :n2ocycle => :e_globalN2Oemissions, :n2oemissions => :e_globalN2Oemissions)
    connect_param!(m, :n2ocycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :n2ocycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :n2oforcing => :c_CH4concentration, :CH4Cycle => :c_CH4concentration)
    connect_param!(m, :n2oforcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    lgemit[:er_LGemissionsgrowth] = scenario[:er_LGemissionsgrowth]
    set_param!(m, :LGcycle, :y_year, [2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.,2250.,2300.])
    set_param!(m, :LGcycle, :y_year_0, 2015.)
    connect_param!(m, :LGcycle => :e_globalLGemissions, :LGemissions => :e_globalLGemissions)
    connect_param!(m, :LGcycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :LGcycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :LGforcing => :c_LGconcentration, :LGcycle => :c_LGconcentration)

    sulfemit[:pse_sulphatevsbase] = scenario[:pse_sulphatevsbase]

    connect_param!(m, :TotalForcing => :f_CO2forcing, :co2forcing => :f_CO2forcing)
    connect_param!(m, :TotalForcing => :f_CH4forcing, :ch4forcing => :f_CH4forcing)
    connect_param!(m, :TotalForcing => :f_N2Oforcing, :n2oforcing => :f_N2Oforcing)
    connect_param!(m, :TotalForcing => :f_lineargasforcing, :LGforcing => :f_LGforcing)
    totalforcing[:exf_excessforcing] = scenario[:exf_excessforcing]
    connect_param!(m, :TotalForcing => :fs_sulfateforcing, :SulphateForcing => :fs_sulphateforcing)

    # next: add vector and panel example
    p = load_parameters(m)
    p["y_year_ann"] = Mimi.dim_keys(m.md, :year) # NEW - for annual version
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)
    m
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
function get_output(coefs,Df,sq,r2,f_terms,Gres,Tot,p,n)
    n_term = length(coefs)
    DF = RegCol("Df")
    R2 = RegCol("R²")
    sumsq = RegCol("SumOfSqs")
    F = RegCol("F")
    P = RegCol("P")

    for i in 1:n_term
        setcoef!(R2, coefs[i]=>r2[i])
        setcoef!(DF, coefs[i]=>Df[i])
        setcoef!(sumsq, coefs[i]=>sq[i])
        setcoef!(F, coefs[i]=>f_terms[i])
        setcoef!(P, coefs[i]=>p[i])
    end
    coefArray = hcat(Df,sq,r2,f_terms,p)
    coefArray= vcat(coefArray,[n-1-sum(Df) Gres 1-sum(r2) NaN NaN ])
    coefArray= vcat(coefArray,[n-1 Tot 1 NaN NaN ])

    setcoef!(R2, "Residual"=>1-sum(r2))
    setcoef!(DF, "Residual"=>n-1-sum(Df))
    setcoef!(sumsq, "Residual"=>Gres)
    setcoef!(R2, "Total"=>1)
    setcoef!(DF, "Total"=>n-1)
    setcoef!(sumsq, "Total"=>Tot)
    return NamedArray(coefArray,(vcat(coefs...,"Residual","Total"),["Df","SumOfSqs","R2","F","P"])), hcat(DF,sumsq,R2,F,P)
end


text:
using JSON
using LinearAlgebra: normInf

"""
    δ(x, y)::Float64

Compute the difference between the two values.

If the two values not numerically comparable (_e.g._ strings) then `Inf` or `0`
is returned (depending on whether or not the inputs are equal).
"""
function δ(x, y)::Float64
    if applicable(-, x, y)
        diff = x - y
        if applicable(normInf, diff)
            return normInf(diff)
        end
    end
    return (x == y) ? 0 : Inf64
end

function δ(x::TestCaseResult, y::TestCaseResult)::Float64
    if x.exception != nothing || y.exception != nothing
        error("Cannot compare results with exceptions.")
    end
    return δ(x.result, y.result)
end

# function output_test_results(solution::TestSuiteResult, submission::TestSuiteResult)
#     outfile = get(ENV, "COMPARATIVE_AUTOGRADER_OUTFILE", nothing)
#     if outfile == nothing
#         output_test_results(Base.stdout, solution, submission)
#     else
#         open(outfile, "w") do f
#             output_test_results(f, solution, submission)
#         end
#     end
# end
#
# function output_test_results(
#     out::IO,
#     solution::TestSuiteResult,
#     submission::TestSuiteResult,
# )
#     output = nothing
#     try
#         output = generate_test_results_dict(solution, submission)
#     catch e
#         output = Dict(
#             "passed" => false,
#             "error" => smallrepr(e),
#             "backtrace" => hasmethod(showerror, (typeof(e), )) ? truncatestring(sprint(showerror(e))) : "",
#         )
#     end
#     if isa(out, AbstractString)
#         error("writing to file not implemented") # todo
#     end
#     log(output["passed"] ? "Submission passed." : "Submission failed.")
#     JSON.print(out, output)
# end
#
# function output_error_dict(error::String, backtrace::Union{Nothing,String}="")
#     output_error_dict(Base.stdout, error, backtrace)
# end
#
# function output_error_dict(out::IO, error::String, backtrace::Union{Nothing,String}="")
#     JSON.print(out, generate_error_dict(error, backtrace))
# end
#
# generate_error_dict(error::String, backtrace::Union{Nothing,String}="") = Dict(
#     "passed" => false,
#     "error" => error,
#     "backtrace" => backtrace == nothing ? "" : backtrace,
# )
#
# function generate_test_results_dict(
#     solution::TestSuiteResult,
#     submission::TestSuiteResult,
# )::Dict
#     testcase_array = []
#     if solution.exception != nothing
#         log("Solution threw an exception: $(solution.exception)\n$(chomp(solution.backtrace))")
#         return generate_error_dict("Solution threw an exception: $(solution.exception)", solution.backtrace)
#     end
#     if submission.exception != nothing
#         log("Submission threw an exception: $(submission.exception)\n$(chomp(submission.backtrace))")
#         return generate_error_dict("Submission threw an exception: $(submission.exception)", submission.backtrace)
#     end
#     if length(solution.results) != length(submission.results)
#         error("Test case dimensions didn't match (solution returned $(length(solution.results)) test cases, submission returned $(length(submission.results))).")
#     end
#     for (sol, sub) in zip(solution.results, submission.results)
#         if sol.exception != nothing
#             log("Solution threw an exception on a test case: $(sol.exception)\n$(chomp(sol.backtrace))")
#             error("Solution threw an exception on a test case.")
#         end
#         if sol.result == nothing
#             log("Solution returned nothing.")
#             error("Solution returned nothing.")
#         end
#         testcase_dict = Dict(
#             "timeElapsed" => sub.time,
#             "stdout" => truncatestring(sub.stdout),
#             "stderr" => truncatestring(sub.stderr),
#         )
#         if sub.exception != nothing
#             testcase_dict["passed"] = false
#             testcase_dict["exception"] = sub.exception
#             testcase_dict["backtrace"] = sub.backtrace
#         else
#             if sub.result == nothing
#                 log("ERROR(?): result is $(sub.result) and exception is $(sub.exception).")
#             end
#             error = calculate_error(sol.result, sub.result)
#             testcase_dict["delta"] = error
#             if error < 1e-5 # todo: no hardcoded constants
#                 testcase_dict["passed"] = true
#             else
#                 testcase_dict["passed"] = false
#             end
#             testcase_dict["result"] = smallrepr(sub.result)
#         end
#         push!(testcase_array, testcase_dict)
#     end
#     final_output = Dict(
#         "passed" => all(testcase["passed"] for testcase in testcase_array),
#         "testcases" => testcase_array,
#     )
#     return final_output
#     JSON.print(out, final_output)
# end
#
# function calculate_error(a::Number, b::Number)
#     return abs(a - b)
# end
#
# # Calculate the error between two arrays.
# function calculate_error(a::Union{AbstractArray,Tuple}, b::Union{AbstractArray,Tuple})
#     return maximum(calculate_error.(a, b))
# end


text:
# only package loading for now
using ClusterManagers


text:
##
# Represent a binary hierarchical Domain, Space, and Fun
##

mutable struct HierarchicalDomain{S,T,HS} <: ApproxFun.UnivariateDomain{T}
    data::HS
    HierarchicalDomain{S,T,HS}(data::HS) where {S,T,HS} = new{S,T,HS}(data)
end

mutable struct HierarchicalSpace{S,T,HS} <: ApproxFun.Space{ApproxFun.AnyDomain,T}
    data::HS
    HierarchicalSpace{S,T,HS}(data::HS) where {S,T,HS} = new{S,T,HS}(data)
end

mutable struct HierarchicalFun{S,T,HS}
    data::HS
    HierarchicalFun{S,T,HS}(data::HS) where {S,T,HS} = new{S,T,HS}(data)
end

for (HDSF,etyp) in ((:HierarchicalDomain,:eltype),(:HierarchicalSpace,:rangetype),(:HierarchicalFun,:eltype))
    @eval begin
        export $HDSF

        $HDSF(data::Tuple{S1,S2}) where {S1,S2} = $HDSF{promote_type(S1,S2),promote_type($etyp(S1),$etyp(S2)),Tuple{S1,S2}}(data)
        $HDSF(data::Tuple{$HDSF{S1,T1,HS1},$HDSF{S2,T2,HS2}}) where {S1,S2,T1,T2,HS1,HS2} = $HDSF{promote_type(S1,S2),promote_type($etyp(S1),$etyp(S2),T1,T2),Tuple{$HDSF{S1,T1,HS1},$HDSF{S2,T2,HS2}}}(data)
        $HDSF(data::Tuple{S1,$HDSF{S,T,HS}}) where {S,S1,T,HS} = $HDSF{promote_type(S,S1),promote_type($etyp(S),$etyp(S1),T),Tuple{S1,$HDSF{S,T,HS}}}(data)
        $HDSF(data::Tuple{$HDSF{S,T,HS},S1}) where {S,S1,T,HS} = $HDSF{promote_type(S,S1),promote_type($etyp(S),$etyp(S1),T),Tuple{$HDSF{S,T,HS},S1}}(data)
        collectdata(H::$HDSF{S,T,NTuple{2,S}}) where {S,T} = collect(data(H))
        collectdata(H::$HDSF{S,T,Tuple{S,$HDSF{S,T,HS}}}) where {S,T,HS} = vcat(H.data[1],collectdata(H.data[2]))
        collectdata(H::$HDSF{S,T,Tuple{$HDSF{S,T,HS},S}}) where {S,T,HS} = vcat(collectdata(H.data[1]),H.data[2])

        $HDSF(data::Vector) = $HDSF(data,round(Int,log2(length(data))))

        function $HDSF(data::Vector,n::Int)
            @assert length(data) == 2^n
            if n == 1
                return $HDSF(tuple(data...))
            elseif n ≥ 2
                return $HDSF(($HDSF(data[1:2^(n-1)],n-1),$HDSF(data[1+2^(n-1):end],n-1)))
            end
        end

        Base.similar(H::$HDSF) = $HDSF(map(similar,data(H)))
        Base.similar(H::$HDSF{SS,V,T}, S) where {SS,V,T} = $HDSF(map(A->similar(A,S),data(H)))

        data(H::$HDSF) = H.data

        nlevels(H::$HDSF) = 1+mapreduce(nlevels,max,data(H))

        partition(H::$HDSF) = data(H)

        function partition(VH::Vector{H}) where H<:$HDSF
            n = length(VH)
            H11,H12 = partition(VH[1])
            H1,H2 = fill(H11,n),fill(H12,n)
            for i=1:n
                H1[i],H2[i] = partition(VH[i])
            end
            H1,H2
        end

        function collectdata(H::$HDSF{S}) where S
            ret = S[]
            append!(ret,mapreduce(collectdata,vcat,data(H)))
            ret
        end

        convert(::Type{$HDSF{S,T,HS}},M::$HDSF{S,T,HS}) where {S,T,HS} =
            M
        convert(::Type{$HDSF{S,T,HS}},M::$HDSF) where {S,T,HS} =
            $HDSF{S,T,HS}(convert(Vector{S},collectdata(M)))
        Base.promote_rule(::Type{$HDSF{S,T,HS}},::Type{$HDSF{SS,TT,HSS}}) where {S,T,HS,SS,TT,HSS} =
            $HDSF{promote_type(S,SS),promote_type(T,TT),promote_type(HS,HSS)}
    end
end

domain(H::HierarchicalDomain) = H
domain(H::HierarchicalSpace) = HierarchicalDomain(map(domain,data(H)))
space(H::HierarchicalFun) = HierarchicalSpace(map(space,data(H)))

ApproxFun.domaindimension(H::HierarchicalSpace) = 1
Space(H::HierarchicalDomain) = HierarchicalSpace(map(Space,data(H)))

PiecewiseSpace(H::HierarchicalSpace) = PiecewiseSpace(map(PiecewiseSpace,data(H))...)

Fun(f::Function,H::HierarchicalDomain) = HierarchicalFun((Fun(f,H.data[1]),Fun(f,H.data[2])))
Fun(f::Function,H::HierarchicalSpace) = HierarchicalFun((Fun(f,H.data[1]),Fun(f,H.data[2])))

# algebra

for op in (:+,:-)
    @eval begin
        $op(H::HierarchicalFun) = HierarchicalFun(map($op,data(H)))
        $op(H::HierarchicalFun,a::Number) = HierarchicalFun(($op(H.data[1],a),$op(H.data[2],a)))
        $op(a::Number,H::HierarchicalFun) = HierarchicalFun(($op(a,H.data[1]),$op(a,H.data[2])))

        $op(H::HierarchicalFun,J::HierarchicalFun) = HierarchicalFun(map($op,data(H),data(J)))
        $op(H::HierarchicalFun,J::Fun) = $op(J,full(H))
    end
end

*(H::HierarchicalFun,a::Number) = HierarchicalFun((H.data[1]*a,H.data[2]*a))
*(a::Number,H::HierarchicalFun) = H*a

for op in (:linebilinearform,:bilinearform)
    @eval begin
        $op(H::HierarchicalFun,J::HierarchicalFun) = $op(H.data[1],J.data[1])+$op(H.data[2],J.data[2])
        $op(H::HierarchicalFun,J::Fun{S}) where {S<:PiecewiseSpace} =
            sum(map($op,collectdata(H),components(J)))
        $op(J::Fun{S},H::HierarchicalFun) where {S<:PiecewiseSpace} =
            sum(map($op,components(J),collectdata(H)))
    end
end

Base.cumsum(H::HierarchicalFun) = HierarchicalFun((cumsum(H.data[1]),sum(H.data[1])+cumsum(H.data[2])))
Base.conj!(H::HierarchicalFun) = (map(conj!,data(H));H)
Base.copy!(H::HierarchicalFun,J::HierarchicalFun) = (map(copy!,data(H),data(J));H)

for op in (:(Base.zero),:(Base.ones),:(Base.abs),:(Base.abs2),:(Base.conj),:(Base.copy))
    @eval begin
        $op(H::HierarchicalFun) = HierarchicalFun(map($op,data(H)))
    end
end

# Formatted operations

for (op,opformatted) in ((:+,:⊕),(:-,:⊖))
    @eval begin
        function $opformatted(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P}) where {S,T,HS1<:HierarchicalFun,HS2<:HierarchicalFun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($opformatted(H1,Fun(components(J)[1:n1],PiecewiseSpace)),
                                $opformatted(H2,Fun(components(J)[1+n1:n1+n2],PiecewiseSpace))))
        end
        function $opformatted(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P}) where {S,T,HS1<:Fun,HS2<:HierarchicalFun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = 1,length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($op(H1,Fun(components(J)[1:n1],PiecewiseSpace)),
                                    $opformatted(H2,Fun(components(J)[1+n1:n1+n2],PiecewiseSpace))))
        end
        function $opformatted(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P}) where {S,T,HS1<:HierarchicalFun,HS2<:Fun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),1
            return HierarchicalFun(($opformatted(H1,Fun(components(J)[1:n1],PiecewiseSpace)),
                                    $op(H2,Fun(components(J)[1+n1:n1+n2],PiecewiseSpace))))
        end
        function $opformatted(H::HierarchicalFun{S,T,Tuple{HS1,HS2}},J::Fun{P}) where {S,T,HS1<:Fun,HS2<:Fun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = 1,1
            return HierarchicalFun(($op(H1,Fun(components(J)[1:n1],PiecewiseSpace)),
                                    $op(H2,Fun(components(J)[1+n1:n1+n2],PiecewiseSpace))))
        end

        function $opformatted(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}}) where {S,T,HS1<:HierarchicalFun,HS2<:HierarchicalFun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($opformatted(Fun(components(J)[1:n1],PiecewiseSpace),H1),
                                    $opformatted(Fun(components(J)[1+n1:n1+n2],PiecewiseSpace),H2)))
        end
        function $opformatted(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}}) where {S,T,HS1<:Fun,HS2<:HierarchicalFun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = 1,length(PiecewiseSpace(space(H2)))
            return HierarchicalFun(($op(Fun(components(J)[1:n1],PiecewiseSpace),H1),
                                    $opformatted(Fun(components(J)[1+n1:n1+n2],PiecewiseSpace),H2)))
        end
        function $opformatted(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}}) where {S,T,HS1<:HierarchicalFun,HS2<:Fun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = length(PiecewiseSpace(space(H1))),1
            return HierarchicalFun(($opformatted(Fun(components(J)[1:n1],PiecewiseSpace),H1),
                                    $op(Fun(components(J)[1+n1:n1+n2],PiecewiseSpace),H2)))
        end
        function $opformatted(J::Fun{P},H::HierarchicalFun{S,T,Tuple{HS1,HS2}}) where {S,T,HS1<:Fun,HS2<:Fun,P<:PiecewiseSpace}
            H1,H2 = partition(H)
            n1,n2 = 1,1
            return HierarchicalFun(($op(Fun(components(J)[1:n1],PiecewiseSpace),H1),
                                    $op(Fun(components(J)[1+n1:n1+n2],PiecewiseSpace),H2)))
        end
    end
end


text:
# set power limit to max value
# SPM
export psp_setplimitmax

const PLIMIT_MAX = [0x53, 0x50, 0x4D, 0x0D]

"""
    psp_setplimitmax(io)

Set power limit to maximum value.
"""
function psp_setplimitmax(io_psp::IO)
  write(io_psp, PLIMIT_MAX)
  return nothing
end


text:
"""
    KissSmoothing: Easily smooth your data!

exports a single function `denoise`, for further help look at it's docstring.
"""
module KissSmoothing
using FFTW: dct, idct
using Statistics: mean

"""
    denoise(V::Array; factor=1.0, rtol=1e-12, dims=ndims(V), verbose = false)

smooth data in `V` according to:

    `V` : data to smooth

    `factor` : smoothing intensity

    `rtol` : relative tolerance on how precise the smoothing intensity is determined

    `dims` : array dimension being smoothed along

    `verbose` : enables some printing of internal info

returns a tuple (S,N) where:

    `S` : is the smoothed data

    `N` : is the extracted noise

in particular `S + N` reconstructs the original data `V`.
"""
function denoise(
    V::AbstractArray{Float64,N};
    factor::Float64 = 1.0,
    rtol::Float64 = 1e-12,
    dims::Int64 = N,
    verbose::Bool = false,
) where {N}
    buf = IOBuffer()
    K1 = sqrt(2 / pi)
    K2 = sqrt(2) * K1
    #K3 = sqrt(6)*K1
    #K4 = sqrt(20)*K1
    lV = size(V, dims)
    if lV < 4
        return copy(V), zero(V)
    end
    iV = dct(V)
    stri = map(i -> ifelse(i == dims, lV, 1), 1:ndims(V))
    X = map(abs2, reshape(LinRange(0.0, 1.0, lV), stri...))
    d = factor * mean(abs, diff(V, dims = dims)) * (K1 / K2)
    σt = 0.5
    σd = 0.25
    f = zero(V)
    for iter = 1:60
        σ = sqrt(lV) * σt / (1 - σt)
        if !isfinite(σ)
            break
        end
        f .= idct(iV .* exp.(-X .* σ))
        c = mapreduce((x, y) -> abs(x - y), +, f, V) / length(V)
        Δ = d - c
        σt += σd * sign(Δ)
        σd /= 2
        if verbose
            println(buf, iter, "  ", σ, "  ", Δ)
        end
        if abs(d - c) < rtol * d
            break
        end
    end
    if verbose
        print(String(take!(buf)))
    end
    f, V .- f
end

function tps(x::AbstractArray{Float64},y::AbstractArray{Float64})
    ϵ = nextfloat(0.0)
    r = ϵ+mapreduce(+, x, y) do a, b
        abs2(a - b)
    end
    r * log(r)
end

struct RBF{G<:AbstractArray{Float64},C<:AbstractArray{Float64}}
    Γ::G
    C::C
end

function evalPhi(xs::AbstractArray{Float64}, cp::AbstractArray{Float64})
    Phi = zeros(size(xs, 1), size(cp, 1)+1)
    for i = 1:size(xs, 1)
        mu = 0.0
        for j = 1:size(cp, 1)
            k = tps(xs[i, :], cp[j, :])
            Phi[i, j] = k
            mu += k
        end
        mu /= size(cp, 1)
        for j = 1:size(cp, 1)
            Phi[i, j] -= mu
        end
        Phi[i, size(cp, 1)+1] = 1.0
    end
    Phi
end

function (net::RBF)(X::AbstractArray{Float64})
    evalPhi(X, net.C) * net.Γ
end

"""
    fit_rbf(xv::Array, yv::Array, cp::Array)

fit thin-plate radial basis function according to:

    `xv` : array NxP, N number of training points, P number of input variables

    `yv` : array NxQ, N number of training points, Q number of output variables

    `cp` : array KxP, K number of control points, P number of input variables

returns a callable RBF object.
"""
function fit_rbf(
    xv::AbstractArray{Float64},
    yv::AbstractArray{Float64},
    cp::AbstractArray{Float64},
)
    RBF(evalPhi(xv, cp) \ yv, collect(cp))
end

export denoise, fit_rbf, RBF
end # module


text:
"""
    AlgoBaseModule

Module for base  algorithms.
"""
module AlgoBaseModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FEMMBaseModule: integratefieldfunction, transferfield!
import LinearAlgebra: norm, dot

function _keymatch(key::String, allowed_keys::Array{String})
    matched_key = nothing
    for  j = 1:length(allowed_keys)
        m = match(Regex("^$key"), allowed_keys[j])
        if (m != nothing)
            matched_key = allowed_keys[j]
            break
        end
    end
    return matched_key
end

function dcheck!(d::FDataDict, recognized_keys::Array{String})
    notmatched = fill("", 0)
    for  key in keys(d)
        matched_key = _keymatch(key, recognized_keys)
        if matched_key == nothing
            push!(notmatched, "Key \"$key\" not matched")
        else
            if key !=  matched_key
                push!(notmatched, "Key \"$key\" not fully matched (partial match \"$matched_key\")")
            end
        end
    end
    return notmatched
end

"""
    richextrapol(solns::AbstractArray{Real}, params::AbstractArray{Real})

Richardson extrapolation.

`solns` =  array of solution values
`params` = array of values of parameters for which the solns have been obtained

This function is applicable only to fixed ratio between the mesh sizes,
  `params[1]/params[2) = params[2)/params[3)`.
Output:
solnestim= estimate of the asymptotic solution from the data points in the solns array
beta= convergence rate
c = constant in the estimate "error=c*h^beta"
residual = residual after equations from which the above quantities were
     solved (this is a measure of how accurately was the system solved).
"""
function richextrapol(solns::FFltVec, params::FFltVec)
    @assert abs(params[1]/params[2] - params[2]/params[3]) < 1e-6 "Parameter pair ratio not fixed"
    nsolns = solns./solns[1];
    solnestim = ((-(nsolns[1]*nsolns[3]-nsolns[2]^2)/(2*nsolns[2]-nsolns[1]-nsolns[3])))*solns[1];
    if (solnestim-solns[1]) <= 0
        beta = log((solnestim-solns[2])/(solnestim-solns[3]))/log(params[2]/params[3]);
    else
        beta = log((solnestim-solns[1])/(solnestim-solns[3]))/log(params[1]/params[3]);
    end
    # just to check things, calculate the residual
    c = (solnestim-solns[1])/params[1]^beta;
    residual = zeros(FFlt,3)
    for I =1:3
        residual[I] = (solnestim-solns[I])-c*params[I]^beta; # this should be close to zero
    end
    return solnestim, beta, c, residual
end

"""
    bisect(fun, xl, xu, tolx, tolf)

Implementation of the bisection method.

Tolerance both on `x` and on `f(x)` is used.
`fun` = function,
`xl`= lower value of the bracket,
`xu`= upper Value of the bracket,
`tolx`= tolerance on the location of the root,
`tolf`= tolerance on the function value
"""
function bisect(fun, xl, xu, tolx, tolf)
    iter = 0;
    if (xl > xu)
        temp = xl; xl = xu; xu = temp;
    end
    fl = fun(xl);
    fu = fun(xu);
    @assert fl*fu < 0.0 "Need to get a bracket"
    if fl == 0.0
        xu = xl;
        return xl, xu;
    end
    if fu == 0.0
        xl = xu;
        return xl, xu;
    end
    while true
        xr = (xu + xl) / 2.0; # bisect interval
        fr = fun(xr); # value at the midpoint
        if (fr*fl < 0.0)
            xu = xr; fu = fr;# upper --> midpoint
        elseif (fr == 0.0)
            xl = xr; xu = xr;# exactly at the root
        else
            xl = xr; fl = fr;# lower --> midpoint
        end
        if (abs(xu-xl) < tolx) || (abs(fr) < tolf)
            break; # We are done
        end
        iter = iter+1;
    end
    return xl, xu;
end

"""
    fieldnorm(modeldata)

Compute norm of the target field.

`modeldata` = data dictionary, mandatory keys:
    - "fens" = finite element node set
    - "regions" = array of regions
    - "targetfields" = array of fields, one for each region
    - "geom" = geometry field
    - "elementsize" = representative element size,
"""
function fieldnorm(modeldata)
    fens = modeldata["fens"]
    regions = modeldata["regions"]
    targetfields = modeldata["targetfields"]
    geom = modeldata["geom"]

    @assert length(regions) == length(targetfields)

    fnorm = 0.0 # Initialize the norm of the difference
    for i = 1:length(regions)
        # Compute the addition to the norm of the field
        fnorm += integratefieldfunction(regions[i]["femm"], geom, targetfields[i], (x, v) -> norm(v)^2, 0.0)
    end

    return sqrt(fnorm)
end

"""
    fielddiffnorm(modeldatacoarse, modeldatafine)

Compute norm of the difference of the target fields.

For both the "coarse"- and "fine"-mesh `modeldata` the data dictionaries need to contain the mandatory keys:
- "fens" = finite element node set
- "regions" = array of regions
- "targetfields" = array of fields, one for each region
- "geom" = geometry field
- "elementsize" = representative element size,
- "geometricaltolerance" = geometrical tolerance (used in field transfer; refer to the documentation of `transferfield!`)
- "parametrictolerance" = parametric tolerance (used in field transfer; refer to the documentation of `transferfield!`)
"""
function fielddiffnorm(modeldatacoarse, modeldatafine)
    # Load coarse-mesh data
    fenscoarse = modeldatacoarse["fens"]
    regionscoarse = modeldatacoarse["regions"]
    targetfieldscoarse = modeldatacoarse["targetfields"]

    # Load fine-mesh data
    fensfine = modeldatafine["fens"]
    regionsfine = modeldatafine["regions"]
    targetfieldsfine = modeldatafine["targetfields"]

    # Tolerances
    geometricaltolerance = modeldatacoarse["geometricaltolerance"]
    parametrictolerance = get(modeldatacoarse, "parametrictolerance", 0.01);

    geom = modeldatafine["geom"]

    @assert length(regionscoarse) == length(regionsfine)
    @assert length(regionscoarse) == length(targetfieldscoarse)
    @assert length(regionsfine) == length(targetfieldsfine)

    diffnorm = 0.0 # Initialize the norm of the difference
    for i = 1:length(regionscoarse)
        # Transfer the result from the coarse mesh to the fine-mesh
        ffine = targetfieldsfine[i]
        fcoarse = targetfieldscoarse[i]
        fcoarsetransferred = deepcopy(ffine)
        fcoarsetransferred = transferfield!(fcoarsetransferred,
            fensfine, regionsfine[i]["femm"].integdata.fes, fcoarse,
            fenscoarse, regionscoarse[i]["femm"].integdata.fes, geometricaltolerance; parametrictolerance = parametrictolerance)
        # Form the difference  field
        diffff = deepcopy(fcoarsetransferred)
        diffff.values[:] = ffine.values - fcoarsetransferred.values
        # Compute the addition to the norm of the difference
        diffnorm += integratefieldfunction(regionsfine[i]["femm"], geom, diffff, (x, v) -> norm(v)^2, 0.0)
    end

    return sqrt(diffnorm)
end

"""
    evalconvergencestudy(modeldatasequence, File)

Evaluate a convergence study from a model-data sequence.

`modeldatasequence` = array of `modeldata` dictionaries.

Refer to methods `fieldnorm` and `fielddiffnorm` for details on the required keys in the dictionaries.
"""
function evalconvergencestudy(modeldatasequence)
    # Find the element sizes
    elementsizes = [md["elementsize"] for md in modeldatasequence]

    # The  finest solution will provide the norm used in the normalization of the errors
    finestsolnorm = fieldnorm(modeldatasequence[end])
    # Compute the approximate errors  as the differences of successive solutions
    diffnorms = FFlt[]
    for i = 1:(length(modeldatasequence) - 1)
        push!(diffnorms, fielddiffnorm(modeldatasequence[i], modeldatasequence[i + 1]))
    end
    # Normalize the errors
    errornorms = diffnorms ./ finestsolnorm
    # Compute the convergence rate
    f = log.(vec(errornorms))
    A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
    p = A \ f
    convergencerate = p[1]

    return elementsizes, errornorms, convergencerate
end

"""
    conjugategradient(A, b, x0, maxiter)

Compute one or more iterations of the conjugate gradient process.
"""
function conjugategradient(A::T, b::FFltVec, x0::FFltVec, maxiter::FInt) where {T}
    x = deepcopy(x0);
    gt = deepcopy(x0);
    d = deepcopy(x0);
    gt .= A*x - b;# transpose of gradient: g = x'*A-b';
    @. d = gt
    for iter = 1:maxiter
        Ad = A*d;
        rho = dot(d, Ad);
        alpha = -dot(gt, d)/rho; # alpha =(-g*d)/(d'*A*d);
        @. x += alpha * d;
        gt .= A*x - b;# negative transpose of gradient
        beta = dot(gt, Ad)/rho; # beta =(g*A*d)/(d'*A*d);
        @. d = beta*d - gt;
    end
    return x
end

"""
    qtrap(ps, xs)

Compute the area under the curve given by a set of parameters along
an interval and the values of the 'function' at those parameter values.
The parameter values need not be uniformly distributed.

Trapezoidal rule is used to evaluate the integral. The 'function' is
assumed to vary linearly inbetween the given points.
"""
function qtrap(ps, xs)
    @assert length(ps) == length(xs)
    num = 0.0
    for i = 1:length(ps)-1
        num += (ps[i+1]-ps[i]) * (xs[i+1] + xs[i])/2.0
    end
    return num
end

"""
    qcovariance(ps, xs, ys; ws = nothing)

Compute the covariance for two functions given by the arrays `xs` and `ys`
at the values of the parameter `ps`. `ws` is the optional weights vector;
if it is not supplied, uniformly distributed default weights are assumed.

Notes:
– The mean is subtracted from both functions.
– This function is not particularly efficient: it computes the mean of
both functions and it allocates arrays instead of overwriting the
contents of the arguments.
"""
function qcovariance(ps, xs, ys; ws = nothing)
    @assert length(ps) == length(xs) == length(ys)
    if (ws == nothing)
        ws = ones(FFlt, length(ps))
    end
    @assert length(ws) == length(xs)
    ws[:] = ws ./ qtrap(ps, ws) # Make sure the integral of the weight is equal to 1.0
    xmean = qtrap(ps, xs .* ws)
    ymean = qtrap(ps, ys .* ws)
    zxs = deepcopy(xs) .- xmean
    zys = deepcopy(ys) .- ymean
    return qtrap(ps, zxs .* zys .* ws)
end

"""
    qvariance(ps, xs; ws = nothing)

Compute the variance of a function given by the array `xs` at
the values of the parameter `ps`. `ws` is the optional weights vector
with unit default weights.
"""
qvariance(ps, xs; ws = nothing) = qcovariance(ps, xs, xs; ws = ws)

end


text:
include("Geo3x3.jl")

code = Geo3x3.encode(35.65858, 139.745433, 14)
println(code)

pos = Geo3x3.decode("E9139659937288")
println(pos)


text:
push!(LOAD_PATH, "./src")

# "Packages"
using Revise # reduce need for recompile
using Plots
using Documenter
using LinearAlgebra
using SparseArrays

# "User defined modules"
using CommonUtils
using Basis1D
using Basis3DHex
using UniformHexMesh

using SetupDG
using UnPack

N = 3
K1D = 8
CFL = .75
T = .25 # endtime

VX,VY,VZ,EToV = uniform_hex_mesh(K1D,K1D,K1D)

rd = init_reference_hex(N,gauss_quad(0,0,N))
md = init_mesh((VX,VY,VZ),EToV,rd)

# make domain periodic
@unpack Nfaces = rd
@unpack xf,yf,zf,K,mapM,mapP,mapB = md
LX = 2; LY = 2; LZ = 2
mapPB = build_periodic_boundary_maps(xf,yf,zf,LX,LY,LZ,
                                    Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

"Time integration"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N+1)*(N+2)*3/2  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

function rhs(u,rd::RefElemData,md::MeshData)

    @unpack Dr,Ds,Dt,LIFT,Vf = rd
    @unpack rxJ,sxJ,txJ,ryJ,syJ,tyJ,rzJ,szJ,tzJ,J = md
    @unpack nxJ,nyJ,nzJ,sJ = md
    @unpack mapP,mapB = md

    uf = Vf*u
    du = uf[mapP]-uf
    uflux = @. .5*(du*nxJ - 0*abs(nxJ)*du)

    ur,us,ut = (A->A*u).((Dr,Ds,Dt))
    uxJ      = @. rxJ*ur + sxJ*us + txJ*ut
    rhsu     = uxJ + LIFT*uflux

    return @. -rhsu/J
end

# set initial conditions
@unpack x,y,z = md
u = @. exp(-25*(x^2+y^2))

resQ = zeros(size(x))
for i = 1:Nsteps
    for INTRK = 1:5
        rhsQ    = rhs(u,rd,md)
        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. u   += rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Number of time steps: $i out of $Nsteps")
    end
end

# "plotting nodes"
gr(aspect_ratio=1,legend=false,
   markerstrokewidth=0,markersize=2)

@unpack Vp = rd
xp,yp,zp = (x->Vp*x).((x,y,z))
vv = Vp*u
ids = map(x->x[1],findall(@. abs(zp[:]-1)<1e-10))
(xp,yp,zp,vv) = (x->x[ids]).((xp,yp,zp,vv))
scatter(xp,yp,vv,zcolor=vv,camera=(0,90))

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
using Yields
using Test

include("generic.jl")
include("Rate.jl")
include("bootstrap.jl")
include("RateCombination.jl")
include("SmithWilson.jl")

include("ActuaryUtilities.jl")
include("misc.jl")
#TODO EconomicScenarioGenerators.jl integration tests


text:
using Test
using MOSLab
SemiUstrip()
S = Silicon()

F = BoltzmanDist()
T = 300.0*Tunit

niB(T,S)
n(T,F,S)
@test isapprox(Eg(T,S.Eg),1.05; atol=0.1)
@test isapprox(log10(ni(T,F,S)),10.0; atol=0.5)
@test n(T,F,S) ≈ ni(T,F,S)
@test p(T,F,S) ≈ ni(T,F,S)
@test Nv(T,S) == 1.04e19
@test Nc(T,S) == 2.8e19
@test p(T,F,S) ≈ n(T,F,S)
@test Ef(T,F,S) ≈ Ei(T,F,S)
@test isapprox(Ei(T,F,S), 0.5*(Ec(T,S)+Ev(T,S)) + kb*T*0.5*(Nv(T,S)/Nc(T,S));atol=0.1)

PD = PDoppingC(1e15)
@test Na(PD) == 1e15
@test Nd(PD) == 0
ND = NDoppingC(1e15)
@test Nd(ND) == 1e15
@test Na(ND) == 0

PS = PSiliconNU(1e15)
@test Na(PS.Dopping) == 1e15
@test P(T,F,PS) ≈ N(T,F,PS)
@test p(T,F,PS) ≈ Na(PS.Dopping)
@test Ef(T,F,PS) < Ei(T,F,PS)
P(T,F,PS)
n(T,F,PS)
NS = NSiliconNU(1e15)
@test Nd(NS.Dopping) == 1e15
@test P(T,F,NS)≈N(T,F,NS)
@test n(T,F,NS) ≈ Nd(NS.Dopping)
@test Ef(T,F,NS) > Ei(T,F,NS)

SDI = SemiconductorData(T,F,Silicon();psi=0.0)
@test isapprox(Eg(SDI),1.05; atol=0.1)
@test isapprox(log10(ni(SDI)),10.0; atol=0.5)
@test n(SDI) ≈ ni(SDI)
@test p(SDI) ≈ ni(SDI)
@test p(SDI) ≈ n(SDI)
@test Ef(SDI) ≈ Ei(SDI)
@test isapprox(Ei(SDI), 0.5*(Ec(SDI)+Ev(SDI)) + kb*T*0.5*(Nv(T,S)/Nc(T,S));atol=0.1)


text:
function sfMeas(sf::Matrix{T}, idx::Vector{Int64}) where T<:Complex

  meas = zeros(ComplexF64, length(idx), size(sf,2))
  for k=1:size(sf,2)
    meas[:,k] .= sf[:,k][idx]
  end

  return meas
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries, IntervalArithmetic

using Test
eeuler = Base.MathConstants.e

@testset "Tests Taylor1 and TaylorN expansions over Intervals" begin
    a = 1..2
    b = -1 .. 1
    p4(x, a) = x^4 + 4*a*x^3 + 6*a^2*x^2 + 4*a^3*x + a^4
    p5(x, a) = x^5 + 5*a*x^4 + 10*a^2*x^3 + 10*a^3*x^2 + 5*a^4*x + a^5

    ti = Taylor1(Interval{Float64}, 10)
    x, y = set_variables(Interval{Float64}, "x y")

    @test eltype(ti) == Interval{Float64}
    @test eltype(x) == Interval{Float64}

    @test p4(ti,-a) == (ti-a)^4
    @test p5(ti,-a) == (ti-a)^5
    @test p4(ti,-b) == (ti-b)^4
    @test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs)

    @test p4(x,-y) == (x-y)^4
    @test p5(x,-y) == (x-y)^5
    @test p4(x,-a) == (x-a)^4
    @test p5(x,-a) == (x-a)^5
    @test p4(x,-b) == (x-b)^4
    for ind in eachindex(p5(x,-b))
        @test all((p5(x,-b)[ind]).coeffs .⊆ (((x-b)^5)[ind]).coeffs)
    end
    @test evaluate(p4(x,-b), IntervalBox(a,b)) == p4(a, b)
    @test (p5(x,-b))(IntervalBox(a,b)) == p5(a, b)

    # Tests `evaluate`
    @test (a-b)^4 ⊆ ((x-y)^4)(a × b)
    @test (((x-y)^4)[4])(a × b) == -39 .. 81
end


text:
using JuliaChem
using JSON
using Plots
include("full_repl.jl")

function start_calculations()
    JuliaChem.initialize()
end

"""
run a hartree fock calculation using JuliaChem
JuliaChem must be initialized before this method is called

Parameters: hydrogen_radius :: AbstractFloat  - The radius between the two hydrogens in Atomic Units
for the hartree fock run.

Returns
energy of the system in hartrees (Atomic Units)
"""
function run_h2_hartree_fock(hydrogen_radius:: AbstractFloat) :: AbstractFloat
    file_path = create_input_file(hydrogen_radius/2)
    rhf_result = full_rhf(file_path)
    return rhf_result[1]["Energy"]
end

function create_input_file(hydrogen_radius) :: String
    input_file_JSON = JSON.parsefile(joinpath(@__DIR__, "hydrogen_template.json"))
    input_file_JSON["molecule"]["geometry"][6] = hydrogen_radius
    file_path = joinpath(@__DIR__, "input_files/H2-radius$(hydrogen_radius).json")
    io = open(file_path, "w+")
    JSON.print(io, input_file_JSON)
    close(io)
    return file_path
end

function end_calculations()
    JuliaChem.finalize()
end


text:
module ZeroIndex
# I really (REALLY) want zero indexd arrays for some DP problems, here's
# a quick hack to get them

export vget
function vget{T<:Any, N<:Integer}(a::Array{T, 1}, i::N)::T
    return a[i + 1]
end

export vset!
function vset!{T<:Any, N<:Integer}(a::Array{T, 1}, i::N, val::T)
    a[i + 1] = val
end

export aget
function aget{T<:Any, N<:Integer}(a::Array{T, 2}, i::N, j::N)::T
    return a[i + 1, j + 1]
end

export aset!
function aset!{T<:Any, N<:Integer}(a::Array{T, 2}, i::N,  j::N, val::T)
    a[i + 1, j + 1] = val
end

end # module ZeroIndex