text:
"""
This module implements various methods of locally approximating a function, such as multi-linear, simplex
and k-nearest-neighbor approximation. An example use case is for locally approximating value
functions in value iteration
"""
module LocalFunctionApproximation

using GridInterpolations
using NearestNeighbors
import NearestNeighbors: NNTree

using Distances

export
	LocalFunctionApproximator,
	LocalGIFunctionApproximator,
	LocalNNFunctionApproximator,
	n_interpolating_points,
	get_all_interpolating_points,
	get_all_interpolating_values,
	get_interpolating_nbrs_idxs_wts,
	compute_value,
	set_all_interpolating_values,
	finite_horizon_extension

abstract type LocalFunctionApproximator end

"""
	n_interpolating_points(lfa::LocalFunctionApproximator)

Return the number of interpolanting points that the approximator is using
"""
function n_interpolating_points end

"""
	get_all_interpolating_points(lfa::LocalFunctionApproximator)

Return the vector of points (in a specific order) that are used to interpolate
"""
function get_all_interpolating_points end

"""
	get_all_interpolating_values(lfa::LocalFunctionApproximator)

Return the vector of all interpolating values (in the same order as the interpolating points)
"""
function get_all_interpolating_values end

"""
	get_interpolating_nbrs_idxs_wts(lfa::LocalFunctionApproximator, v::AbstractVector)

Return a tuple of (indices, weights) for the interpolants for a specific query v
"""
function get_interpolating_nbrs_idxs_wts end

"""
	compute_value(lfa::LocalFunctionApproximator, v::AbstractVector)

Return the value of the function at some query point v, based on the local function approximator

	compute_value(lfa::LocalFunctionApproximator, v_list::AbstractVector{V}) where V <: AbstractVector{Float64}

Return the value of the function for a list of query points, based on the local function approximator
"""
function compute_value end

"""
	set_all_interpolating_values(lfa::LocalFunctionApproximator, vals::AbstractVector)

Set the values of all interpolating points to the input vector
"""
function set_all_interpolating_values end

"""
	finite_horizon_extension(lfa::LocalFunctionApproximator, hor::Int64)

Extend the LFA appropriately along a new dimension to allow for finite-horizon approximations
"""
function finite_horizon_extension end

include("local_gi_fa.jl")
include("local_nn_fa.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libjulia_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libjulia")
JLLWrappers.@generate_main_file("libjulia", UUID("5ad3ddd2-0711-543a-b040-befd59781bbf"))
end  # module libjulia_jll

ext: jl


text:
struct JacobianCache{CacheType1,CacheType2,CacheType3,fdtype,returntype,inplace}
    x1  :: CacheType1
    fx  :: CacheType2
    fx1 :: CacheType3
end

function JacobianCache(
    x,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x),
    inplace    :: Type{Val{T3}} = Val{true}) where {T1,T2,T3}

    if eltype(x) <: Real && fdtype==Val{:complex}
        x1 = fill(zero(Complex{eltype(x)}), size(x))
        _fx = fill(zero(Complex{eltype(x)}), size(x))
    else
        x1 = similar(x)
        _fx = similar(x)
    end

    if fdtype==Val{:complex}
        _fx1  = nothing
    else
        _fx1 = similar(x)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype,inplace)
end

function JacobianCache(
    x ,
    fx,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(x),
    inplace    :: Type{Val{T3}} = Val{true}) where {T1,T2,T3}

    if eltype(x) <: Real && fdtype==Val{:complex}
        x1 = fill(zero(Complex{eltype(x)}), size(x))
    else
        x1 = similar(x)
    end

    if eltype(fx) <: Real && fdtype==Val{:complex}
        _fx = fill(zero(Complex{eltype(x)}), size(fx))
    else
        _fx = similar(fx)
    end

    if fdtype==Val{:complex}
        _fx1  = nothing
    else
        _fx1 = similar(fx)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype,inplace)
end

function JacobianCache(
    x1 ,
    fx ,
    fx1,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(fx),
    inplace    :: Type{Val{T3}} = Val{true}) where {T1,T2,T3}

    if fdtype==Val{:complex}
        !(returntype<:Real) && fdtype_error(returntype)

        if eltype(fx) <: Real
            _fx = fill(zero(Complex{eltype(x1)}), size(fx))
        else
            _fx = fx
        end
        if eltype(x1) <: Real
            _x1 = fill(zero(Complex{eltype(x1)}), size(x1))
        else
            _x1 = x1
        end
    else
        _x1 = x1
        @assert eltype(fx) == T2
        @assert eltype(fx1) == T2
        _fx = fx
    end
    JacobianCache{typeof(_x1),typeof(_fx),typeof(fx1),fdtype,returntype,inplace}(_x1,_fx,fx1)
end

function finite_difference_jacobian(f, x::AbstractArray{<:Number},
    fdtype     :: Type{T1}=Val{:central},
    returntype :: Type{T2}=eltype(x),
    inplace    :: Type{Val{T3}}=Val{true},
    f_in       :: Union{T2,Nothing}=nothing) where {T1,T2,T3}

    cache = JacobianCache(x,fdtype,returntype,inplace)
    finite_difference_jacobian(f,x,cache,f_in)
end

function finite_difference_jacobian(f,x,cache::JacobianCache,f_in=nothing)
    J = fill(zero(eltype(x)), length(x), length(x))
    finite_difference_jacobian!(J,f,x,cache,f_in)
    J
end

function finite_difference_jacobian!(J::AbstractMatrix{<:Number},
    f,x::AbstractArray{<:Number},
    cache::JacobianCache{T1,T2,T3,fdtype,returntype,inplace},
    f_in::Union{T2,Nothing}=nothing) where {T1,T2,T3,fdtype,returntype,inplace}

    m, n = size(J)
    x1, fx, fx1 = cache.x1, cache.fx, cache.fx1
    copyto!(x1, x)
    vfx = vec(fx)
    if fdtype == Val{:forward}
        vfx1 = vec(fx1)
        epsilon_factor = compute_epsilon_factor(Val{:forward}, eltype(x))
        @inbounds for i ∈ 1:n
            epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
            x1_save = x1[i]
            x1[i] += epsilon
            if inplace == Val{true}
                f(fx1, x1)
                if f_in isa Nothing
                    f(fx, x)
                else
                    vfx = vec(f_in)
                end
                J[:,i] = (vfx1 - vfx) / epsilon
            else
                fx1 .= f(x1)
                if f_in isa Nothing
                    fx .= f(x)
                else
                    vfx = vec(f_in)
                end
                J[:,i] = (vfx1 - vfx) / epsilon
            end
            x1[i] = x1_save
        end
    elseif fdtype == Val{:central}
        vfx1 = vec(fx1)
        epsilon_factor = compute_epsilon_factor(Val{:central}, eltype(x))
        @inbounds for i ∈ 1:n
            epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
            x1_save = x1[i]
            x_save = x[i]
            x1[i] += epsilon
            x[i]  -= epsilon
            if inplace == Val{true}
                f(fx1, x1)
                f(fx, x)
                @. J[:,i] = (vfx1 - vfx) / (2*epsilon)
            else
                fx1 .= f(x1)
                fx .= f(x)
                J[:,i] = (vfx1 - vfx) / (2*epsilon)
            end
            x1[i] = x1_save
            x[i]  = x_save
        end
    elseif fdtype==Val{:complex} && returntype<:Real
        epsilon = eps(eltype(x))
        @inbounds for i ∈ 1:n
            x1_save = x1[i]
            x1[i] += im * epsilon
            if inplace == Val{true}
                f(fx,x1)
                @. J[:,i] = imag(vfx) / epsilon
            else
                fx .= f(x1)
                J[:,i] = imag(vfx) / epsilon
            end
            x1[i] = x1_save
        end
    else
        fdtype_error(returntype)
    end
    J
end


text:
"""
    setup_dump_state_diagnostics(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create the "DumpState" `DiagnosticsGroup` which contains all the
prognostic state variables. These are output on `x`, `y`, `z` or `lat`,
`long`, `level` dimensions of an interpolated grid (`interpol` _must_
be specified) as well as a `time` dimension at the specified `interval`.
"""
function setup_dump_state_diagnostics(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "DumpState",
        Diagnostics.dump_state_init,
        Diagnostics.dump_state_fini,
        Diagnostics.dump_state_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

dump_state_init(dgngrp, currtime) = dump_init(dgngrp, currtime, Prognostic())

function dump_state_collect(dgngrp, currtime)
    interpol = dgngrp.interpol
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q.data)
    bl = dg.balance_law
    mpirank = MPI.Comm_rank(mpicomm)

    istate = similar(Q.data, interpol.Npl, number_states(bl, Prognostic()))
    interpolate_local!(interpol, Q.data, istate)

    if interpol isa InterpolationCubedSphere
        # TODO: get indices here without hard-coding them
        _ρu, _ρv, _ρw = 2, 3, 4
        project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))
    end

    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)

    if mpirank == 0
        statenames = flattenednames(vars_state(bl, Prognostic(), FT))
        varvals = OrderedDict()
        for (vari, varname) in enumerate(statenames)
            varvals[varname] = all_state_data[:, :, :, vari]
        end
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function dump_state_fini(dgngrp, currtime) end


text:
"""
This module represents an explicit, direct sparse 2D finite-difference Poisson solver for heterogeneous media,
i.e. media having spatially varying (space-dependent) medium parameters.
The following functionality is currently available in this module:
* 'Forward' problem: given the source, and the conductivity distribution, solve for the electrical potential ψ:
	 Example problem: solve Aψ=∇⋅j; A=∇⋅(σ(x,z)∇), for ψ.
* 'Inverse' problem: given wave propagation related pore pressure, and mechanical medium properties, calculate the source term:
	 Example problem: solve Aψ=-∇⋅j; A=∇⋅([Q*k/η](x,z)∇B*P), for ∇⋅j;
* The Boundary Value Problem that is currently implemented assumes Neumann boundary conditions at all boundaries.
Developed by:\\\n
Niels Grobbe, Massachusetts Institute of Technology, USA.\\\n
In collaboration with: Aimé Fournier & Laurent Demanet, Massachusetts Institute of Technology, USA.\\\n
Date: October, 2017 \\\n
Contact: ngrobbe@gmail.com
"""
module Poisson

using SparseArrays
using StatsBase

"""
Poisson.solve(field,mpars,solflag) solves the forward or inverse Poisson problem in a heterogeneous medium.
# Function input arguments:
* `field` : either source term used in 'forward' mode to get field, or the field (e.g. pore pressure) used in 'inverse' mode to calculate source-term
* `mpars` : medium parameters for forward or inverse mode
* `solflag` : Solution flag determining whether to use forward (solflag=1) or inverse mode (solflag=-1) of Poisson solver.
"""
function solve(field,mpars,solflag::Int64) # fields: either field or source, mpars=medium pars

if(solflag==-1)
	pp=field; #
elseif(solflag==1)
	src=field;
else
	error("In Poisson.solve: change solflag to solflag=1 (forward problem) or solflag=-1 (inverse problem)")
end

nz=size(mpars,1);
nx=size(mpars,2);
dz= inv(nz - 1);
dx= inv(nx - 1);
z=(0 : nz - 1)*dz;
x=(0:nx-1)*dx;

#--------------------  Build explicit operator matrix A=∇⋅(mpars(x,z)∇) for direct Poisson solver --------------------
# Some notes:
# We build the operator matrix explicitly like this to properly deal with the Neumann boundary conditions at all boundaries.
# The standard Kronecker product approach for building sparse operator matrices seemed to generate some errors at the boundaries.
# Medium properties calculated on staggered grid.
# Directions: z-direction row-wise, x-direction column-wise.

# Allocate non-zero operator entries
k = 5*(nz - 2)*(nx - 2)+ # pentadiagonal operator entry allocation: 5, since using 5-point finite-difference stencil (star-shaped)
    2*(2*(nz - 2) + 2*(nx - 2)+  # bidiagonal non-corner 4-boundary allocation # for boundaries (2 points only)
    4); # bidiagonal corner boundary allocation
# The sparse() function is often a handy way to construct sparse matrices.
# It takes as its input a vector I of row indices, a vector J of column indices, and a vector V of nonzero values.
# Sparse(I,J,V) constructs a sparse matrix such that S[I[k], J[k]] = V[k].
iz=Array{Float64,1};
ix=Array{Float64,1};
iz=collect(1:nz*nx); # same row: row index=1 --> vectors must have same length [k]
ix=collect(1:nz*nx); # column indices
val=zeros(nz*nx); # actual values

global A=sparse(ix,iz,val,nz*nx,nz*nx); # pre-allocate sparse operator matrix A

dz = dz.^2; # for simplicity below
dx = dx.^2;
println("fffffff", dz*dx)

#-------------- Create operator A for the interior points (without boundary)---------------------------------------
for j = 2 : nx - 1 # interior x (column) loop (for amount of internal nodes)
    k = 2 + nz*(j - 1); # initial interior z for this x. To jump from j (new column), need nz*... to proceed to next column vector k
    for i = 2 : nz - 1  # interior z (row) loop:
        A[k,k - nz] = .5*(mpars[i,j] + mpars[i,j - 1])/dx; # moving left (and placing actual staggered 2D grid medium value there)
        A[k,k - 1] = .5*(mpars[i,j] + mpars[i - 1,j])/dz; # moving up the grid
        A[k,k] = -(mpars[i,j] + .5*(mpars[i,j - 1] + mpars[i,j + 1]))/dx -
        (mpars[i,j] + .5*(mpars[i - 1,j] + mpars[i + 1,j]))/dz; # diagonal elements: ∇⋅\medpars(x,z)∇ calculated using 5 mpars points
        A[k,k + 1] = .5*(mpars[i,j] + mpars[i + 1,j])/dz; # moving down the grid
        A[k,k + nz] = .5*(mpars[i,j]+ mpars[i,j + 1])/dx; # moving right on grid (need nz values to go to next column)
        k = k + 1; # to next z
    end # The foregoing A terms are operations by hand ...
