text:
GBScalar{T}() where {T} = GBScalar{T}(libgb.GxB_Scalar_new(toGrB_Type(T)))

Base.unsafe_convert(::Type{libgb.GxB_Scalar}, s::GBScalar) = s.p

function GBScalar(v::T) where {T <: valid_union}
    x = GBScalar{T}()
    x[] = v
    return x
end

function Base.deepcopy(s::GBScalar{T}) where {T}
    return GBScalar{T}(libgb.GxB_Scalar_dup(s))
end

clear(s::GBScalar) = libgb.GxB_Scalar_clear(s)

for T ∈ valid_vec
    func = Symbol(:GxB_Scalar_setElement_, suffix(T))
    @eval begin
        function Base.setindex!(value::GBScalar{$T}, s::$T)
            libgb.$func(value, s)
        end
    end
    func = Symbol(:GxB_Scalar_extractElement_, suffix(T))
    @eval begin
        function Base.getindex(value::GBScalar{$T})
            libgb.$func(value)
        end
    end
end
Base.eltype(::Type{GBScalar{T}}) where{T} = T

function Base.show(io::IO, ::MIME"text/plain", s::GBScalar{T}) where {T}
    print(io, "GBScalar{", string(T), "}: ", string(s[]))
end

#I don't love this for multiple reasons. nnz is a misnomer, it should really be nzvals.
# The second reason is that I'm not sure if I should provide a length or size. Since
# the length can be zero, but it's a scalar... Uncomment if needed somewhere.

#nvals(s::GBScalar) = libgb.GxB_Scalar_nvals(s)
nnz(s::GBScalar) = libgb.GxB_Scalar_nvals(s)
#Base.length(s::GBScalar) = nnz(s)


text:
using Test, GameOfLife
using SharedArrays
using DistributedArrays: distribute

interior(grid) = @view grid[2:end-1,2:end-1]

function generate(world, s=Serial())
    m, n = size(world)
    grid = BitArray(undef, m+2, n+2)
    interior(grid) .= world
    grid
end

function generate(world, ::ProcParallel)
    m, n = size(world)
    grid = SharedArray{Bool,2}(m+2, n+2)
    interior(grid) .= world
    grid
end

function generate(world, ::DistributedTasks)
    grid = generate(world, Serial())
    distribute(grid)
end

world0 = [0 0 0 0 1 1 1;
          0 1 1 1 0 0 0;
          0 0 0 0 1 0 1]
world1 = [1 0 1 0 1 0 1;
          1 0 1 1 0 0 1;
          1 0 1 0 1 0 1]

@testset "general simulation tools" begin
    m, n = size(world0)
    g0 = generate(world0)
    prepare!(g0)
    @test g0 == [1 0 0 0 0 1 0 1 0;
                 1 0 0 0 0 1 1 1 0;
                 0 0 1 1 1 0 0 0 0;
                 1 0 0 0 0 1 0 1 0;
                 1 0 0 0 0 1 1 1 0]
    g1 = similar(g0)
    I = (2:m+1,2:n+1)
    life_rule!(g1, g0, I, I)
    @test interior(g1) == world1
end

const Strategies = [Serial, ThreadParallel, ProcParallel, LocalTasks, DistributedTasks]

@testset "distribution strategies" begin
    @testset "$S" for S in Strategies
        s = S()
        old = generate(world0, s)
        new = similar(old)
        step!(new, old, s)
        @test interior(old) == world0
        @test interior(new) == world1
    end
end

function randomize!(x)
    for i in eachindex(x)
        x[i] = rand(eltype(x))
    end
end

# Props to ffevotte: https://discourse.julialang.org/t/print-debug-info-for-failed-test/22311
onfail(body, x) = error("I might have overlooked something: $x")
onfail(body, _::Test.Pass) = nothing
onfail(body, _::Tuple{Test.Fail,T}) where {T} = body()

gettuples(x) = map(Tuple, x |> interior |> findall)

E_new = Dict()

# If this suite fails, run the script via including the file and plot the grid:
#
# using Plots
# scatter(gettuples(E_ref .!= E_new[S]), markersize=1)
@testset "slightly bigger dataset" begin
    n = 2*GameOfLife.TILESIZE + GameOfLife.TILESIZE÷2
    g = BitArray(undef, n, n)
    randomize!(g)
    ref = similar(g)
    step!(ref, g, Serial())

    @testset "$S" for S in Strategies[2:end]
        s = S()
        old = generate(interior(g), s)
        new = similar(old)
        step!(new, old, s)
        onfail(@test interior(new) == interior(ref)) do
            global E_old = g
            global E_ref = ref
            global E_new[S] = new
        end
    end
end

ext: jl


text:
# NOTE: This file is to be copied across the `runtests.jl` of all Banyan Julia
# projects

using Test
using Banyan
using BanyanArrays
using BanyanDataFrames

using Arrow
using CSV
using DataFrames
using FilePathsBase
using Parquet
using AWSS3
using Statistics
using Random

function include_tests_to_run(args...)
    clear_jobs()
    for arg in args
        include(arg)
    end
end

get_enabled_tests() = lowercase.(ARGS)

# NOTE: For testing, please provide the following:
# - AWS_DEFAULT_PROFILE (if you don't already have the desired default AWS account)
# - BANYAN_USERNAME (optional)
# - BANYAN_USER_ID (optional)
# - BANYAN_API_KEY (optional)
# - BANYAN_CLUSTER_NAME (required)
# - BANYAN_NWORKERS (defaults to 2)
# - BANYAN_NWORKERS_ALL (defaults to "false")
# - BANYAN_SCHEDULING_CONFIG_ALL (defaults to "false")
# - NUM_TRIALS (defaults to 1)

# TODO: Copy the below to other Banyan projects' test suites if changes are made

username = get(ENV, "BANYAN_USERNAME", nothing)
user_id = get(ENV, "BANYAN_USER_ID", nothing)
api_key = get(ENV, "BANYAN_API_KEY", nothing)
cluster_name = get(ENV, "BANYAN_CLUSTER_NAME", nothing)
nworkers = get(ENV, "BANYAN_NWORKERS", "2")
ntrials = parse(Int, get(ENV, "BANYAN_NTRIALS", "1"))

global job = create_job(
    username = username,
    user_id = user_id,
    api_key = api_key,
    cluster_name = cluster_name,
    nworkers = parse(Int32, nworkers),
    banyanfile_path = "file://res/Banyanfile.json",
    return_logs = true,
    sample_rate = get(ENV, "BANYAN_TEST_WITH_STRESS", "0") == "1" ? 1024 : parse(Int32, nworkers)
)

function run_with_job(test_fn, name)
    # This function should be used for tests that need a job to be already
    # created to run. We look at environment variables for a specification for
    # how to authenticate and what cluster to run on

    if isempty(get_enabled_tests()) ||
       any([occursin(t, lowercase(name)) for t in get_enabled_tests()])
        if get(ENV, "BANYAN_NWORKERS_ALL", "false") == "true"
            for nworkers in [16, 8, 4, 2, 1]
                with_job(
                    username = username,
                    api_key = api_key,
                    cluster_name = cluster_name,
                    nworkers = parse(Int32, nworkers),
                    banyanfile_path = "file://res/Banyanfile.json",
                    user_id = user_id,
                    destroy_job_on_exit=false
                ) do j
                    test_fn(j)
                end
            end
        elseif !isnothing(nworkers)
            with_job(job=job, destroy_job_on_exit=false) do j
                for i in 1:ntrials
                    if ntrials > 1
                        @time test_fn(j)
                    else
                        test_fn(j)
                    end
                end
            end
        end
    end
end

function run(test_fn, name)
    # This function should be used for tests that test cluster/job managemnt
    # and so they only need environment variables to dictate how to
    # authenticate. These can be read in from ENV on a per-test basis.

    if isempty(get_enabled_tests()) ||
       any([occursin(t, lowercase(name)) for t in get_enabled_tests()])
        test_fn()
    end
end

function include_all_tests()
    include_tests_to_run("test_io.jl")
    include_tests_to_run("test_small_dataset.jl")
    include_tests_to_run("test_medium_dataset.jl")
    include_tests_to_run("test_stress.jl")
    include_tests_to_run("test_groupby_filter_indexing.jl")
end

# bucket - name of cluster's s3 bucket
# path - path to write file to in bucket
# download_path - either http(s) link to a file or a local Path indicating the source of the file
function verify_file_in_s3(bucket, path, download_path)
     if !s3_exists(Banyan.get_aws_config(), bucket, path)
        if typeof(download_path) == String && (startswith(download_path, "https://") || startswith(download_path, "http://"))
            download(download_path, S3Path("s3://$(bucket)/$(path)", config=Banyan.get_aws_config()))
        else  # upload local file
            cp(Path(download_path), S3Path("s3://$(bucket)/$(path)", config=Banyan.get_aws_config()))
        end
    end
end

with_job(job=job) do j
    configure_scheduling(report_schedule=true)
    if get(ENV, "BANYAN_SCHEDULING_CONFIG_ALL", "false") == "true"
        configure_scheduling(encourage_parallelism=true, encourage_parallelism_with_batches=true)
        include_all_tests()
        configure_scheduling(encourage_parallelism=true, encourage_parallelism_with_batches=false)
        include_all_tests()
        configure_scheduling(encourage_parallelism=false, encourage_parallelism_with_batches=false)
        include_all_tests()
    else
        configure_scheduling(encourage_parallelism_with_batches=true)
        include_all_tests()
    end
end

# Requirements:
# - Name testgroups and testsets and filter and run
# - Create sample data in S3 only if needed
# - Create jobs only if needed
# - Run everything on different data and on different job configs

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
## Inter-Worker Communication Residual Constrained Alternating Minimization
# Function Definitions

function iwcpAltMinChol{ET<:Number}(b::AbstractArray{ET}, r::Int, K::Int)
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Generate L0
    m::Int,n::Int = size(b)

    # May need to micro-manage this distribution
    Ld= distribute(collect(ET, randn(m,r)))
    Rd= distribute(zeros(ET,m,r))

    # Find non-zero inds for each col, reduce data into collections
    indn = (1:n)::UnitRange{Int64}

    # This global declaration is a workaround to open issue #15276
    const indnzd = distribute([find(b[:,i] .!= 0) for i::Int in indn])
    const bsubd = distribute([b[indnzd[i],i] for i in indn])

    # Find non zero inds for each row, reduce data into collections
    indm = (1:m)::UnitRange{Int64}
    const indnz_trand = distribute([find(b[i,:] .!= 0) for i::Int in indm])
    const bsub_trand = distribute([b[i,indnz_trand[i]] for i in indm])

    b = []

    # Begin Alterations
    for k = 1:K

        # Solve for R
        remotes1 = [@spawnat p altupdate!(Rd, Ld, bsubd, indnzd) for p in workers()]

        # Wait for all updates
        [wait(remotes1[w]) for w = 1:length(remotes1)]

        # Solve for L
        remotes2 = [@spawnat p altupdate!(Ld, Rd, bsub_trand, indnz_trand) for p in workers()]

        # Wait for all updates
        [wait(remotes2[w]) for w = 1:length(remotes2)]

    end

   return Ld, Rd

end
export iwcpAltMinChol

function altupdate!{ET<:Number}(Rd::AbstractArray{ET}, Ld::AbstractArray{ET}, bsubd::AbstractArray,
                                indnzd)
# Updates the local portion of R on a single worker

    #Get info about local part of darrays
    m = localindexes(Ld)

    # Get local part
    R_tmp = localpart(Rd)

    # Update the ith row of Rd
    for i::Int in m[1]

        # Solve for what this row or R will be updated to
        R_update = collect(ET,conj(pminFrobRsubChol(convert(Array{ET},Ld[indnzd[i],:]), bsubd[i])))

        # Update this row of R
        setindex!(R_tmp,vec(R_update),i - m[1][1] + 1,:)

    end

end