end

#------------- Create operator A for the top and bottom boundaries including Neumann conditions (without the corners)-------------

if(solflag==-1) # for 'inverse' problem: source plays no role here, so just allocate zeros with proper size.
	src=zeros(nz,nx);

elseif(solflag==1) # for 'forward' problem
	src = src[:]; # vectorize source (which was function input argument)
end

for j = 2 : nx - 1 # interior x (column) loop:
    k = 1 + nz*(j - 1);# initial boundary z for this x (k = i + m*(j - 1))
    for i = [1 nz]  # % north & south boundary z (row) loop:
        o = k - Int(sign(i - .5*nz)); #% sign function determines correct direction to: inward from north or south boundary # to define direction(from boundary down or up)
        A[k,k] = -A[o,k]; #% build in symmetry ... --> trick to get symmetry, the next line is leading.
        A[k,o] = -A[k,k]; #% cancelation implies equation # Neumann: value of psi = same. The operator does not know value. Just knows has to subtract two values --> minus sign
        src[k] = 0; # 0 to equate boundary values in source term ... (?)
        k = k + nz - 1; #% from north to south boundary
    end
end

# Create A for the east and west boundaries without corner points
for j = [1 nx] # west & east boundary x (column) loop:
    k = 2 + nz*(j - 1); #% initial interior z for this x (k = i + m*(j - 1))
    for i = 2 : nz - 1 # interior z (row) loop:
        o = k - Int(sign(j - .5*nx))*nz; # inward from west or east boundary
        A[k,k] = -A[o,k]; #% build in symmetry ...
        A[k,o] = -A[k,k]; #% cancelation implies equation
        src[k]= 0;  #% 0 to equate boundary values
        k = k + 1; #% to next z
    end
 end

# Create A for the corner points
for j = [1 nx] #% west & east boundary x (column) loop:
    k = 1 + nz*(j - 1); #% start on north boundary
    for i = [1 nz] #% north & south boundary y (row) loop:
       o = k - Int(sign(i - .5*nz))- #% inward from north or south boundary
        Int(sign(j - .5*nx))*nz; # inward from west or east boundary # two terms to move diagonal
        A[k,k] = A[o,o]; #% nonzero, alas asymmetric ...
        A[k,o] = -A[k,k]; #% ... cancelation implies equation
        src[k] = 0; #% 0 to equate boundary values
        k = k + nz - 1; #% from north to south boundary
    end
 end

if(solflag == -1) #for inverse problem: fields + medium --> source term
	pp_vec=pp[:]; # vectorize 2D pore pressure field (according to organization operator A).
	divj=-reshape(A*pp_vec,nz,nx);

    	return divj  # function command --> return these variables, in this case electrical source term ∇⋅j (used in inverse problem)

elseif(solflag==1) # forward problem: source+medium --> fields
# To 'fix' the nullspace, and make the source consistent with the nullspace --> not always fast/memory allowable since
# julia backslash operator does not like non-square matrices... (?)
	#ix=Array{Float64,1};
	#ix=ones(nz*nx); # same row: row index=1
	#iz=collect(1:nz*nx); # column indices
	#v=ones(nz*nx); # actual values
	#onemat=sparse(ix,iz,v);
	# A = [onemat;A]; # fix null space
	# psi = reshape(A\[0; src], nz, nx); # make source consistent with null space

	psi_vec=A\src[:]; # direct inverse solve using backslash operator
	psi = reshape(psi_vec, nz, nx); # reshape electrical potential vector psi into 2D electrical potential map (z,x)
	psi = psi .- StatsBase.mean(StatsBase.mean(psi)); # alternative to 'fix' nullspace: make zero-mean...
#	src=reshape(src,nz,nx); # reshape source vector into 2D source map for plotting purposes only.

    	return psi,A  # function command --> return these variables, in this case electrical potential ψ
end # end solflag if statement

end # function command

end # end module Poisson


text:
function as_symbols(them)
    map(Symbol, (them...,))
end

split_keyword(keyword::Expr) =
    if keyword.head === :kw
    Pair(keyword.args[1], keyword.args[2])
else
    error("Cannot split keyword $keyword")
end

# Split a function call into its pieces
# Normalize non-function-like patterns into function calls
function split_call(call_expression::Expr)
    if @capture call_expression call_(arguments__; keywords__)
        (call, arguments...), (; map(split_keyword, keywords)...)
    elseif @capture call_expression call_(arguments__)
        (call, arguments...), ()
    else
        error("$call_expression is not a function call")
    end
end


text:
abstract type Vertebrate end
abstract type Mammal <: Vertebrate end
abstract type Reptile <: Vertebrate end

struct Cat <: Mammal
    name
end

struct Dog <: Mammal
    name
end

struct Crocodile <: Reptile
    name
end

Base.show(io::IO, cat::Cat) = print(io, "Cat ", cat.name)
Base.show(io::IO, dog::Dog) = print(io, "Dog ", dog.name)
Base.show(io::IO, croc::Crocodile) = print(io, "Crocodile ", croc.name)

# adopt new pet
function adopt(m::Mammal)
    println(m, " is now adopted.")
    return m
end

# quick test about Liskov substitution principle

adopt(Cat("Felix"));
adopt(Dog("Clifford"));
#=
julia> adopt(Cat("Felix"))
"Cat Felix is now adopted."

julia> adopt(Dog("Clifford"))
"Dog Clifford is now adopted."
=#

adopt(Crocodile("Solomon"));
#=
julia> adopt(Crocodile("Solomon"))
ERROR: MethodError: no method matching adopt(::Crocodile)
Closest candidates are:
  adopt(::Mammal) at REPL[9]:2
=#

# Covariant?  Does Cat <: Mammal imply Array{Cat} <: Array{Mammal}?
adopt(ms::Array{Mammal,1}) = "adopted " * string(ms)

adopt([Cat("Felix"), Cat("Garfield")])
#=
julia> adopt([Cat("Felix"), Cat("Garfield")])
ERROR: MethodError: no method matching adopt(::Array{Cat,1})
Closest candidates are:
  adopt(::Array{Mammal,1}) at REPL[48]:1
  adopt(::Mammal) at REPL[33]:3
=#

adopt(Mammal[Cat("Felix"), Cat("Garfield")])
#=
julia> adopt(Mammal[Cat("Felix"), Cat("Garfield")])
"adopted Mammal[Cat Felix, Cat Garfield]"
=#

# The answer is no.
# But this this works because an Array{Mammal} was passed.
adopt([Cat("Felix"), Dog("Clifford")])
#=
julia> adopt([Cat("Felix"), Dog("Clifford")])
"accepted Mammal[Cat Felix, Dog Clifford]"
=#

# That's because Array{Mammal} is an array of pointers since
# Mammal isn't concrete.

# What we should have done?

# homongeneous array of objects with the same concrete type
function adopt(ms::Array{T,1}) where {T <: Mammal}
    return "accepted same kind:" * string(ms)
end
#=
julia> methods(adopt)
# 3 methods for generic function "adopt":
[1] adopt(ms::Array{Mammal,1}) in Main at REPL[22]:1
[2] adopt(m::Mammal) in Main at REPL[16]:2
[3] adopt(ms::Array{T,1}) where T<:Mammal in Main at REPL[26]:1
=#

adopt([Cat("Felix"), Cat("Garfield")])
adopt([Dog("Clifford"), Dog("Astro")])
adopt([Cat("Felix"), Dog("Clifford")])

#=
julia> adopt([Cat("Felix"), Cat("Garfield")])
"accepted Cat[Cat Felix, Cat Garfield]"

julia> adopt([Dog ("Clifford"), Dog ("Astro")])
"accepted Dog [Dog  Clifford, Dog  Astro]"

julia> adopt([Cat("Felix"), Dog ("Clifford")])
"accepted Mammal[Cat Felix, Dog  Clifford]"
=#

# This is great because the functional behavior of adoption
# may depend on the kind of pet being adopted.

# ------------------------------------------------------
# Method arguments

friend(m::Mammal, f::Mammal) = "$m and $f become friends."

Tuple{Cat,Cat} <: Tuple{Mammal,Mammal}
Tuple{Cat,Dog} <: Tuple{Mammal,Mammal}
Tuple{Dog,Cat} <: Tuple{Mammal,Mammal}
Tuple{Dog,Dog} <: Tuple{Mammal,Mammal}
#=
julia> Tuple{Cat,Cat} <: Tuple{Mammal,Mammal}
true

julia> Tuple{Cat,Dog} <: Tuple{Mammal,Mammal}
true

julia> Tuple{Dog,Cat} <: Tuple{Mammal,Mammal}
true

julia> Tuple{Dog,Dog} <: Tuple{Mammal,Mammal}
true
=#

Tuple{Cat,Crocodile} <: Tuple{Mammal,Mammal}
#=
julia> Tuple{Cat,Crocodile} <: Tuple{Mammal,Mammal}
false
=#

# ------------------------------------------------------
# Function type