function pminFrobRsubChol{ET<:Number}(L2::AbstractArray{ET},bsub::AbstractArray{ET})
# Solves a single row, meant to be mapped over a DArray

    bsub = transpose(L2)*bsub

    # Chol fact
    L2 = L2'*L2

    # Faster to do in w/o mem alloc, even though it does chol twice
    v = backsub(chol(L2),fwdsub((chol(L2))',bsub))

    return v

end
export pminFrobRsubChol

function backsub{ET<:Number}(A::AbstractArray{ET}, y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)

    x[m] = y[m]/A[m,m]

    for i = m-1:-1:1
        x[i] = (y[i]-dot(A[i,i+1:m], x[i+1:m]))/A[i,i]
    end

    return x
end

function fwdsub{ET<:Number}(A::AbstractArray{ET},y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)
    x[1] = y[1]/A[1,1]

    for i = 2:1:m
        x[i] = (y[i]-dot(A[i,i-1:-1:1], x[i-1:-1:1]))/A[i,i]
    end

    return x
end


text:
version = "0.9.1"
linux_build_id = "hf484d3e_0"
macos_build_id = "h0a44026_0"
windows_build_id = "he025d50_0"

function is_64_bits()
    return Int == Int64
end

function unpack(file, directory)
    if Sys.isunix()
        run(`mkdir -p $directory`)
        run(`tar xjf $file --directory=$directory`)
    elseif Sys.iswindows()
        exe7z = joinpath(Sys.BINDIR, "7z.exe")
        run(pipeline(`$exe7z x $file -y -so`, `$exe7z x -si -y -ttar -o$directory`))
    end
end

if Sys.iswindows()
    build_id = windows_build_id
    unpacked_file = joinpath("Library", "bin", "chemfiles.dll")
    if !is_64_bits()
        @error "There is no prebuilt chemfiles library for 32-bit windows"
    end
    platform = "win-64"
elseif Sys.islinux()
    build_id = linux_build_id
    unpacked_file = joinpath("lib", "libchemfiles.so.$version")
    if !is_64_bits()
        @error "There is no prebuilt chemfiles library for 32-bit Linux"
    end
    platform = "linux-64"
elseif Sys.isapple()
    build_id = macos_build_id
    unpacked_file = joinpath("lib", "libchemfiles.$version.dylib")
    if !is_64_bits()
        @error "There is no prebuilt chemfiles library for 32-bit macOS"
    end
    platform = "osx-64"
else
    error("Could not find a prebuilt chemfiles library for the current platform")
end

URL = "https://anaconda.org/conda-forge/chemfiles-lib/$version/download/$platform/chemfiles-lib-$version-$build_id.tar.bz2"
LOCAL_ARCHIVE = joinpath(@__DIR__, basename(URL))

@info "Downloading pre-build library ..."
download(URL, LOCAL_ARCHIVE)
@info "Unpacking library ..."
unpack(LOCAL_ARCHIVE, joinpath(@__DIR__, "usr"))

LIBPATH = joinpath(@__DIR__, "usr", unpacked_file)

write(joinpath(@__DIR__, "deps.jl"), """
# This is an auto-generated file; do not edit

using Libdl

# Macro to load a library
macro checked_lib(libname, path)
    Libdl.dlopen_e(path) == C_NULL && @error "Unable to load \\n\\n\$libname (\$path)\\n\\nPlease re-run Pkg.build(package), and restart Julia."
    quote const \$(esc(libname)) = \$path end
end

# Load dependencies
@checked_lib libchemfiles "$(escape_string(LIBPATH))"
""")

@info "Done"


text:
#=
Threshold methods

Running thresholds

=#

export ThresholdMeanP, ThresholdMeanN

#running median

#running mean

mutable struct ThresholdMeanP <: Threshold
    m_k::Float64
    k::Int64
    m_l::Float64
    s_k::Float64
    s_l::Float64
    mean::Float64
    stds::Float64
end

ThresholdMeanP()=ThresholdMeanP(0.0,1,0.0,0.0,0.0,0.0,3.0)

ThresholdMeanP(stds::Float64)=ThresholdMeanP(0.0,1,0.0,0.0,0.0,0.0,stds)

function threshold(t::ThresholdMeanP,sort::Sorting,v,i)

    sort.t.k+=1

    sort.t.m_k=sort.t.m_l+(v[i,sort.id]-sort.t.m_l)/sort.t.k
    sort.t.s_k=sort.t.s_l+(v[i,sort.id]-sort.t.m_l)*(v[i,sort.id]-sort.t.m_k)

    myvar=sort.t.s_k/(sort.t.k-1)

    sort.t.mean+=v[i,sort.id]
    mymean=sort.t.mean/sort.t.k

    sort.thres=mymean+sort.t.stds*sqrt(myvar)

    sort.t.m_l=sort.t.m_k
    sort.t.s_l=sort.t.s_k

    nothing
end

mutable struct ThresholdMeanN <: Threshold
    m_k::Float64
    k::Int64
    m_l::Float64
    s_k::Float64
    s_l::Float64
    mean::Float64
    stds::Float64
end

ThresholdMeanN()=ThresholdMeanN(0.0,1,0.0,0.0,0.0,0.0,3.0)

ThresholdMeanN(stds::Float64)=ThresholdMeanN(0.0,1,0.0,0.0,0.0,0.0,stds)

function threshold(t::ThresholdMeanN,sort::Sorting,v,i)

    sort.t.k+=1

    sort.t.m_k=sort.t.m_l+(v[i,sort.id]-sort.t.m_l)/sort.t.k
    sort.t.s_k=sort.t.s_l+(v[i,sort.id]-sort.t.m_l)*(v[i,sort.id]-sort.t.m_k)

    myvar=sort.t.s_k/(sort.t.k-1)

    sort.t.mean+=v[i,sort.id]
    mymean=sort.t.mean/sort.t.k

    sort.thres=mymean-sort.t.stds*sqrt(myvar)

    sort.t.m_l=sort.t.m_k
    sort.t.s_l=sort.t.s_k

    nothing
end


text:
using GrammarExpts
using DecisionTrees


text:
# Structure Learning
using Clustering

"""
    chainruleSPN(cardinalities::AbstractVector,depth::Int64=1,root::Variable=1)::SumProductNetwork

Generates a sum-product network encoding a naive joint distribution factorization according to the application of the Chain Rule for finite-valued variables with given `cardinalities`. The network is rooted at `root` node id, and performs a recursive decomposition starting at variable `depth`.

# Arguments

- `depth::Integer`: id of variable to decompose at.
- `root::Variable`: id of the root node of the SPN.

"""
function chainruleSPN(cardinalities::AbstractVector,depth::Integer=1,root::Integer=1)::SumProductNetwork
    cardinality = cardinalities[depth]
    nodelist = Vector{Node}()
    if depth == length(cardinalities)
        # final layer
        push!(nodelist,SumNode([i+root for i=1:cardinality],ones(cardinality)/cardinality))
        for value = 1:cardinality
            push!(nodelist,IndicatorFunction(depth,value))
        end
        return SumProductNetwork(nodelist)
        # return SumProductNetwork(Vector([CategoricalDistribution(depth,ones(cardinality)/cardinality)]))
    end
    push!(nodelist,SumNode([i+root for i=1:cardinality],ones(cardinality)/cardinality))
    index = root + cardinality + 1
    children = Vector{Node}()
    for value = 1:cardinality
        # Add product node with left child Indicator X=value and Left child the recursion
        push!(nodelist,ProductNode([index,index+1]))
        push!(children,IndicatorFunction(depth,value))
        index += 1
        # Now recur to build subSPNs
        S = chainruleSPN(cardinalities,depth+1,index)
        append!( children, nodes(S) )
        index += length(S)
    end
    append!(nodelist,children)
    SumProductNetwork(nodelist)
end

"""
    fit(node, Data; params)

Learns (smoothed) Maximum Likelihood Estimators of univariate distribution from vector Data of i.i.d. samples.
"""
function fit(node::LeafNode, Data::AbstractVector; params...)
    if isa(node,CategoricalDistribution)
        for datum in Data
           node.values[Int(datum)] += 1
        end
        node.values ./= length(Data)
    else
        @error "Learning of this type of node not implemented"
    end

end

# TODO: Implement correlation matrix and splitting
# This part was largely inspired by https://github.com/trappmartin/SumProductNetworks.jl/blob/master/src/structureUtilities.jl

"""
    splitrows(Data, scope, instances, minsamplesize)

Learns a sum node by clustering instances of columns in scope into two groups by k-means clustering. Returns a list of pairs (cluster,weight) where cluster is a vector of indexes view of a cluster and weights is the proportion of instances in that cluster.
"""
function splitrows(Data::AbstractMatrix, scope, instances, minsamplesize)
    @assert length(instances) > minsamplesize "Sample size too small for clustering"
    data = view(Data, instances, scope)
    results = kmeans(transpose(data), 2)
    cluster1 = instances[results.assignments .== 1] # view(Data, results.assignments .== 1, :)
    cluster2 = instances[results.assignments .!= 1] # view(Data, results.assignments .!= 1, :)
    weight = length(cluster1)/length(instances)

    if length(cluster1) < minsamplesize || length(cluster2) < minsamplesize
        return ((instances,1.0),)
    end
    return ((cluster1, weight), (cluster2, 1.0 - weight))
end

"""
    splitcolumns(Data, scope, instances, minsamplesize)

Learns a product node by clustering columns using the connected components of the association matrix after filtering out low score edges. Returns list of scopes of child nodes.
"""
function splitcolumns(Data::AbstractMatrix, scope, instances, minsamplesize)
    @assert length(scope) > 1 "Scope has size ($length(scope)); Must have at least two variables to split columns"
    if length(instances) < minsamplesize
        # stopping criterion, too small sample size: factorize
        return ((variable,) for variable in scope)
    else
        # generate random binary partition
        leftscope = scope[Random.randperm(length(scope))[1:floor(Integer,length(scope)/2)]]
        rightscope = setdiff(scope,leftscope)
        return (leftscope, rightscope)
    end
end

"""
    learnspn(Data,cardinalities,minsamplesize)

Implement Gen and Domingos 2013 top-down learning algorithm.
"""
function learnspn(Data::AbstractMatrix,cardinalities,minsamplesize=100)::SumProductNetwork

    numrows,numcolumns = size(Data)
    queue = Vector{Tuple}()
    nodes = Vector{Node}([])

    children = []
    to = Int64[]
    from = Int64[]
    weights = Float64[]

    # Start with a sum node as root
    # push!(queue, (1,1.0,collect(1:numcolumns),view(Data,:,:)))
    parent,scope,instances = 1,collect(1:numcolumns),collect(1:numrows)
    push!(nodes,SumNode(scope)) # add node to list
    node = 1 # node id
    push!(children,[])
    for (cluster,weight) in splitrows(Data, scope, instances, minsamplesize)
        push!(queue,(node,weight,scope,cluster)) # enqueue
    end

    while length(queue) > 0

        parent,weight,scope,instances = popfirst!(queue) # dequeue
        #parent,weight,scope,dataset = pop!(queue) # remove from top

        # Stopping criterion, create univariate leaf
        if  length(scope) == 1
            leaf = CategoricalDistribution(scope[1],zeros(cardinalities[scope[1]]))
            fit(leaf,view(Data,instances,scope[1])) # learn node distribution from data
            push!(nodes,leaf) # add node to node list (node id is list size after insertion)
            push!(children[parent],length(nodes)) # add arc parent -> leaf
            push!(children,[]) # add outbound arcs (none)
            if isa(nodes[parent],SumNode)
                push!(from,parent)
                push!(to,node)
                push!(weights,weight) # add weight to arc parent -> node
            end
        elseif isa(nodes[parent],ProductNode)
            # Create sum node
            push!(nodes,SumNode(scope)) # add node to list
            node = length(nodes) # node id is node list size
            push!(children,[])
            push!(children[parent],node)
            for (cluster,weight) in splitrows(Data, scope, instances, minsamplesize)
                push!(queue,(node,weight,scope,cluster)) # enqueue
            end
        else # isa(nodes[parent],SumNode)
            # Create product node
            push!(nodes,ProductNode(scope)) # add node to list
            node = length(nodes) # node id
            push!(children[parent],node)
            push!(children,[])
            push!(from,parent)
            push!(to,node)
            push!(weights,weight) # add weight to arc parent -> node
            for subscope in splitcolumns(Data, scope, instances, minsamplesize)
                push!(queue,(node,NaN,subscope,instances)) #enqueue
            end
        end

    end

    weightmatrix = sparse(to,from,weights)

    spn = SumProductNetwork(nodes, children, weightmatrix)

    return spn
end


text:
# Only run coverage on linux with julia v1.0.
get(ENV, "TRAVIS_OS_NAME", "")       == "linux"   || exit()
get(ENV, "TRAVIS_JULIA_VERSION", "") == "1.0" || exit()

using Pkg; Pkg.instantiate()

using Coverage

cd(joinpath(@__DIR__, "..")) do
    # push coverage results to Coveralls
    # Coveralls.submit(Coveralls.process_folder())'
    # push coverage results to Codecov
    Codecov.submit(Codecov.process_folder())'
end


text:
function remove(collection::GeoJSON.FeatureCollection, key::AbstractString, value::AbstractString)
    fc = GeoJSON.FeatureCollection()
    for feature in collection.features
        if haskey(feature.properties, key) && feature.properties[key] != value
            push!(fc.features, deepcopy(feature))
        end
    end
    fc
end

function remove(collection::GeoJSON.FeatureCollection, key::AbstractString, value::Integer)
    fc = GeoJSON.FeatureCollection()
    for feature in collection.features
        if haskey(feature.properties, key) && feature.properties[key] != value
            push!(fc.features, deepcopy(feature))
        end
    end
    fc
end

function filter(collection::GeoJSON.FeatureCollection, property::AbstractString, value::AbstractString)
    fc = GeoJSON.FeatureCollection()
    for feature in collection.features
        if haskey(feature.properties, key) && feature.properties[key] == value
            push!(fc.features, deepcopy(feature))
        end
    end
    fc
end

function filter(collection::GeoJSON.FeatureCollection, key::AbstractString, value::Integer)
    fc = GeoJSON.FeatureCollection()
    for feature in collection.features
        if haskey(feature.properties, key) && feature.properties[key] == value
            push!(fc.features, deepcopy(feature))
        end
    end
    fc
end


text:
export exp_J2!,
       exp_J1!

"""
    exp_J2!(state::Vector{T},
            basis::Vector{Int},
            α::AbstractMatrix{T}) where {T<:Number}

Apply the exponential operator e^{∑_pq α_pq n_p n_q} to `state`.
"""
function exp_J2!(state::Vector{T},
                 basis::Vector{Int},
                 α::AbstractMatrix{T}) where {T<:Number}
    # Deduce the number of sites from the highest basis state.
    L = 0
    max_bs = basis[end]
    for i=0:31 # Arbitrary very high number.
        if (max_bs>>i)&1 == 1
            L = i+1
        end
    end

    for (ibs, bs) in enumerate(basis)
        c = zero(T)
        for i=0:L-1
            (bs>>i)&1 != 1 && continue
            c += α[i+1, i+1]
            for j=i+1:L-1
                if (bs>>j)&1 == 1
                    c += α[i+1, j+1] + α[j+1, i+1]
                end
            end
        end
        state[ibs] *= exp(c)
    end
    return state
end

"""
    exp_J1!(state::Vector{T},
            basis::Vector{Int},
            α::AbstractMatrix{T}) where {T<:Number}

Apply the exponential operator e^{∑_p α_p n_p} to `state`.
"""
function exp_J1!(state::Vector{T},
                 basis::Vector{Int},
                 α::Vector{T}) where {T<:Number}
    # Deduce the number of sites from the highest basis state.
    L = 0
    max_bs = basis[end]
    for i=0:31 # Arbitrary very high number.
        if (max_bs>>i)&1 == 1
            L = i+1
        end
    end

    for (ibs, bs) in enumerate(basis)
        c = zero(T)
        for i=0:L-1
            if (bs>>i)&1 == 1
                c += α[i+1]
            end
        end
        state[ibs] *= exp(c)
    end
    return state
end


text:
# acc satisfiability script
using OVERTVerify
using LazySets
using Dates

@testset "feasibility" begin

    # ACC CONTROLLER
    # 3x20 controller here: https://github.com/souradeep-111/sherlock/tree/master/systems_with_networks/ARCH_2019/ACC
    controller = "acc_controller"
    controller_filepath = "../nnet_files/jmlr/"*controller*".nnet"
    # println("Controller is: ", controller)
    query = OvertQuery(
        ACC,  # problem
        controller_filepath,    # network file
        Id(),      	# last layer activation layer Id()=linear, or ReLU()=relu
        "MIP",     	# query solver, "MIP" or "ReluPlex"
        2,        	# ntime
        0.1,       	# dt
        -1,        	# N_overt
        )

    # x1,x2,x3 are lead vehicle variables
    # x4,x5,x6 are ego vehicle variables
    x_lead = [90,110]
    v_lead = [32, 32.2]
    gamma_lead = [0,0]
    x_ego = [10,11]
    v_ego = [30, 30.2]
    gamma_ego = [0, 0]
    var_list = [x_lead, v_lead, gamma_lead, x_ego, v_ego, gamma_ego]
    input_set = Hyperrectangle(
        low=[v_range[1] for v_range in var_list],
        high=[v_range[2] for v_range in var_list]
        )
    # test using Constraint on output instead of hyperrectangle!
    # desired property to be proven is:
    # D_rel = x_lead - x_ego >= D_safe = D_default + T_gap * v_ego
    # x_lead + -1*x_ego - T_gap*v_ego >= D_default = 10
    # BUT WE MUST NEGATE THIS (and so we flip >= to <=) (and then want it to always be unsat)
    # THe measurement model gives x_lead + -1*x_ego - T_gap*v_ego  so we just have to pass it through with coeff=1.
    output_constraint = Constraint([1.], :(<=), 10)
    # NOTE: this constraint is written for use with the measurement, so we want to apply the measurement model here, and we set apply_meas to true in the call to symbolic_satisfiability below

    SATus, vals, stats = symbolic_satisfiability(query, input_set, output_constraint, apply_meas=true)

    @test SATus == "unsat"

end


text:
function testemptygroup()
  Group([])
end
testemptygroup()

function testgroupnesting1()
  g1 = Group([Node(1,1,1),
      Node(2,2,2),
      Node(3,3,3)])
  verified =
  """
  N_1 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  N_2 x=2.000000000e+00 y=2.000000000e+00 z=2.000000000e+00
  N_3 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  """
  testelement(g1,verified)
  g2 = Group([Node(4,4,4),
      g1,
      Node(:a,5,5,5)])
  verified =
  """
  N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  N_2 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  N_3 x=2.000000000e+00 y=2.000000000e+00 z=2.000000000e+00
  N_4 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  Na x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  """
  testelement(g2,verified)
  g3 = Group([g2,Node(6,6,6)])
  verified =
  """
  N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  N_2 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  N_3 x=2.000000000e+00 y=2.000000000e+00 z=2.000000000e+00
  N_4 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  Na x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  N_5 x=6.000000000e+00 y=6.000000000e+00 z=6.000000000e+00
  """
  testelement(g3,verified)
  g4 = Group([Node(7,7,7),g3])
  verified =
  """
  N_1 x=7.000000000e+00 y=7.000000000e+00 z=7.000000000e+00
  N_2 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  N_3 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  N_4 x=2.000000000e+00 y=2.000000000e+00 z=2.000000000e+00
  N_5 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  Na x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  N_6 x=6.000000000e+00 y=6.000000000e+00 z=6.000000000e+00
  """
  testelement(g4,verified)
end
testgroupnesting1()

function testgroupunits1()
  g1 = Group([Node(1,1,1),
              Units("in"),
              Node(3,3,3)],
              Dict(),
              Units("mm"))
  g2 = Group([Node(4,4,4)], Dict(), Units("cm"))
  g3 =  Group([Node(5,5,5)])
  g4 = Group([Node(6,6,6),
              Units("um"),
              Node(7,7,7)])
  g5 = Group([g4,g2])
  verified =
  """
  N_1 x=6.000000000e+00 y=6.000000000e+00 z=6.000000000e+00
  .units um
  N_2 x=7.000000000e+00 y=7.000000000e+00 z=7.000000000e+00
  .units cm
  N_3 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units um
  """
  testelement(g5, verified)
  g6 = Group([g2,g4])
  verified =
  """
  .units cm
  N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units mm
  N_2 x=6.000000000e+00 y=6.000000000e+00 z=6.000000000e+00
  .units um
  N_3 x=7.000000000e+00 y=7.000000000e+00 z=7.000000000e+00
  """
  testelement(g6, verified)
  g7 = Group([g3,g2])
  verified =
  """
  N_1 x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  .units cm
  N_2 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units mm
  """
  testelement(g7, verified)
  g8 = Group([g2,g3])
  verified =
  """
  .units cm
  N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units mm
  N_2 x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  """
  testelement(g8, verified)
  g9 = Group([Units("um"),g1,g2,g1,g2])
  verified =
  """
  .units um
  .units mm
  N_1 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  .units in
  N_2 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  .units um
  .units cm
  N_3 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units um
  .units mm
  N_1 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  .units in
  N_2 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  .units um
  .units cm
  N_3 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units um
  """
  testelement(g9, verified)
  g10 = Group([Group([g1,g3]),Group([g4,g2])])
  verified =
  """
  .units mm
  N_1 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  .units in
  N_2 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  .units mm
  N_3 x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  N_4 x=6.000000000e+00 y=6.000000000e+00 z=6.000000000e+00
  .units um
  N_5 x=7.000000000e+00 y=7.000000000e+00 z=7.000000000e+00
  .units cm
  N_6 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units um
  """
  testelement(g10, verified)
  g11 = Group([Node(4,4,4), g1, Node(:a,5,5,5)])
  verified =
  """
  N_1 x=4.000000000e+00 y=4.000000000e+00 z=4.000000000e+00
  .units mm
  N_2 x=1.000000000e+00 y=1.000000000e+00 z=1.000000000e+00
  .units in
  N_3 x=3.000000000e+00 y=3.000000000e+00 z=3.000000000e+00
  .units mm
  Na x=5.000000000e+00 y=5.000000000e+00 z=5.000000000e+00
  """
  testelement(g11, verified)
end
testgroupunits1()


text:
import Stochy: rewrite_address, rewrite_param_addresses

@test rewrite_address(list(1,5), list(), list()) == list(1,5)
@test rewrite_address(list(2,5), list(), list(1)) == list(1,2,5)
@test rewrite_address(list(1,5), list(1), list()) == list(5)
@test rewrite_address(list(1,5), list(1), list(2)) == list(2,5)
@test rewrite_address(list(1,2,4), list(1,2), list(3)) == list(3,4)
@test rewrite_address(list(1,4,5), list(1), list(2,3)) == list(2,3,4,5)

@test rewrite_param_addresses({list(3,1,0) => :param}, list(1,0), list(2)) == {list(3,2) => :param}


text:
let
    for (name, prob) in Optim.UnivariateProblems.examples
        results = optimize(prob.f, prob.bounds..., method = GoldenSection())

        @assert Optim.converged(results)
        @assert norm(Optim.minimizer(results) - prob.minimizers) < 1e-7
    end
end


text:
@enum BrainflowExitCodes begin

    STATUS_OK = 0
    PORT_ALREADY_OPEN_ERROR = 1
    UNABLE_TO_OPEN_PORT_ERROR = 2
    SER_PORT_ERROR = 3
    BOARD_WRITE_ERROR = 4
    INCOMMING_MSG_ERROR = 5
    INITIAL_MSG_ERROR = 6
    BOARD_NOT_READY_ERROR = 7
    STREAM_ALREADY_RUN_ERROR = 8
    INVALID_BUFFER_SIZE_ERROR = 9
    STREAM_THREAD_ERROR = 10
    STREAM_THREAD_IS_NOT_RUNNING = 11
    EMPTY_BUFFER_ERROR = 12
    INVALID_ARGUMENTS_ERROR = 13
    UNSUPPORTED_BOARD_ERROR = 14
    BOARD_NOT_CREATED_ERROR = 15
    ANOTHER_BOARD_IS_CREATED_ERROR = 16
    GENERAL_ERROR = 17
    SYNC_TIMEOUT_ERROR = 18
    JSON_NOT_FOUND_ERROR = 19
    NO_SUCH_DATA_IN_JSON_ERROR = 20
    CLASSIFIER_IS_NOT_PREPARED_ERROR = 21
    ANOTHER_CLASSIFIER_IS_PREPARED_ERROR = 22
    UNSUPPORTED_CLASSIFIER_AND_METRIC_COMBINATION_ERROR = 23
end

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow


text:
# packages and code loaded on all processes
using DelimitedFiles, Statistics
using StaticArrays, Parameters, Dierckx, SimulatedAnnealing, Setfield, Interpolations
using OrdinaryDiffEq, Plots

include("musclemodel/torque_generator.jl")
include("model/parameters.jl")
include("model/eom.jl")
include("model/functions.jl")
include("model/plotting.jl")
include("model/functions_ext.jl")
include("setup.jl")
include("callbacks.jl")
include("optimisation.jl")

# set up initial model
inputs = load_inputs()#=const =#
inputs.initial_conditions[:vcmx] = 9.6
p, u₀ = set_values(inputs)#=const =#
matching_data = inputs.matching_data#=const =#

# time span
tspan = (0.0, 0.317)

# initialise problem
prob = ODEProblem(eom, u₀, tspan, p)

# swing time
TSW = 0.374


text:
@testset "fsm_ack_rst_win_relaxed.jl" begin
    base_seq = WrappingInt32(1 << 31)
    DEFAULT_CAPACITY = 64000
    @testset "in ESTABLISHED, send unacceptable segments and ACKs" begin
        conn = fsm_in_established(base_seq - 1, base_seq - 1)

        send_ack!(conn, base_seq, base_seq - 1)
        expect_no_seg(conn)
        send_ack!(conn, base_seq, base_seq - 1)
        expect_no_seg(conn)

        send_byte!(conn, base_seq - 1, base_seq, Char(1))
        @test conn |> unassembled_bytes == 0
        expect_one_seg(conn, ackno=base_seq)

        send_byte!(conn, base_seq + DEFAULT_CAPACITY - 1, base_seq, Char(1))
        @test conn |> unassembled_bytes == 1
        expect_one_seg(conn, ackno=base_seq)
        expect_no_data(conn)

        send_byte!(conn, base_seq, base_seq, Char(1))
        @test conn |> unassembled_bytes == 1
        expect_one_seg(conn, ackno=base_seq+1)
        expect_data!(conn)
        send_rst!(conn, base_seq + 1)
        expect_state(conn, JLSponge.RESET)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
function getradec(obs)::Tuple{Float64, Float64}
    return (float(obs.ra), float(obs.dec))
end

function getmjd(obs)::Int
    return Int(obs.mjd)
end

function getsource(obs)::Symbol
    return Symbol(obs.source)
end


text:
function batch_reactor(; store_results=false::Bool)
    datadir = nothing
    if store_results
        datadir = generate_data_subdir(@__FILE__)
    end

    function system!(du, u, p, t, controller)
        # fixed parameters
        α, β, γ, δ = 0.5f0, 1.0f0, 1.0f0, 1.0f0

        # neural network outputs controls taken by the system
        y1, y2 = u
        c1, c2 = controller(u, p)

        # dynamics of the controlled system
        y1_prime = -(c1 + α * c1^2) * y1 + δ * c2
        y2_prime = (β * c1 - γ * c2) * y1

        # update in-place
        @inbounds begin
            du[1] = y1_prime
            du[2] = y2_prime
        end
    end

    # define objective function to optimize
    function loss(controlODE, params)
        sol = solve(controlODE, params)
        return -Array(sol)[2, end]  # second variable, last value, maximize
    end

    # initial conditions and timepoints
    u0 = [1.0f0, 0.0f0]
    tspan = (0.0f0, 1.0f0)
    Δt = 1f-2

    # set arquitecture of neural network controller
    controller = FastChain(
        FastDense(2, 12, tanh_fast),
        FastDense(12, 12, tanh_fast),
        FastDense(12, 2),
        (x, p) -> 5 * sigmoid_fast.(x),  # controllers ∈ (0, 5)
    )

    controlODE = ControlODE(controller, system!, u0, tspan; Δt)

    θ = initial_params(controller)

    # variables for streamplots
    phase_time = 0.0f0
    xlims, ylims = (0.0f0, 1.5f0), (0.0f0, 1.5f0)
    coord_lims = [xlims, ylims]
    # xwidth = xlims[end] - xlims[1]
    # ywidth = ylims[end] - ylims[1]
    # start_points_x = range(u0[1] - 1e-4, u0[1] - xwidth/5; length=3)
    # start_points_y = range(u0[2] + 1e-4, u0[2] + ywidth/5; length=3)

    _, states_raw, _ = run_simulation(controlODE, θ)
    marker = FinalState(; points=states_raw[:, end], fmt="m*", markersize=20)
    phase_portrait(
        controlODE,
        θ,
        phase_time,
        coord_lims;
        markers=[marker],
        # start_points_x, start_points_y,
        start_points=reshape(u0 .+ (-1e-4, 0), 1, 2),
        title="Initial policy",
    )

    # closures to comply with required interface
    loss(params) = loss(controlODE, params)
    function plotting_callback(params, loss)
        return plot_simulation(controlODE, params; only=:controls, show=loss)
    end

    @info "Training..."
    result = sciml_train(
        loss,
        θ,
        LBFGS(; linesearch=BackTracking(; iterations=10));
        # iterations=100,
        allow_f_increases=true,
        # cb=plotting_callback,
    )

    store_simulation(
        "optimized",
        controlODE,
        result.minimizer;
        metadata=Dict(:loss => loss(result.minimizer)),
        datadir,
    )
    plot_simulation(controlODE, result.minimizer; only=:controls)

    θ_opt = result.minimizer

    _, states_opt, _ = run_simulation(controlODE, θ_opt)
    marker = FinalState(; points=states_opt[:, end], fmt="m*", markersize=20)
    return phase_portrait(
        controlODE,
        θ_opt,
        phase_time,
        coord_lims;
        markers=[marker],
        # start_points_x, start_points_y,
        start_points=reshape(u0 .+ (-1e-4, 0), 1, 2),
        title="Optimized policy",
    )
end  # script wrapper


text:
################################################################################
#
#   ADDING SINGLE BONDS
#       - to unitcells
#       - to lattices
#
################################################################################

# Add a single bond to a unitcell and return its index in the list
function addBond!(
            unitcell :: U,
            from  :: Integer,
            to    :: Integer,
            label :: L = getDefaultLabel(L),
            wrap  :: NTuple{N,Int64} = NTuple{N,Int64}(zeros(N)),
            returning_bond :: Bool = true
        ) :: Int64 where {L,N,B<:AbstractBond{L,N},S,U<:AbstractUnitcell{S,B}}

    # create a new site
    bond = newBond(B, from, to, label, wrap)

    # obtain the list of bonds in the unitcell
    bond_list = bonds(unitcell)
    # maybe push the returning bond into the list first
    if returning_bond
        # build return wrap
        return_wrap = NTuple{N,Int64}([-w for w in wrap])
        # push to list
        push!(bond_list, newBond(B, to, from, label, return_wrap))
    end
    # push to the list of bonds
    push!(bond_list, bond)
    # set the list of bonds in the unitcell
    bonds!(unitcell, bond_list)

    # return the new bond index
    return length(bond_list)
end

# Add a single bond to a lattice and return its index in the list
function addBond!(
            lattice  :: LA,
            from  :: Integer,
            to    :: Integer,
            label :: L = getDefaultLabel(L),
            wrap  :: NTuple{N,Int64} = NTuple{N,Int64}(zeros(N)),
            returning_bond :: Bool = true
        ) :: Int64 where {L,N,B<:AbstractBond{L,N},S,U,LA<:AbstractLattice{S,B,U}}

    # create a new site
    bond = newBond(B, from, to, label, wrap)

    # obtain the list of bonds in the lattice
    bond_list = bonds(lattice)
    # maybe push the returning bond into the list first
    if returning_bond
        # build return wrap
        return_wrap = NTuple{N,Int64}([-w for w in wrap])
        # push to list
        push!(bond_list, newBond(B, to, from, label, return_wrap))
    end
    # push to the list of bonds
    push!(bond_list, bond)
    # set the list of bonds in the lattice
    bonds!(lattice, bond_list)

    # return the new bond index
    return length(bond_list)
end

# export the function
export addBond!


text:
using Gridap
using GridapGmsh

using GridapMHD
using GridapMHD: analytical_shercliff_u
using GridapMHD: analytical_shercliff_j
using GridapMHD: transient_driver_inductionless_MHD
using GridapMHD: driver_inductionless_MHD

# Problem setting
ρ = 1.0
ν = 1.0
σ = 1.0

U0 = 10.0
B0 = 10.0
L = 1.0
Re = U0 * L / ν
Ha = B0 * L * sqrt(σ/(ρ*ν))

# Mesh load
#model = GmshDiscreteModel("driver/duct.msh")
partition = (6,6,2)
domain = (-1,1,-1,1,0,0.1)
model0 = CartesianDiscreteModel(domain,partition)
model = simplexify(model0)

# Boundary conditions location

function hartmann_walls(coord)
  y = coord[2]
  tol = 1e-6
  if (1.0-tol < y ) | ( y < -1.0+tol )
    return true
  end
  false
end

function side_walls(coord)
  x = coord[1]
  tol = 1e-6

  if (1.0-tol < x ) | ( x < -1.0+tol )
    return true
  end
  false
end

function inlet(coord)
  z = coord[3]
  tol = 1e-6
  if (z < tol )
    return true
  end
  false
end

function outlet(coord)
  z = coord[3]
  tol = 1e-6

  if (0.1 - tol < z )
    return true
  end
  false
end

# Name inlet, outlet, and walls element sets.
add_entity!(model,hartmann_walls,"hartmann_walls")
add_entity!(model,side_walls,"side_walls")
add_entity!(model,inlet,"inlet")
add_entity!(model,outlet,"outlet")

side_wall_semilength = 1.0
hartmann_wall_semilength = 1.0
number_fourier_sumands = 10
wall_thickness = 1.0
wall_conductivity = 0.0
K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3
u0(x) = analytical_shercliff_u(side_wall_semilength, hartmann_wall_semilength,
  wall_thickness, wall_conductivity, σ, ν*ρ, ∂p∂z, Ha, number_fourier_sumands, x)/U0

j0(x) = analytical_shercliff_j(side_wall_semilength, hartmann_wall_semilength,
  wall_thickness, wall_conductivity, σ, ν*ρ, ∂p∂z, Ha, number_fourier_sumands, x)/(σ * B0 * U0)

# Boundary conditions values
function g_u(coord,t)
  if (abs(coord[3]) < 1e-6) | (abs(coord[3]-0.1) < 1e-6)
    return u0(coord)
  else
    return VectorValue(0.0, 0.0, 0.0)
  end
end
g_u(t::Real) = x -> g_u(x,t)
g_u(x::VectorValue) =  g_u(x,0.0)

function g_j(x,t)
  #return j0(x)
  return VectorValue(0.0,0.0,0.0)
end
g_j(t::Real) = x -> g_j(x,t)
g_j(x::VectorValue) = g_j(x,0.0)

# Body force
function f_u(x)
#  VectorValue(0.0,0.0, -∂p∂z) * L/U0^2
  VectorValue(0.0,0.0, 0.0)
end

# Initial conditions
#function u0(coord)
#  x = coord[1]; y = coord[2]
#  return VectorValue(0.0, 0.0, (x-1)*(x+1) * (y-1)*(y+1) )
#end

xh_t, trian, dΩ = driver_inductionless_MHD(;
  #t0=0.0,
  #tF=500.0,
  #Δt=500.0,
  #θ=1.0,
  Re=Re,
  Ha=Ha,
  model=model,
  fluid_dirichlet_tags = ["inlet","side_walls","hartmann_walls"],
  fluid_neumann_tags = ["outlet"],
  magnetic_dirichlet_tags = ["inlet","side_walls","hartmann_walls","outlet"],
  magnetic_neumann_tags = [],
  c_w = 0.0,
  α = 10.0,
  fluid_dirichlet_conditions = g_u ,
  magnetic_dirichlet_conditions = g_j,
  fluid_body_force = f_u,
  #u0 = u0,
  #j0 = j0,
  constraint_presures = (false,false),
  usegmres = false,
  resultsfile = nothing
)

uh, ph, jh, φh = xh_t
divj = ∇⋅jh
writevtk(trian, "TransientDuct", order=2,
  cellfields=["uh"=>uh, "ph"=>ph, "jh"=>jh, "φh"=>φh, "divj"=>divj])

# xh, tf = xh_t[end]
# uh, ph, jh, φh = xh
# divj = (∇⋅jh)
#
# # Scale unknowns
# uh = uh * U0
# jh = jh * σ * B0 * U0
#
# # Analyical solutions
# side_wall_semilength = 1.0
# hartmann_wall_semilength = 1.0
# number_fourier_sumands = 10
# wall_thickness = 1.0
# wall_conductivity = 0.0
# K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
# ∂p∂z = -Re * K / L^3


text:
@testset "Shift register" begin
    registers = 8191
    for i = 1:8191
        output_xb, registers = @inferred GNSSSignals.shift_register(
            registers, [1, 3, 4, 6, 7, 8, 12, 13]
        )
        results = [2788, 2056 , 3322, 2087, 6431]
        if (i in [266, 804, 1559, 3471, 5343])
            @test registers in results
        end
    end
    @test registers == 8191
end

@testset "GPS L5" begin

    gpsl5 = GPSL5(use_gpu = Val(true))
    @test @inferred(get_center_frequency(gpsl5)) == 1.17645e9Hz
    @test @inferred(get_code_length(gpsl5)) == 10230
    @test @inferred(get_secondary_code_length(gpsl5)) == 10
    CUDA.@allowscalar @test @inferred(get_code(gpsl5, 0, 1)) == 1
    CUDA.@allowscalar @test @inferred(get_code(gpsl5, 0.0, 1)) == 1
    CUDA.@allowscalar @test @inferred(get_code_unsafe(gpsl5, 0.0, 1)) == 1
    @test @inferred(get_data_frequency(gpsl5)) == 100Hz
    @test @inferred(get_code_frequency(gpsl5)) == 10230e3Hz
    CUDA.@allowscalar gpsl5.codes[:,1] == L5_SAT1_CODE_GPU

end

@testset "Neuman sequence" begin
    gpsl5 = GPSL5(use_gpu = Val(true))
    CUDA.@allowscalar begin
    code = get_code.(gpsl5, 0:103199, 1)
    satellite_code = code[1:10230]
    neuman_hofman_code = [0,0,0,0,1,1,0,1,0,1]
    for i = 1:10
        @test code[1+10230*(i-1):10230*i] == (
            satellite_code .* (Int8(-1)^neuman_hofman_code[i])
        )
    end
    @test code[1:10230] == code[10231:20460]
    end
end


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This routine solves Poisson's equation using piecewise linear finite element methods.
The matrix is assembled locally, and boundary nodes are not removed."

m = 50 # number of elements
dt = .001
T = .2 # final time

α(t) = 1.0 * cos(10*pi*t)
β(t) = 2.0 * sin(10*pi*t)
f(x) = 0
κ(x) = 1.0 + 2.0*(x > 0)
u0(x) = exp(-25*x^2)

# define spatial grid
x = LinRange(-1,1,m+2)

# construct local FEM matrix
A = spzeros(m+2,m+2)
M = spzeros(m+2,m+2)
b = zeros(m+2)
A_local = [1 -1;-1 1]
M_local = [2 1; 1 2] ./ 3
for e = 1:m+1 # loop over the intervals
    ids = e:e+1
    h_e = x[e+1]-x[e]
    x_midpt = .5*(x[e+1]+x[e])

    A[ids,ids] .+= κ(x_midpt) * A_local / h_e
    b[ids] .+= .5 * h_e * f(x_midpt) # f(.5*(x_i+x_i+1))

    M[ids,ids] .+= M_local * .5* h_e
end

# modify for boundary conditions
function impose_Dirichlet_BC!(A,b,i,val)
    b .-= A[:,i]*val
    b[i] = val
    A[i,:] .= 0
    A[:,i] .= 0
    A[i,i] = 1.0
end
function impose_Neumann_BC!(A,b,i,val)
    b[i] += val
end

Nsteps = ceil(Int,T/dt)
dt = T/Nsteps
u = u0.(x)
interval = 1
@gif for i = 1:Nsteps
    B = (M + .5*dt*A) # Crank-Nicolson matrix
    b = M*u - .5*dt*A*u # Crank-Nicolson RHS for f(x,t) = 0
    tk = (i - 1/2)*dt # in between tk, tk+1
    impose_Dirichlet_BC!(B,b,1,α(tk))
    impose_Dirichlet_BC!(B,b,m+2,β(tk))
    u .= B\b
    if i%interval==0
        println("on step $i out of $Nsteps")
        plot(x,u,legend=false,ylims=(-2,3))
    end
end every interval


text:
plt = Observable{Any}(plot_signal(raw_ismrmrd; darkmode))
ui = dom"div"(plt)
map!(p->plot_signal(p; darkmode), plt, sig_obs)
content!(w, "div#content", ui)


text:
function _separator()
    if Sys.iswindows()
        return ';'
    end
    return ':'
end


text:
using Revise
using NNFEM
using ADCME
using ADCMEKit
using PyPlot
using ProgressMeter

m = 40
n = 20
h = 1/n
NT = 100
Δt = 1/NT
tau = 0.5
c = 1.0

elements = []
prop = Dict("name"=> "PlaneStrainViscoelasticityProny", "rho"=> 2000.0, "E"=> 37844812616.2842, "nu"=> 0.307699122884735, "tau"=>tau, "c"=>c)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 2
        nodes = [
        (i-1)*h (j-1)*h
        i*h (j-1)*h
        i*h j*h
        (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)

for i = 1:m+1
    idx = n*(m+1)+i
    EBC[idx,:] .= -1
end

# for j = 1:n
#     idx = (j-1)*(m+1) + m + 1
#     FBC[idx, :] .= -1
#     f[idx,:] = [-1.0;0.0]
# end

Edge_traction_data = []
IdDict = Dict(
    (1,2)=>1,
    (2,3)=>2,
    (3,4)=>3,
    (1,4)=>4
)
for k = 1:length(elements)
    e = elements[k]
    if sum(e.coords[:,1] .> 2.0 - 1e-3)==2
        i = findall(e.coords[:,1] .> 2.0 - 1e-3)
        push!(Edge_traction_data, [k IdDict[(minimum(i), maximum(i))] 0])
    end
end
Edge_traction_data = Int64.(vcat(Edge_traction_data...))

function Edge_func(x, y, time, id)
    return [-1e6 * ones(length(x)) zeros(length(x))]
end

ndims = 2
domain = Domain(coords, elements, ndims, EBC, g, FBC, f, Edge_traction_data)
globaldata = GlobalData(missing, missing, missing, missing, domain.neqs, nothing, nothing, nothing, Edge_func)

# LinearGeneralizedAlphaSolverInit!(globaldata, domain, Δt)
@showprogress for i = 1:NT
    # global globaldata, domain = LinearGeneralizedAlphaSolverStep(globaldata, domain, Δt)
    global globaldata, domain = GeneralizedAlphaSolverStep(globaldata, domain, Δt)
end

# visualize_von_mises_stress_on_scoped_body(getStateHistory(domain), domain; scale_factor=9.0e7, vmin=0.012, vmax=27.992)
figure()
visualize_von_mises_stress_on_scoped_body(getStateHistory(domain)[end:end,:], domain;
                frames=1, scale_factor=8.3e2, vmin=2339, vmax=3.56e6)
title("von Mises Stress")
savefig("prony_vm.png")
figure()
visualize_total_deformation_on_scoped_body(getStateHistory(domain)[end:end,:], domain;
                frames=1)
title("Total Displacement")
savefig("prony_td.png")

figure()
visualize_x_deformation_on_scoped_body(getStateHistory(domain)[end:end,:], domain;
                frames=1, scale_factor=8.3e2)
title("x-directional Displacement")
savefig("prony_x.png")

figure()
visualize_y_deformation_on_scoped_body(getStateHistory(domain)[end:end,:], domain;
                frames=1, scale_factor=8.3e2)
title("y-directional Displacement")
savefig("prony_y.png")


text:
#################################################################
## expression to graph testing
#################################################################

### for loop #1
ex = quote
    a=zeros(10)
    for i in 1:10
        t = x+z
        a[i] = b[i]+t
    end
end

g = m.tograph(ex);
@test sort(collect(values(g.exti.kv))) == [:b, :x, :z]
@test sort(collect(values(g.seti.kv))) == [:a]
@test length(g.seto.kv) == 0
@test length(g.exto.kv) == 0

exout = striplinenumbers(quote
    _tmp1 = zeros(10)
    for i = 1:10
        t = x + z
        _tmp1[i] = b[i] + t
    end
    a = _tmp1
end)
m.resetvar()
@test m.tocode(g) == exout

### for loop #2
ex = quote
    a=zeros(10)
    z = 12
    for i in 1:10
        t = x+z
        for j in 1:10
            u = t+z+v
            a[i] = b[i]+u
        end
    end
end
g = m.tograph(ex)

@test sort(collect(values(g.exti.kv))) == [:b, :v, :x]
@test sort(collect(values(g.seti.kv))) == [:a, :z]
@test length(g.seto.kv) == 0
@test length(g.exto.kv) == 0

g2 = g.nodes[7].main[2]  # first level loop
@test sort(collect(values(g2.exti.kv))) == [:a, :b, :i, :v, :x, :z]
@test sort(collect(values(g2.seti.kv))) == [:a, :t]
@test sort(collect(values(g2.exto.kv))) == [:a, :b, :v, :x, :z]
@test sort(collect(values(g2.seto.kv))) == [:a]

g3 = g2.nodes[7].main[2] # second level loop
@test sort(collect(values(g3.exti.kv))) == [:a, :b, :i, :t, :v, :z]
@test sort(collect(values(g3.seti.kv))) == [:a, :u]
@test sort(collect(values(g3.exto.kv))) == [:a, :b, :i, :t, :v, :z]
@test sort(collect(values(g3.seto.kv))) == [:a]

m.resetvar()
exout = striplinenumbers(quote
        z = 12
        _tmp1 = zeros(10)
        for i = 1:10
            t = x + z
            for j = 1:10
                u = t + z + v
                _tmp1[i] = b[i] + u
            end
        end
        a = _tmp1
    end)
exout2 = striplinenumbers(quote  # alternate possible result
        _tmp1 = zeros(10)
        z = 12
        for i = 1:10
            t = x + z
            for j = 1:10
                u = t + z + v
                _tmp1[i] = b[i] + u
            end
        end
        a = _tmp1
    end)
res = m.tocode(g)
@test (res == exout) | (res == exout2)

#################################################################
## full cycle  : tograph -> splitnary -> simplify -> prune -> tocode
#################################################################
function fullcycle(ex) # ex = :( y = a[2] ; x = 2 ; y )
    g = m.tograph(ex)

    length(g.seti.kv) == 0 && error("nothing defined here")

    # isolate a single variable of interest
    if haskey( g.seti.vk, nothing) # last statement has priority
        lastnode = g.seti.vk[nothing]
    else # find last evaluated otherwise
        m.evalsort!(g)
        nvars = collect(keys(g.seti))
        lastnode = g.nodes[ maximum(indexin(nvars, g.nodes)) ]
    end
    # lastnode = (g.nodes[1], :out)
    # for n in reverse(g.nodes)
    #     haskey(g.seti, n) || continue
    #     lastnode = n
    #     break
    # end

    sym = g.seti[lastnode]
    # g.seti = m.BiDict{m.ExNode,Any}([lastnode], [ sym==nothing ? :out : sym ])
    g.seti = m.BiDict{m.ExNode,Any}([lastnode], [ sym ])

    g |> m.splitnary! |> m.simplify! |> m.prune!
    m.simplify!(g)
    m.prune!(g)

    m.resetvar()
    m.tocode(g)
end

@test fullcycle(:(a = b+6 ))        == Expr(:block, :(a = b+6) )
@test fullcycle(:(sin(y);a=3))      == Expr(:block, :(a = 3) )
@test fullcycle(:(a += b+6))        == Expr(:block, :(a = a + (b+6)) )
@test fullcycle(:(a -= b+6))        == Expr(:block, :(a = a - (b+6)) )
@test fullcycle(:(a *= b+6))        == Expr(:block, :(a = a * (b+6)) )
@test fullcycle(:(a = b'))          == Expr(:block, :(a = b') )
@test fullcycle(:(a = [1,2]))       == Expr(:block, :(a = [1,2]) )
@test fullcycle(:(a = 4:5 ))        == Expr(:block, :(a = 4:5) )

@test fullcycle(:(a = b+4+5))       == Expr(:block, :(a = b+9) )
@test fullcycle(:(a = b+0))         == Expr(:block, :(a = b) )
@test fullcycle(:(a = b*0))         == Expr(:block, :(a = 0) )
@test fullcycle(:(a = b*1))         == Expr(:block, :(a = b) )
@test fullcycle(:(a = b*(0.5+0.5))) == Expr(:block, :(a = b) )
@test fullcycle(:(a = b/1))         == Expr(:block, :(a = b) )

@test fullcycle(:(5))                          == Expr(:block, :(5) )
@test fullcycle(:(a = 2 ; b = 2 ; a:b))        == Expr(:block, :(2:2) )
@test fullcycle(:(a = 2 ; a))                  == Expr(:block, :(2) )
@test fullcycle(:(a = x ; b = a ))             == Expr(:block, :( b = x ) )
@test fullcycle(:(a = x ; b = a ; a + b))      == Expr(:block, :( x+x ) )
@test fullcycle(:(a = x ; b = a ; c = b ; b))  == Expr(:block, :( x ) )

@test fullcycle(:( a[2] ))                        == Expr(:block, :( a[2]) )
@test fullcycle(:( y = a[2] ; y ))                == Expr(:block, :( a[2]) )
@test fullcycle(:( y = a[2] ; y[1] ))             == Expr(:block, :( a[2][1]) )
@test fullcycle(:( y[1] = a[2] ; y[1] ))          == :(y[1] = a[2]; y[1])
@test fullcycle(:( y = a+1 ; y[2]+y[1] ))         == :(_tmp1 = a+1 ; _tmp1[2]+_tmp1[1])
@test fullcycle(:( a[2] = x ; a[3] ))             == :(a[2] = x ; a[3])
@test fullcycle(:( a[2] = x ; y=a[3] ; y ))       == :(a[2] = x ; a[3])
@test fullcycle(:( b = a ; b[2] = x; 1 + b[2] ))  == :(a[2] = x ; 1+a[2])
@test fullcycle(:( b = a ; b[2] = x; 1 + b[1] ))  == :(a[2] = x ; 1+a[1])
@test fullcycle(:( a[1] + a[2] ))                 == Expr(:block, :( a[1] + a[2]) )
@test fullcycle(:( a[1:2] ))                      == Expr(:block, :( a[1:2]) )
@test fullcycle(:( a = x ; a[1:2] = a[1:2] ))     == Expr(:block, :( a = x ) )
@test fullcycle(:( a = x ; a[1:2,3] = a[1:2,3] )) == Expr(:block, :( a = x ) )
@test fullcycle(:( a = x ; a[1:2,3] = a[1:2,4] )) == :( _tmp1 = 1:2 ; x[_tmp1,3] = x[_tmp1,4])

@test fullcycle(:( x[:] ))               == Expr(:block, :( x[1:length(x)] ) )
@test fullcycle(:( x[a+b, c:d] ))        == Expr(:block, :( x[a + b,c:d] ) )
@test fullcycle(:( x[1:4] ))             == Expr(:block, :( x[1:4] ) )
@test fullcycle(:( x[1:end] ))           == Expr(:block, :( x[1:length(x)]) )
@test fullcycle(:( a[1:end, :, 10:15] )) == Expr(:block, :( a[1:size(a,1),1:size(a,2),10:15]) )

@test fullcycle(:( a.x ))                     == Expr(:block, :( a.x) )
@test fullcycle(:( y = a.x ))                 == Expr(:block, :(y = a.x) )
@test fullcycle(:( y = a.x + 1 ; y.b + y.c )) == :(_tmp1 = +(a.x,1) ; _tmp1.b+_tmp1.c)
@test fullcycle(:( a.x = x ; a[3] ))          == :(a.x = x; a[3])
@test fullcycle(:( a.x = x ; y = a.y ; y ))   == :(a.x = x ; a.y )
@test fullcycle(:( b = a; b.x = x; 1 + b.y )) == :(a.x = x ; 1+a.y )
@test fullcycle(:( a.x + a.y ))               == Expr(:block, :( a.x+a.y ) )

@test fullcycle(:( a = b.f[i]))        == Expr(:block, :(a = b.f[i]) )
@test fullcycle(:( a = b[j].f[i]))     == Expr(:block, :(a = b[j].f[i]) )

ex = quote
    for i in 1:10
        a[i] = b[i]+2
    end
    c=3
end
exout = Expr(:block, :( c=3 ))
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a=0
    for i in 1:10
        a += x
    end
end
exout = quote
    _tmp1=0
    for i in 1:10
        _tmp1 = _tmp1+x
    end
    a=_tmp1
end
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a=zeros(10)
    for i in 1:10
        a[i] = b[i]+2
    end
end
exout = quote
    _tmp1 = zeros(10)
    for i = 1:10
        _tmp1[i] = b[i] + 2
    end
    a = _tmp1
end
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a=zeros(10)
    z=sum(a)
    for i in 1:10
        a[i] = b[i]+2
    end
end
exout = quote
    _tmp1 = zeros(10)
    for i = 1:10
        _tmp1[i] = b[i] + 2
    end
    a = _tmp1
end
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a=zeros(10)
    for i in 1:10
        a[i] = b[i]+2
    end
    z=sum(a)
end
exout = quote
    _tmp1 = zeros(10)
    for i = 1:10
        _tmp1[i] = b[i] + 2
    end
    z=sum(_tmp1)
end
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a=zeros(10+6)
    for i in 1:10
        t = 4+3+2
        a[i] += b[i]+t
    end
    z=sum(a)
end
exout = quote
    _tmp1 = zeros(16)
    for i = 1:10
        _tmp1[i] = _tmp1[i] + (b[i] + 9)
    end
    z=sum(_tmp1)
end
@test fullcycle(ex) == striplinenumbers(exout)

###  test evalconstants, simplify
ex = quote
    y = x * a * 1
    y2 = (x * a) + 0 + 3
    x += 1
    y3 = x * a
    y + y2 + y3 + 12
end
exout = :( _tmp1 = *(x,a) ; +(_tmp1,+(+(_tmp1,3),+(*(+(x,1),a),12))) )
@test fullcycle(ex) == exout

###  test respect of allocations
ex = quote
    a=zeros(2)
    a[2] = x
    sum(a)
end
exout = quote
    _tmp1 = zeros(2)
    _tmp1[2] = x
    sum(_tmp1)
end
@test fullcycle(ex) == striplinenumbers(exout)

ex = quote
    a = zeros(5)
    x = sum(a)
    a[2] = 1
    y = sum(a)
    x + y
end
exout = quote
    _tmp1 = zeros(5)
    _tmp2 = sum(_tmp1)
    _tmp1[2] = 1
    _tmp2 + sum(_tmp1)
end

@test fullcycle(ex) == striplinenumbers(exout)


text:
############################################################
# linked list; last elt points to itself
############################################################

mutable struct Orbit5{T, S} <: AbstractOrbit{T,S}
    pt::T # starting point
    last::T
    elts::Dict{T, Tuple{S, T}} # idx -> (content, nextptr)
end

function Orbit5(pt::T, val::S) where {T,S}
    return Orbit5(pt, pt, Dict{T, Tuple{S, T}}(pt=>(val,pt)))
end

Orbit5(pt::T) where {T} = Orbit5(pt, nothing)

@inline Base.in(x, orb::Orbit5) = haskey(orb.elts, x)

@inline function Base.push!(orb::Orbit5{T,S}, tup::Tuple{T,S}) where {T,S}
    new, val = tup
    orb.elts[orb.last] = (orb.elts[orb.last][1], new)
    orb.last = new
    orb.elts[new] = (val, new)
    return orb
end

@inline Base.push!(orb::Orbit5{T,Nothing}, n::T) where T = push!(orb, (n,nothing))

# o is a point in orbit, return (value, nextptr)
@inline Base.getindex(orb::Orbit5, o) = orb.elts[o][1]

@inline Base.iterate(orb::Orbit5) = orb.pt, orb.pt

@inline islast(orb::Orbit5, s) = s == orb.last
@inline next(orb::Orbit5, s) = (last(orb.elts[s]), last(orb.elts[s]))

# s is the previous elt
@inline Base.iterate(orb::Orbit5, s) = (islast(orb, s) ? nothing : next(orb, s))

@inline Base.length(orb::Orbit5) = length(orb.elts)

@inline function Base.:(==)(o1::Orbit5, o2::Orbit5)
    o1.pt == o2.pt && o1.last == o2.last || return false
    return o1.elts == o2.elts
end

@inline Base.first(orb::Orbit5) = orb.pt


text:
###########################
#### Surface roughness ####
###########################

#' Roughness Reynolds Number
#'
#' calculates the Roughness Reynolds Number.
#'
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - ustar     Friction velocity (m s-1)
#' - z0m       Roughness length (m)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  pressure0 - reference atmospheric pressure at sea level (Pa)
#'                  Tair0 - reference air temperature (K)
#'
#' # Details
 The Roughness Reynolds Number is calculated as in Massman 1999a:
#'
#'            ``Re = z0m * ustar / v``
#'
#'          where `v` is the kinematic viscosity (m2 s-1).
#'
#' # Value
 - Re -: Roughness Reynolds Number (-)
#'
#' #References
#' Massman, W_J., 1999a: A model study of kB H- 1 for vegetated surfaces using
#'            'localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.
#'
#' ```@example; output = false
#' ```
#' Reynolds_Number(25,100,0.5,z0m=0.5)
#'
"""
"""
function Reynolds_Number(Tair,pressure,ustar,z0m,constants=BigleafConstants())

  v  = kinematic_viscosity(Tair,pressure,constants)
  Re = z0m*ustar/v

  return(Re)
end

#' Roughness Parameters
#'
#' A simple approximation of the two roughness parameters displacement height (d)
#'              and roughness length for momentum (z0m).
#'
#' - method    Method to use, one of `Val(:canopy_height),Val(:canopy_height_LAI),Val(:wind_profile)`
#'                  NOTE: if `method = Val(:canopy_height)`, only the following three arguments
#'                  are used. If `method = Val(:canopy_height_LAI)`, only `zh, LAI, cd`,
#'                  and `hs` are required.
#' - zh        Vegetation height (m)
#' - frac_d    Fraction of displacement height on canopy height (-)
#' - frac_{z0m}  Fraction of roughness length on canopy height (-)
#' - LAI       Leaf area index (-)
#' - zr        Instrument (reference) height (m)
#' - cd        Mean drag coefficient for individual leaves. Defaults to 0.2.
#'                  Only needed if `method = Val(:canopy_height_LAI)`.
#' - hs        roughness length of the soil surface (m). Only needed if `method = Val(:canopy_height_LAI)`
#'                  The following arguments are only needed if `method = Val(:wind_profile)`!
#' - data      DataFrame or matrix containing all required variables
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - wind      Wind speed at height zr (m s-1)
#' - ustar     Friction velocity (m s-1)
#' - H         Sensible heat flux (W m-2)
#' - d         Zero-plane displacement height (m); optional
#' - z0m       Roughness length for momentum (m); optional
#' - stab_roughness   Should stability correction be considered? Default is `true`.
#' - stab_formulation Stability correction function used (If `stab_correction = true`).
#'                         Either `Val(:Dyer_1970)` or `Val(:Businger_1971)`.
#' - constants k - von-Karman constant (-)
#'                  Kelvin - conversion degree Celsius to Kelvin
#'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  g - gravitational acceleration (m s-2)
#'                  se_median - conversion standard error (SE) of the mean to SE of the median
#'
#'
#' # Details
 The two main roughness parameters, the displacement height (d)
#'          and the roughness length for momentum (z0m) can be estimated from simple
#'          empirical relationships with canopy height (zh). If `method = Val(:canopy_height)`,
#'          the following formulas are used:
#'
#'            ``d = frac_d * zh``
#'
#'            ``z0m = frac_{z0m} * zh``
#'
#'          where frac_d defaults to 0.7 and frac_{z0m} to 0.1.
#'
#'          Alternatively, d and z0m can be estimated from both canopy height and LAI
#'          (If `method = Val(:canopy_height_LAI)`).
#'          Based on data from Shaw & Pereira 1982, Choudhury & Monteith 1988 proposed
#'          the following semi-empirical relations:
#'
#'            ``X = cd * LAI``
#'
#'            ``d = 1.1 * zh * ln(1 + X^(1/4))``
#'
#'            ``z0m = hs + 0.3 * zh * X^(1/2)   for 0 <= X <= 0.2``
#'
#'            ``z0m = hs * zh * (1 - d/zh)   for 0.2 < X``
#'
#'          If `method = Val(:wind_profile)`, z0m is estimated by solving
#'          the wind speed profile for z0m:
#'
#'            ``z0m = median((zr - d) * exp(-k*wind / ustar - psi_m)``
#'
#'          By default, d in this equation is fixed to 0.7*zh, but can be set to any
#'          other value. psi_m is 0 if `stab_roughness = false`.
#'
#' # Value
 a DataFrame with the following columns:
#'         - d: Zero-plane displacement height (m)
#'         - z0m: Roughness length for momentum (m)
#'         - z0m_se: Only if `method = wind_profile`: Standard Error of the median for z0m (m)
#'
#' #References
#' Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat
#'             budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
#'
#'             Shaw, R. H., Pereira, A., 1982: Aerodynamic roughness of a plant canopy:
#'             a numerical experiment. Agricultural Meteorology, 26, 51-65.
#'
#' #See also
#' [`wind_profile`](@ref)
#'
#' ```@example; output = false
#' ```
#' # estimate d and z0m from canopy height for a dense (LAI=5) and open (LAI=2) canopy
#' roughness_parameters(method=Val(:canopy_height_LAI),zh=25,LAI=5)
#' roughness_parameters(method=Val(:canopy_height_LAI),zh=25,LAI=2)
#'
#' # fix d to 0.7*zh and estimate z0m from the wind profile
#' df = DataFrame(Tair=c(25,25,25),pressure=100,wind=c(3,4,5),ustar=c(0.5,0.6,0.65),H=200)
#' roughness_parameters(method=Val(:wind_profile),zh=25,zr=40,frac_d=0.7,data=df)
#'
#' # assume d = 0.8*zh
#' roughness_parameters(method=Val(:wind_profile),zh=25,zr=40,frac_d=0.8,data=df)
#'
#' @importFrom stats median sd complete_cases
#' @export
function roughness_parameters(method=c(Val(:canopy_height),Val(:canopy_height_LAI),Val(:wind_profile)),zh,
                                 frac_d=0.7,frac_{z0m}=0.1,LAI,zr,cd=0.2,hs=0.01,data,Tair="Tair",pressure="pressure",
                                 wind="wind",ustar="ustar",H="H",d=nothing,z0m=nothing,
                                 stab_roughness=true,stab_formulation=c(Val(:Dyer_1970),Val(:Businger_1971)),
                                 constants=BigleafConstants())

  method           = match_arg(method)
  stab_formulation = match_arg(stab_formulation)

  if (method == Val(:canopy_height))

    d      = frac_d*zh
    z0m    = frac_{z0m}*zh
    z0m_se = missing

elseif (method == Val(:canopy_height_LAI))

    X = cd * LAI
    d = 1.1 * zh * log(1 + X^(1/4))

    if (X >= 0 & X <= 0.2)
      z0m = hs + 0.3 * X^(1/2)
else
      z0m = 0.3 * zh * (1 - d/zh)
end
    z0m_se = missing

elseif (method == Val(:wind_profile))

    check_input(data,Tair,pressure,wind,ustar,H)

    if (isnothing(d))

      d = frac_d * zh

end

    if (stab_roughness)

      zeta  = stability_parameter(data=data,Tair=Tair,pressure=pressure,ustar=ustar,H=H,
                                   zr=zr,d=d,constants=constants)
      psi_m = stability_correction(zeta,formulation=stab_formulation)[,"psi_m"]
      z0m_all = (zr - d) * exp(-constants.k*wind / ustar - psi_m)

else

      z0m_all = (zr - d) * exp(-constants.k*wind / ustar)

end

    z0m_all[z0m_all > zh] = missing

    z0m    = median(z0m_all,na_rm=true)
    z0m_se = constants.se_median] * (sd(z0m_all,na_rm=true) / sqrt(length(z0m_all[complete_cases(z0m_all))))

end

  return(DataFrame(d,z0m,z0m_se))
end

#' Wind Speed at Given Heights in the Surface Layer
#'
#' Wind speed at a given height above the canopy estimated from single-level
#'              measurements of wind speed.
#'
#' - data      DataFrame or matrix containing all required variables
#' - z         Heights above ground for which wind speed is calculated.
#'                  Needs to be of same length as `data` or of length 1
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - ustar     Friction velocity (m s-1)
#' - H         Sensible heat flux (W m-2)
#' - wind      Wind speed at height zr (m s-1); only used if `stab_correction = true`
#' - zr        Instrument (reference) height (m)
#' - zh        Canopy height (m)
#' - d         Zero-plane displacement height (-)
#' - frac_d    Fraction of displacement height on canopy height (-);
#'                  only used if `d` is not available
#' - z0m       Roughness length (m), optional; only used if `stab_correction = false` (default=0.1)
#' - frac_{z0m}  Fraction of roughness length on canopy height (-), optional; only used if `z0m` is not provided.
#'                  Default is 0.1.
#' - estimate_z0m Should `z0m` be estimated from the logarithmic wind profile? If `true` (the default),
#'                     arguments `z0m` and `frac_{z0m}` are ignored.
#'                     See [`roughness_parameters`](@ref) for details.
#' - stab_correction Should stability correction be applied? Defaults to `true`
#' - stab_formulation Stability correction function used (If `stab_correction = true`).
#'                         Either `Val(:Dyer_1970)` or `Val(:Businger_1971)`.
#' - constants k - von-Karman constant (-)
#'                  Kelvin - conversion degree Celsius to Kelvin
#'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  g - gravitational acceleration (m s-2)
#'
#' # Details
 The underlying assumption is the existence of a logarithmic wind profile
#'          above the height d + z0m (the height at which wind speed mathematically reaches zero
#'          according to the Monin-Obhukov similarity theory).
#'          In this case, the wind speed at a given height z is given by:
#'
#'            ``u(z) = (ustar/k) * (ln((z - d) / z0m) - \\psi_m``
#'
#'          The roughness parameters zero-plane displacement height (d) and roughness length (z0m)
#'          can be approximated from [`roughness_parameters`](@ref). ``\\psi_m`` is omitted
#'          if `stab_correction = false` (not recommended). If `estimate_z0m = true`,
#'          z0m is first estimated from the wind profile equation and then used in the equation
#'          above for the calculation of `u(z)` (see e.g. Newman & Klein 2014).
#'
#' # Note
#' Note that this equation is only valid for z >= d + z0m, and it is not
#'       meaningful to calculate values closely above d + z0m. All values in `heights`
#'       smaller than d + z0m will return 0.
#'
#' # Value
 A vector of wind speed at heights `z`.
#'
#' #References
#' Monteith, J_L., Unsworth, M_H., 2008: Principles of Environmental Physics.
#'             3rd edition. Academic Press, London.
#'
#'             Newman, J_F., Klein, P_M., 2014: The impacts of atmospheric stability on
#'             the accuracy of wind speed extrapolation methods. Resources 3, 81-105.
#'
#' #See also
#' [`roughness_parameters`](@ref)
#'
#' ```@example; output = false
#' ```
#' heights = seq(18,40,2)  # heights above ground for which to calculate wind speed
#' df = DataFrame(Tair=25,pressure=100,wind=c(3,4,5),ustar=c(0.5,0.6,0.65),H=c(200,230,250))
#' ws = DataFrame(matrix(missing,ncol=length(heights),nrow=nrow(df)))
#' colnames(ws) = paste0(heights,"m")
#' for (i in seq_along(heights))
#'   ws[,i] = wind_profile(df,z=heights[i],zr=40,zh=25,d=16)
#' }
#'
#' @export
function wind_profile(data,z,Tair="Tair",pressure="pressure",ustar="ustar",H="H",wind="wind",
                         zr,zh,d=nothing,frac_d=0.7,z0m=nothing,frac_{z0m}=nothing,estimate_z0m=true,
                         stab_correction=true,stab_formulation=c(Val(:Dyer_1970),Val(:Businger_1971)),
                         constants=BigleafConstants())

  stab_formulation = match_arg(stab_formulation)

  check_input(data,ustar)

  ## determine roughness parameters
  if (isnothing(d))
    if (isnothing(frac_d))
      stop("Either 'd' or 'frac_d' must be specified")
end
    d = frac_d * zh
end

  if (isnothing(z0m) & !estimate_z0m)
    if (isnothing(frac_{z0m}))
      stop("Either 'z0m' or 'frac_{z0m}' must be specified if 'estimate_z0m' = false")
end
    z0m = frac_{z0m} * zh
end

  if (estimate_z0m)

    if (!isnothing(z0m) | !isnothing(frac_{z0m}))
      cat("Note that arguments 'z0m' and 'frac_{z0m}' are ignored if 'estimate_z0m' = true. z0m is
           calculated from the logarithmic wind_profile equation.",fill=true)
end

    check_input(data,Tair,pressure,wind,ustar,H)

    z0m = roughness_parameters(method=Val(:wind_profile),zh=zh,zr=zr,d=d,data=data,
                                Tair=Tair,pressure=pressure,wind=wind,ustar=ustar,H=H,
                                stab_roughness=true,stab_formulation=stab_formulation,
                                constants=constants)[,"z0m"]
end

  if ( any(z < (d + z0m) & !ismissing(d + z0m)) )
    @warn"function is only valid for heights above d + z0m! Wind speed for heights below d + z0m will return 0!")
end

  ## calculate wind speeds at given heights z
  if (stab_correction)

    zeta  = stability_parameter(data=data,Tair=Tair,pressure=pressure,ustar=ustar,H=H,
                                 zr=z,d=d,constants=constants)
    psi_m = stability_correction(zeta,formulation=stab_formulation)[,"psi_m"]
    wind_heights = pmax(0,(ustar / constants.k) * (log(pmax(0,(z - d)) / z0m) - psi_m))

else

    wind_heights = pmax(0,(ustar / constants.k) * (log(pmax(0,(z - d)) / z0m)))

end

  return(wind_heights)
end


text:
# Pairwise-Alignment Result
# =========================
#
# Result of pairwise alignment.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Result of pairwise alignment
"""
type PairwiseAlignmentResult{T,S1,S2}
    # alignment score/distance
    value::T
    isscore::Bool
    aln::Nullable{PairwiseAlignment{S1,S2}}
end

function PairwiseAlignmentResult(value, isscore, seq, ref)
    return PairwiseAlignmentResult(value, isscore,
                                   Nullable(PairwiseAlignment(seq, ref)))
end

@compat function (::Type{PairwiseAlignmentResult{S1,S2}}){S1,S2}(value, isscore)
    return PairwiseAlignmentResult(value, isscore,
                                   Nullable{PairwiseAlignment{S1,S2}}())
end

# Accessors
# ---------

"""
    score(alignment_result)

Return score of alignment.
"""
score(aln::PairwiseAlignmentResult) = aln.value

"""
    distance(alignment_result)

Retrun distance of alignment.
"""
distance(aln::PairwiseAlignmentResult) = aln.value

"""
    hasalignment(alignment_result)

Check if alignment is stored or not.
"""
hasalignment(aln::PairwiseAlignmentResult) = !isnull(aln.aln)

"""
    alignment(alignment_result)

Return alignment if any.

See also: `hasalignment`
"""
function alignment(aln::PairwiseAlignmentResult)
    if !hasalignment(aln)
        throw(ArgumentError("alignment is not stored"))
    end
    return get(aln.aln)
end

# Printer
# -------

function Base.show{T,S1,S2}(io::IO, aln::PairwiseAlignmentResult{T,S1,S2})
    println(io, summary(aln), ':')
    if aln.isscore
        print(io, "  score: ", aln.value)
    else
        print(io, "  distance: ", aln.value)
    end
    if hasalignment(aln)
        println(io)
        print(io, alignment(aln))
    end
end


text:
abstract Move

# note that we don't have to keep track of which player made a move on
# the move itself, since we can always figure out how to take the last
# move back from the global game state: whoever's turn it *isn't* made
# the move we're trying to take back

immutable WallPlacement <: Move
    wall::Wall
end

show(io::IO, wp::WallPlacement) = write(io, "WallPlacement$(wp.wall)")

immutable PlayerMove <: Move
    from::Loc
    to::Loc
end

show(io::IO, pm::PlayerMove) = write(io, "PlayerMove($(pm.from)->$(pm.to))")

isoutofbounds(pm::PlayerMove) = isoutofbounds(pm.to)

blocks(wall::Wall, pm::PlayerMove) = separates(wall, pm.from. pm.to)
blocks(wall::Wall, wp::WallPlacement) = center(wall) == center(move.wall)

function any_blocks(walls::Vector{Wall}, move::Move)
    # TODO it's a shame this has to duplicate any_separates
    for wall in walls
        if blocks(wall, move)
            return true
        end
    end
    return false
end


text:
# Autogenerated wrapper script for jemalloc_jll for aarch64-apple-darwin
export libjemalloc

JLLWrappers.@generate_wrapper_header("jemalloc")
JLLWrappers.@declare_library_product(libjemalloc, "@rpath/libjemalloc.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libjemalloc,
        "lib/libjemalloc.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Core type of `NCDatasets`
the `Attributes` and `Group` parts as well.
and the `Attributes` part of it.

High-level interface is at the "High-level" section about actually
loading/reading/making datasets.
=#

# Exception type for error thrown by the NetCDF library
mutable struct NetCDFError <: Exception
    code::Cint
    msg::String
end

# base type of attributes list
# concrete types are Attributes (single NetCDF file) and
# MFAttributes (multiple NetCDF files)

abstract type BaseAttributes
end

abstract type AbstractDataset
end

abstract type AbstractDimensions
end

abstract type AbstractGroups
end

############################################################
# Types and subtypes
############################################################

# List of attributes (for a single NetCDF file)
# all ids should be Cint

mutable struct Attributes{TDS<:AbstractDataset} <: BaseAttributes
    ds::TDS
    varid::Cint
end

mutable struct Groups{TDS<:AbstractDataset} <: AbstractGroups
    ds::TDS
end

mutable struct Dimensions{TDS<:AbstractDataset} <: AbstractDimensions
    ds::TDS
end

abstract type AbstractVariable{T,N} <: AbstractArray{T,N} end

# Variable (as stored in NetCDF file, without using
# add_offset, scale_factor and _FillValue)
mutable struct Variable{NetCDFType,N,TDS<:AbstractDataset} <: AbstractVariable{NetCDFType, N}
    ds::TDS
    varid::Cint
    dimids::NTuple{N,Cint}
    attrib::Attributes{TDS}
end

mutable struct NCDataset{TDS} <: AbstractDataset where TDS <: Union{AbstractDataset,Nothing}
    # parent_dataset is nothing for the root dataset
    parentdataset::TDS
    ncid::Cint
    iswritable::Bool
    # true of the NetCDF is in define mode (i.e. metadata can be added, but not data)
    # need to be a reference, so that remains syncronised when copied
    isdefmode::Ref{Bool}
    attrib::Attributes{NCDataset{TDS}}
    dim::Dimensions{NCDataset{TDS}}
    group::Groups{NCDataset{TDS}}
    # mapping between variables related via the bounds attribute
    # It is only used for read-only datasets to improve performance
    _boundsmap::Dict{String,String}
    function NCDataset(ncid::Integer,
                       iswritable::Bool,
                       isdefmode::Ref{Bool};
                       parentdataset = nothing,
                       )

        function _finalize(ds)
            @debug begin
                ccall(:jl_, Cvoid, (Any,), "finalize $ncid $timeid \n")
            end
            # only close open root group
            if (ds.ncid != -1) && (ds.parentdataset == nothing)
                close(ds)
            end
        end
        ds = new{typeof(parentdataset)}()
        ds.parentdataset = parentdataset
        ds.ncid = ncid
        ds.iswritable = iswritable
        ds.isdefmode = isdefmode
        ds.attrib = Attributes(ds,NC_GLOBAL)
        ds.dim = Dimensions(ds)
        ds.group = Groups(ds)
        ds._boundsmap = Dict{String,String}()
        if !iswritable
            initboundsmap!(ds)
        end
        timeid = Dates.now()
        @debug "add finalizer $ncid $(timeid)"
        finalizer(_finalize, ds)
        return ds
    end
end

"Alias to `NCDataset`"
const Dataset = NCDataset

# Variable (with applied transformations following the CF convention)
mutable struct CFVariable{T,N,TV,TA,TSA}  <: AbstractVariable{T, N}
    # this var is generally a `Variable` type
    var::TV
    # Dict-like object for all attributes read from disk
    attrib::TA
    # a named tuple with fill value, scale factor, offset,...
    # immutable for type-stability
    _storage_attrib::TSA
end

# Multi-file related type definitions

mutable struct MFAttributes{T} <: BaseAttributes where T <: BaseAttributes
    as::Vector{T}
end

function Base.getindex(a::MFAttributes,name::AbstractString)
    return a.as[1][name]
end

function Base.setindex!(a::MFAttributes,data,name::AbstractString)
    for a in a.as
        a[name] = data
    end
    return data
end

mutable struct MFVariable{T,N,M,TA,A} <: AbstractVariable{T,N}
    var::CatArrays.CatArray{T,N,M,TA}
    attrib::MFAttributes{A}
    dimnames::NTuple{N,String}
    varname::String
end

mutable struct MFDimensions{T} <: AbstractDimensions where T <: AbstractDimensions
    as::Vector{T}
    aggdim::String
end

mutable struct MFGroups{T} <: AbstractGroups where T <: AbstractGroups
    as::Vector{T}
    aggdim::String
end

mutable struct MFDataset{T,N,S<:AbstractString,TA,TD,TG} <: AbstractDataset where T <: AbstractDataset
    ds::Array{T,N}
    aggdim::S
    attrib::MFAttributes{TA}
    dim::MFDimensions{TD}
    group::MFGroups{TG}
    _boundsmap::Dict{String,String}
end

# DeferDataset are Dataset which are open only when there are accessed and
# closed directly after. This is necessary to work with a large number
# of NetCDF files (e.g. more than 1000).

struct Resource
    filename::String
    mode::String
    metadata::OrderedDict
end

mutable struct DeferAttributes <: BaseAttributes
    r::Resource
    varname::String # "/" for global attributes
    data::OrderedDict
end

mutable struct DeferDimensions <: AbstractDimensions
    r::Resource
    data::OrderedDict
end

mutable struct DeferGroups <: AbstractGroups
    r::Resource
    data::OrderedDict
end

mutable struct DeferDataset <: AbstractDataset
    r::Resource
    groupname::String
    attrib::DeferAttributes
    dim::DeferDimensions
    group::DeferGroups
    data::OrderedDict
    _boundsmap::Union{Nothing,Dict{String,String}}
end

mutable struct DeferVariable{T,N} <: AbstractVariable{T,N}
    r::Resource
    varname::String
    attrib::DeferAttributes
    data::OrderedDict
end


text:
function asinh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && return x
    signbit(x) && return - asinh(-x)
    result = abs(x)
    result = result + sqrt(square(result) + 1.0)
    result = log(result)
    return copysign(result, x)
end

function acosh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && !signbit(x) && return x
    (x < 1.0 || isinf(x)) && throw(DomainError("$x"))
    result = x + sqrt(square(x) - 1.0)
    result = log(result)
    return result
end

function atanh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    (abs(x) > 1.0 || isinf(x)) && throw(DomainError("$x"))
    twox = DoubleFloat{T}(x.hi+x.hi, x.lo+x.lo)
    result = 1.0 + twox / (1.0 - x)
    result = log(result)
    result = DoubleFloat{T}(result.hi*0.5, result.lo*0.5)
    return result
end

function acsch(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && return zero(DoubleFloat{T})
    iszero(x) && return signbit(x) ? -DoubleFloat{T}(Inf) : DoubleFloat{T}(Inf)
    invx = inv(x)
    invx2 = inv(x*x)
    result = log(sqrt(1 + invx2) + invx)
    return result
end

function asech(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    (x < 0.0 || x > 1.0 || isinf(x)) && throw(DomainError("$x"))
    iszero(x) && return inf(DoubleFloat{T})
    isone(x) && return zero(DoubleFloat{T})
    invx = inv(x)
    result = fma(sqrt(invx + 1), sqrt(invx - 1), invx)
    return log(result)
end

function acoth(x::DoubleFloat{T}) where {T<:IEEEFloat}
    isnan(x) && return x
    isinf(x) && return copysign(zero(DoubleFloat{T}), x)
    isone(x) && return inf(DoubleFloat{T})
    isone(-x) && return -inf(DoubleFloat{T})
    (-1.0 < x < 1.0) && throw(DomainError("$x"))
    invx = inv(x)
    result = (log(1+invx) - log(1-invx))/2
    return result
end


text:
using ClearswiApp
@time msg = clearswi_main(ARGS)
println(msg)


text:
using Documenter
using POMDPSimulators

makedocs(
    modules = [POMDPSimulators],
    format = Documenter.HTML(),
    sitename = "POMDPSimulators.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPSimulators.jl.git",
)


text:
# We have one global JIT and TM
const orc = Ref{LLVM.OrcJIT}()
const tm = Ref{LLVM.TargetMachine}()

function __init__()
    CACHE[NativeInterpreter] = CodeCache(cpu_invalidate)

    opt_level = Base.JLOptions().opt_level
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end

    tm[] = LLVM.JITTargetMachine(; optlevel=optlevel)
    LLVM.asm_verbosity!(tm[], true)

    orc[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
    atexit() do
        return LLVM.dispose(orc[])
    end
end


text:
abstract type Node end
abstract type Layer end

mutable struct NodeGene <: Node
    #= A node gene encodes the basic artificial neuron model.
    nodetype should be "INPUT", "HIDDEN", or "OUTPUT" =#
    id::Int64
    ntype::Symbol
    bias::Float64
    activation::Symbol
    ninputs::Int64
    weights::Vector{Float64}
    enable_bit::Bool
    function NodeGene(id::Int64, nodetype::Symbol, ninputs::Int64, bias::Float64=0.0,
                      activation::Symbol=:sigm)
        @assert  activation in [:none, :sigm, :tanh, :relu, :identity]
        weights = randn(ninputs)
        new(id, nodetype, bias, activation, ninputs, weights, true)
    end
    function NodeGene(id::Int64, nodetype::Symbol, weights::Vector{Float64},
                      bias::Float64=0.0, activation::Symbol=:sigm)
        @assert  activation in [:none, :sigm, :tanh, :relu, :identity]
        new(id, nodetype, bias, activation, length(weights), weights, true)
    end
end

function Base.show(io::IO, ng::NodeGene)
    @printf(io, "Node %2d -> type=%6s, bias=%+2.10s, inputs=%+2.10s, %s()",
            ng.id, ng.ntype, ng.bias, ng.ninputs, ng.activation)
end

get_weights(ng::NodeGene, enabled_weights::Vector{Bool}) = ng.weights[enabled_weights]
get_bias(ng::NodeGene) = ng.bias
is_same_innov(n1::NodeGene, n2::NodeGene) = n1.id == n2.id

function get_child(ng::NodeGene, other::NodeGene, fittest_parent::Int64)
    # Creates a new NodeGene ramdonly inheriting its attributes from parents
    assert(ng.id == other.id)
    if fittest_parent == 1
        ng = NodeGene(ng.id, ng.ntype, ng.ninputs, ng.bias, ng.activation)
    elseif fittest_parent == 2
        ng = NodeGene(other.id, other.ntype, other.ninputs, other.bias, other.activation)
    else
        ng = NodeGene(ng.id, ng.ntype,
                      randbool() ? ng.bias : other.bias,
                      randbool() ? ng.ninputs : other.ninputs,
                      ng.activation)
    end
    return ng
end

function mutate_enable_bit!(ng::NodeGene)
    ng.enable_bit = xor(ng.enable_bit, true)
end

function mutate_bias!(ng::NodeGene, cf::Config)
    ng.bias += randn() * cf.bias_mutation_power
    #check_weight_value!(ng.bias, cf)
    return
end

function mutate_weights!(ng::NodeGene, cf::Config)
    ng.weights .+= randn() * cf.weight_mutation_power
    #check_weight_value!.(ng.weights, [cf])
    return
end

function check_weight_value!(weight::Float64, cf::Config)
    if weight > cf.max_weight
        weight = cf.max_weight
    elseif weight < cf.min_weight
        weight = cf.min_weight
    end
end

function mutate_node!(ng::NodeGene, cf::Config)
    if rand() < cf.prob_mutatebias mutate_bias!(ng, cf) end
    if rand() < cf.prob_mutate_weight mutate_weights!(ng, cf) end
end

get_new_innov_number!(g::Global) =  g.innov_number += 1
get_layer_id!(g::Global) = g.layerCnt += 1

mutable struct MLP_layer <: Layer
    ntype::Symbol
    position::Int64
    activation::Symbol
    nodes_number::Int64
    nodes::Vector{Node}
    layer_id::Int64
    function MLP_layer(g::Global, ntype::Symbol, position::Int64, activation::Symbol, ninputs::Int64,
                       nodes_number::Int64)
        @assert position >= 1
        nodes = Vector{Node}(undef, nodes_number)
        for i = 1:nodes_number
            nodes[i] = NodeGene(get_new_innov_number!(g), ntype, ninputs, 0.0, activation)
        end
        new(ntype, position, activation, nodes_number, nodes, get_layer_id!(g))
    end
    function MLP_layer(g::Global, ntype::Symbol, position::Int64, activation::Symbol, nodes::Vector{<:Node})
        @assert position >= 1
        new(ntype, position, activation, length(nodes), nodes, get_layer_id!(g))
    end
end
#get_layer_id!(g)
function get_layer_params(l::MLP_layer, previous_enabled_nodes::Vector{Bool}, concatenate::Bool=false)
    number_enabled_nodes = count((n) -> n.enable_bit, l.nodes)
    (get_weight_matrix(l, previous_enabled_nodes, concatenate=concatenate)..., l.activation, number_enabled_nodes)
end

get_enabled_nodes(l::MLP_layer) = filter((n) -> n.enable_bit, l.nodes)
get_not_enabled_nodes(l::MLP_layer) = filter((n) -> !(n.enable_bit), l.nodes)

maxInnov_layer(l::MLP_layer) = map((n) -> n.id , l.nodes) |> maximum

function get_layer_dictionary(layer::MLP_layer)
    dict_keys = [n.id for n in layer.nodes]
    dict_values = map((p) -> (layer.position, p), 1:length(layer.nodes))
    zip(dict_keys, dict_values) |> Dict
end

function get_enabled_bit_position(layer::MLP_layer)
    positions = findall((n) -> n.enable_bit, layer.nodes)
    return rand(positions)
end

function mutate_node_weight!(g::Global, layer::MLP_layer)
    node_pos = get_enabled_bit_position(layer)
    mutate_node!(layer.nodes[node_pos], g.cf)
end

function get_weight_matrix(l::MLP_layer, enabled_weights::Vector{Bool}; concatenate::Bool=false)
    enabled_nodes = get_enabled_nodes(l)
    weights = get_weights.(enabled_nodes, [enabled_weights])
    biases = get_bias.(enabled_nodes)
    weights_mtx = vcat(map((a) -> a', weights)...)
    if !concatenate
        return (weights_mtx, biases)
    else
        return hcat(weights_mtx, biases)
    end
end

function add_node_weight!(l::MLP_layer)
    for node in l.nodes
        push!(node.weights, randn())
        node.ninputs += 1
    end
end

function remove_node_weight!(l::MLP_layer, pos::Int=-1)
    for node in l.nodes
        pos == -1 ? deleteat!(node.weights, length(node.weights)) : deleteat!(node.weights, pos)
        node.ninputs -= 1
    end
end

function swapcols!(X::AbstractMatrix, i::Integer, j::Integer)
    @inbounds for k = 1:size(X,1)
        X[k,i], X[k,j] = X[k,j], X[k,i]
    end
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
    Compute the log probability of a sample `z` under the given DPP.
"""
function Distributions.logpdf(dpp::DeterminantalPointProcess, z::AbstractVector{<:Int})
    L_z_eigvals = eigvals(dpp.L[z, z])
    return sum(log.(L_z_eigvals)) - sum(log.(dpp.Lfact.values .+ 1))
end

"""
    Compute the log probability of a sample `z` under the given k-DPP.
"""
function Distributions.logpdf(kdpp::kDeterminantalPointProcess, z::AbstractArray{<:Int})
    L_z_eigvals = eigvals(kdpp.dpp.L[z, z])
    return sum(log.(L_z_eigvals)) .-
           log(elem_symm_poly(kdpp.dpp.Lfact.values, kdpp.k)[end, end])
end

"""
    Compute the probability of a sample `z` under the given DPP.
"""
function Distributions.pdf(dpp::DeterminantalPointProcess, z::AbstractArray{<:Int})
    return exp(logpdf(dpp, z))
end

"""
    Compute the probability of a sample `z` under the given k-DPP.
"""
function Distributions.pdf(kdpp::kDeterminantalPointProcess, z::AbstractArray{<:Int})
    return exp(logpdf(kdpp, z))
end