#=
julia> all(
all(x::Tuple{Bool,Bool,Bool}) in Base at tuple.jl:390
all(x::Tuple{Bool,Bool}) in Base at tuple.jl:389
all(x::Tuple{Bool}) in Base at tuple.jl:388
all(x::Tuple{}) in Base at tuple.jl:387
all(B::BitArray) in Base at bitarray.jl:1627
all(a::AbstractArray; dims) in Base at reducedim.jl:664
all(f::Function, a::AbstractArray; dims) in Base at reducedim.jl:665
all(itr) in Base at reduce.jl:642
all(f, itr) in Base at reduce.jl:724
=#

all(isodd, [1, 2, 3, 4, 5])
#=
julia> all(isodd, [1, 2, 3, 4, 5])
false
=#

typeof(isodd) <: Function
#=
julia> typeof(isodd) <: Function
true
=#

typeof(isodd)
typeof(isodd) |> supertype
isabstracttype(Function)
#=
julia> typeof(isodd)
typeof(isodd)

julia> typeof(isodd) |> supertype
Function

julia> isabstracttype(Function)
true
=#

typeof(println) <: Function
all(println, [1, 2, 3, 4, 5])
#=
julia> typeof(println) <: Function
true

julia> all(println, [1, 2, 3, 4, 5])
1
ERROR: TypeError: non-boolean (Nothing) used in boolean context
=#

# specific function types
const SignFunctions = Union{typeof(isodd),typeof(iseven)};
myall(f::SignFunctions, a::AbstractArray) = all(f, a);
myall(isodd, [1, 3, 5])
myall(iseven, [2, 4, 6])
myall(println, [2, 4, 6])
#=
julia> const SignFunctions = Union{typeof(isodd),typeof(iseven)};

julia> myall(f::SignFunctions, a::AbstractArray) = all(f, a);

julia> myall(isodd, [1, 3, 5])
true

julia> myall(iseven, [2, 4, 6])
true

julia> myall(println, [2, 4, 6])
ERROR: MethodError: no method matching myall(::typeof(println), ::Array{Int64,1})
=#

# ------------------------------------------------------
# Function subtyping

female_dogs = [Dog("Pinky"), Dog("Pinny"), Dog("Moonie")]
female_cats = [Cat("Minnie"), Cat("Queenie"), Cat("Kittie")]

select(::Type{Dog}) = rand(female_dogs)
select(::Type{Cat}) = rand(female_cats)

# Function: Mammal -> Union{Dog,Cat}
match(m::Mammal) = select(typeof(m))

match(Dog("Astro"))
match(Cat("Garfield"))
#=
julia> match(Dog("Astro"))
Dog Moonie

julia> match(Cat("Garfield"))
Cat Kittie
=#

# It's ok to kiss mammals :-)
kiss(m::Mammal) = "$m kissed!"

# Meet a partner
function meet_partner(finder::Function, self::Mammal)
    partner = finder(self)
    kiss(partner)
end

meet_partner(match, Cat("Felix"))
#=
julia> meet_partner(match, Cat("Felix"))
"Cat Kittie kissed!"
=#

# How about Mammal -> Vertebrate?
neighbor(m::Mammal) = Crocodile("Solomon")

meet_partner(neighbor, Cat("Felix"))
#=
julia> meet_partner(neighbor, Cat("Felix"))
ERROR: MethodError: no method matching excite(::Crocodile)
Closest candidates are:
  excite(::Mammal) at REPL[28]:2
=#

# So, function return type needs to be the same or subtype of the
# expected type.  Returning super-type is unsafe.

# What about function arguments?

# Function: Cat -> Mammal
buddy(cat::Cat) = rand([Dog("Astro"), Dog("Goofy"), Cat("Lucifer")])

meet_partner(buddy, Cat("Felix"))
meet_partner(buddy, Dog("Chef"))
#=
julia> meet_partner(buddy, Cat("Felix"))
"Cat Lucifer kissed!"

julia> meet_partner(buddy, Dog("Chef"))
ERROR: MethodError: no method matching buddy(::Dog)
=#

# "Be liberal in what you accept and conservative in what you produce."

# ----------------------------------------------------------------
# Revisting Base.all function.

# Wrap a function in a parametric type that captures the argument types and return type
struct PredicateFunction{T,S}
    f::Function
end

# Call the underlying function
(pred::PredicateFunction{T,S})(x::T; kwargs...) where {T,S} =
    pred.f(x; kwargs...)

# quick experiment
PredicateFunction{Number,Bool}(iseven)(1)
PredicateFunction{Number,Bool}(iseven)(2)
#=
julia> PredicateFunction{Number,Bool}(iseven)(1)
false

julia> PredicateFunction{Number,Bool}(iseven)(2)
true
=#

# our safe version of `all`
function safe_all(pred::PredicateFunction{T,S}, a::AbstractArray) where
        {T <: Any, S <: Bool}
    all(pred, a)
end

#=
julia> safe_all(PredicateFunction{Number,Bool}(iseven), [1,2,3])
false

julia> safe_all(PredicateFunction{Number,Bool}(iseven), [2,4,6])
true
=#

# ----------------------------------------------------------------
# OPTIONAL.  Animal kingdom example continues.

# Technically, what does `meet_partner` need?

# The signature doesn't tell us much.
methods(meet_partner)
#=
julia> methods(meet_partner)
# 1 method for generic function "meet_partner":
[1] meet_partner(finder::Function, self::Mammal) in Main at REPL[29]:3
=#

# But, we know that needs a function of this type: Mammal -> Mammal
# Given that arguments are contravariant and return type is convariant,
# We can infer that any function with the following type would work:
# T -> S where {T >: Mammal, S <: Mammal}

# Until Julia supports typed functions, we can work around by
# formulating as such.

# One argument with type T, having return type S.
struct TypedFunction1{T,S}
    f::Function
end

(tf::TypedFunction1{T,S})(x::T; kwargs...) where {T,S} =
    tf.f(x; kwargs...)

# Let's redefine `meet_partner` to require that.
function meet_partner(finder::TypedFunction1{T,S}, self::Mammal) where
        {T >: Mammal, S <: Mammal}
    partner = finder(self)
    kiss(partner)
end

# Create a callable struct that takes Vertebrate and returns Cat
kissy_match(v::Vertebrate) = Cat("Kissy")
best_cat_matcher = TypedFunction1{Vertebrate,Cat}(kissy_match)

# Now we can pass this typed matcher
meet_partner(best_cat_matcher, Dog("Goofy"))

# Can we make a mistsake?  Let's try a Cat -> Cat.
cat_cat_match(cat::Cat) = Cat("Brownie")
brownie_matcher = TypedFunction1{Cat,Cat}(cat_cat_match)

meet_partner(brownie_matcher, Dog("Goofy"))
#=
julia> meet_partner(brownie_matcher, Dog("Goofy"))
ERROR: MethodError: no method matching meet_partner(::TypedFunction1{Cat,Cat}, ::Dog)
=#

# Cool, so we're completely type-safe with function arguments.

# Perhaps the syntax can be improved with metaprogramming? e.g.
# @tf F = Mammal ▶ Mammal function meet_partner(finder::F, self::Mammal)
#     partner = finder(self)
#     excite(partner)
# end


text:
KRED="\x1B[31m"
KGRN="\x1B[32m"
KRES="\033[0m"

using PdbTool
using Printf

macro spath()
	return dirname(Base.source_path()) * "/"
end

function TEST_parsePdb()
	@printf("Testing parsePdb ...")
	try
		pdb=PdbTool.parsePdb(@spath() * "5PTI.pdb");
		@printf("%s SUCCES%s\n",KGRN,KRES)
	catch
		@printf("%s FAIL (function exited with error)%s\n", KRED,KRES)
		return false
	end
	return true
end

function TEST_mapChainToHmm()
	@printf("Testing mapChainToHmm ...")
	try
		pdb=PdbTool.parsePdb(@spath() * "5PTI.pdb");
		PdbTool.mapChainToHmm(pdb.chain["A"],@spath() * "Kunitz_BPTI.hmm")
		@printf("%s SUCCESS %s\n",KGRN,KRES)
	catch
		@printf("%s FAIL (function exited with error) %s\n", KRED, KRES);
		return false
	end
	return true
end

function TEST_mapChainToHmmLegacy()
	@printf("Testing mapChainToHmmLegacy ...")
	try
		pdb=PdbTool.parsePdb(@spath() * "5PTI.pdb");
		PdbTool.mapChainToHmmLegacy(pdb.chain["A"],@spath() * "Kunitz_BPTI.hmm")
		@printf("%s SUCCESS\n%s",KGRN,KRES)
	catch
		@printf("%s FAIL (function exited with error) %s\n", KRED, KRES);
		return false
	end
	return true
end

# This function returns always true - the user is just informed that things
# don't work. This is because some tools (cmsearch e.g.) are not necessary for
# most functions and IF necessary the functions will check the external
# functions explicitly
function TEST_externals()

	@printf("Testing access to external functions ... \n")
	@printf("\thmmalign ... ")
	if !PdbTool.EXT_TEST.EXT_TEST_hmmalign()
		@printf("%s NOT AVAILABLE (mapChainToHmm will not work)%s\n", KRED, KRES);
	else
		@printf("%s AVAILABLE\n%s",KGRN,KRES)
	end
        # FIXME!!!
	# @printf("\thmmsearch ... ")
	# if !PdbTool.EXT_TEST.EXT_TEST_hmmsearch()
	# 	@printf("%s NOT AVAILABLE (mapChainToHmmLegacy will not work)%s\n", KRED, KRES);
	# else
	# 	@printf("%s AVAILABLE\n%s",KGRN,KRES)
	# end

	@printf("\tcmsearch ... ")
	if !PdbTool.EXT_TEST.EXT_TEST_hmmsearch()
		@printf("%s NOT AVAILABLE (mapping RNAs will not work)%s\n", KRED, KRES);
	else
		@printf("%s AVAILABLE\n%s",KGRN,KRES)
	end

	return true
end

function TEST_all()

	if !TEST_externals()
		return false
	end

	if !TEST_parsePdb()
		return false
	end

	if !TEST_mapChainToHmm()
		return false
	end

	# if !TEST_mapChainToHmmLegacy() #FIXME!
	# 	return false
	# end

	return true
end

function testall()
	if !TEST_all()
		@printf("%s Stopped - if you cannot fix the problem please consider reporting to https://github.com/christophfeinauer/PdbTool/issues %s\n", KRED,KRES);
	else
		@printf("%sAll tests finished with success. %s", KGRN,KRES)
	end
end

testall()

max_stars_repo_name: christophfeinauer/PdbTool

max_issues_repo_name: christophfeinauer/PdbTool

max_forks_repo_name: christophfeinauer/PdbTool


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8aac8df3-1551-4c9f-a8bd-a62751a29b2a
md"""
## Path 1
"""

# ╔═╡ 03307e43-cb61-4321-95ac-7bbb16e0cfc6
@bind x html"<input type=range max=10000>"

# ╔═╡ 692746e0-7c96-47ac-b1ee-5ff34ee66751
@bind y html"<input type=range max=10000>"

# ╔═╡ b18c2329-18d7-4041-962c-0ef98f8aa591
(x,y)

# ╔═╡ 399ef117-2085-4c9d-9d8d-d03a03baf5ef
md"""
## Path 2
"""

# ╔═╡ a3a04d5f-b0f0-4740-9b37-92570864f142
high_res = true

# ╔═╡ a822ac82-5691-4f8e-a60a-1a4582cf59e7
dog_file = if high_res
		download("https://upload.wikimedia.org/wikipedia/commons/e/ef/Pluto_in_True_Color_-_High-Res.jpg")
else
	download("https://upload.wikimedia.org/wikipedia/commons/thumb/e/ef/Pluto_in_True_Color_-_High-Res.jpg/240px-Pluto_in_True_Color_-_High-Res.jpg")
end

# ╔═╡ 627461ce-9e80-4707-b0ba-ddc6bb9b4269
begin
	struct Dog end
	function Base.show(io::IO, ::MIME"image/jpg", ::Dog)
		write(io, read(dog_file))
	end
end

# ╔═╡ 5ce8ebc6-b509-42f0-acd5-8008673b04ab
md"Downloaded image is $(filesize(dog_file) / 1000) kB"

# ╔═╡ 1f48fe19-3ee8-44ac-a591-7b4df2d2f93a
md"""
This cell will have very large Uint8Arrays in the output body
"""

# ╔═╡ 5539db10-b0d2-48b6-8985-ef437b8ae0b5
@bind show_dogs html"<input type=checkbox>"

# ╔═╡ 74329553-ab9b-4b6c-a77b-9c24ac48490b
show_dogs === true && Dog()

# ╔═╡ 8811b8bd-7d44-4f4a-b52d-0948dad39a51
md"""
## Path 3
"""

# ╔═╡ f0575125-d7dd-4cf5-bfd3-3275d6bdd0ce
a = 100

# ╔═╡ 802be3da-6f24-4693-b48e-0479aec4a02c
begin
	a
	@bind b html"<input>"
end

# ╔═╡ 997f714f-1df4-4b2c-ab0a-50dd52cb4a82
md"""
## Path 4
"""

# ╔═╡ a504bc18-fff4-4cd7-b74c-173abcf1ef68
begin
	a
	@bind c html"<input>"
end

# ╔═╡ 00a0169e-fd93-4dfe-b45b-f088312e24ab
md"""
## Path 5
"""

# ╔═╡ e7b5f5ec-d673-4397-ac0c-7a4bfc364fde
f(x) = x

# ╔═╡ 0a61c092-e61a-4219-a57f-172a8c8c4117
begin
	f(1)
	@bind five1 html"<input>"
end

# ╔═╡ 31695eb5-7d77-4d6b-a32a-ff31cf1fd8e9
md"""
## Path 6
"""

# ╔═╡ 25774e63-94fb-4b03-a43a-72a62f08f0c5
begin
	f(1)
	@bind five2 html"<input>"
end

# ╔═╡ a04672a7-3d3f-4b0b-8c8f-f8b73f208de4
md"""
## Path 7
"""

# ╔═╡ 2d0f4354-20e0-47a0-8a86-99cd50bca80f
@bind six1 html"<input>"

# ╔═╡ fee1f75a-5f3e-4213-8ba3-872871cf7f68
@bind six2 html"<input>"

# ╔═╡ 20ebc38b-24c1-4e39-97b0-dcd53a9a5bf7
@bind six3 html"<input>"

# ╔═╡ 5545ac33-82e9-4994-be49-5776d512e2c1
(six1, six2)

# ╔═╡ b776eba5-60bc-4d0f-8e93-e2baf2f695bc
md"""
## Path 8
"""

# ╔═╡ a56b8b24-bf38-49bf-b19e-d8feadd55db3
(six2, six3)

# ╔═╡ 6784f0e5-5108-48ca-99ac-9d154b1d3c55
md"""
## Path 9
"""

# ╔═╡ aacc0632-797f-49d6-b9ce-3e728fade6c3
begin
	@bind cool1 let
		@bind cool2 html"nothing"

		html"<input type=range>"
	end
end

# ╔═╡ 6357cfaf-e6d6-49f2-b3bb-5a27b28cf0fa
cool1

# ╔═╡ 4fd43c26-dd05-4520-93e9-901760ef49b4
cool2

# ╔═╡ 357762fc-52fe-4727-b0a8-af55eea466b7
md"""
## Path 10
"""

# ╔═╡ 0024288a-9ee3-42b5-82c6-d996f45be9ed
let
	md"""

	Hello $(@bind world html"<input value=world>")
	"""
end

# ╔═╡ ef732ddf-b034-42cf-815f-1c8b53da6401
world

# ╔═╡ c143b2de-78c5-46ad-852f-3c2a9115cb72
md"""
## Path 11
"""

# ╔═╡ cf628a57-933b-4984-a317-63360c345534
@bind boring html"<input>"

# ╔═╡ 22659c85-700f-4dad-a22a-7aafa71225c0
# boring is never referenced

# ╔═╡ Cell order:
# ╟─8aac8df3-1551-4c9f-a8bd-a62751a29b2a
# ╠═03307e43-cb61-4321-95ac-7bbb16e0cfc6
# ╠═692746e0-7c96-47ac-b1ee-5ff34ee66751
# ╠═b18c2329-18d7-4041-962c-0ef98f8aa591
# ╟─399ef117-2085-4c9d-9d8d-d03a03baf5ef
# ╠═627461ce-9e80-4707-b0ba-ddc6bb9b4269
# ╠═a3a04d5f-b0f0-4740-9b37-92570864f142
# ╠═a822ac82-5691-4f8e-a60a-1a4582cf59e7
# ╟─5ce8ebc6-b509-42f0-acd5-8008673b04ab
# ╟─1f48fe19-3ee8-44ac-a591-7b4df2d2f93a
# ╠═5539db10-b0d2-48b6-8985-ef437b8ae0b5
# ╠═74329553-ab9b-4b6c-a77b-9c24ac48490b
# ╟─8811b8bd-7d44-4f4a-b52d-0948dad39a51
# ╠═f0575125-d7dd-4cf5-bfd3-3275d6bdd0ce
# ╠═802be3da-6f24-4693-b48e-0479aec4a02c
# ╟─997f714f-1df4-4b2c-ab0a-50dd52cb4a82
# ╠═a504bc18-fff4-4cd7-b74c-173abcf1ef68
# ╟─00a0169e-fd93-4dfe-b45b-f088312e24ab
# ╠═e7b5f5ec-d673-4397-ac0c-7a4bfc364fde
# ╠═0a61c092-e61a-4219-a57f-172a8c8c4117
# ╟─31695eb5-7d77-4d6b-a32a-ff31cf1fd8e9
# ╠═25774e63-94fb-4b03-a43a-72a62f08f0c5
# ╟─a04672a7-3d3f-4b0b-8c8f-f8b73f208de4
# ╠═2d0f4354-20e0-47a0-8a86-99cd50bca80f
# ╠═fee1f75a-5f3e-4213-8ba3-872871cf7f68
# ╠═20ebc38b-24c1-4e39-97b0-dcd53a9a5bf7
# ╠═5545ac33-82e9-4994-be49-5776d512e2c1
# ╟─b776eba5-60bc-4d0f-8e93-e2baf2f695bc
# ╠═a56b8b24-bf38-49bf-b19e-d8feadd55db3
# ╟─6784f0e5-5108-48ca-99ac-9d154b1d3c55
# ╠═aacc0632-797f-49d6-b9ce-3e728fade6c3
# ╠═6357cfaf-e6d6-49f2-b3bb-5a27b28cf0fa
# ╠═4fd43c26-dd05-4520-93e9-901760ef49b4
# ╟─357762fc-52fe-4727-b0a8-af55eea466b7
# ╠═0024288a-9ee3-42b5-82c6-d996f45be9ed
# ╠═ef732ddf-b034-42cf-815f-1c8b53da6401
# ╟─c143b2de-78c5-46ad-852f-3c2a9115cb72
# ╠═cf628a57-933b-4984-a317-63360c345534
# ╠═22659c85-700f-4dad-a22a-7aafa71225c0


text:
using Earth2014
using Test

# For CI, make sure the downloads do not hang
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

@testset "Earth2014.jl" begin
    @testset "$res resolution" for res in ["1min", "5min"]
        x,y,z = Earth2014.load(res=res)
        @test x isa Vector
        @test y isa Vector
        @test z isa Array
        @test size(z) == (length(y), length(x))
    end
end


text:
# Destabilizing polarization as a function of number of attributes.
# Figure 1
using DrWatson
quickactivate(@__DIR__)

using PolarizationFramework

ns = [5, 9]#[9,15,25]
reps = [10, 10]
reps_dict = Dict(zip(ns, reps))
gs = [1:2:21..., 25, 29, 33, 37, 41, 45, 49, 55, 61, 67, 73, 81, 89, 97]

threshold = 0.5;
vs = [4, @onlyif("attr_types" == "OA", 1000)] #includes CA

attr_types = ["BA", "UA", "OA", "UPA"]

# sg = [-1.5, -1.1, -1.01, -0.99, -0.5, -0.1, -0.01]
# gammas=sort([sg..., 0, (-sg)...]);
# gammas = [2, 2.5, 3, 3.5]
gammas = [0.5, 1.5, 2, 2.5, 3, 3.5, 4, 6]
# gammas = [0.5]
# gammas = [6]

larger_sizes = [@onlyif(i <= "ns" < 2 * i, i) for i = 1:maximum(ns)]

all_params = @strdict(ns, gs, threshold, vs, attr_types, larger_sizes)
dicts = dict_list(all_params)
[d["reps"] = reps_dict[d["ns"]] for d in dicts] #adding reps
# [d["gammas"] = gammas for d in dicts] #adding gammas

for params in dicts
    n, g, attr_type, v, rep, larger_size = let
        @unpack ns, threshold, reps, gs, attr_types, vs, larger_sizes = params
        ns, gs, attr_types, vs, reps, larger_sizes
    end
    println("Started ", @ntuple(n, g, attr_type, v, larger_size))

    if attr_type == "UA"
        attr = UnorderedAttributes(g, threshold, v)
    elseif attr_type == "BA"
        attr = BinaryAttributes(g)
    elseif attr_type == "OA"
        attr = OrderedAttributes(g, threshold, v)
    elseif attr_type == "UPA"
        attr = UnorderedPositiveAttributes(g, threshold, v)
    else
        throw(attr_type)
    end

    r = using_heider_attr_destab(
        n,
        attr,
        gammas,
        larger_size,
        rep,
        3000.0,
        "Heider7!",
        disp_each = 0,
        disp_more_every = 600,
        save_each = 600,
        files_folder = ["data", "sims"],
        filename_prefix = "DestabFig1",
    )

end

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
## features

"""
    get_front_neighbor(env::MergingEnvironment, scene::Scene, egoid::Int64)
returns the front neighbor of `egoid` in its lane. It returns an object of type `NeighborLongitudinalResult`
"""
function get_front_neighbor(env::MergingEnvironment, scene::Scene, egoid::Int64)
    ego_ind = findfirst(egoid, scene)
    ego = scene[ego_ind]
    # merge neighbor
    ego_lane = get_lane(env.roadway, ego)
    if ego_lane == main_lane(env)
        fore_res = get_neighbor_fore_along_lane(scene, ego_ind, env.roadway)
    else
        lane = main_lane(env)
        # s_base = lane[env.merge_proj.curveproj.ind, env.roadway].s
        s_base = env.main_lane_length
        fore_res = get_neighbor_fore_along_lane(scene, env.roadway, lane.tag, s_base, index_to_ignore=ego_ind)
        # fore_res = NeighborLongitudinalResult(fore_res.ind, fore_res.Δs - dist_to_merge(env, ego))
    end
end

"""
    get_neighbors(env::MergingEnvironment, scene::Scene, egoid::Int64)
returns the following neighbors id and relative distance (if they exist)
    - the front neighbor of vehicle `egoid`
    - the vehicle right behind the merge point (if `egoid` is on the main lane)
    - the front neighbor of the projection of `egoid` on the main lane
    - the rear neighbor of the projection of `egoid` on the merge lane
"""
function get_neighbors(env::MergingEnvironment, scene::Scene, egoid::Int64)
    ego_ind = findfirst(egoid, scene)
    ego = scene[ego_ind]

    #front neighbor
    front = get_front_neighbor(env, scene, egoid)

    # merge neighbor
    ego_lane = get_lane(env.roadway, ego)
    if ego_lane == main_lane(env)
        fore_res = get_neighbor_rear_along_lane(scene, ego_ind, env.roadway)
    else
        lane = main_lane(env)
        # s_base = lane[env.merge_proj.curveproj.ind, env.roadway].s
        s_base = env.main_lane_length
        fore_res = get_neighbor_rear_along_lane(scene, env.roadway, lane.tag, s_base, index_to_ignore=ego_ind)
        # fore_res = NeighborLongitudinalResult(fore_res.ind, fore_res.Δs + dist_to_merge(env, ego))
    end
    # two closest car in main lane
    proj_lane = main_lane(env)
    main_lane_proj = proj(ego.state.posG, proj_lane, env.roadway)
    s_main = proj_lane[main_lane_proj.curveproj.ind, env.roadway].s
    fore_main = get_neighbor_fore_along_lane(scene, env.roadway, proj_lane.tag, s_main, index_to_ignore=ego_ind)
    rear_main = get_neighbor_rear_along_lane(scene, env.roadway, proj_lane.tag, s_main, index_to_ignore=ego_ind)
    return front, fore_res, fore_main, rear_main
end

"""
    dist_to_merge(env::MergingEnvironment, veh::Vehicle)
returns the distance to the merge point.
"""
function dist_to_merge(env::MergingEnvironment, veh::Vehicle)
    lane = get_lane(env.roadway, veh)
    if lane == main_lane(env)
        frenet_merge = get_frenet_relative_position(veh.state.posG, env.merge_index, env.roadway)
        dist = frenet_merge.Δs
    else
        dist = veh.state.posF.s - get_end(lane)
    end
    return dist
end

"""
    time_to_merge(env::MergingEnvironment, veh::Vehicle, a::Float64 = 0.0)
return the time to reach the merge point using constant acceleration prediction.
If the acceleration, `a` is not specified, it performs a constant velocity prediction.
"""
function time_to_merge(env::MergingEnvironment, veh::Vehicle, a::Float64 = 0.0)
    d = -dist_to_merge(env, veh)
    v = veh.state.v
    t = Inf
    if isapprox(a, 0)
        t =  d/veh.state.v
    else
        delta = v^2 + 2.0*a*d
        if delta < 0.0
            t = Inf
        else
            t = (-v + sqrt(delta)) / a
        end
        if t < 0.0
            t = Inf
        end
    end
    return t
end

"""
    find_merge_vehicle(env::MergingEnvironment, scene::Scene)
returns the id of the merging vehicle if there is a vehicle on the merge lane.
"""
function find_merge_vehicle(env::MergingEnvironment, scene::Scene)
    for veh in scene
        lane = get_lane(env.roadway, veh)
        if lane == merge_lane(env)
            return veh
        end
    end
    return nothing
end

"""
    constant_acceleration_prediction(env::MergingEnvironment, veh::Vehicle, acc::Float64, time::Float64, v_des::Float64)
returns the state of vehicle `veh` after time `time` using a constant acceleration prediction.

# inputs
- `env::MergingEnvironment` the environment
- `veh::Vehicle` the initial state of the vehicle
- `acc::Float64` the current acceleration of the vehicle
- `time::Float64` the prediction horizon
- `v_des::Float64` the desired speed of the vehicle (assumes that the vehicle will not exceed that speed)
"""
function constant_acceleration_prediction(env::MergingEnvironment,
                                          veh::Vehicle,
                                          acc::Float64,
                                          time::Float64,
                                          v_des::Float64)
        # act = LaneFollowingAccel(acc)
        # vehp = propagate(veh, act, env.roadway, time, true)
        v1 = veh.state.v
        v2 = clamp(veh.state.v + acc*time , 0.0, v_des)
        if acc ≈ 0.0
            Δs = v1*time
        else
            Δs = (v2^2 - v1^2) / (2*acc)
        end
        Δs = max(0., Δs)
        sp = veh.state.posF.s + Δs
        lane = get_lane(env.roadway, veh)
        vehp = vehicle_state(sp, lane, v2, env.roadway)
        return Vehicle(vehp, veh.def, veh.id)
end

"""
    distance_projection(env::MergingEnvironment, veh::Vehicle)
Performs a projection of `veh` onto the main lane. It returns the longitudinal position of the projection of `veh` on the main lane.
The projection is computing by conserving the distance to the merge point.
"""
function distance_projection(env::MergingEnvironment, veh::Vehicle)
    if get_lane(env.roadway, veh) == main_lane(env)
        return veh.state.posF.s
    else
        dm = -dist_to_merge(env, veh)
        return env.roadway[env.merge_index].s - dm
    end
end

"""
    collision_time(env::MergingEnvironment, veh::Vehicle, mergeveh::Vehicle, acc_merge::Float64, acc_min::Float64)
compute the collision time between two vehicles assuming constant acceleration.
"""
function collision_time(env::MergingEnvironment,
                        veh::Vehicle,
                        mergeveh::Vehicle,
                        acc_merge::Float64,
                        acc_min::Float64)
    rel_vel = mergeveh.state.v - veh.state.v
    rel_pos = distance_projection(env, mergeveh) - distance_projection(env, veh)
    rel_acc = acc_merge - acc_min
    delta = rel_vel^2 - 2*rel_acc*rel_pos
    if delta < 0.0
        return nothing
    elseif rel_acc != 0.0
        t_coll = (-rel_vel + sqrt(delta))/rel_acc
        return t_coll
    elseif rel_vel != 0.0
        t_coll = - rel_pos / rel_vel
        return t_coll
    else
        t_coll = nothing
        return t_coll
    end
end

"""
    braking_distance(v::Float64, t_coll::Float64, acc::Float64)
computes the distance to reach a velocity of 0. at constant acceleration `acc` in time `t_coll` with initial velocity `v`
"""
function braking_distance(v::Float64, t_coll::Float64, acc::Float64)
    brake_dist = v*t_coll + 0.5*acc*t_coll^2
    return brake_dist
end

"""
    AutomotiveDrivingModels.propagate(veh::Entity{VehicleState,D,I}, action::LaneFollowingAccel, roadway::Roadway, ΔT::Float64, nobackup::Bool) where {D,I}
A propagate method for `LaneFollowingAccel` that prevents the car from backing up in `nobackup = true`.
"""
function AutomotiveDrivingModels.propagate(veh::Entity{VehicleState,D,I}, action::LaneFollowingAccel, roadway::Roadway, ΔT::Float64, nobackup::Bool) where {D,I}

    a_lon = action.a

    ds = veh.state.v

    ΔT² = ΔT*ΔT

    Δs = ds*ΔT + 0.5*a_lon*ΔT²
    v₂ = ds + a_lon*ΔT

    if nobackup
        Δs = max(Δs, 0.0)
        v₂ = max(v₂, 0.0)
    end

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    posG = roadway[roadind].pos
    posF = Frenet(roadind, roadway, t=veh.state.posF.t, ϕ=veh.state.posF.ϕ)
    VehicleState(posG, posF, v₂)
end

function AutomotiveDrivingModels.tick!(
    scene::EntityFrame{S,D,I},
    roadway::R,
    actions::Vector{A},
    Δt::Float64,
    nobackup::Bool
    ) where {S,D,I,R,A}

    for i in 1 : length(scene)
        veh = scene[i]
        state′ = propagate(veh, actions[i], roadway, Δt, nobackup)
        scene[i] = Entity(state′, veh.def, veh.id)
    end

    return scene
end


text:
# TODO: Rotations.jl

# General axis
function rotation_matrix(axis::NTuple{3, <: Real}, α::Real; out=Matrix{Float64}(undef, 3, 3))
	# Arbitrary axis
	# Apply vector from the right
	(sinα, cosα) = sincos(α)
	one_minus_cosα = (1 - cosα)
	# First column
	out[1, 1] = cosα + axis[1]^2 * one_minus_cosα
	out[2, 1] = axis[1] * axis[2] * one_minus_cosα + axis[3] * sinα
	out[3, 1] = axis[1] * axis[3] * one_minus_cosα - axis[2] * sinα
	# Second column
	out[1, 2] = axis[1] * axis[2] * one_minus_cosα - axis[3] * sinα
	out[2, 2] = cosα + axis[2]^2 * one_minus_cosα
	out[3, 2] = axis[2] * axis[3] * one_minus_cosα + axis[1] * sinα
	# Third column
	out[1, 3] = axis[1] * axis[3] * one_minus_cosα + axis[2] * sinα
	out[2, 3] = axis[2] * axis[3] * one_minus_cosα - axis[1] * sinα
	out[3, 3] = cosα + axis[3]^2 * one_minus_cosα
	return out
end

# Around axis in plane alinged with coordinate system
function rotation_matrix(plane::Val{:xy}, axis::NTuple{2, <: Real}, α::Real; out=Matrix{Float64}(undef, 3, 3))
	# axis in xy plane, axis must be normalised!
	# Apply vector from the right
	(sinα, cosα) = sincos(α)
	one_minus_cosα = (1 - cosα)
	# First column
	out[1, 1] = cosα + axis[1]^2 * one_minus_cosα
	out[2, 1] = axis[1] * axis[2] * one_minus_cosα
	out[3, 1] = -axis[2] * sinα
	# Second column
	out[1, 2] = axis[1] * axis[2] * one_minus_cosα
	out[2, 2] = cosα + axis[2]^2 * one_minus_cosα
	out[3, 2] = axis[1] * sinα
	# Third column
	out[1, 3] = axis[2] * sinα
	out[2, 3] = -axis[1] * sinα
	out[3, 3] = cosα
	return out
end

# Around axis in plane alinged with coordinate system
function rotation_matrix(plane::Val{:xz}, axis::NTuple{2, <: Real}, α::Real; out=Matrix{Float64}(undef, 3, 3))
	# Arbitrary axis
	# Apply vector from the right
	(sinα, cosα) = sincos(α)
	one_minus_cosα = (1 - cosα)
	# First column
	out[1, 1] = cosα + axis[1]^2 * one_minus_cosα
	out[2, 1] = axis[2] * sinα
	out[3, 1] = axis[1] * axis[2] * one_minus_cosα
	# Second column
	out[1, 2] = -axis[2] * sinα
	out[2, 2] = cosα
	out[3, 2] = axis[1] * sinα
	# Third column
	out[1, 3] = axis[1] * axis[2] * one_minus_cosα
	out[2, 3] = -axis[1] * sinα
	out[3, 3] = cosα + axis[2]^2 * one_minus_cosα
	return out
end

# Axis defined by phase
function rotation_matrix(plane::Val{:xy}, θ::Real, α::Real; out=Matrix{Float64}(undef, 3, 3))
	# Rotate by α, around axis specified by θ
	# Right handed coordinate system, x,y,z
	# This is for B1 in xy plane
	# Apply vector from the right
	(sinα, cosα) = sincos(α)
	(sinθ, cosθ) = sincos(θ)
	# First column
	out[1, 1] = sinθ^2 + cosα * cosθ^2
	out[2, 1] = sinθ * cosθ - cosα * sinθ * cosθ
	out[3, 1] = -sinα * cosθ
	# Second column
	out[1, 2] = sinθ * cosθ - cosα * sinθ * cosθ
	out[2, 2] = cosθ^2 + cosα * sinθ^2
	out[3, 2] = sinα * sinθ
	# Third column
	out[1, 3] = sinα * cosθ
	out[2, 3] = -sinα * sinθ
	out[3, 3] = cosα
	return out
end

function rotation_matrix(plane::Val{:xz}, θ::Real, α::Real; out=Matrix{Float64}(undef, 3, 3))
	# Rotate by α, around axis specified by θ
	# Right handed coordinate system, x,y,z
	# This is for B1 in xz plane
	# Apply vector from the right
	(sinα, cosα) = sincos(α)
	(sinθ, cosθ) = sincos(θ)
	# First column
	out[1, 1] = sinθ^2 + cosα * cosθ^2
	out[2, 1] = -sinα * cosθ
	out[3, 1] = sinθ * cosθ - cosα * sinθ * cosθ
	# Second column
	out[1, 2] = sinα * cosθ
	out[2, 2] = cosα
	out[3, 2] = -sinα * sinθ
	# Third column
	out[1, 3] = sinθ * cosθ - cosα * sinθ * cosθ
	out[2, 3] = sinα * sinθ
	out[3, 3] = cosθ^2 + cosα * sinθ^2

	return out
end


text:
@testset "bit_flipping_env" begin
    rng = StableRNG(123)
    env = BitFlippingEnv(; N = 7, T = 2, rng = rng)
    test_state = state(env, GoalState())
    env(1)
    @test is_terminated(env) == false
    env(1)
    @test is_terminated(env) == true
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
"""
The Boid.jl module defines types and functions to support autonomous,
distributed computing/processing.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
module Boid

include("AbstractProcessingCore.jl")

include("AbstractChannelData.jl")
include("InputChannel.jl")
include("OutputChannel.jl")

include("AbstractControlLogicCore.jl")
include("ControlUnit.jl")

include("Node.jl")

# Methods
include("utils.jl")

end  # End of Boid.jl module


text:
using QuantEcon: tauchen

model = :CapitalAdjustModel
dipar_1 = Dict(:nK => 100, :nz => 20, :β => 0.9, :ρ => 0.6, :σ => 0.3)
p_1_all = createmodel(model; dipar_1..., intdim = :All)
p_1_sep = createmodel(model; dipar_1..., intdim = :Separable)
p_1_sep_states = createmodel(model; dipar_1..., intdim = :Separable_States)
p_1_sep_exogstates = createmodel(model; dipar_1..., intdim = :Separable_ExogStates)

T_1_all = transitionmatrix(p_1_all)
T_1_sep = transitionmatrix(p_1_sep)
T_1_sep_states = transitionmatrix(p_1_sep_states)
T_1_sep_exogstates = transitionmatrix(p_1_sep_exogstates)

model = :CapitalAdjustModel2
dipar_2 = Dict(:nK => 20, :nN=>15, :nz => 3, :β => 0.9, :ρ => 0.6, :σ => 0.3)
p_2_all = createmodel(model; dipar_2..., intdim = :All)
p_2_sep = createmodel(model; dipar_2..., intdim = :Separable)
p_2_sep_states = createmodel(model; dipar_2..., intdim = :Separable_States)
p_2_sep_exogstates = createmodel(model; dipar_2..., intdim = :Separable_ExogStates)

T_2_all = transitionmatrix(p_2_all)
T_2_sep = transitionmatrix(p_2_sep)
T_2_sep_states = transitionmatrix(p_2_sep_states)
T_2_sep_exogstates = transitionmatrix(p_2_sep_exogstates)

@testset "transition matrix" begin

	@testset "compare" begin

		@testset "CapitalAdjustModel" begin

			T_all = T_1_all
			T_sep = T_1_sep
			T_sep_states = T_1_sep_states
			T_sep_exogstates = T_1_sep_exogstates

			nK = dipar_1[:nK]
			nz = dipar_1[:nz]

			# compare All and Separable
			# to do

			# compare Separable and Separable_States
			# to do

			# compare Separable_States and Separable_ExogStates
			tT_test = [reshape(T_sep_states,(nK, nz, nz))[j,:,:] for j = 1 : 5]
			@test all([T_sep_exogstates == T_test for T_test in tT_test])

		end

		@testset "CapitalAdjustModel2" begin

			T_all = T_2_all
			T_sep = T_2_sep
			T_sep_states = T_2_sep_states
			T_sep_exogstates = T_2_sep_exogstates

			nK = dipar_2[:nK]
			nN = dipar_2[:nN]
			nz = dipar_2[:nz]

			# compare All and Separable
			# to do

			# compare Separable and Separable_States
			# to do

			# compare Separable_States and Separable_ExogStates
			tT_test = [reshape(T_sep_states, (nK,nN,nz,nz))[j,j,:,:] for j = 1 : 5]
			@test all([T_sep_exogstates == T_test for T_test in tT_test])

		end

	end # compare

	@testset "external transitionmatrix" begin

		nz = dipar_1[:nz]
		ρ = dipar_1[:ρ]
		σ = dipar_1[:σ]

		# have same nodes as in tauchen?
		@test isapprox(collect(tauchen(nz, ρ, σ).state_values),
			p_1_sep_exogstates.tStateVectors[2], rtol = mytol)
		solve(p_1_sep_exogstates, mTransition = tauchen(nz, ρ, σ).p)

	end

end # transition matrix


text:
module UtilsTest

using Test
using ProximalBase
using Random
using LinearAlgebra
using SparseArrays
using Statistics

@testset "shrink" begin

  @test shrink(3.0, 1.0) == 2.0
  @test shrink(-2.0, 1.0) == -1.0
  @test shrink(0.5, 1.0) == 0.

  Random.seed!(1)
  p = 10
  x = randn(p)
  c = norm(x) + 0.1
  @test shrinkL2!(similar(x), x, c) == zeros(p)
  c = c - 0.15
  @test shrinkL2!(similar(x), x, c) == max(1. - c / norm(x), 0,) * x

  x = randn(p, 2*p)
  c = norm(x) + 0.1
  @test shrinkL2!(similar(x), x, c) == zeros(p, 2*p)
  c = c - 0.15
  @test shrinkL2!(similar(x), x, c) == max(1. - c / norm(x), 0,) * x
end

@testset "multiply" begin

  Random.seed!(1)
  n, m = 100, 100
  A = randn(n, m)
  b = randn(m)
  x = A * b
  xt = A' * b
  for i=1:n
    @test A_mul_B_row(A, b, i) ≈ x[i]
    @test At_mul_B_row(A, b, i) ≈ xt[i]
  end

  bs = sprand(m, 0.3)
  x = A * bs
  xt = A' * bs
  bi = SparseIterate(bs)
  xi = A * bi
  xti = A' * bi
  for i=1:n
    @test A_mul_B_row(A, bs, i) ≈ x[i]
    @test At_mul_B_row(A, bs, i) ≈ xt[i]

    @test A_mul_B_row(A, bi, i) ≈ x[i]
    @test At_mul_B_row(A, bi, i) ≈ xt[i]

    @test xi[i] ≈ x[i]
    @test xti[i] ≈ xt[i]
  end

  # symmetric A * X * B
  X = randn(100, 50)
  Σx = Symmetric(cov(X))
  Y = randn(100, 50)
  Σy = Symmetric(cov(Y))
  Δ = sprand(50, 50, 0.1)
  Δ = (Δ + Δ') / 2.

  A = Σx * Δ * Σy
  Δs = SymmetricSparseIterate(Δ)
  @test A_mul_X_mul_B(Σx, Δs, Σy) ≈ A

  # symmetric A * U * U' * B
  X = randn(100, 50)
  Σx = Symmetric(cov(X))
  Y = randn(100, 50)
  Σy = Symmetric(cov(Y))
  U = randn(50, 5)

  A = (Σx * U)*(U' * Σy)
  @test A_mul_UUt_mul_B(Σx, U, Σy) ≈ A

  # asymmetric A * X * B
  X = randn(100, 50)
  Σx = Symmetric(cov(X))
  Y = randn(100, 50)
  Σy = Symmetric(cov(Y))
  Δ = sprand(50, 50, 0.1)
  A = Σx * Δ * Σy

  Δs = SparseIterate( Δ[:] )
  R = zeros(50, 50)
  for col=1:50
      for row=1:50
          R[row, col] = A_mul_X_mul_B_rc(Σx, Δs, Σy, row, col)
      end
  end

  @test R ≈ A
end

end


text:
module MathLibSignatures
using ..Concepts

export MathematicalObject,
  Cone,
  Interval,
  SemidefiniteCone,
  ClosedInterval

abstract type MathematicalObject end
abstract type Cone     <: MathematicalObject  end
abstract type Interval <: MathematicalObject  end

struct SemidefiniteCone <: Cone
  rank::Optional{Int64}
  function SemidefiniteCone(; rank::Optional{Int64} = nothing)
    if rank == nothing
      return new(nothing)
    end
    if rank <= 0
      @warn("Rank should never be below zero. Check again.")
      throw(MethodError())
    end
    return new(rank)
  end
end

struct ClosedInterval{T} <: Interval where T<:Real
  ll::Optional{T}
  rr::Optional{T}

  function ClosedInterval{T}(ll::T, rr::T) where T<:Real
    if ll > rr
      @warn("Right end point of the interval shoud be bigger than left, return an empty interval instead.")
      return new{T}(nothing, nothing)
    end
    return new{T}(ll, rr)
  end
end

function project() end

function project!() end

function scope_test()
  @info("inside MathLib.Concepts")
end

end


text:
using Roots

H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

function YfindNP(Q,P,H)
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+( 1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))

        Y_find1(y)=Hamil(0,y,Q,P)-H
    try
        Y=find_zeros(Y_find1,0,10, maxeval=100,maxfnevals=300,tol=1e-15)
     catch
        Y=zeros(0)
    end
end
# function YfindNP(Q,P,H)
# Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+( 1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
#
#         Y_find1(y)=Hamil(0,y,Q,P)-H
#     try
#         Y=find_zero(Y_find1,0.01, maxeval=100,maxfnevals=300,tol=1e-15)
#      catch
#         Y=zeros(0)
#     end
# end


text:
function _commits_behind_remote(repodir::String, url::String, remote_branch::String; deep = 10)
    _fetch(repodir, url; verbose = false) || return -1
    clist = _list_commits(repodir, remote_branch; count = deep)
    ch = GitLinks._HEAD_hash(repodir)
    chidx = findfirst(isequal(ch), clist)
    return isnothing(chidx) ? -1 : chidx - 1
end
_commits_behind_remote(gl::GitLink, remote_branch::String; kwargs...) =
    _commits_behind_remote(repo_dir(gl), remote_url(gl), remote_branch; kwargs...)


text:
function test_generator(process)
    gen = DataGenerator(
        process,
        batch_size=10,
        x_context=Uniform(-1, 3),
        x_target=Uniform(5, 10),
        num_context=DiscreteUniform(0, 5),
        num_target=DiscreteUniform(8, 8),
        σ²=1e-2
    )

    epoch = gen(500)  # `num_batches = 500`.

    # Check `num_batches`.
    @test length(epoch) == 500

    # Check minimum numbers of context points.
    @test minimum([size(batch[1], 1) for batch in epoch]) == 0
    @test minimum([size(batch[2], 1) for batch in epoch]) == 0

    # Check maximum numbers of context points.
    @test maximum([size(batch[1], 1) for batch in epoch]) == 5
    @test maximum([size(batch[2], 1) for batch in epoch]) == 5

    # Check the numbers of target points.
    @test all([size(batch[3], 1) == 8 for batch in epoch])
    @test all([size(batch[4], 1) == 8 for batch in epoch])

    for batch in epoch
        # Check `x_context`.
        if !isempty(batch[1])
            @test minimum(batch[1]) >= -1
            @test maximum(batch[1]) <= 3
        end

        # Check `x_target`.
        if !isempty(batch[3])
            @test minimum(batch[3]) >= 5
            @test maximum(batch[3]) <= 10
        end

        # Check `batch_size`.
        for x in batch
            @test size(x, 3) == 10
        end

        # Check consistency of context and target set.
        @test size(batch[1], 1) == size(batch[2], 1)
        @test size(batch[3], 1) == size(batch[4], 1)
    end
end

@testset "data.jl" begin
    @testset "UniformUnion" begin
        n = 10000
        dist = UniformUnion(Uniform(1, 2), Uniform(3, 5))
        x = rand(dist, n)

        # Check that the probabilities are calculated correctly.
        @test sum(x .<= 2) ≈ n / 3 rtol=5e-2
        @test sum(x .>= 3) ≈ 2n / 3 rtol=5e-2
    end
    @testset "Generate from sawtooth" begin
        test_generator(Sawtooth())
    end
    @testset "Generate from a Bayesian ConvNP" begin
        test_generator(BayesianConvNP())
    end
    @testset "Generate from a mixture" begin
        test_generator(Mixture(Sawtooth(), BayesianConvNP()))
    end
end


text:
@testset "ChainedVector" begin

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    @test x == 1:10
    @test length(x) == 10
    @test Base.IndexStyle(x) == Base.IndexLinear()

    x[1] = 0
    x[end] = 11
    @test x[1] == 0
    @test x[end] == 11

    @test copy(x) == x
    empty!(x)
    @test length(x) == 0
    @test copy(x) == x

    @test_throws ArgumentError resize!(x, -1)
    resize!(x, 10)
    @test length(x) == 10
    resize!(x, 20)
    @test length(x) == 20
    resize!(x, 15)
    @test length(x) == 15
    resize!(x, 5)
    @test length(x) == 5

    push!(x, 1)
    @test x[end] == 1
    empty!(x)
    push!(x, 1)
    @test x[1] == x[end] == 1

    pushfirst!(x, 2)
    @test x[1] == 2
    empty!(x)
    pushfirst!(x, 2)
    @test x[1] == x[end] == 2
    pushfirst!(x, 3)
    @test x[1] == 3

    @test pop!(x) == 2
    @test popfirst!(x) == 3
    @test isempty(x)
    @test_throws ArgumentError pop!(x)
    @test_throws ArgumentError popfirst!(x)

    @test_throws BoundsError insert!(x, 0, 1)
    @test_throws BoundsError insert!(x, 2, 1)
    insert!(x, 1, 1)
    @test x[1] == 1
    insert!(x, 1, 2)
    @test x[1] == 2

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = ChainedVector([[11,12,13], [14,15,16], [17,18,19,20]])

    z = vcat(x, y)
    @test length(z) == 20
    @test z == 1:20

    z = ChainedVector([[21,22,23], [24,25,26], [27,28,29,30]])
    a = vcat(x, y, z)
    @test length(a) == 30
    @test a == 1:30

    empty!(x)
    z = vcat(x, y)
    @test z == y

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    append!(x, y)
    @test length(x) == 20
    @test x == 1:20

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    append!(x, collect(11:20))
    @test length(x) == 20
    @test x == 1:20

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    append!(x, 11:20)
    @test length(x) == 20
    @test x == 1:20

    empty!(x)
    append!(x, y)
    @test x == y

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = ChainedVector([[11,12,13], [14,15,16], [17,18,19,20]])

    prepend!(y, x)
    @test length(y) == 20
    @test y == 1:20

    y = ChainedVector([[11,12,13], [14,15,16], [17,18,19,20]])
    prepend!(y, collect(1:10))
    @test length(y) == 20
    @test y == 1:20

    y = ChainedVector([[11,12,13], [14,15,16], [17,18,19,20]])
    prepend!(y, 1:10)
    @test length(y) == 20
    @test y == 1:20

    empty!(y)
    prepend!(y, x)
    @test y == x

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    deleteat!(x, 1)
    @test x[1] == 2
    deleteat!(x, 1:4)
    @test x == 6:10
    deleteat!(x, [2, 4])
    @test x == [6, 8, 10]

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    deleteat!(x, 1)
    deleteat!(x, 1)
    deleteat!(x, 1)
    @test length(x.arrays) == 2

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    deleteat!(x, 4)
    deleteat!(x, 4)
    deleteat!(x, 4)
    @test length(x.arrays) == 2

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = 0
    for i in x
        y += i
    end
    @test y == 55

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    b = [true, false, false, false, false, false, false, false, false, false]
    deleteat!(x, b)
    @test x[1] == 2
    @test length(x) == 9

    deleteat!(x, Int[])
    @test length(x) == 9

    #30
    x = ChainedVector([Vector{String}(undef, 3), ["hey", "ho"]])
    @test !isassigned(x, 1)
    @test isassigned(x, 4)

    #36
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9], [10,11,12]])
    deleteat!(x, [1, 2, 10, 11])
    @test x == [3, 4, 5, 6, 7, 8, 9, 12]

    #38
    A = ChainedVector([collect(((i - 1) * 153 + 1):(i * 153 + 1)) for i = 1:16])
    inds = collect(1:1883)
    deleteat!(A, inds)
    @test length(A) == 581

    # similar
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = similar(x, length(x))
    @test length(x) == length(y)
    @test eltype(x) == eltype(y)
    y = similar(x, Union{Missing, eltype(x)}, 2 * length(x))
    @test 2 * length(x) == length(y)
    @test eltype(y) == Union{Missing, eltype(x)}
    y = similar(x, 1)
    @test length(y) == 1
    y = similar(x, 0)
    @test length(y) == 0
    @test eltype(x) == eltype(y)
    y2 = similar(y, 10)
    @test length(y2) == 10
    @test eltype(y2) == eltype(x)

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    @test 10 in x
    @test !(11 in x)
    cnt = Ref(0)
    foreach(x -> cnt[] += 1, x)
    @test cnt[] == length(x)
    y = map(identity, x)
    @test x == y
    y = map(x -> x + 1, x)
    @test all(x -> x[1] == x[2] + 1, zip(y, x))

    # https://github.com/JuliaData/SentinelArrays.jl/issues/57
    cx1 = ChainedVector([[1, 2], [3]])
    cx2 = ChainedVector([[1.1, 2.2], [3.3]])
    @test ChainedVector([cx1, cx2]) isa ChainedVector{Float64, <:ChainedVector{Float64}}

    x = ChainedVector([[1], [2], [3]])
    y = map(v -> v == 1 ? missing : v, x)
    @test y isa ChainedVector{Union{Missing,Int}}
    @test isequal(y, ChainedVector([[missing], [2], [3]]))

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = map(v -> v == 1 ? missing : v, x)
    @test y isa ChainedVector{Union{Missing,Int}}
    @test isequal(y, ChainedVector([[missing,2,3], [4,5,6], [7,8,9,10]]))

    x = ChainedVector(Vector{Float64}[])
    y = map(v -> v > 1, x)
    @test y isa ChainedVector{Bool}

    # map!
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = copy(x)
    map!(x -> x + 1, y, x)
    @test all(x -> x[1] + 1 == x[2], zip(x, y))
    map!(x -> x + 1, x, y)
    @test all(x -> x[1] + 1 == x[2], zip(y, x))
    map!(x -> 1, x, x)
    @test all(x -> x == 1, x)
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = ChainedVector([[1,2,3], [4,5,6,7], [8,9,10]])
    map!(x -> x + 1, x, y)
    @test all(x -> x[1] + 1 == x[2], zip(y, x))

    # reductions
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    @test any(x -> iseven(x), x)
    @test any(map(x -> iseven(x), x))
    @test !all(x -> iseven(x), x)
    @test !all(map(x -> iseven(x), x))
    @test reduce(+, x) == 55
    @test foldl(+, x) == 55
    @test foldr(+, x) == 55
    @test mapreduce(x -> x + 1, +, x) == 65
    @test mapfoldl(x -> x + 1, +, x) == 65
    @test mapfoldr(x -> x + 1, +, x) == 65
    @test count(x -> iseven(x), x) == 5
    @test count(map(x -> iseven(x), x)) == 5
    @test extrema(x) == (1, 10)
    @test extrema(x -> x + 1, x) == (2, 11)
    @test minimum(x) == 1
    @test minimum(x -> x + 1, x) == 2
    @test maximum(x) == 10
    @test maximum(x -> x + 1, x) == 11

    # finds
    @test findmax(x) == (10, 10)
    @test findmax(x -> x + 1, x) == (11, 10)
    @test findmin(x) == (1, 1)
    @test findmin(x -> x + 1, x) == (2, 1)
    @test argmax(x) == 10
    @test argmax(x -> x + 1, x) == 10
    @test argmin(x) == 1
    @test argmin(x -> x + 1, x) == 1
    @test argmax(abs, ChainedVector([-10:5])) == -10
    @test findfirst(iseven, x) == 2
    @test findfirst(x -> x == 10, x) == 10
    @test findfirst(x -> x == 11, x) === nothing
    @test findfirst(map(iseven, x)) == 2
    @test findlast(iseven, x) == 10
    @test findlast(x -> x == 1, x) == 1
    @test findlast(x -> x == 11, x) === nothing
    @test findlast(map(iseven, x)) == 10
    @test findnext(iseven, x, 2) == 2
    @test findnext(map(iseven, x), 2) == 2
    @test findnext(iseven, x, 3) == 4
    @test findnext(x -> x == 11, x, 3) === nothing
    @test findprev(iseven, x, 10) == 10
    @test findprev(map(iseven, x), 10) == 10
    @test findprev(iseven, x, 9) == 8
    @test findprev(x -> x == 11, x, 9) === nothing
    @test findall(map(iseven, x)) == [2, 4, 6, 8, 10]
    @test findall(iseven, x) == [2, 4, 6, 8, 10]
    @test filter(iseven, x) == [2, 4, 6, 8, 10]
    filter!(iseven, x)
    @test x == [2, 4, 6, 8, 10]
    @test length(x) == 5
    @test replace(iseven, x) == map(iseven, x)
    @test replace!(iseven, copy(x)) == replace!(iseven, x)
    @test replace(x, 1 => 2) == map(x -> x == 1 ? 2 : x, x)
    @test map(x -> x == 1 ? 2 : x, x) == replace!(x, 1 => 2)
    @test all(==(2), x)
    @test length(x) == 5

    x = ChainedVector(Vector{Float64}[])
    @test !any(x -> iseven(x), x)
    @test !any(map(x -> iseven(x), x))
    @test all(x -> iseven(x), x)
    @test all(map(x -> iseven(x), x))
    @test_throws Union{ArgumentError,MethodError} reduce(+, x)
    @test_throws Union{ArgumentError,MethodError} foldl(+, x)
    @test_throws Union{ArgumentError,MethodError} foldr(+, x)
    @test_throws Union{ArgumentError,MethodError} mapreduce(x -> x + 1, +, x)
    @test_throws Union{ArgumentError,MethodError} mapfoldl(x -> x + 1, +, x)
    @test_throws Union{ArgumentError,MethodError} mapfoldr(x -> x + 1, +, x)
    @test count(x -> iseven(x), x) == 0
    @test count(map(x -> iseven(x), x)) == 0
    @test_throws ArgumentError extrema(x)
    @test_throws ArgumentError extrema(x -> x + 1, x)
    @test_throws Union{ArgumentError,MethodError} minimum(x)
    @test_throws Union{ArgumentError,MethodError} minimum(x -> x + 1, x)
    @test_throws Union{ArgumentError,MethodError} maximum(x)
    @test_throws Union{ArgumentError,MethodError} maximum(x -> x + 1, x)

    @test_throws ArgumentError findmax(x)
    @test_throws ArgumentError findmax(x -> x + 1, x)
    @test_throws ArgumentError findmin(x)
    @test_throws ArgumentError findmin(x -> x + 1, x)
    @test_throws ArgumentError argmax(x)
    @test_throws ArgumentError argmax(x -> x + 1, x)
    @test_throws ArgumentError argmin(x)
    @test_throws ArgumentError argmin(x -> x + 1, x)
    @test findfirst(iseven, x) === nothing
    @test findfirst(x -> x == 10, x) === nothing
    @test findfirst(x -> x == 11, x) === nothing
    @test findfirst(map(iseven, x)) === nothing
    @test findlast(iseven, x) === nothing
    @test findlast(x -> x == 1, x) === nothing
    @test findlast(x -> x == 11, x) === nothing
    @test findlast(map(iseven, x)) === nothing
    @test findnext(iseven, x, 2) === nothing
    @test findnext(iseven, x, 3) === nothing
    @test findnext(x -> x == 11, x, 3) === nothing
    @test findprev(iseven, x, 10) === nothing
    @test findprev(iseven, x, 9) === nothing
    @test findprev(x -> x == 11, x, 9) === nothing
    @test findall(map(iseven, x)) == Union{}[]
    @test findall(iseven, x) == Union{}[]
    @test filter(iseven, x) == Float64[]
    filter!(iseven, x)
    @test isempty(x)
    @test replace(iseven, x) == map(iseven, x)
    @test replace(x, 1 => 2) == map(x -> x == 1 ? 2.0 : x, x)
    @test replace!(iseven, copy(x)) == replace!(iseven, x)
    @test map(x -> x == 1 ? 2.0 : x, x) == replace!(x, 1 => 2)
    @test isempty(x)

    # copyto!
    # ChainedVector dest: doffs, soffs, n
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = copy(map(x -> x + 1, x))
    x2 = similar(x)
    copyto!(x2, y)
    @test x2 == y
    copyto!(x2, 2, y[1:end-1])
    @test x2[1] == y[1]
    @test x2[2:end] == y[1:end-1]
    @test_throws ArgumentError copyto!(x2, 1, y, 1, -1)
    @test_throws ArgumentError copyto!(x2, 1, y, 1, 30)
    copyto!(x2, 2, y, 9, 1)
    @test x2[2] == y[9]
    y2 = copy(x2);
    copyto!(x2, 1, copy(x), 1, 0)
    @test x2 == y2
    copyto!(x2, 2, y2, 2)
    @test x2[2:end] == y2[2:end]
    # ChainedVector src: doffs, soffs, n
    y2 = copy(y);
    copyto!(y2, x)
    @test y2 == x
    copyto!(y2, 2, view(x, 1:(length(x) - 1)))
    @test y2[2:end] == x[1:end-1]
    y3 = copy(y2)
    copyto!(y2, 1, x, 1, 0)
    @test y2 == y3
    @test_throws ArgumentError copyto!(y2, 1, x, 1, 30)
    @test_throws ArgumentError copyto!(y2, 1, x, 1, -1)
    copyto!(y2, 2, x, 1, 1)
    @test y2[2] == x[1]
    y2 = copy(y);
    copyto!(y2, 2, x, 2)
    @test y2[2:end] == x[2:end]
    # ChainedVector dest & src: doffs, soffs, n
    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    x2 = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    y = map(x -> x + 1, x)
    copyto!(x2, y)
    @test x2 == y
    x2 = map(identity, x)
    copyto!(x2, 2, view(y, 2:length(y)))
    @test x2[2:end] == y[2:end]
    x2 = map(identity, x)
    copyto!(x2, 2, y, 2)
    @test x2[2:end] == y[2:end]

    # https://github.com/JuliaData/SentinelArrays.jl/issues/45
    a = []
    append!(a, ChainedVector([[1, 2, 3]]))
    @test length(a) == 3

    # https://github.com/JuliaData/CSV.jl/issues/842
    @test size(eachindex(ChainedVector([["a"]]))) == (1,)

    x = [1, 2, 3]
    append!(x, ChainedVector([[1, 2, 3]]))
    @test length(x) == 6
    @test x == [1, 2, 3, 1, 2, 3]

    x = ChainedVector([[1,2,3], [4,5,6], [7,8,9,10]])
    for (i, idx) in enumerate(eachindex(x))
       @test i == idx
    end

    i1 = first(eachindex(x))
    @test convert(Int,  i1) isa Int
    @test convert(Int8,  i1) isa Int8
    @test i1 == 1
    @test i1 < 2
    @test isless(i1, 2)
end

@testset "iteration protocol on ChainedVector" begin
    for len in 0:6
        cv = ChainedVector([1:len])
        @test length(cv) == len
        c = 0
        for (i, v) in enumerate(cv)
            c += 1
            @test i == v
        end
        @test c == len
        for j in 0:len
            cv = ChainedVector([1:j, j+1:len])
            @test length(cv) == len
            c = 0
            for (i, v) in enumerate(cv)
                c += 1
                @test i == v
            end
            @test c == len
            for k in j:len
                cv = ChainedVector([1:j, j+1:k, k+1:len])
                @test length(cv) == len
                c = 0
                for (i, v) in enumerate(cv)
                    c += 1
                    @test i == v
                end
                @test c == len

                for l in k:len
                    cv = ChainedVector([1:j, j+1:k, k+1:l, l+1:len])
                    @test length(cv) == len
                    c = 0
                    for (i, v) in enumerate(cv)
                        c += 1
                        @test i == v
                    end
                    @test c == len
                end
            end
        end
    end
    x = ChainedVector([collect(1:i) for i = 10:100])
    y = copy(x)
    for (a, b) in zip(x, y)
        @test a == b
    end
    for (aidx, bidx) in zip(eachindex(x), eachindex(y))
        # getindex w/ custom ChainedVectorIndex
        @test x[aidx] == y[bidx]
        @test hash(aidx) == hash(bidx)
    end
    for (i, idx) in enumerate(eachindex(x))
        # setindex! w/ custom ChainedVectorIndex
        x[idx] = i
    end
    @test x == 1:length(x)
    ei = eachindex(x)
    @test length(ei) == length(x)
    @test eltype(ei) <: SentinelArrays.ChainedVectorIndex
    idx = collect(Iterators.take(ei, 20))
    @test x[1:20] == x[idx]
    x = ChainedVector(Vector{Int}[])
    @test isempty(x)
    @test length(x) == 0
    cnt = 0
    for y in x
        cnt += 1
    end
    @test cnt == 0
    for i in eachindex(x)
        cnt += x[i]
    end
    @test cnt == 0
end

@testset "empty ChainedVector" begin
    x = ChainedVector([[]])
    @test isempty(x)
    @test length(x) == 0
    @test_throws BoundsError x[1]
    @test_throws BoundsError x[end]

    push!(x, 1)
    @test x[end] == 1

    empty!(x)
    append!(x, [2])
    @test x[end] == 2
end


text:
using Documenter, Turf

makedocs(sitename="Turf.jl", doctest = false, pages=["Home" => ["index.md"],
    "Manual" => ["Getting Started" => "getting-started.md",
    "Examples" => "examples.md",
    hide("points-inside-polygon.md"), hide("simplify.md"),
    hide("nearest-point.md"), hide("grids.md"),
    "Methods" => "methods.md"]], format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"))

deploydocs(
    repo = "github.com/philoez98/Turf.jl.git",
)


text:
include("src/main.jl")

main(ARGS)


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Asynchronously reads the IO and buffers the read content. When end marker
# (GastonDone, which is stored in `gmarker` to account for Unix/Windows
# variability in line endings) is found in the content, sends everything
# between start (GastonBegin) and end markers to the returned channel.
# In case of timeout sends :timeout; in case of end of file, sends :eof.

function async_reader(io::IO, timeout_sec)::Channel
    ch = Channel(1)
    task = @async begin
        reader_task = current_task()
        function timeout_cb(timer)
            put!(ch, :timeout)
            Base.throwto(reader_task, InterruptException())
        end

        buf = ""
        while (match_done = findfirst(gmarker_done, buf)) == nothing
            timeout = Timer(timeout_cb, timeout_sec)
            data = String(readavailable(io))
            if data == ""; put!(ch, :eof); return; end
            timeout_sec > 0 && close(timeout) # Cancel the timeout
            buf *= data
        end
        match_begin = findfirst(gmarker_start, buf)
        start = (match_begin != nothing) ? last(match_begin)+1 : 1
        put!(ch, buf[start:first(match_done)-1])
    end
    bind(ch, task)
    return ch
end

# llplot() is our workhorse plotting function
function llplot(fig::Figure;print=false)
    global gnuplot_state

    # if figure has no data, stop here
    if isempty(fig)
        return
    end

    # In order to capture all output produced by our plotting commands
    # (error messages and figure text in case of text terminals), we
    # send "marks" to the stdout and stderr streams before the first
    # and after the last command. Everything between these marks will
    # be returned by our async_readers.
    gnuplot_send("\nreset session\n")
    gnuplot_send("set print \"-\"") # Redirect print to stdout
    gnuplot_send("print \"GastonBegin\"")
    gnuplot_send("printerr \"GastonBegin\"")

    # Build terminal setup string and send it to gnuplot
    gnuplot_send(termstring(fig,print))

    # Datafile filename. This is where we store the coordinates to plot.
    # This file is then read by gnuplot to do the actual plotting. One file
    # per figure handle is used; this avoids polutting /tmp with too many files.
    filename = joinpath(tempdir(),"gaston-$(tmpprefix)-$(fig.handle)")
    f = open(filename,"w")

    # Send appropriate coordinates and data to gnuplot, depending on
    # whether we are doing 2-d, 3-d or image plots.

    # 2-d plot: Z is empty or plostyle is {,rgb}image
    if isempty(fig.curves[1].Z) ||
        fig.curves[1].conf.plotstyle == "image" ||
        fig.curves[1].conf.plotstyle == "rgbimage"
        # create data file
        for i in fig.curves
            ps = i.conf.plotstyle
            if ps == "errorbars" || ps == "errorlines"
                if isempty(i.E.yhigh)
                    # ydelta (single error coordinate)
                    writedlm(f,[i.x i.y i.E.ylow],' ')
                else
                    # ylow, yhigh (double error coordinate)
                    writedlm(f,[i.x i.y i.E.ylow i.E.yhigh],' ')
                end
            elseif ps == "financebars"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x),5)
                # output matrix
                for col = 1:length(i.x)
                    tmparr[col,1] = i.x[col]
                    tmparr[col,2] = i.F.open[col]
                    tmparr[col,3] = i.F.low[col]
                    tmparr[col,4] = i.F.high[col]
                    tmparr[col,5] = i.F.close[col]
                end
                writedlm(f,tmparr,' ')
            elseif ps == "image"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x)*length(i.y),3)
                tmparr_row_index = 1  # index into tmparr row
                # output matrix
                for row = 1:length(i.y)
                    x = length(i.x)
                    for col = 1:length(i.x)
                        tmparr[tmparr_row_index,1] = i.x[col]
                        tmparr[tmparr_row_index,2] = i.y[row]
                        tmparr[tmparr_row_index,3] = i.Z[row,col]
                        tmparr_row_index = tmparr_row_index+1
                        x = x-1
                    end
                end
                writedlm(f,tmparr,' ')
            elseif ps == "rgbimage"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x)*length(i.y), 5)
                tmparr_row_index = 1
                # output matrix
                for col = 1:length(i.x)
                    y = length(i.y)
                    for row = 1:length(i.y)
                        tmparr[tmparr_row_index,1] = i.x[col]
                        tmparr[tmparr_row_index,2] = i.y[row]
                        tmparr[tmparr_row_index,3] = i.Z[y,col,1]
                        tmparr[tmparr_row_index,4] = i.Z[y,col,2]
                        tmparr[tmparr_row_index,5] = i.Z[y,col,3]
                        tmparr_row_index = tmparr_row_index+1
                        y = y-1
                    end
                end
                writedlm(f,tmparr,' ')
            else
                writedlm(f,[i.x i.y],' ')
            end
            write(f,"\n\n")
        end
        close(f)

        # Send gnuplot commands.
        # Build figure configuration to gnuplot
        gnuplot_send_fig_config(fig.axes)
        # Send user command to gnuplot
        !isempty(fig.gpcom) && gnuplot_send(fig.gpcom)
        # send plot command to gnuplot
        gnuplot_send(linestr(fig.curves, "plot", filename))

    # 3-d plot: Z is not empty and plotstyle is not {,rgb}image
    elseif !isempty(fig.curves[1].Z) &&
            fig.curves[1].conf.plotstyle != "image" &&
            fig.curves[1].conf.plotstyle != "rgbimage"
        # create data file
        if isa(fig.curves[1].Z, Matrix)  # surface plot
            for i in fig.curves
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(1, 3)
                tmparr_row_index = 1
                for row in 1:length(i.x)
                    for col in 1:length(i.y)
                        tmparr[1,1] = i.x[row]
                        tmparr[1,2] = i.y[col]
                        tmparr[1,3] = i.Z[row,col]
                        writedlm(f,tmparr,' ')
                    end
                    write(f,"\n")
                end
                write(f,"\n\n")
            end
        else # scatter plot
            for i in fig.curves
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(1, 3)
                tmparr_row_index = 1
                for k in 1:length(i.x)
                    tmparr[1,1] = i.x[k]
                    tmparr[1,2] = i.y[k]
                    tmparr[1,3] = i.Z[k]
                    writedlm(f,tmparr,' ')
                end
                write(f,"\n")
            end
            write(f,"\n\n")
        end

        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(fig.axes)
        # Send user command to gnuplot
        !isempty(fig.gpcom) && gnuplot_send(fig.gpcom)
        # send command to gnuplot
        gnuplot_send(linestr(fig.curves, "splot",filename))
    end
    # Close output files, if any
    gnuplot_send("set output")

    # Make sure gnuplot is done; if terminal is text, read data
    # reset error handling
    err = ""
    gnuplot_state.gp_lasterror = err
    gnuplot_state.gp_error = false

    gnuplot_send("print \"GastonDone\"")
    gnuplot_send("printerr \"GastonDone\"")

    # Start reading gnuplot's streams in "background"
    ch_out = async_reader(P.gstdout, config[:timeouts][:stdout_timeout])
    ch_err = async_reader(P.gstderr, config[:timeouts][:stderr_timeout])

    out = take!(ch_out)
    err = take!(ch_err)

    out === :timeout && error("Gnuplot is taking too long to respond.")
    out === :eof     && error("Gnuplot crashed")

    # We don't care about stderr timeouts.
    err === :timeout && (err = "")
    err === :eof     && (err = "")

    # check for errors while plotting
    if err != ""
        gnuplot_state.gp_lasterror = err
        gnuplot_state.gp_error = true
        @warn("Gnuplot returned an error message:\n  $err)")
    end

    # if there was no error and text terminal, read all data from stdout
    if err == ""
        if (config[:term][:terminal] ∈ term_text) || (config[:mode] == "ijulia")
            fig.svg = out
        end
    end

    return nothing

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

using IntervalSets

struct apb_test <: BAT.AbstractParamBounds
    nparams::Integer
end

BAT.nparams(a::apb_test) = a.nparams
BAT.unsafe_intersect(a::apb_test, b::apb_test) = true

@testset "parameter bounds" begin
    @testset "BAT.oob" begin
        @test BAT.isoob(BAT.oob(1.))
        @test BAT.isoob(1.) == false
        @test BAT.isoob(BAT.oob(1))
        @test BAT.isoob(1) == false
        @test BAT.isoob(BAT.oob(ones(Float64,2,2)))
        @test BAT.isoob(ones(Float64,2,2)) == false
    end

    @testset "BAT.apply_bounds" begin
        @test BAT.apply_bounds(+0.3, -1, 2, hard_bounds) ≈ +0.3
        @test BAT.apply_bounds(-0.3, -1, 2, reflective_bounds) ≈ -0.3
        @test BAT.apply_bounds(+1.3, -1, 2, cyclic_bounds) ≈ +1.3

        @test isnan(BAT.apply_bounds(-1.3, -1, 2, hard_bounds))
        @test isnan(BAT.apply_bounds(2.3, -1, 2, hard_bounds))
        @test BAT.apply_bounds(-1.3, -1, 2, hard_bounds, 5) ≈ 5
        @test BAT.apply_bounds(2.3, -1, 2, hard_bounds, 5) ≈ 5

        @test BAT.apply_bounds(-1.3, -1, 2, reflective_bounds) ≈ -0.7
        @test BAT.apply_bounds(-4.3, -1, 2, reflective_bounds) ≈ +1.7
        @test BAT.apply_bounds(-7.3, -1, 2, reflective_bounds) ≈ -0.7
        @test BAT.apply_bounds(+2.3, -1, 2, reflective_bounds) ≈ +1.7
        @test BAT.apply_bounds(+5.3, -1, 2, reflective_bounds) ≈ -0.7
        @test BAT.apply_bounds(+8.3, -1, 2, reflective_bounds) ≈ +1.7

        @test BAT.apply_bounds(-1.3, -1, 2, cyclic_bounds) ≈ +1.7
        @test BAT.apply_bounds(-4.3, -1, 2, cyclic_bounds) ≈ +1.7
        @test BAT.apply_bounds(-7.3, -1, 2, cyclic_bounds) ≈ +1.7
        @test BAT.apply_bounds(+2.3, -1, 2, cyclic_bounds) ≈ -0.7
        @test BAT.apply_bounds(+5.3, -1, 2, cyclic_bounds) ≈ -0.7
        @test BAT.apply_bounds(+8.3, -1, 2, cyclic_bounds) ≈ -0.7

        @test typeof(@inferred BAT.apply_bounds(+0.3, -1, 2, hard_bounds)) == Float64
        @test typeof(@inferred BAT.apply_bounds(+0.3f0, -1, 2, hard_bounds)) == Float32
        @test typeof(@inferred BAT.apply_bounds(+0.3f0, -1, 2.0, hard_bounds)) == Float64

        @test BAT.apply_bounds(+5.3, ClosedInterval(-1, 2), reflective_bounds) ≈ -0.7
    end

    @testset "BAT.NoParamBounds" begin
        n = 2
        @test typeof(@inferred BAT.NoParamBounds(n)) == BAT.NoParamBounds

        params = [-1000., 1000]
        uparams = BAT.NoParamBounds(n)
        @test nparams(uparams) == n
        @test params in uparams
        @test in( hcat(params, params), uparams, 1)

        @test BAT.apply_bounds!(params, uparams) == params
    end

    @testset "BAT.HyperRectBounds" begin
        @test typeof(@inferred BAT.HyperRectBounds([-1., 0.5], [2.,1], BAT.hard_bounds)) <: ParamVolumeBounds{Float64, HyperRectVolume{Float64}}
        @test_throws ArgumentError BAT.HyperRectBounds([-1.], [2.,1], [hard_bounds, reflective_bounds])
        @test_throws ArgumentError BAT.HyperRectBounds([-1., 2.], [2.,1], [hard_bounds, reflective_bounds])

        hyperRectBounds = @inferred BAT.HyperRectBounds([-1., -1.], [0.5,1], [hard_bounds, hard_bounds])
        @test nparams(hyperRectBounds) == 2
        @test [0.0, 0.0] in hyperRectBounds
        @test ([0.5, 2] in hyperRectBounds) == false

        hyperRectBounds = @inferred BAT.HyperRectBounds([ClosedInterval(-1.,0.5),
            ClosedInterval(-1.,1.)], [hard_bounds, hard_bounds])
        @test nparams(hyperRectBounds) == 2
        @test [0.0, 0.0] in hyperRectBounds
        @test ([0.5, 2] in hyperRectBounds) == false

        @test BAT.apply_bounds!([0.3, -4.3, -7.3], BAT.HyperRectBounds([-1.,-1,-1], [2.,2,2],
            [hard_bounds, reflective_bounds, cyclic_bounds])) ≈ [+0.3, 1.7, 1.7]

        @test BAT.apply_bounds!([0.3 0.3 0.3; 0.3 -7.3 +8.3; 0.3 -7.3 +8.3],
            BAT.HyperRectBounds([-1., -1., -1], [2., 2.,2.], [hard_bounds, reflective_bounds, cyclic_bounds])) ≈ [+0.3 0.3 0.3;0.3 -0.7 1.7;0.3 1.7 -0.7]

        @test BAT.isoob(BAT.apply_bounds!(rand!(MersenneTwister(7002), zeros(Float64, 2, 2)), hyperRectBounds))
        @test BAT.isoob(BAT.apply_bounds!(rand!(MersenneTwister(7002),
            BAT.spatialvolume(hyperRectBounds), zeros(Float64, 2, 2)), hyperRectBounds)) == false

        @test BAT.isoob(BAT.apply_bounds!([0,2.],hyperRectBounds))
        @test BAT.apply_bounds!([0,2.],hyperRectBounds,false) ≈ [0,2.]
    end

    @testset "similar" begin
        a = @inferred BAT.HyperRectBounds([-1., 0.5], [2.,1], BAT.hard_bounds)
        c = @inferred similar(a)
        @test typeof(c) <:HyperRectBounds{Float64}
        @test typeof(c.vol) <: BAT.HyperRectVolume{Float64}
        @test typeof(c.bt) <: Array{BoundsType,1}
    end

    @testset "BAT.intersect" begin
        a = @inferred apb_test(3)
        b = @inferred apb_test(4)
        @test_throws ArgumentError intersect(a,b)
        b = @inferred apb_test(3)
        @test intersect(a,b)
        b = @inferred NoParamBounds(3)
        @test BAT.unsafe_intersect(b, NoParamBounds(3)) == b
        @test BAT.unsafe_intersect(a, NoParamBounds(3)) == a
        @test BAT.unsafe_intersect(NoParamBounds(3), a) == a

        hb = hard_bounds
        rb = reflective_bounds
        cb = cyclic_bounds

        @test intersect(hb, rb) == hard_bounds
        @test intersect(cb, hb) == hard_bounds
        @test intersect(rb, cb) == reflective_bounds
        @test intersect(cb, rb) == reflective_bounds
        @test intersect(cb, cb) == cyclic_bounds

        hyperRectBounds_a = @inferred BAT.HyperRectBounds([-1., -1.], [2.,3.],
            [hard_bounds, cyclic_bounds])
        hyperRectBounds_b = @inferred BAT.HyperRectBounds([-0.5, -0.5], [1.,2.],
            [cyclic_bounds, cyclic_bounds])

        res = @inferred intersect(hyperRectBounds_a, hyperRectBounds_b)
        @test res.vol.lo ≈ hyperRectBounds_b.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_b.vol.hi
        @test res.bt == hyperRectBounds_b.bt

        res = @inferred intersect(hyperRectBounds_b, hyperRectBounds_a)
        @test res.vol.lo ≈ hyperRectBounds_b.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_b.vol.hi
        @test res.bt == hyperRectBounds_b.bt

        hyperRectBounds_c = @inferred BAT.HyperRectBounds([-0.5, -0.5], [4.,4.],
            [cyclic_bounds, reflective_bounds])

        res = @inferred intersect(hyperRectBounds_a, hyperRectBounds_c)
        @test res.vol.lo ≈ hyperRectBounds_c.vol.lo
        @test res.vol.hi ≈ hyperRectBounds_a.vol.hi
        @test res.bt == [hard_bounds, reflective_bounds]

    end
end


text:
function mra(wt::WaveletTransform)
  D = fill(0.0, wt.L, wt.level, wt.series)
  S = fill(0.0, wt.L, wt.level, wt.series)

  # compute the details and smooths
  for j = 1:wt.level
    DWj = wt.W
    SVj = wt.V
    SWj = fill(0.0, wt.L, wt.level, wt.series)
    DVj = fill(0.0, wt.L, wt.level, wt.series)

    for k=j:-1:1
      for i=1:wt.series
        DVj[:,j,i] = eval(Expr(:call, symbol(string(typeof(wt),"Backward")), DWj[:,j,i], DVj[:,j,i], wt.filter, k))
        SVj[:,j,i] = eval(Expr(:call, symbol(string(typeof(wt),"Backward")), SWj[:,j,i], SVj[:,j,i], wt.filter, k))
      end

      DWj = fill(0.0, wt.L, wt.level, wt.series)
      SWj = DWj
    end

    D[:,j,:] = DVj[:,j,:]
    S[:,j,:] = SVj[:,j,:]
  end

  return (D, S)
end


text:
# import of Rmath functions

module RFunctions

using Compat
import Rmath: libRmath

### import macro

function _import_rmath(rname::Symbol, jname::Symbol, pargs)
    # C function names
    if rname == :norm
        dfun = Expr(:quote, "dnorm4")
        pfun = Expr(:quote, "pnorm5")
        qfun = Expr(:quote, "qnorm5")
        rfun = Expr(:quote, "rnorm")
    else
        dfun = Expr(:quote, string('d', rname))    # density
        pfun = Expr(:quote, string('p', rname))    # cumulative probability
        qfun = Expr(:quote, string('q', rname))    # quantile
        rfun = Expr(:quote, string('r', rname))    # random sampling
    end

    # Julia function names
    pdf = @compat Symbol(jname, "pdf")
    cdf = @compat Symbol(jname, "cdf")
    ccdf = @compat Symbol(jname, "ccdf")

    logpdf = @compat Symbol(jname, "logpdf")
    logcdf = @compat Symbol(jname, "logcdf")
    logccdf = @compat Symbol(jname, "logccdf")

    invcdf = @compat Symbol(jname, "invcdf")
    invccdf = @compat Symbol(jname, "invccdf")
    invlogcdf = @compat Symbol(jname, "invlogcdf")
    invlogccdf = @compat Symbol(jname, "invlogccdf")

    rand = @compat Symbol(jname, "rand")
    has_rand = true
    if rname == :nbeta || rname == :nf || rname == :nt
        has_rand = false
    end

    # arguments & argument types
    _pts = fill(:Cdouble, length(pargs))

    dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
    ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
    qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
    rtypes = Expr(:tuple, _pts...)

    pdecls = [Expr(:(::), ps, :Real) for ps in pargs] # [:(p1::Real), :(p2::Real), ...]

    # function implementation
    quote
        $pdf($(pdecls...), x::Real) =
            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)

        $logpdf($(pdecls...), x::Real) =
            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)

        $cdf($(pdecls...), x::Real) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)

        $ccdf($(pdecls...), x::Real) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 0)

        $logcdf($(pdecls...), x::Real) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 1)

        $logccdf($(pdecls...), x::Real) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 1)

        $invcdf($(pdecls...), q::Real) =
            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 1, 0)

        $invccdf($(pdecls...), q::Real) =
            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 0, 0)

        $invlogcdf($(pdecls...), lq::Real) =
            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 1, 1)

        $invlogccdf($(pdecls...), lq::Real) =
            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 0, 1)

        if $has_rand
            $rand($(pdecls...)) =
                ccall(($rfun, libRmath), Float64, $rtypes, $(pargs...))
        end
    end
end

macro import_rmath(rname, jname, pargs...)
    esc(_import_rmath(rname, jname, pargs))
end

### Import specific functions

@import_rmath beta beta α β
@import_rmath binom binom n p
@import_rmath chisq chisq k
@import_rmath f fdist d1 d2
@import_rmath gamma gamma α β
@import_rmath hyper hyper ms mf n
@import_rmath norm norm μ σ
@import_rmath nbinom nbinom r p
@import_rmath pois pois λ
@import_rmath t tdist k

@import_rmath nbeta nbeta α β λ
@import_rmath nchisq nchisq k λ
@import_rmath nf nfdist k1 k2 λ
@import_rmath nt ntdist k λ

end


text:
#=
sios = HTSLIB.open("/home/guo/haplox/Github/HTSLIB/test/data/100_sort.bam","rb","bam")

@show sios.handle
@show convert(Ptr{HTSLIB.HTSFile},sios.handle)
htsfl = unsafe_load(convert(Ptr{HTSLIB.HTSFile},sios.handle))
@show bytestring(htsfl.fn)

pidx = HTSLIB.sam_index_load(sios,"/home/guo/haplox/Github/HTSLIB/test/data/100_sort.bam")
phdr = sios.phdr
@show pidx,phdr

iters  = HTSLIB.sam_itr_querys(pidx,phdr,"chr1:65432520-65432550")
@show iters

b = HTSLIB.bam_init1()
i = 0
while HTSLIB.sam_itr_next!(sios.handle,iters,b)
    i+=1

end
=#
bios = HTSLIB.open("data/100_sort.bam","rb","bam")
htsfl = unsafe_load(convert(Ptr{HTSLIB.HTSFile},bios.handle))
pidx = HTSLIB.sam_index_load(bios,"data/100_sort.bam")
iters = HTSLIB.sam_itr_querys(pidx,bios.phdr,"chr1:65432520-65432550")
precord = HTSLIB.bam_init1()
while HTSLIB.sam_itr_next!(bios.handle,iters,precord)
    HTSLIB.sam_format!(bios.phdr,precord,bios.pkstr)
    kstr = HTSLIB.strptr(bios.pkstr)
    println(kstr)
end