text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related with the IERS (International Earth Orientation and
#   Reference Systems Service) EOP (Earth Orientation Parameters) data.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] IERS (2010). Transformation between the International Terrestrial
#       Reference System and the Geocentric Celestial Reference System. IERS
#       Technical Note No. 36, Chapter 5.
#
#   [2] ftp://hpiers.obspm.fr/eop-pc/models/uai2000.package
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export get_iers_eop, get_iers_eop_iau_1980, get_iers_eop_iau_2000A
export read_iers_eop
export deps_dpsi

################################################################################
#                                  Functions
################################################################################

"""
    get_iers_eop([data_type]; force_download = false)

Download and parse the IERS EOP C04 data. The data type is specified by
`data_type`. Supported values are:

- `Val(:IAU1980)`: Get IERS EOP C04 IAU1980 data.
- `Val(:IAU2000A)`: Get IERS EOP C04 IAU2000A data.

If `data_type` is omitted, then it defaults to `Val(:IAU1980)`.

The files are downloaded using the `RemoteFile` package with daily updates.
Hence, if one desires to force a download before the scheduled time, then set
the keyword `force_download` to `true`.

!!! note
    The files will be downloaded from the default URL. If the user want to use
    another one, then use the specialized functions
    [`get_iers_eop_iau_1980`](@ref), [`get_iers_eop_iau_2000A`](@ref)

See also: [`get_iers_eop_iau_1980`](@ref), [`get_iers_eop_iau_2000A`](@ref)

# Returns

A structure ([`EOPData_IAU1980`](@ref) or [`EOPData_IAU2000A`](@ref), depending
on `data_type`) with the interpolations of the EOP parameters. Notice that the
interpolation indexing is set to the Julian Day.

# Examples

```julia-repl
julia> eop = get_iers_eop()
[ Info: Downloading file 'EOP_IAU1980.TXT' from 'https://datacenter.iers.org/data/csv/finals.all.csv' with cURL.
EOPData_IAU1980:
     Data │ Timespan
 ─────────┼──────────────────────────────────────────────
        x │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
        y │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
  UT1-UTC │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
      LOD │ 1973-01-02T00:00:00 -- 2021-05-19T00:00:00
       dψ │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
       dϵ │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00

julia> eop = get_iers_eop(Val(:IAU2000A))
[ Info: Downloading file 'EOP_IAU2000A.TXT' from 'https://datacenter.iers.org/data/csv/finals2000A.all.csv' with cURL.
EOPData_IAU2000A:
     Data │ Timespan
 ─────────┼──────────────────────────────────────────────
        x │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
        y │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
  UT1-UTC │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
      LOD │ 1973-01-02T00:00:00 -- 2021-05-19T00:00:00
       dX │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
       dY │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
```
"""
get_iers_eop(;kwargs...) = get_iers_eop(Val(:IAU1980); kwargs...)

function get_iers_eop(::Val{:IAU1980}; force_download = false)
    return get_iers_eop_iau_1980(force_download = force_download)
end

function get_iers_eop(::Val{:IAU2000A}; force_download = false)
    return get_iers_eop_iau_2000A(force_download = force_download)
end

"""
    get_iers_eop_iau_1980(url::String = "https://datacenter.iers.org/data/csv/finals.all.csv"; force_download = false)

Get the IERS EOP C04 IAU1980 data from the URL `url`.

If `url` is omitted, then it defaults to
https://datacenter.iers.org/data/csv/finals.all.csv

The file is downloaded using the `RemoteFile` package with daily updates. Hence,
if one desires to force a download before the scheduled time, then set the
keyword `force_download` to `true`.

!!! note
    The interpolation of every field in [`EOPData_IAU1980`](@ref) between two
    points in the grid is linear. If extrapolation is needed, then if will use
    the nearest value (flat extrapolation).

See also: [`get_iers_eop`](@ref), [`get_iers_eop_iau_2000A`](@ref)

# Returns

The structure [`EOPData_IAU1980`](@ref) with the interpolations of the EOP
parameters. Notice that the interpolation indexing is set to the Julian Day.

# Examples

```julia-repl
julia> eop = get_iers_eop_iau1980()
[ Info: Downloading file 'EOP_IAU1980.TXT' from 'https://datacenter.iers.org/data/csv/finals.all.csv' with cURL.
EOPData_IAU1980:
     Data │ Timespan
 ─────────┼──────────────────────────────────────────────
        x │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
        y │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
  UT1-UTC │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
      LOD │ 1973-01-02T00:00:00 -- 2021-05-19T00:00:00
       dψ │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
       dϵ │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
```
"""
function get_iers_eop_iau_1980(
    url::String = "https://datacenter.iers.org/data/csv/finals.all.csv";
    force_download = false
)
    _eop_iau1980 = @RemoteFile(
        @eval($url),
        file="EOP_IAU1980.TXT",
        updates=:daily
    )

    download(_eop_iau1980; force = force_download, force_update = true)

    # Parse the data removing the header.
    eop, ~ = readdlm(path(_eop_iau1980), ';'; header = true)

    # Return the parsed EOP data.
    return _parse_iers_eop_iau_1980(eop)
end

"""
    get_iers_eop_iau_2000A(url::String = "https://datacenter.iers.org/data/latestVersion/224_EOP_C04_14.62-NOW.IAU2000A224.txt"; force_download = false)

Get the IERS EOP C04 IAU2000A data from the URL `url`.

If `url` is omitted, then it defaults to
https://datacenter.iers.org/data/csv/finals2000A.all.csv.

The file is downloaded using the `RemoteFile` package with daily updates. Hence,
if one desires to force a download before the scheduled time, then set the
keyword `force_download` to `true`.

!!! note
    The interpolation of every field in [`EOPData_IAU2000A`](@ref) between two
    points in the grid is linear. If extrapolation is needed, then if will use
    the nearest value (flat extrapolation).

See also: [`get_iers_eop`](@ref), [`get_iers_eop_iau_1980`](@ref)

# Returns

The structure `EOPData_IAU2000A` with the interpolations of the EOP parameters.
Notice that the interpolation indexing is set to the Julian Day.

# Examples

```julia-repl
julia> eop = get_iers_eop(Val(:IAU2000A))
[ Info: Downloading file 'EOP_IAU2000A.TXT' from 'https://datacenter.iers.org/data/csv/finals2000A.all.csv' with cURL.
EOPData_IAU2000A:
     Data │ Timespan
 ─────────┼──────────────────────────────────────────────
        x │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
        y │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
  UT1-UTC │ 1973-01-02T00:00:00 -- 2022-05-28T00:00:00
      LOD │ 1973-01-02T00:00:00 -- 2021-05-19T00:00:00
       dX │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
       dY │ 1973-01-02T00:00:00 -- 2021-08-05T00:00:00
```
"""
function get_iers_eop_iau_2000A(
    url::String = "https://datacenter.iers.org/data/csv/finals2000A.all.csv";
    force_download = false
)
    _eop_iau2000A = @RemoteFile(
        @eval($url),
        file="EOP_IAU2000A.TXT",
        updates=:daily
    )

    download(_eop_iau2000A; force = force_download, force_update = true)

    # Parse the data removing the header.
    eop, ~ = readdlm(path(_eop_iau2000A), ';'; header = true)

    # Return the parsed EOP data.
    return _parse_iers_eop_iau_2000A(eop)
end

"""
    read_iers_eop(filename::String[, data_type])

Read IERS EOP Data from the file `filename`. The user must specify if the data
is related to the model IAU 1980 (`data_type = Val(:IAU1980)`), which is the
default, or to the model IAU 2000A (`data_type = Val(:IAU2000A)`).

!!! note
    The input file **must be exactly the same** as provided by IERS in CSV format.
    One can download it using the following commands:

    * IAU 1980

        curl -O https://datacenter.iers.org/data/csv/finals.all.csv
        wget https://datacenter.iers.org/data/csv/finals.all.csv

    * IAU 2000A

        curl -O https://datacenter.iers.org/data/csv/finals2000A.all.csv
        wget https://datacenter.iers.org/data/csv/finals2000A.all.csv

# Returns

A structure ([`EOPData_IAU1980`](@ref) or [`EOPData_IAU2000A`](@ref), depending
on `data_type`) with the interpolations of the EOP parameters. Notice that the
interpolation indexing is set to the Julian Day.
"""
read_iers_eop(filename::String) = read_iers_eop(filename, Val(:IAU1980))

function read_iers_eop(filename::String, ::Val{:IAU1980})
    eop, ~ = readdlm(filename, ';'; header = true)
    return _parse_iers_eop_iau_1980(eop)
end

function read_iers_eop(filename::String, ::Val{:IAU2000A})
    eop, ~ = readdlm(filename, ';'; header = true)
    return _parse_iers_eop_iau_2000A(eop)
end

################################################################################
#                                      IO
################################################################################

function show(io::IO, eop::EOPData_IAU1980)
    print(io, "EOPData_IAU1980")
    return nothing
end

function show(io::IO, mime::MIME"text/plain", eop::EOPData_IAU1980)
    # Check if IO has support for colors.
    color = get(io, :color, false)::Bool

    if color
        _b = string(Crayon(bold = true))
        _g = string(crayon"dark_gray")
        _r = string(Crayon(reset = true))
    else
        _b = ""
        _g = ""
        _r = ""
    end

    println(io, "EOPData_IAU1980:")
    println(io, _b, "     Data ", _g, "│ ", _r, _b, "Timespan", _r)
    println(io, _g, " ─────────┼──────────────────────────────────────────────", _r)
    println(io, _b, "        x ", _g, "│ ", _r, _get_iers_eop_timespan(eop.x))
    println(io, _b, "        y ", _g, "│ ", _r, _get_iers_eop_timespan(eop.y))
    println(io, _b, "  UT1-UTC ", _g, "│ ", _r, _get_iers_eop_timespan(eop.UT1_UTC))
    println(io, _b, "      LOD ", _g, "│ ", _r, _get_iers_eop_timespan(eop.LOD))
    println(io, _b, "       dψ ", _g, "│ ", _r, _get_iers_eop_timespan(eop.dPsi))
    print(io,   _b, "       dϵ ", _g, "│ ", _r, _get_iers_eop_timespan(eop.dEps))

    return nothing
end

function show(io::IO, eop::EOPData_IAU2000A)
    print(io, "EOPData_IAU2000A")
    return nothing
end

function show(io::IO, mime::MIME"text/plain", eop::EOPData_IAU2000A)
    # Check if IO has support for colors.
    color = get(io, :color, false)::Bool

    if color
        _b = string(Crayon(bold = true))
        _g = string(crayon"dark_gray")
        _r = string(Crayon(reset = true))
    else
        _b = ""
        _g = ""
        _r = ""
    end

    println(io, "EOPData_IAU2000A:")
    println(io, _b, "     Data ", _g, "│ ", _r, _b, "Timespan", _r)
    println(io, _g, " ─────────┼──────────────────────────────────────────────", _r)
    println(io, _b, "        x ", _g, "│ ", _r, _get_iers_eop_timespan(eop.x))
    println(io, _b, "        y ", _g, "│ ", _r, _get_iers_eop_timespan(eop.y))
    println(io, _b, "  UT1-UTC ", _g, "│ ", _r, _get_iers_eop_timespan(eop.UT1_UTC))
    println(io, _b, "      LOD ", _g, "│ ", _r, _get_iers_eop_timespan(eop.LOD))
    println(io, _b, "       dX ", _g, "│ ", _r, _get_iers_eop_timespan(eop.dX))
    print(io,   _b, "       dY ", _g, "│ ", _r, _get_iers_eop_timespan(eop.dY))

    return nothing
end

################################################################################
#                              Private Functions
################################################################################

# Parse the IERS EOP (IAU 1980) data.
function _parse_iers_eop_iau_1980(eop::Matrix)
    # Create the EOP Data structure by creating the interpolations.
    #
    # The interpolation will be linear between two points in the grid. The
    # extrapolation will be flat, considering the nearest point.
    knots::Vector{Float64} = Vector{Float64}(eop[:, 1] .+ 2400000.5)

    return EOPData_IAU1980(
        _create_iers_eop_interpolation(knots, eop[:, 6]),
        _create_iers_eop_interpolation(knots, eop[:, 8]),
        _create_iers_eop_interpolation(knots, eop[:, 11]),
        _create_iers_eop_interpolation(knots, eop[:, 13]),
        _create_iers_eop_interpolation(knots, eop[:, 16]),
        _create_iers_eop_interpolation(knots, eop[:, 18]),
        _create_iers_eop_interpolation(knots, eop[:, 7]),
        _create_iers_eop_interpolation(knots, eop[:, 9]),
        _create_iers_eop_interpolation(knots, eop[:, 12]),
        _create_iers_eop_interpolation(knots, eop[:, 14]),
        _create_iers_eop_interpolation(knots, eop[:, 17]),
        _create_iers_eop_interpolation(knots, eop[:, 19]),
    )
end

# Parse the IERS EOP (IAU 2000A) data.
function _parse_iers_eop_iau_2000A(eop::Matrix)
    # Create the EOP Data structure by creating the interpolations.
    #
    # The interpolation will be linear between two points in the grid. The
    # extrapolation will be flat, considering the nearest point.
    knots::Vector{Float64} = Vector{Float64}(eop[:, 1] .+ 2400000.5)

    EOPData_IAU2000A(
        _create_iers_eop_interpolation(knots, eop[:, 6]),
        _create_iers_eop_interpolation(knots, eop[:, 8]),
        _create_iers_eop_interpolation(knots, eop[:, 11]),
        _create_iers_eop_interpolation(knots, eop[:, 13]),
        _create_iers_eop_interpolation(knots, eop[:, 20]),
        _create_iers_eop_interpolation(knots, eop[:, 22]),
        _create_iers_eop_interpolation(knots, eop[:, 7]),
        _create_iers_eop_interpolation(knots, eop[:, 9]),
        _create_iers_eop_interpolation(knots, eop[:, 12]),
        _create_iers_eop_interpolation(knots, eop[:, 14]),
        _create_iers_eop_interpolation(knots, eop[:, 21]),
        _create_iers_eop_interpolation(knots, eop[:, 23]),
    )
end

# Create the interpolation object for the `knots` and `field` from IERS.
function _create_iers_eop_interpolation(
    knots::AbstractVector,
    field::AbstractVector
)
    # Obtain the last available index of the field.
    last_id = findlast(!isempty, field)
    last_id === nothing && (last_id = length(field))

    # Convert the field to a `Vector{Float64}`.
    field_float::Vector{Float64} = Vector{Float64}(field[1:last_id])

    # Create the interpolation object.
    interp = extrapolate(interpolate(
        (knots[1:last_id],),
        field_float,
        Gridded(Linear())),
        Flat()
    )

    return interp
end

function _get_iers_eop_timespan(itp::AbstractInterpolation)
    str = string(jd_to_date(DateTime, first(first(itp.itp.knots)))) * " -- " *
          string(jd_to_date(DateTime, last(first(itp.itp.knots))))
    return str
end

################################################################################
#                                   Helpers
################################################################################

"""
    deps_dpsi(eop_iau2000a::EOPData_IAU2000A, JD::Number)

Returns the celestial pole offsets in obliquity (`δϵ_2000`) and longitude
(`δΨ_2000`) [arcsec]. This function obtains those values by converting the
celestial pole offsets with respect to the GCRS (`dX` and `dY`). These values
are necessary in the equinox-based IAU-2006 theory.

The algorithm was obtained from **[1]**(eq. 5.25) and
**[2]**(`DPSIDEPS2000_DXDY2000`).

# References

- **[1]**: IERS (2010). Transformation between the International Terrestrial
    Reference System and the Geocentric Celestial Reference System. IERS
    Technical Note No. 36, Chapter 5.
- **[2]**: ftp://hpiers.obspm.fr/eop-pc/models/uai2000.package

"""
function deps_dpsi(eop_iau2000a::EOPData_IAU2000A, JD::Number)
    # Constants.
    d2r = π / 180
    a2d = 1 / 3600
    a2r = a2d * d2r

    # Obtain the parameters `dX` and `dY` [arcseg].
    dX = eop_iau2000a.dX(JD)
    dY = eop_iau2000a.dY(JD)

    # Compute the Julian century.
    T_TT = (JD - JD_J2000) / 36525

    # Luni-solar precession [rad].
    Ψ_a = @evalpoly(T_TT, 0, +5038.47875, -1.07259, -0.001147) * a2r

    # Planetary precession [rad].
    χ_a = @evalpoly(T_TT, 0, +10.5526, -2.38064, -0.001125) * a2r

    sϵ₀, cϵ₀ = sincos(84381.406 * a2r)

    aux = Ψ_a * cϵ₀ - χ_a
    den = aux^2 * sϵ₀ - sϵ₀
    δϵ  = (aux * sϵ₀ * dX - sϵ₀ * dY) / den
    δΨ  = (dX - aux * dY) / den

    return δϵ, δΨ
end


text:
__precompile__()
module MyStats

export hist, hist_indices, histND, histND_indices, min_max, min_max_ind, condmean, Bins, dbkhist_indices, dbkBins, threaded_sum, tmean, tstd, histstd

function hist(field::AbstractArray, min::Real, max::Real, nbins::Integer=100, pdf::Bool=true)
    dx = max - min
    indices = Vector{Float64}(nbins)

    @inbounds for i in linearindices(field)
        n = trunc(Int, ((field[i] - min)/dx - eps())*nbins) + 1
        indices[n] += 1.
    end

    x = Bins(min,max,nbins)
    if pdf
        area = sum(indices)*step(x)
        indices ./= area
    end

    return [x indices]
end

hist(field::AbstractArray, nbins::Integer=100, pdf::Bool=true) = hist(field, min_max(field)..., nbins, pdf)

hist(field::AbstractArray, indices::AbstractVector{<:Integer}, nbins::Integer=100, pdf::Bool=true) = hist(view(field,indices), nbins, pdf)

# return a Vector of Vectors of indices corresponding to the bin.
function hist_indices(field::AbstractArray,min::Real,max::Real,nbins::Integer=30)
    dx = max - min
    indices = Vector{Vector{Int}}(nbins)
    hint = length(field)÷nbins
    @inbounds for i in linearindices(indices)
        indices[i] = Vector{Int}()
        sizehint!(indices[i],hint)
    end

    @inbounds for i in linearindices(field)
        n = trunc(Int, ((field[i] - min)/dx - eps())*nbins) + 1
        (1 <= n <= nbins) && push!(indices[n],i)
    end

    return indices
end

hist_indices(field::AbstractArray,nbins::Integer=30) = hist_indices(field,min_max(field)...,nbins)

function histND_indices(field::NTuple{N,AbstractArray},min::NTuple{N,Real},max::NTuple{N,Real},nbins::NTuple{N,Integer}) where {N}
    dx = max .- min
    indices = Array{Vector{Int},N}(nbins)

    @inbounds for i in linearindices(indices)
        indices[i] = Vector{Int}()
    end

    @inbounds for i in linearindices(field[1])
        n = @. trunc(Int, ((getindex(field,i) - min)/dx - eps())*nbins) + 1
        all(1 .<= n .<= nbins) && push!(indices[n...],i)
    end

    return indices
end

histND_indices(field::NTuple{N,AbstractArray},nbins::NTuple{N,Integer}) where {N} = histND_indices(field,min_max(field)...,nbins)

histND_indices(field::NTuple{3,AbstractArray}) = histND_indices(field,(5,5,5))

histND_indices(field::NTuple{2,AbstractArray}) = histND_indices(field,(12,12))

###################################

function histND(field::NTuple{N,AbstractArray},min::NTuple{N,Real},max::NTuple{N,Real},nbins::NTuple{N,Integer},pdf::Bool=true) where {N}
    dx = max .- min
    indices = zeros(Float64,nbins)

    @inbounds for i in linearindices(field[1])
        n = @. trunc(Int, ((getindex(field,i) - min)/dx - eps())*nbins) + 1
        indices[n...] += 1.
    end
    x = Bins.(min,max,nbins)
    barea = prod((max .- min) ./ nbins)
    if pdf
        area = sum(indices)*barea
        indices ./= area
    end

    return x, indices
end

histND(field::NTuple{N,AbstractArray},nbins::NTuple{N,Integer}) where {N} = histND(field,min_max(field)...,nbins)

histND(field::NTuple{3,AbstractArray}) = histND(field,(5,5,5))

histND(field::NTuple{2,AbstractArray}) = histND(field,(12,12))

function min_max(f::NTuple{N,AbstractArray}) where N
    aux = min_max.(f)
    return (getindex.(aux,1), getindex.(aux,2))
end

function min_max(f::AbstractArray)
    nt::Int = Threads.nthreads()
    const pmin = zeros(eltype(f),nt)
    const pmax = zeros(eltype(f),nt)

    Threads.@threads for i in 1:nt
        get_min_max(f,pmin,pmax,i)
    end

    return minimum(pmin), maximum(pmax)
end

function get_min_max(f,pmi,pma,i)
    nt::Int = Threads.nthreads()

    r = divide_range(length(f), nt)[i]

    minv = f[r[1]]
    maxv = f[r[1]]

    @inbounds for j in r
        if f[j] > maxv
            maxv = f[j]
        elseif f[j] < minv
            minv = f[j]
        end
    end

    pmi[i] = minv
    pma[i] = maxv
    return nothing
end

function divide_range(l::Integer,np::Integer)
    p = l÷np

    m = p
    a = Vector{UnitRange{Int}}(np)
    a[1] = 1:m
    @inbounds for i=2:(np-1)
        a[i] = (m+1):(m+p)
        m = m+p
    end
    a[np] = (m+1):l
    return a
end

function min_max_ind(f::AbstractArray)
    nt::Int = Threads.nthreads()
    const imin = zeros(Int,nt)
    const imax = zeros(Int,nt)

    Threads.@threads for i in 1:nt
        get_min_max_ind(f,imin,imax,i)
    end

    return compare_ind(f,imin,Base.:<), compare_ind(f,imax,Base.:>)
end

function get_min_max_ind(f,pmi,pma,i)
    nt::Int = Threads.nthreads()

    r = divide_range(length(f), nt)[i]

    minv = r[1]
    maxv = r[1]

    @inbounds for j in r
        if f[j] > f[maxv]
            maxv = j
        elseif f[j] < f[minv]
            minv = j
        end
    end

    pmi[i] = minv
    pma[i] = maxv
    return nothing
end

function compare_ind(v::AbstractArray,ind,f::Function)
    tv = ind[1]

    @inbounds for j in ind
        if f(v[j],v[tv])
            tv = j
        end
    end

    return tv
end

function condmean(field::AbstractArray,condindices)
    const result = zeros(length(condindices))
    const err = zeros(length(condindices))
    Threads.@threads for i in linearindices(condindices)
        @inbounds begin
            ind = condindices[i]
            l = length(ind)
            for j in ind
                result[i] += field[j]
            end
            m = result[i]/l
            result[i] = m
            for j in ind
                err[i] += (field[j] - m)^2
            end
            err[i] = sqrt(err[i]/(l*(l-1)))
        end
    end
    return result, err
end

struct Bins <: AbstractVector{Float64}
    minv::Float64
    maxv::Float64
    n::Int
    dx::Float64
    Bins(minv,maxv,n) = new(minv,maxv,n,(maxv-minv)/n)
end

Base.length(a::Bins) = a.n
Base.size(a::Bins) = (a.n,)
Base.step(a::Bins) = a.dx
Base.IndexStyle(::Type{Bins}) = IndexLinear()

@inline function Base.getindex(a::Bins,i::Integer)
    @boundscheck((1 <= i <= a.n) || throw(BoundsError(a,i)))
    dx = a.dx
    return (a.minv + 0.5*dx) + (i-1)*dx
end

function dbkhist_indices(field::AbstractArray, nbins::Integer=30)

    indicesp = sortperm(vec(field))

    indices = Vector{Vector{Int}}(nbins)

    ranges = divide_range(length(field),nbins)
    @inbounds for i in 1:nbins
        indices[i] = indicesp[ranges[i]]
    end

    return indices

end

struct dbkBins <: AbstractVector{Float64}
    edges::Vector{Float64}
    n::Int
    dbkBins(edges) = new(edges,length(endges)-1)
end

Base.length(a::dbkBins) = a.n
Base.size(a::dbkBins) = (a.n,)
Base.IndexStyle(::Type{dbkBins}) = IndexLinear()

@inline function Base.getindex(a::dbkBins,i::Integer)
    @boundscheck((1 <= i <= a.n) || throw(BoundsError(a,i)))
    ed = a.edges
    return 0.5*(ed[i] + ed[i+1])
end

function threaded_sum(v,func::Function=Base.identity)
    nt::Int = Threads.nthreads()
    const result = zeros(eltype(v),nt)
    ranges = divide_range(length(v),nt)
    Threads.@threads for i in 1:nt
        get_sum!(result,v,func,ranges[i],i)
    end

    return sum(result)
end

function get_sum!(result::AbstractVector,v::AbstractArray,f::Function,ran,j::Integer)
    @inbounds begin
        r = zero(eltype(v))
        @simd for i in ran
            r += f(v[i])
        end
        result[j] = r
    end
    return nothing
end

tmean(v::AbstractArray,func::Function=Base.identity) = threaded_sum(v,func)/length(v)

function tstd(v::AbstractArray,func::Function=Base.identity,m::Number=tmean(v,func))
    f = x->((func(x)-m)^2)
    return sqrt(threaded_sum(v,f)/(length(v)-1))
end

tstd(v::AbstractArray,m::Number) = tstd(v,Base.identity,m)

function histstd(field::AbstractArray, min::Real, max::Real, nbins::Integer=100, pdf::Bool=true)
    dx = max - min
    indices = Vector{Float64}(nbins)
    @inbounds for i in linearindices(field)
        n = trunc(Int, ((field[i] - min)/dx - eps())*nbins) + 1
        (1 <= n <= nbins) && (indices[n] += 1.)
    end

    x = Bins(min,max,nbins)
    if pdf
        area = sum(indices)*step(x)
        indices ./= area
    end

    return [x indices]
end

histstd(field::AbstractArray,m::Real,stand::Real,ntimes::Real,nbins::Integer=100,pdf::Bool=true) = histstd(field,m-ntimes*stand,m+ntimes*stand,nbins,pdf)

function histstd(field::AbstractArray,ntimes::Real,nbins::Integer=100,pdf::Bool=true)
    m = tmean(field)
    stand = tstd(field,m)
    histstd(field,m-ntimes*stand,m+ntimes*stand,nbins,pdf)
end

end # module

ext: jl


text:
using OffsetArrays

const program = OffsetVector(parse.(Int, split(readchomp("input"), ',')), -1)

retrieve(program::AbstractVector, x::Integer, mode::Integer) =
    mode == 0 ? program[x] :
    mode == 1 ? x :
    error("unrecognized mode $mode")

function run!(program::AbstractVector, input::AbstractVector)
    output = Int[]
    pos = 0
    while pos < length(program) - 1 && program[pos] != 99
        code = digits(program[pos], pad=5)
        opcode, _, mode1, mode2, mode3 = code
        @assert mode3 == 0
        if opcode == 1
            a, b, c = program[pos+1:pos+3]
            program[c] = retrieve(program, a, mode1) + retrieve(program, b, mode2)
            pos += 4
        elseif opcode == 2
            a, b, c = program[pos+1:pos+3]
            program[c] = retrieve(program, a, mode1) * retrieve(program, b, mode2)
            pos += 4
        elseif opcode == 3
            program[program[pos+1]] = pop!(input)
            pos += 2
        elseif opcode == 4
            push!(output, retrieve(program, program[pos+1], mode1))
            pos += 2
        elseif opcode == 5
            if retrieve(program, program[pos+1], mode1) > 0
                pos = retrieve(program, program[pos+2], mode2)
            else
                pos += 3
            end
        elseif opcode == 6
            if retrieve(program, program[pos+1], mode1) == 0
                pos = retrieve(program, program[pos+2], mode2)
            else
                pos += 3
            end
        elseif opcode == 7
            a, b, c = program[pos+1:pos+3]
            program[c] = Int(retrieve(program, a, mode1) < retrieve(program, b, mode2))
            pos += 4
        elseif opcode == 8
            a, b, c = program[pos+1:pos+3]
            program[c] = Int(retrieve(program, a, mode1) == retrieve(program, b, mode2))
            pos += 4
        else
            error("unexpected opcode $opcode as position $pos")
        end
    end
    return output
end

println("Part 1 answer: ", last(run!(copy(program), [1])))
println("Part 2 answer: ", last(run!(copy(program), [5])))

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
@testset "fill" begin
    s = """
        @def x = 5
        {{fill x}}
        """ |> fd2html_td
    @test isapproxstr(s, "5")
    @test_throws F.HTMLFunctionError ("{{fill x y z}}" |> fd2html)
end

@testset "href" begin
    s = "{{href EQR foo}}" |> fd2html_td
    @test isapproxstr(s, "<b>??</b>")
end

@testset "redirect" begin
    # only one arg
    @test_throws F.HTMLFunctionError ("{{redirect foo bar}}" |> fd2html)
    # must end with html
    @test_throws F.HTMLFunctionError ("{{redirect foo}}" |> fd2html)
end


text:
module TestVine

using Copulas
using Base.Test

#######################
## Tree functions ##
#######################

## getindex
##---------

@test tP[1] == [1]
@test tP[3] == [2, 1]
@test_throws BoundsError tP[8]

@test tP[3, 2] == 1
@test tP[1, 1] == 1
@test_throws BoundsError tP[1, 3]

@test tP[[2, 2], [1, 2]] == [1, 2]
@test tP[[4, 4], [2, 3]] == [3, 4]
@test tP[[4, 1], [3, 1]] == [1, 4]
@test tP[[1, 2], [1, 1]] == [1]
@test_throws Exception tP[[1, 2], [3]]

## equality
##---------

paths2 = Array{Int, 1}[[1], [1, 2, 3], [2, 3, 4], [2, 1, 3]]
tP2 = Copulas.Tree(5, paths2)

@test tP == tP
@test !(tP == tP2)

## conditional set
##----------------

tP = Copulas.Tree(6, [1 2 3], [1 2 4], [5 3])
@test Copulas.condSetChk([2, 1], tP)
@test !(Copulas.condSetChk([2, 4], tP))

end


text:
abstract type Model{T} end

floating_base_dim(::Model) = 3
nominal_configuration(model::Model) = zeros(model.nq)


text:
@testset "HSA Status Error" begin
    errorcode = AMDGPU.HSAError(HSA.STATUS_SUCCESS)
    @test AMDGPU.description(errorcode) == "HSA_STATUS_SUCCESS: The function has been executed successfully."
end

@testset "HSA Async Queue Error" begin
    kernel() = (AMDGPU.trap(); nothing)
    @test_throws AMDGPU.QueueError wait(@roc kernel())
end


text:
"""
    insert_parents!(node::Node, template, <keyword arguments>)
    insert_generations!(node::Node, template, <keyword arguments>)
    insert_children!(node::Node, template, <keyword arguments>)
    insert_siblings!(node::Node, template, <keyword arguments>)

Insert new nodes in the mtg following filters rules. It is important to note the function
always return the root node, whether it is the old one or a new inserted one, so the user is
encouraged to assign the results to an object.

Insert nodes programmatically in an MTG as:
- new parents of the filtered nodes: `insert_parents!`
- new children of the filtered nodes: `insert_children!`
- new siblings of the filtered node: `insert_siblings!`
- new children of the filtered nodes, but the previous children of the filtered node become
the children of the inserted node: `insert_generations!`

# Arguments

## Mandatory arguments

- `node::Node`: The node to start at.
- `template`:
    - A template [`NodeMTG`](@ref) or [`MutableNodeMTG`](@ref) used for the inserted node,
    - A NamedTuple with values for link, symbol, index, and scale
    - Or a function taking the node as input and returning said template
- `attr`: Attributes for the node. Similarly to `template`, can be:
    - An attribute of the same type as of node attributes (*e.g.* a Dict or a NamedTuple)
    - A function to compute new attributes (should also return same type for the attributes)

## Keyword Arguments (filters)

- `scale = nothing`: The scale at which to insert. Usually a Tuple-alike of integers.
- `symbol = nothing`: The symbol at which to insert. Usually a Tuple-alike of Strings.
- `link = nothing`: The link with at which to insert. Usually a Tuple-alike of Char.
- `all::Bool = true`: Continue after the first insertion (`true`), or stop.
- `filter_fun = nothing`: Any function taking a node as input, e.g. [`isleaf`](@ref) to decide
on which node the insertion will be based on.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","A1B1.mtg")
mtg = read_mtg(file)

# Insert new Shoot nodes before all scale 2 nodes:
mtg = insert_parents!(mtg, MultiScaleTreeGraph.MutableNodeMTG("/", "Shoot", 0, 1), scale = 2)

mtg
```
"""
insert_parents!, insert_generations!, insert_children!, insert_siblings!

function insert_parents!(
    node,
    template,
    attr = typeof(node.attributes)();
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true,
    filter_fun = nothing
)

    insert_nodes!(
        node, template, insert_parent!, attr;
        scale = scale, symbol = symbol, link = link, all = all, filter_fun = filter_fun
    )
end

function insert_generations!(
    node,
    template,
    attr = typeof(node.attributes)();
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true,
    filter_fun = nothing
)
    insert_nodes!(
        node, template, insert_generation!, attr;
        scale = scale, symbol = symbol, link = link, all = all, filter_fun = filter_fun
    )
end

function insert_children!(
    node,
    template,
    attr = typeof(node.attributes)();
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true,
    filter_fun = nothing
)
    insert_nodes!(
        node, template, insert_child!, attr;
        scale = scale, symbol = symbol, link = link, all = all, filter_fun = filter_fun
    )
end

function insert_siblings!(
    node,
    template,
    attr = typeof(node.attributes)();
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true,
    filter_fun = nothing
)
    insert_nodes!(
        node, template, insert_sibling!, attr;
        scale = scale, symbol = symbol, link = link, all = all, filter_fun = filter_fun
    )
end

"""
Actual workhorse of insert_parents!, insert_generations!, insert_children!, insert_siblings!
"""
function insert_nodes!(
    node,
    template,
    fn,
    attr = typeof(node.attributes)();
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true, # like continue in the R package, but actually the opposite
    filter_fun = nothing
)

    max_node_id = [max_id(node)]
    # # Check the filters once, and then compute the descendants recursively using `descendants_`
    check_filters(node, scale = scale, symbol = symbol, link = link)
    keep = is_filtered(node, scale, symbol, link, filter_fun) # true is inserted²²

    if isempty(methods(attr)) && !isa(attr, Function)
        # Attr is not given as a function, making it a function
        attr_fun = x -> attr
    else
        attr_fun = attr
    end

    # Only go to the children if we keep the current node and don't want all values:
    if all || !keep
        insert_nodes!_(node, template, fn, attr_fun, max_node_id, scale, symbol, link, all, filter_fun)
    end

    # We apply the function *after* visiting the children to be sure we don't add nodes indefinitely:
    if keep
        node = fn(node, template, attr_fun, max_node_id)
    end

    # Always return the root, whether it is the same one or a new one
    return get_root(node)
end

function insert_nodes!_(node, template, fn, attr_fun, max_node_id, scale, symbol, link, all, filter_fun)

    # Is there any filter happening for the current node? (true is inserted):
    keep = is_filtered(node, scale, symbol, link, filter_fun)

    # Only go to the children if we keep the current node and don't want all values:
    if !isleaf(node) && (all || !keep)
        # First we apply the algorithm recursively on the children:
        for chnode in ordered_children(node)
            insert_nodes!_(chnode, template, fn, attr_fun, max_node_id, scale, symbol, link, all, filter_fun)
        end
    end

    # We apply the function *after* visiting the children to be sure we don't add nodes indefinitely:
    if keep
        node = fn(node, template, attr_fun, max_node_id)
    end

    return node
end

"""
    new_node_MTG(node, template<:Union{NodeMTG,MutableNodeMTG,NamedTuple,MutableNamedTuple})
    new_node_MTG(node, fn)

Returns a new NodeMTG matching the one used in node (either `NodeMTG` or `MutableNodeMTG`)
based on a template, or on a function that takes a node as input and return said template.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

# using a NodeMTG as a template:
MultiScaleTreeGraph.new_node_MTG(mtg, NodeMTG("/", "Leaf", 1, 2))
# Note that it returns a MutableNodeMTG because `mtg` is using this type instead of a `NodeMTG`

# using a NamedTuple as a template:
MultiScaleTreeGraph.new_node_MTG(mtg, (link = "/", symbol = "Leaf", index = 1, scale = 2))

# using a function that returns a template based on the first child of the node:
MultiScaleTreeGraph.new_node_MTG(
    mtg,
    x -> (
            link = x[1].MTG.link,
            symbol = x[1].MTG.symbol,
            index = x[1].MTG.index,
            scale = x[1].MTG.scale)
        )
```
"""
function new_node_MTG(node, fn)
    template = fn(node)
    new_node_MTG(node, template)
end

function new_node_MTG(node, template::T) where {T<:Union{NodeMTG,MutableNodeMTG,NamedTuple,MutableNamedTuple}}
    typeof(node.MTG)(template.link, template.symbol, template.index, template.scale)
end

"""
    insert_parent!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])
    insert_generation!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])
    insert_child!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])
    insert_sibling!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])

Insert a node in an MTG as:

- a new parent of node: `insert_parent!`
- a new child of node: `insert_child!`
- a new sibling of node: `insert_sibling!`
- a new child of node, but the children of node become the children of the inserted node:
`insert_generation!`

# Arguments

- `node::Node`: The node from which to insert a node (as its parent, child or sibling).
- `template`:
    - A template [`NodeMTG`](@ref) or [`MutableNodeMTG`](@ref) used for the inserted node,
    - A NamedTuple with values for link, symbol, index, and scale
    - Or a function taking the node as input and returning said template
- `attr_fun`: A function to compute new attributes based on the filtered node. Must return
attribute values of the same type as the one used in other nodes from the MTG (*e.g.* Dict or
NamedTuple). If you just need to pass attributes values to a node use `x -> your_values`.
- `max_id::Vector{Int64}`: The maximum id of the nodes in the MTG as a vector of length one.
Used to compute the name of the inserted node. It is incremented in the function, and use by
default the value from [`max_id`](@ref).

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","A1B1.mtg")
mtg = read_mtg(file)

template = MultiScaleTreeGraph.MutableNodeMTG("/", "Shoot", 0, 1)
mtg = insert_parent!(mtg[1][1], template)
mtg

# The template can be a function that returns the template. For example a dummy example would
# be a function that uses the NodeMTG of the first child of the node:

insert_parent!(
    mtg[1][1],
    node -> (
        link = node[1].MTG.link,
        symbol = node[1].MTG.symbol,
        index = node[1].MTG.index,
        scale = node[1].MTG.scale)
    )
)
```
"""
insert_parent!, insert_generation!, insert_child!, insert_sibling!

function insert_parent!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])

    max_id[1] += 1

    if isroot(node)

        new_node = Node(
            join(["node_", max_id[1]]),
            max_id[1],
            nothing,
            Dict{Int,Node}(node.id => node),
            nothing,
            new_node_MTG(node, template),
            attr_fun(node)
        )

        # Add to the new root the mandatory root attributes:
        root_attrs = Dict(
            :symbols => node[:symbols],
            :scales => node[:scales],
            :description => node[:description]
        )

        append!(new_node, root_attrs)

        # Add the new root node as the parent of the previous one:
        node.parent = new_node
    else
        new_node = Node(
            join(["node_", max_id[1]]),
            max_id[1],
            node.parent,
            Dict{Int,Node}(node.id => node),
            nothing,
            new_node_MTG(node, template),
            attr_fun(node)
        )

        # Add the new node to the parent:
        pop!(node.parent.children, node.id)
        push!(node.parent.children, new_node.id => new_node)

        # Add the new node as the parent of the previous one:
        node.parent = new_node
    end

    return node
end

function insert_child!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])

    max_id[1] += 1

    addchild!(node, max_id[1], new_node_MTG(node, template), attr_fun(node))

    return node
end

function insert_sibling!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])

    max_id[1] += 1

    new_node = Node(
        join(["node_", max_id[1]]),
        max_id[1],
        parent(node),
        nothing,
        parent(node).children,
        new_node_MTG(node, template),
        attr_fun(node)
    )

    # Add the new node to the children of the parent node:
    push!(parent(node).children, new_node.id => new_node)

    return node
end

function insert_generation!(node, template, attr_fun = node -> typeof(node.attributes)(), max_id = [max_id(node)])

    max_id[1] += 1

    new_node = Node(
        join(["node_", max_id[1]]),
        max_id[1],
        node,
        node.children,
        nothing,
        new_node_MTG(node, template),
        attr_fun(node)
    )

    # Add the new node as the only child of the node:
    node.children = Dict{Int,Node}(new_node.id => new_node)

    return node
end

@deprecate insert_node!(node, template, max_id) insert_parent!(node, template, max_id)


text:
__precompile__()

module Materials

using LinearAlgebra
using ..adiff

# Material types
struct Hooke
  E     ::Float64
  ν     ::Float64
  small ::Bool
  Hooke(E,ν;small=false) = new(Float64(E),Float64(ν), small)
end
struct MooneyRivlin
  C1  ::Float64
  C2  ::Float64
  K   ::Float64
  MooneyRivlin(C1, C2)    = new(Float64(C1), Float64(C2), NaN)
  MooneyRivlin(C1, C2, K) = new(Float64(C1), Float64(C2), Float64(K))
end
struct NeoHooke
  C1   ::Float64
  K    ::Float64
  NeoHooke(C1)    = new(Float64(C1), NaN)
  NeoHooke(C1, K) = new(Float64(C1), Float64(K))
end
Material = Union{Hooke,MooneyRivlin,NeoHooke}
HyperEla = Union{MooneyRivlin,NeoHooke}
# default convergence tolerance for 2D stress
dTol     = 1e-7
maxiter  = 30
function setmaxiter(x)
  global maxiter = Int64(x)
end
function setdTol(x)
  global dTol = Float64(x)
end
# elastic energy evaluation functions for materials
function getϕ(F::Array{N,2}, mat::M) where {N<:Number, M<:HyperEla}
  C = transpose(F)F
  if length(C) == 9
    (I1,I2,I3) = getInvariants(C)
  else
    L3 = isnan(mat.K) ? (F[1]F[4]-F[2]F[3])^-2 : 1
    (I1,I2,I3) = getInvariants(C, L3)
  end
  Materials.getϕ(I1,I2,I3,mat)
end
function getϕ(I1, I2, I3, mat::MooneyRivlin)

  C1, C2, K = mat.C1, mat.C2, mat.K
  if isnan(K)
    ϕ  = C1*(I1-3) + C2*(I2-3)
  else
    J  = sqrt(I3)
    I1 = I1*J^(-2/3)
    I2 = I2*J^(-4/3)
    # ϕ  = C1*(I1-3) + C2*(I2-3) + K*log(J)^2
    ϕ  = C1*(I1-3) + C2*(I2-3) + K*(J-1)^2
  end
  return ϕ
end
function getϕ(I1, I2, I3, mat::NeoHooke)

  C1, K  = mat.C1, mat.K
  if isnan(K)
    ϕ  = C1*(I1-3)
  else
    J  = sqrt(I3)
    I1 = I1*J^(-2/3)
    # ϕ  = C1*(I1-3) + K*log(J)^2
    ϕ  = C1*(I1-3) + K*(J-1)^2
  end
  return ϕ
end
function getϕ(C11::N where N<:Number, mat::Hooke)
  ϕ = 0.5mat.E*(C11-1)^2
end
function getϕ(F::Array{N,2} where N<:Number, mat::Hooke)

  if mat.small
    E = 0.5(F+transpose(F))-I   # the symmetric part of G
  else
    E = 0.5(transpose(F)F-I)    # the Green-Lagrange strain tensor
  end

  ν, Es = mat.ν, mat.E
  λ = Es*ν/(1+ν)/(1-2ν)
  μ = Es/2/(1+ν)

  if size(F) == (3,3)
    ϕ =  (μ+λ/2) * (E[1]^2   + E[5]^2   + E[9]^2)
    ϕ += λ       * (E[1]E[5] + E[5]E[9] + E[9]E[1])
    ϕ += 2μ      * (E[2]^2   + E[3]^2   + E[6]^2)
  else
    # 2D is plain strain
    ϕ = (μ+λ/2)*(E[1]^2+E[4]^2) + λ*E[1]E[4] + 2μ*E[2]^2
  end

  return ϕ
end
# status retrieving functions
function getP(F::Array{Float64,2}, mat) # 1st PK tensor from F
  ϕ = getϕ(adiff.D2(F), mat)
  reshape(adiff.grad(ϕ), (3,3))
end
function getσ(F::Array{Float64,2}, mat) # Cauchy stress
  J = det(F)
  P = getP(F,mat)         # 1st PK stress
  (P*transpose(F))/J
end
function getinfo(F::Array{Float64,2}, mat; info = :ϕ)
  if info == :ϕ
    getϕ(F, mat)
  elseif info == :detF
    det(F)
  elseif info == :G
    F-I
  elseif info == :P
    getP(F, mat)
  elseif info == :σ
    getσ(F, mat)
  elseif info == :σVM
    σ = getσ(F, mat)
    sqrt(0.5)*sqrt((σ[1]-σ[5])^2+(σ[5]-σ[9])^2+(σ[9]-σ[1])^2+6(σ[4]^2+σ[7]^2+σ[8]^2))
  elseif info == :S
    F\getP(F, mat)
  elseif info == :I
    getInvariants(transpose(F)F)
  elseif info == :LE
    svdF = svd(F)
    svdF.V*diagm(0=>log.(svdF.S))*transpose(svdF.V)
  else
    F
  end
end
# helper functions for 2D and 1D stress
function getϕ(F11::N, mat::M; binc=true) where {N<:Number, M<:HyperEla}

  C11 = F11^2
  L1  = Real(C11)
  L3  = binc ? getL3(L3->getϕ(L1,L3,L3,mat), L1) :
  getL3(L3->getϕ(L1,L3,L3,mat), √L1)

  if isa(C11, adiff.D2)
    ϕ1    = getϕ(adiff.D2([L1, L3, L3])..., mat)
    dL21  = -ϕ1.h[2]/(ϕ1.h[3]+ϕ1.h[5])
    L2    = adiff.D2(L3, dL21*C11.g)
    ϕ     = getϕ(C11,L2,L2,mat)
  else
    ϕ  = getϕ(C11,L3,L3,mat)
  end

  return ϕ
end
function getL3(func, L3)
  iter   = 0
  bdone  = false
  while !bdone
    ϕ = func(adiff.D2(L3))
    r = adiff.grad(ϕ)[1]
    if abs(r) < dTol
      bdone = true
    else
      updt = -r/adiff.hess(ϕ)[1]
      L3   = L3+updt<0 ? L3/2 : L3+updt
    end
    iter +=1
    if iter > maxiter
      println("failed in getL3, with L3: ",L3," r: ", r, " and iter: ",iter)
      error("failed in getL3")
      bdone = true
    end
  end
  return  L3
end
function getInvariants(C, C33)
  I1 = C[1]+C[4]+C33
  I2 = C[1]C[4]+C[4]C33+C[1]C33-C[2]^2
  I3 = C[1]C[4]C33-C33*C[2]^2

  (I1,I2,I3)
end
function getInvariants(C)
  I1 = C[1]+C[5]+C[9]
  I2 = C[1]C[5]+C[5]C[9]+C[1]C[9]-C[2]^2-C[3]^2-C[6]^2
  I3 = C[1]C[5]C[9]+2C[2]C[3]C[6]-C[1]C[6]^2-C[5]C[3]^2-C[9]C[2]^2

  (I1,I2,I3)
end
end


text:
#%%
function getParam(iLine)
    return strip( split(lines[iLine], "=")[2] , ' ')
end

function getParam2(n, iLine, typeNumber)
    v = zeros(n)
    l = split( strip( getParam(iLine), ['{', '}'] ), ',' )
    for i in 1:n
        if typeNumber == "Int"
            v[i] = parse(Int, l[i])
        elseif typeNumber == "Float64"
            v[i] = parse(Float64, l[i])
        end
    end
    return v
end

function getParam3(n, iLine)
    v = zeros(Int, (n, 4))
    l = split( strip( getParam(iLine), ['{', '}'] ), "}, {"  )
    for i in 1:n
            temp = split( l[i], ',' )
            for j in 1:4
                v[i, j] = parse(Int, temp[j])
            end
    end
    return v
end

function load(filename="workdir/pot.mtp")
    lines = readlines(filename)

    @assert lines[1] == "MTP" "Can read only MTP format potentials"

    version = getParam(2)
    @assert version == "1.1.0" "MTP file must have version \"1.1.0\""

    potential_name     = getParam(3)
    scaling            = parse(Float64, getParam(4))
    species_count      = parse(Int, getParam(5))
    potential_tag      = getParam(6)
    radial_basis_type  = getParam(7)
    min_dist           = parse(Float64, getParam(8))
    max_dist           = parse(Float64, getParam(9))
    radial_basis_size  = parse(Int, getParam(10))
    radial_funcs_count = parse(Int, getParam(11))

    iline = 13
    regression_coeffs = zeros( (species_count, species_count, radial_funcs_count, radial_basis_size) )
    for s1 in 1:species_count
        for s2 in 1:species_count
            iline += 1
            for i in 1:radial_funcs_count
                list_t = split( strip( lines[15], ['\t', '{', '}'] ), ',' )
                iline += 1
                for j in 1:radial_basis_size
                    regression_coeffs[s1, s2, i, j] = parse(Float64, list_t[j])
                end
            end
        end
    end
    #

    alpha_moments_count = parse(Int, getParam(iline))

    alpha_index_basic_count = parse(Int, getParam(iline + 1))
    alpha_index_basic       = getParam3(alpha_index_basic_count, iline + 2)

    radial_func_max = maximum([ alpha_index_basic[i, 1] for i in 1:alpha_index_basic_count ])
    @assert radial_func_max == (radial_funcs_count - 1) "Wrong number of radial functions specified"

    alpha_index_times_count = parse(Int, getParam(iline + 3))
    alpha_index_times       = getParam3(alpha_index_times_count, iline + 4)

    alpha_scalar_moments = parse(Int, getParam(iline + 5))
    alpha_moment_mapping = getParam2(alpha_scalar_moments, iline + 6, "Int")

    alpha_count = alpha_scalar_moments + 1

    species_coeffs       = getParam2(species_count, iline + 7, "Float64")
    moment_coeffs        = getParam2(alpha_scalar_moments, iline + 8, "Float64")
    #
    params = Dict(
            "scaling" => scaling,
            "species_count" => species_count,
            "rbasis_type" => radial_basis_type,
            "min_dist" => min_dist,
            "max_dist" => max_dist,
            "rb_size" => radial_basis_size,
            "radial_func_count" => radial_funcs_count,
            "regression_coeffs" => regression_coeffs,
            "alpha_moments_count" => alpha_moments_count,
            "alpha_index_basic_count" => alpha_index_basic_count,
            "alpha_index_basic" => alpha_index_basic,
            "alpha_index_times" => alpha_index_times,
            "alpha_scalar_moments" => alpha_scalar_moments,
            "alpha_count" => alpha_count,
            # "linear_coeffs":linear_coeffs,
            "species_coeffs" => species_coeffs,
            "moment_coeffs" => moment_coeffs,
            "alpha_index_times_count" => alpha_index_times_count,
            "alpha_moment_mapping" => alpha_moment_mapping,
            )
    #
    return params
end

function init_vecs(parameters)
    rb_vals     = zeros(parameters["rb_size"])
    rb_ders     = zeros(parameters["rb_size"])
    moment_vals = zeros(parameters["alpha_moments_count"])
    basis_vals  = zeros(parameters["alpha_count"])
    site_energy_ders_wrt_moments_ = zeros(parameters["alpha_moments_count"])

    mult        = 2.0 / (parameters["max_dist"] - parameters["min_dist"])

    max_alpha_index_basic = maximum(parameters["alpha_index_basic"]) + 1
    inv_dist_powers_ = zeros(max_alpha_index_basic)
    coords_powers_   = zeros( (max_alpha_index_basic, 3) )

    linear_mults = ones(parameters["alpha_scalar_moments"])
    max_linear = 1e-10 * ones(parameters["alpha_scalar_moments"])

    # # calculate moment_vals[i] for i only in the following list:
    # alpha_index_basic_count = parameters["alpha_index_basic_count"]
    # alpha_index_times_count = parameters["alpha_index_times_count"]
    # alpha_index_times = parameters["alpha_index_times"]
    # alpha_moment_mapping = parameters["alpha_moment_mapping"]
    # l1 = [ i for i in range(alpha_index_basic_count) if i in alpha_moment_mapping ]
    # l2 = [ alpha_index_times[i, 0] for i in range(alpha_index_times_count) if alpha_index_times[i, 0] < alpha_index_basic_count ]
    # l3 = [ alpha_index_times[i, 1] for i in range(alpha_index_times_count) if alpha_index_times[i, 1] < alpha_index_basic_count ]
    # l4 = [ alpha_index_times[i, 3] for i in range(alpha_index_times_count) if alpha_index_times[i, 3] < alpha_index_basic_count ]
    # l1 = l1 + list( set(l2) - set(l1) )
    # l1 = l1 + list( set(l3) - set(l1) )
    # i_moment_vals = l1 + list( set(l4) - set(l1) )

    initializedVecs = Dict(
                "rb_vals" => rb_vals,
                "rb_ders" => rb_ders,
                "moment_vals" => moment_vals,
                "basis_vals" => basis_vals,
                "site_energy_ders_wrt_moments_" => site_energy_ders_wrt_moments_,
                "max_alpha_index_basic" => max_alpha_index_basic,
                "inv_dist_powers_" => inv_dist_powers_,
                "coords_powers_" => coords_powers_,
                "linear_mults" => linear_mults,
                "max_linear" => max_linear,
                "mult" => mult
    )
    #
    return initializedVecs
end

function correctIndexesInParameters!(parameters)
    # mutates `parameters`
    # correcting to index starting in 1 in Julia, instead of 0 as in C++ (MTP code)
    for i in 1:parameters["alpha_index_basic_count"]
        parameters["alpha_index_basic"][i, 1] += 1 # = mu
    end
    for i in 1:parameters["alpha_index_times_count"]
        parameters["alpha_index_times"][i, 1] += 1 # val0 = moment_vals[ alpha_index_times[i, 0] ] in finish_moment_vals()
        parameters["alpha_index_times"][i, 2] += 1 # val1 = moment_vals[ alpha_index_times[i, 1] ]
        parameters["alpha_index_times"][i, 4] += 1 # moment_vals[alpha_index_times[i, 3]] += val2 * val0 * val1
    end
    for i in 1:parameters["alpha_scalar_moments"]
        parameters["alpha_moment_mapping"][i] += 1
    end
    parameters["alpha_index_basic"] = alpha_index_basic
    parameters["alpha_index_times"] = alpha_index_times
    parameters["alpha_moment_mapping"] = alpha_moment_mapping
end

function belongs(x, xmin, xmax)
    return (xmin <= x) && (x <= xmax)
end

# rb_Calc!(r, min_dist, max_dist, mult, rb_vals, rb_ders, scaling, rb_size)
function rb_Calc!(r, min_dist, max_dist, mult, rb_vals, rb_ders, scaling, rb_size)
    # mutates rb_vals, rb_ders
    # from src/radial_basis.cpp: void RadialBasis_Chebyshev::RB_Calc(double r)

    # if parameters["rbasis_type"] == "RBChebyshev":
    @assert belongs(r, min_dist, max_dist) "r does not belong to
                [min_dist, max_dist]: " + string(r) + " min_dist=" +
                string(min_dist) + " max_dist=" + string(max_dist)
    #

    ksi = ( (2 * r) - (min_dist + max_dist)) / (max_dist - min_dist)
    R   = r - max_dist
    R2  = R ^ 2

    # scaling is not read in RadialBasis_Chebyshev::RB_Calc(), so it keeps its default value = 1
    rb_vals[1] = R2 # instead of `scaling * (1 * R2)``
	# rb_ders[0] = scaling * (0 * R2 + 2 * R) <<< why 0 * something??? I have to simplify below
    rb_ders[1] = 2.0 * R

    rb_vals[2] = ksi * R2
    rb_ders[2] = (mult * R2) + (2.0 * ksi * R)

    for i in 3:rb_size
        v1 = rb_vals[i - 1]
        v2 = rb_vals[i - 2]
        d1 = rb_ders[i - 1]
        d2 = rb_ders[i - 2]
        rb_vals[i] = (2.0 * ksi * v1) - v2
        rb_ders[i] = 2.0 * ( (mult * v1) + (ksi * d1) ) - d2
    end
    # return rb_vals, rb_ders
end

# _pows(r, inv_dist_powers_, coords_powers_, max_alpha_index_basic, NeighbVect_j)
function _pows!(r, inv_dist_powers_, coords_powers_, max_alpha_index_basic, NeighbVect_j)
    # mutates inv_dist_powers_, coords_powers_
    inv_dist_powers_[1] = 1.0
    coords_powers_[1]   = [1.0, 1.0, 1.0]
    for k in 2:max_alpha_index_basic
        inv_dist_powers_[k] = inv_dist_powers_[k - 1] / r
        for i in 1:3
            coords_powers_[k][i] = coords_powers_[k - 1][i] * NeighbVect_j[i]
        end
    end
    return inv_dist_powers_, coords_powers_
end

# Next: calculating non-elementary b_i
function finish_moment_vals!(
                moment_vals,
                alpha_index_times_count,
                alpha_index_times,
                )
    # mutates moment_vals
    for i in 1:alpha_index_times_count
        val0 = moment_vals[ alpha_index_times[i, 1] ] # float
        val1 = moment_vals[ alpha_index_times[i, 2] ] # float
        val2 = alpha_index_times[i, 3]                # integer
        moment_vals[alpha_index_times[i, 4]] += val2 * val0 * val1
	end
end

# val, der = get_val_der( mu, regression_coeffs, type_central, type_outer, rb_vals, rb_ders )
function get_val_der(    mu,
                    regression_coeffs,
                    type_central,
                    type_outer,
                    rb_vals,
                    rb_ders
                )
    #
    val = regression_coeffs[type_central, type_outer, mu, :] .* rb_vals
    der = regression_coeffs[type_central, type_outer, mu, :] .* rb_ders
    return val, der
end

function update_moment_vals!(
                    i_moment_vals,
                    moment_vals,
                    alpha_index_basic_count,
                    alpha_index_basic,
                    inv_dist_powers_,
                    regression_coeffs,
                    type_central,
                    type_outer,
                    rb_vals,
                    rb_ders,
                    coords_powers_
                    )
    # update_moment_vals
    #
    for i in 1:alpha_index_basic_count
    # for i in i_moment_vals:
        # print(i)
        val = 0.0
        # der = 0.0
        mu = alpha_index_basic[i, 1]

        # for l in range(rb_size):
        #     m = (type_central*C + type_outer)*K*R + mu * R + l
        #     val += regression_coeffs[m] * rb_vals[l]
        #     der += regression_coeffs[m] * rb_ders[l]

        # for l in range(rb_size):
        #     coef = regression_coeffs[type_central, type_outer, mu, l]
        #     val += coef * rb_vals[l]
        #     # der += coef * rb_ders[l]
        # #
        # depend on `r`:
        val, der = get_val_der( mu, regression_coeffs,
                                type_central, type_outer,
                                rb_vals, rb_ders )
        #

        a1 = alpha_index_basic[i, 2]
        a2 = alpha_index_basic[i, 3]
        a3 = alpha_index_basic[i, 4]

        k = a1 + a2 + a3
        inv_powk = inv_dist_powers_[k]
        val *= inv_powk
        # der = (der * inv_powk) - (k * val / r)

        pow0 = coords_powers_[a1, 1]
        pow1 = coords_powers_[a2, 2]
        pow2 = coords_powers_[a3, 3]

        mult0 = pow0 * pow1 * pow2

        #
        moment_vals[i] += val * mult0
        # mult0 *= der / r
        # moment_jacobian_[i, j, 0] += mult0 * NeighbVect_j[0]
        # moment_jacobian_[i, j, 1] += mult0 * NeighbVect_j[1]
        # moment_jacobian_[i, j, 2] += mult0 * NeighbVect_j[2]

        # if a1 != 0:
        #     prod = val * a1 * coords_powers_[a1 - 1, 0] * pow1 * pow2
        #     moment_jacobian_[i, j, 0] += prod
        # #
        # if a2 != 0:
        #     prod = val * a2 * pow0 * coords_powers_[a2 - 1, 1] * pow2
        #     moment_jacobian_[i, j, 1] += prod
        # #
        # if a3 != 0:
        #     prod = val * a3 * pow0 * pow1 * coords_powers_[a3 - 1, 2]
        #     moment_jacobian_[i, j, 2] += prod
        # #
    end
end

function calcSiteEnergyDers(
                        nbh,
                        type_central, # type of central atom at `i`
                        species_coeffs,
                        moment_coeffs,
                        lenNbh,
                        alpha_index_basic,
                        alpha_index_basic_count,
                        max_alpha_index_basic,
                        alpha_index_times,
                        alpha_index_times_count,
                        alpha_scalar_moments,
                        alpha_moment_mapping,
                        regression_coeffs,
                        moment_vals,
                        inv_dist_powers_,
                        coords_powers_,
                        species_count,
                        radial_func_count,
                        rb_size,
                        min_dist,
                        max_dist,
                        mult,
                        rb_vals,
                        rb_ders,
                        scaling,
                        linear_mults,
                        max_linear,
                        i_moment_vals
                        )
    #
    # from dev_src/mtpr.cpp: void MLMTPR::CalcSiteEnergyDers(const Neighborhood& nbh)
    #

    buff_site_energy_ = 0.0

    # initialize
    moment_vals *= 0.0

    # lenNbh = len(nbh)

    # dicTypes = {"C":0, "O": 1}
    # types = [0,0,0,0, 1,1] #just an example

    # C = species_count   		#number of different species in current potential
    # K = radial_func_count		#number of radial functions in current potential
    # R = rb_size                 #number of Chebyshev polynomials constituting one radial function

    # moment_jacobian_ = np.zeros((alpha_index_basic_count, lenNbh, 3))

    @assert type_central < species_count "Too few species count in the MTP potential!"

    for j in 1:lenNbh
        # NeighbVect_j = nbh.vecs[j] #<<<<<<
        # r = nbh.dists[j]  #<<<<<<

        NeighbVect_j = nbh[j]["D"]
        r = nbh[j]["d"]

        rb_vals, rb_ders = rb_Calc(r, min_dist, max_dist, mult, rb_vals, rb_ders, scaling, rb_size)
        #

        rb_vals *= scaling # rb_vals is numpy array, so we can directly multyply by a float if array is float
        rb_ders *= scaling

        inv_dist_powers_, coords_powers_ = _pows(r, inv_dist_powers_, coords_powers_, max_alpha_index_basic, NeighbVect_j)

        # type_outer = nbh.types[j] #<<<<<<
        type_outer = nbh[j]["type"]

        moment_vals = update_moment_vals(
                            i_moment_vals,
                            moment_vals,
                            alpha_index_basic_count,
                            alpha_index_basic,
                            inv_dist_powers_,
                            regression_coeffs,
                            type_central,
                            type_outer,
                            rb_vals,
                            rb_ders,
                            coords_powers_
                            )
        #

        ## Repulsive term
        ## I think it was not implemented in the C++ MTP code)
		## if (p_RadialBasis->GetRBTypeString() == "RBChebyshev_repuls")
        ## this seems arbitrary, I removed it:
        ## if r < min_dist:
		## 	multiplier = 10000;
		## 	buff_site_energy_ += multiplier*(exp(-10*(r-1)) - exp(-10*(min_dist-1)))
		## 	for (int a = 0; a < 3; a++)
		## 		buff_site_energy_ders_[j][a] += -10 * multiplier*(exp(-10 * (r - 1))/ nbh.dists[j])*nbh.vecs[j][a];
		# #
    end

    # # Next: calculating non-elementary b_i
    # for i in range(alpha_index_times_count):
    #     val0 = moment_vals[ alpha_index_times[i, 0] ] # float
    #     val1 = moment_vals[ alpha_index_times[i, 1] ] # float
    #     val2 = alpha_index_times[i, 2]                # integer
    #     moment_vals[alpha_index_times[i, 3]] += val2 * val0 * val1
	# #
    moment_vals = finish_moment_vals(moment_vals, alpha_index_times_count, alpha_index_times)

    # renewing maximum absolute values
    # for i in range(alpha_scalar_moments):
        # max_linear[i] = max(
        #     max_linear[i],
        #     abs(linear_coeffs[species_count + i] * moment_vals[alpha_moment_mapping[i]])
        #     )
        # #
    #

    # convolving with coefficients
    buff_site_energy_ += species_coeffs[type_central]

    for i in 1:alpha_scalar_moments
        # buff_site_energy_ += linear_coeffs[species_count + i] * linear_mults[i] * moment_vals[alpha_moment_mapping[i]]
        # I simplified because linear_mults[i] = 1
        buff_site_energy_ += moment_coeffs[i] * moment_vals[alpha_moment_mapping[i]]
    end

    return buff_site_energy_
end

function CalcEFS(atoms, neighborhoods, type_centrals, params, vecs)
    # from src/basic_mlip.cpp:  void AnyLocalMLIP::CalcEFS(Configuration& cfg)
    energy = 0.0
    for i in 1:atoms.get_global_number_of_atoms()
        nbh = neighborhoods[i]
        lenNbh = len(nbh)
        type_central = type_centrals[i]

        # energy += calcSiteEnergyDers(nbh, type_central)
        #
        energy += calcSiteEnergyDers(   nbh,
                                        type_central, # type of central atom at `i`
                                        params["species_coeffs"],
                                        params["moment_coeffs"],
                                        lenNbh,
                                        params["alpha_index_basic"],
                                        params["alpha_index_basic_count"],
                                        vecs["max_alpha_index_basic"],
                                        params["alpha_index_times"],
                                        params["alpha_index_times_count"],
                                        params["alpha_scalar_moments"],
                                        params["alpha_moment_mapping"],
                                        params["regression_coeffs"],
                                        vecs["moment_vals"],
                                        vecs["inv_dist_powers_"],
                                        vecs["coords_powers_"],
                                        params["species_count"],
                                        params["radial_func_count"],
                                        params["rb_size"],
                                        params["min_dist"],
                                        params["max_dist"],
                                        vecs["mult"],
                                        vecs["rb_vals"],
                                        vecs["rb_ders"],
                                        params["scaling"],
                                        vecs["linear_mults"],
                                        vecs["max_linear"],
                                        vecs["i_moment_vals"],
                                    )
        #
        # print(energy)
    end
    return energy
end

function get_neighborhoods(atoms, listCutOffs, dictionaryTypes)
    neighborList = neighborlist.NeighborList(listCutOffs, self_interaction=False, bothways=True)
    neighborList.update(atoms)
    li, lj, ld, lD = neighborlist.neighbor_list("ijdD", atoms, listCutOffs)

    neighborhoods = Dict( i => [] for i in range(len(atoms)) )
    # neighborhoods = [ [] for i in range(len(atoms)) ]
    for k in 1:len(li)
        i = li[k]
        j = lj[k]
        d = ld[k]
        D = lD[k]
        t = dictionaryTypes[ atoms[j].symbol ]
        neighborhoods[i].append( Dict( "j"=>j, "d"=>d, "D"=>D, "type"=>t ) )
        # neighborhoods[i].append( [ j, d, D, t ] )
    end
    return neighborhoods
end

# type_centrals = get_type_centrals(atoms, dictionaryTypes)
function get_type_centrals(atoms, dictionaryTypes)
    # dictionaryTypes = { "C":0, "O":1 }
    type_centrals = [dictionaryTypes[i] for i in atoms.get_chemical_symbols()]
    return type_centrals
end

#%%

# x = 3

#%%
using ASE
atoms = bulk("Cu", cubic=true) * 2        # generate periodic Cu supercell
# deleteat!(at, 1)                       # vacancy defect
# emt = pyimport("ase.calculators.emt")  # import the EMT model
# calc = ASECalculator(emt.EMT())        # wrap it into a Julia Object
# @show energy(calc, at)

# atoms.get_global_number_of_atoms()
# len(atoms)

#%%
using JuLIP
at = bulk(:Si, cubic=true)

max_stars_repo_name: chinchay/mtpy

max_issues_repo_name: chinchay/mtpy

max_forks_repo_name: chinchay/mtpy


text:
export Symbolic, AbstractVariable, SymbolicVariable, BasicVariable, processExpr, @sexpr
export SymbolParameter, simplify
import Base.show, Base.(==)

#################################################################
#
# Symbolic types
#   Symbolic - top of the tree
#   AbstractVariable - inherit from this for custom symbolic
#                      variable types
#   SymbolicVariable - use this in method argument typing
#
#################################################################

Compat.@compat abstract type Symbolic end
Compat.@compat abstract type AbstractVariable <: Symbolic end
const SymbolicVariable = Union{Symbol, AbstractVariable}

#################################################################
#
# BasicVariable - an example to use during testing
#
#################################################################

type BasicVariable <: AbstractVariable
    sym::Symbol
end
# The following is probably too plain.
show(io::IO, x::BasicVariable) = print(io, x.sym)
(==)(x::BasicVariable, y::BasicVariable) = x.sym == y.sym

#################################################################
#
# @sexpr - return an Expr with variables spliced in
# processExpr - do the Expr splicing
#
#################################################################

function processExpr(x::Expr)
    if x.head == :call
        quoted = Expr(:quote,x.args[1])
        code = :(Expr(:call,$quoted))
        for y in x.args[2:end]
            push!(code.args,processExpr(y))
        end
        return code
    else
        return x
    end
end

processExpr(x::Any) = x

macro sexpr(x)
    esc(processExpr(x))
end

#################################################################
#
# SymbolParameter
#   used to be able to dispatch on the symbol representing a
#   function
#
#################################################################

type SymbolParameter{T}
end
SymbolParameter(s::Symbol) = SymbolParameter{s}()

#################################################################
#
# simplify()
#
#################################################################

isnumber(a) = isa(a, Number)

# Numbers and symbols can't be simplified further
simplify(x) = x
simplify(n::Number) = n
simplify(s::SymbolicVariable) = s

# The default is just to simplify arguments.
simplify{T}(x::SymbolParameter{T}, args) = Expr(:call, T, map(simplify, args)...)

function simplify(ex::Expr)
    if ex.head != :call
        return ex
    end
    if all(isnumber, ex.args[2:end]) && length(ex.args) > 1
        return eval(current_module(), ex)
    end
    new_ex = simplify(SymbolParameter(ex.args[1]), ex.args[2:end])
    while new_ex != ex
        new_ex, ex = simplify(new_ex), new_ex
    end
    return new_ex
end

function sum_numeric_args(args)
    sum = 0
    sym_args = Any[]
    for arg in args
        if isnumber(arg)
            sum += arg
        else
            push!(sym_args, arg)
        end
    end
    (sum, sym_args)
end

function mul_numeric_args(args)
    prod = 1
    sym_args = Any[]
    for arg in args
        if isnumber(arg)
            prod *= arg
        else
            push!(sym_args, arg)
        end
    end
    (prod, sym_args)
end

# Handle `args` of all lengths
function simplify(::SymbolParameter{:+}, args)
    # Remove any 0's in a sum
    args = map(simplify, filter(x -> x != 0, args))
    if length(args) == 0
        return 0
    # Special Case: simplify(:(+x)) == x
    elseif length(args) == 1
        return args[1]
    else
        (sum, sym_args) = sum_numeric_args(args)
        args = sum==0 ? sym_args : [sum; sym_args]
        return Expr(:call, :+, args...)
    end
end

isminus(ex::Expr) = ex.head == :call && ex.args[1] == :- && length(ex.args) == 2
isminus(ex) = false

# Assume length(args) == 3
function simplify(::SymbolParameter{:-}, args)
    # Remove any 0's in a subtraction
    args = map(simplify, filter(x -> x != 0, args))
    if length(args) == 0
        return 0
    # Special Case: simplify(:(x - x)) == 0
    elseif length(args) == 2 && args[1] == args[2]
        return 0
    # Special Case: simplify(:(x - (-y))) == x + y
    elseif length(args) == 2 && isminus(args[2])
        return Expr(:call, :+, args[1], args[2].args[2])
    else
        return Expr(:call, :-, args...)
    end
end

# Handle `args` of all lengths
function simplify(::SymbolParameter{:*}, args)
    # Remove any 1's in a product
    args = map(simplify, filter(x -> x != 1, args))
    if length(args) == 0
        return 1
    # Special Case: simplify(:(*(x))) == x
    elseif length(args) == 1
        return args[1]
    # Special Case: simplify(:(x * y * z * 0)) == 0
    elseif any(args .== 0)
        return 0
    # Special Case: simplify(:(*(-1,x))) == -x
    elseif length(args) == 2 && args[1] == -1
        return Expr(:call, :-, args[2])
    else
        (prod, sym_args) = mul_numeric_args(args)
        args = prod==1 ? sym_args : [prod; sym_args]
        return Expr(:call, :*, args...)
    end
end

# Assume length(args) == 3
function simplify(::SymbolParameter{:/}, args)
    args = map(simplify, args)
    # Special Case: simplify(:(x / 1)) == x
    if args[2] == 1
        return args[1]
    # Special Case: simplify(:(0 / x)) == 0
    elseif args[1] == 0
        return 0
    else
        return Expr(:call, :/, args...)
    end
end

# Assume length(args) == 3
function simplify(::SymbolParameter{:^}, args)
    args = map(simplify, args)
    # Special Case: simplify(:(x ^ 0)) == 1
    if args[2] == 0
        return 1
    # Special Case: simplify(:(x ^ 1)) == x
    elseif args[2] == 1
        return args[1]
    # Special Case: simplify(:(0 ^ x)) == 0
    elseif args[1] == 0
        return 0
    # Special Case: simplify(:(1 ^ x)) == 1
    elseif args[1] == 1
        return 1
    else
        return Expr(:call, :^, args...)
    end
end


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0,1,20))
end

# GLAbstraction.const_lift is defined as
# const_lift(F, args...) = Reactive.map(F, map(Signal, args)...)
# and allows to also use constant arguments without manually wrapping
# them into signals
π_signal = const_lift(*, timesignal, pi)

# Reactive.map takes a signal and transforms it by applying F into a new signal
# it's similar to Base.map, just over the time dimension (like registering a
# callback in a more traditional event system)
s = """The quick brown fox jumped over
some lazy text sample.
He wasn't really into numbers, but it's
really important to try out number rendering:
This number goes from 0 to π in no time!
And then back to 0 again... Wow!
This is real crazy stuff,
but it gets even more ludicrous:
∮ E⋅da = Q,  n → ∞, ∑ f(i) = ∏ g(i),
∀x∈ℝ: ⌈x⌉ = −⌊−x⌋, α ∧ ¬β = ¬(¬α ∨ β),
ℕ ⊆ ℕ₀ ⊂ ℤ ⊂ ℚ ⊂ ℝ ⊂ ℂ,
⊥ < a ≠ b ≡ c ≤ d ≪ ⊤ ⇒ (A ⇔ B),
2H₂ + O₂ ⇌ 2H₂O, R = 4.7 kΩ, ⌀
I can't even...
"""

# _view and visualize it!
_view(visualize(s,
    model=translationmatrix(Vec3f0(0,600,0)), # move this up, since the text starts at 0 and goes down from there
), window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using NumCME
using Test

TestFixture = NamedTuple{(:stoich_mat, :x0)}

## Bursting gene model
bursting_gene = TestFixture(
    (
        [[-1, 1, 0] [1, -1, 0] [0, 0, 1] [0, 0, -1]],
        [1,0,0]
    )
)
# Test constructors
@test stateset = StateSpaceSparse(bursting_gene.stoich_mat, bursting_gene.x0) isa Any
@test stateset = StateSpaceSparse(bursting_gene.stoich_mat, [bursting_gene.x0]) isa Any

# Test that state exploration routines return state spaces with correct sizes
bursting_statespace = StateSpaceSparse(bursting_gene.stoich_mat, bursting_gene.x0)
@test get_state_count(bursting_statespace) == 1
@test expand!(bursting_statespace, 0) isa Any
@test get_state_count(bursting_statespace) == 1

expand!(bursting_statespace, 1)
@test get_state_count(bursting_statespace) == 3

bursting_statespace = StateSpaceSparse(bursting_gene.stoich_mat, bursting_gene.x0)
expand!(bursting_statespace, 3)
@test get_state_count(bursting_statespace) == 7

bursting_statespace = StateSpaceSparse(bursting_gene.stoich_mat, bursting_gene.x0)
expand!(bursting_statespace, 4)
bursting_expected_states = sort([[1, 0, 0], [0, 1, 0], [1, 0, 1], [0, 1, 1], [1, 0, 2], [0, 1, 2], [1, 0, 3], [0, 1, 3], [1, 0, 4]])
bursting_generated_states = sort(bursting_statespace.states)
@test bursting_expected_states == bursting_generated_states

## Toggle-switch model
toggle_switch = TestFixture(
    (
    [[1; 0] [-1; 0] [0; 1] [0; -1]],
    [0, 0]
)
)

# Test constructors
@test stateset = StateSpaceSparse(toggle_switch.stoich_mat, toggle_switch.x0) isa Any
@test stateset = StateSpaceSparse(toggle_switch.stoich_mat, [toggle_switch.x0]) isa Any

# Test that state exploration routines return state spaces with correct sizes
toggle_statespace = StateSpaceSparse(toggle_switch.stoich_mat, toggle_switch.x0)
@test get_state_count(toggle_statespace) == 1
@test expand!(toggle_statespace, 0) isa Any
@test get_state_count(toggle_statespace) == 1

expand!(toggle_statespace, 1)
@test get_state_count(toggle_statespace) == 3

toggle_statespace = StateSpaceSparse(toggle_switch.stoich_mat, toggle_switch.x0)
expand!(toggle_statespace, 3)
@test get_state_count(toggle_statespace) == 10
toggle_expected_states = sort([[0, 0], [1, 0], [2, 0], [3, 0], [0, 1], [1,1], [2, 1], [0, 2], [1, 2], [0, 3]])
toggle_generated_states = sort(get_states(toggle_statespace))
@test toggle_expected_states == toggle_generated_states


text:
@doc raw"""
```
Z,p,q = retrieve_surface(algorithm::DiscreteShape, img::AbstractArray, iterations::Int=2000)
```
Attempts to produce a heightmap from a grayscale image by minimization of a set
of Euler-Lagrange equations. This is done discretely at each point in the image
utilizing the second derivates of the surface normals and their Fourier Transforms.
# Output
Returns an M by N array of Float `Z` that represents the reconstructed height at
the point and the gradients in M by N arrays of Float 'p' and 'q'.
# Details
The algorithm attempts to minimise the brightness deviation ``ϵ`` through the
brightness constraint ``\epsilon_1`` and smoothness constraint ``\epsilon_2``
as defined bellow:
```math
\epsilon=\iint((\epsilon_1+\lambda\epsilon_2))dxdy=\iint((E(x,y)-R(p,q))^2+
\lambda(p_x^2+q_x^2+p_y^2+q_y^2))
```
This is minimised using the Euler-Lagrange equations defined as;
```math
\dfrac{\delta\epsilon}{\delta p}-\dfrac{\delta}{\delta x}\dfrac{\delta\epsilon}{
\delta p_x}-\dfrac{\delta}{\delta y}\dfrac{\delta\epsilon}{\delta p_y}=0
```
and
```math
\dfrac{\delta
\epsilon}{\delta q}-\dfrac{\delta}{\delta x}\dfrac{\delta\epsilon}{\delta q_x}-
\dfrac{\delta}{\delta y}\dfrac{\delta\epsilon}{\delta q_y}=0
```
which become:
```math
\begin{gathered}
-2(E-R)\dfrac{\delta R}{\delta p}-2\lambda p_{xx}-2\lambda p_{yy}=0\\
-2(E-R)\dfrac{\delta R}{\delta p}-2\lambda q_{xx}-2\lambda q_{yy}=0
\end{gathered}
```
which can be further simplified to give:
```math
\begin{gathered}
\nabla^2p=\dfrac{1}{\lambda}(R-E)\dfrac{\delta R}{\delta p}\\
\nabla^2q=\dfrac{1}{\lambda}(R-E)\dfrac{\delta R}{\delta q}
\end{gathered}
```
where ``\nabla^2p=p_{xx}+p_{yy}`` and ``\nabla^2q=q_{xx}+q_{yy}`` are Laplacians
of p and q.

However, for computation we are dealing with a discreate case of these equations
which can be defined as below:
```math
\begin{gathered}
p_{i,j}=\bar{p}_{i,j}+\dfrac{1}{4\lambda}(E-R)\dfrac{\delta R}{\delta p}\\q_{i,j}=
\bar{q}_{i,j}+\dfrac{1}{4\lambda}(E-R)\dfrac{\delta R}{\delta q}\\
\end{gathered}
```
where
```math
\bar{p}_{i,j}=\dfrac{p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j-1}}{4}
```
and
```math
\bar{q}_{i,j}=\dfrac{q_{i+1,j}+q_{i-1,j}+q_{i,j+1}+q_{i,j-1}}{4}
```
Finally, the algorithm needs to enforce integrability on p and q and retrieve
the surface Z. This can be done by taking the Fast Fourier Transform of p and q
()``c_p(\omega_x,\omega_y)`` and ``c_q(\omega_x,\omega_y)``) and then using the
Inverse Fast Fourier Transform to recover Z and update p and q as per bellow:
```math
\begin{gathered}
p=\sum c_p(\omega_x,\omega_y)e^{j(\omega_xx+\omega_yy)}\\
q=\sum c_q(\omega_x,\omega_y)e^{j(\omega_xx+\omega_yy)}\\Z=\sum c(\omega_x,\omega_y)
e^{j(\omega_xx+\omega_yy)}\\
\end{gathered}
```
where
```math
c(\omega_x,\omega_y)=\dfrac{-j(\omega_xc_p(\omega_x,\omega_y)+\omega_yc_q(\omega_x
,\omega_y))}{\omega_x^2+\omega_y^2}
```

The `albedo` and `illumination_direction` can be manually defined using the
function signature:
```
Z,p,q = retrieve_surface(algorithm::DiscreteShape, img::AbstractArray, albedo::Real, illumination_direction::Vector{T} where T <: Real, iterations::Int=2000)
```
Note: if `albedo` and `illumination_direction` are not supplied they will be
calculated at runtime using `estimate_img_properties`.
# Arguments
The function arguments are described in more detail below.
##  `img`
An `AbstractArray` storing the grayscale value of each pixel within
the range [0,1].
## `albedo`
A `Real` that specifies the albedo (amount of light reflected) of the image. If
`albedo` is specified to must the `illumination_direction`.
## `illumination_direction`
A `Vector{T} where T <: Real` that specifies the tilt value to be used by the
algorithm. The `illumination_direction` should be a vector of the form [x,y,z]
where x,y,z are int he range [0,1]. If `illumination_direction` is specified
to must the `albedo`.
## `iterations`
An `Int` that specifies the number of iterations the algorithm is to perform. If
left unspecified a default value of 2000 is used.
## `smoothness`
An `Int` that specifies the strength of the smoothness constraint in the minimised
function.
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(), 1, [0.2,0,0.9], radius = 5)

#calculate the heightmap (using 500 iterations)
Z,p,q = retrieve_surface(DiscreteShape(), img, 500)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.
"""
function retrieve_surface(algorithm::DiscreteShape, img::AbstractArray, iterations::Int=2000; smoothness::Int=1000)
    ρ,I,σ,τ = estimate_img_properties(img)
    λ=smoothness
    return retrieve_surface(DiscreteShape(), img, ρ, I, iterations, smoothness=λ)
end

function retrieve_surface(algorithm::DiscreteShape, img::AbstractArray, albedo::Real, illumination_direction::Vector{T} where T <: Real, iterations::Int=2000; smoothness::Int=1000)
    ρ = albedo
    I = illumination_direction
    E = Array{Float64}(img)

    #initialize variables
    p = zeros(Complex{Float64},axes(E))
    q = zeros(Complex{Float64},axes(E))
    R = zeros(Complex{Float64},axes(E))
    Z = zeros(axes(E))
    return solve_EulerLagrange(ρ, I, iterations, p, q, R, smoothness, E, Z)
end


text:
#Finding PROPERTIES WITH NUMERICAL VALUES in Trivial inputs
tpropwithnumval = Set();
@info "Finding trivials with numerical value"
for p in coolproptrivialparameters
    missed = Set();
    for fluid in coolpropfluids
        try
            res = ("$(PropsSI(p, Compat.String(fluid)))");
            push!(tpropwithnumval, p);
        catch err
            push!(missed, fluid);
        end
    end
    if (length(missed) == length(coolpropfluids))
    # if (findfirst(fails_any_props_trivals, p) == 0)
        if findfirst(x->x==p, fails_any_props_trivals) == nothing
            @warn "None of the defined fluids return a number for $p (expected to have numerical value.)";
        else
            println("$p is trivial without any numericalvalue");
        end
    else
        length(missed) > 0 && println("Only $(collect(setdiff(coolpropfluids, missed))) return number for $p");
    end
end
if ((tpropwithnumval == Set(trivalwithnumval)) && (Set(setdiff(coolproptrivialparameters, tpropwithnumval)) == Set(fails_any_props_trivals)))
    @info "Trivial inputs with numerical values are as expected.";
else
    @warn "Trivial inputs with numerical values are not as expected.";
end


text:
using Distributions
a = 10
N = 1000
x = ones(N)
y = ones(N)
for t = 2:N
    u1 = rand() * exp(-abs(x[t-1]))
    u2 = rand() * exp(-abs(y[t-1]))
    u3 = rand() * exp(-a*abs(x[t-1]-y[t-1]))
    x[t] = rand(Uniform(log(u1), -log(u1)))
    while a*abs(x[t]-y[t-1]) > -log(u3)
        x[t] = rand(Uniform(log(u1), -log(u1)))
    end
    y[t] = rand(Uniform(log(u2), -log(u2)))
    while a*abs(x[t]-y[t]) > -log(u3)
        y[t] = rand(Uniform(log(u2), -log(u2)))
    end
end
using Plots
plot(x, y, seriestype=:scatter, legend=false)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# rangebars

# Makie
using CairoMakie

f = Figure()
Axis(f[1, 1])

vals = -1:0.1:1
lows = zeros(length(vals))
highs = LinRange(0.1, 0.4, length(vals))

rangebars!(vals, lows, highs, color = LinRange(0, 1, length(vals)),
    whiskerwidth = 10, direction = :x, linewidth=2)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
import Pkg;
Pkg.activate(".")
# %%
path = pwd() * "/functions"
push!(LOAD_PATH,path)
# %%
using Revise
using Calibration
using JLD2, FileIO
# %%
cal_y = calibrate_y()
# %%
cal_m = calibrate_m()
# %%
cal_c = calibrate_c()
# %%
cal_ssξ = calibrate_ssξ()
# %%
cal_u = calibrate_u()
# %%
JLD2.@save "output/calibrations.jld2" cal_y cal_c cal_m cal_ssξ cal_u


text:
#=
    The tree plotting functions are adapted from the tree plotting functionality
    of the EcoJulia Phylo package.

    Further documentation and the original functions can be found here:
        https://github.com/EcoJulia/Phylo.jl
=#
"""
    f(root::T; <keyword arguments>) where T<:GerneralNode

Recipe that handles plotting of MCPhylo Trees. Takes the root node as input.

# Arguments
- 'treetype::Symbol=:dendrogram' :  specifies how the tree is visualized.
                                    Options :fan and :dendrogram
- 'showtips::Bool=:true' : controls if leaf labels are shown or not
- 'tipfont::Union{Tuple{Int64}, Tuple{Int64,Symbol}}=(7,)' :
                specify the font and font colour of the annotations
- 'marker_group=nothing' : currently not supported
- 'line_group=nothing' : currently not supported
"""
@recipe function f(root::T; treetype=:dendrogram, marker_group=nothing,
                   line_group=nothing, showtips=true, tipfont=(7,)
                   ) where T<:AbstractNode

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

    lz = get(plotattributes, :line_z, nothing)
    mz = get(plotattributes, :marker_z, nothing)
    isnothing(lz) || (line_z := _handlez(lz, root))
    isnothing(mz) || (marker_z := _handlez(mz, root))
    mg = _handlez(marker_group, root)
    lg = _handlez(line_group, root)

    d, h, n = _findxy(root)
    adj = 0.03maximum(values(d))
    tipannotations = map(x->(d[x] + adj, h[x], x.name), get_leaves(root))
    x, y = Float64[], Float64[]
    for node ∈ pre_order(root)
        if !node.root
            m = get_mother(node)
            push!(x, d[m], d[m], d[node], NaN)
            push!(y, h[m], h[node], h[node], NaN)
        end
    end

    marker_x, marker_y = _handlemarkers(plotattributes, mg, root, d, h)

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

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

@recipe function f(dend::Dendrogram)
    ex = extrema(filter(isfinite, dend.x))
    xlims --> (ex[1] - 0.05 * ex[2], ex[2] * 1.15)

    sa = get(plotattributes, :series_annotations, nothing)
    dend.showtips && (annotations := map(x -> (x[1], x[2], (x[3], :left, dend.tipfont...)), dend.tipannotations))
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing
        label --> ""
    #    primary := false

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

        dend.x, dend.y
    end
    if !isempty(dend.marker_x) || sa !== nothing
        if isnothing(dend.marker_group)
            @series begin
                seriestype := :scatter
                sa !== nothing && (series_annotations := sa)
                label --> ""
                dend.marker_x, dend.marker_y
            end
        else
            groups = sort(unique(dend.marker_group))
            for group in groups
                idxs = findall(==(group), dend.marker_group)
                @series begin
                    seriestype := :scatter
                    sa !== nothing && (series_annotations := sa[idxs])
                    label --> string(group)
                    dend.marker_x[idxs], dend.marker_y[idxs]
                end
            end
        end
    end
    primary = false
    label = ""
    nothing
end

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

    sa = get(plotattributes, :series_annotations, nothing)
    aspect_ratio := 1
    mx = maximum(filter(isfinite, fan.x))
    if fan.showtips
        xlims --> (1.5 .* (-mx, mx))
        ylims --> (1.5 .* (-mx, mx))
        annotations := map(x -> (_tocirc(x[1], adjust(x[2]))..., (x[3], :left,
            rad2deg(adjust(x[2])), fan.tipfont...)), fan.tipannotations)
    end
    @series begin
        seriestype := :path
        markersize := 0
        markershape := :none
        series_annotations := nothing
        label := ""

        x, y = _circle_transform_segments(fan.x, adjust(fan.y))
        lc = _extend(get(plotattributes, :linecolor, nothing), x)
        lc !== nothing && (linecolor := lc)
        la = _extend(get(plotattributes, :linealpha, nothing), x)
        la !== nothing && (linealpha := la)
        lz = _extend(get(plotattributes, :line_z, nothing), x)
        lz !== nothing && (line_z := lz)
        x, y
    end
    if !isempty(fan.marker_x) || sa !== nothing
        if isnothing(fan.marker_group)
            @series begin
                seriestype := :scatter
                sa !== nothing && (series_annotations := sa)
                label --> ""
                _xcirc.(adjust(fan.marker_y), fan.marker_x), _ycirc.(adjust(fan.marker_y), fan.marker_x)
            end
        else
            groups = sort(unique(fan.marker_group))
            for group in groups
                idxs = findall(==(group), fan.marker_group)
                @series begin
                    seriestype := :scatter
                    sa !== nothing && (series_annotations := sa[idxs])
                    label --> string(group)
                    _xcirc.(adjust(fan.marker_y[idxs]), fan.marker_x[idxs]), _ycirc.(adjust(fan.marker_y[idxs]), fan.marker_x[idxs])
                end
            end
        end
    end
    nothing
end

_handlez(x, root) = x
# need to add this back if i.e. gradiently coloured plots are needed
# _handlez(x::Union{String, Symbol}, tree) = getnodedata.(tree, traversal(tree, preorder), x)
_mylength(x) = 1
_mylength(x::AbstractVector) = length(x)

function _handlemarkers(plotattributes, marker_group, root, d, h)
    marker_x, marker_y = Float64[], Float64[]
    markerfields = filter(x->occursin(r"marker", String(x)), keys(plotattributes))
    isempty(markerfields) && isnothing(marker_group) && return(marker_x, marker_y)
    maxlengthfields = isempty(markerfields) ? 1 : maximum([_mylength(plotattributes[k]) for k in markerfields])
    maxlengthgroup = isnothing(marker_group) ? 1 : length(marker_group)
    maxlength = max(maxlengthfields, maxlengthgroup)
    f = maxlength ∈ (1, count(x-> root.nchild != 0, pre_order(root))) ? filter(x -> x.nchild != 0, pre_order(root)) : pre_order(root)
    append!(marker_x, getindex.(Ref(d), f))
    append!(marker_y, getindex.(Ref(h), f))
    marker_x, marker_y
end

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

function _findxy(root::T)::Tuple{Dict{T, Float64}, Dict{T, Float64}, Vector{String}} where T<:AbstractNode

    # two convenience recursive functions using captured variables
    function findheights!(node::T) where T<:AbstractNode
        if !in(node, keys(height))
            for child in node.children
                findheights!(child)
            end
        end
        if node.nchild != 0
            ch_heights = [height[child] for child in node.children]
            height[node] = (maximum(ch_heights) + minimum(ch_heights)) / 2.
        end
    end

    function finddepths!(node::T, parentdepth::Float64 = 0.0) where T<:AbstractNode
        mydepth = parentdepth
        push!(names, node.name)
        mydepth += node.inc_length
        depth[node] = mydepth
        for child in node.children
            finddepths!(child, mydepth)
        end
    end

    # root_name = root.name / num
    height = Dict(tip => float(i) for (i, tip) in enumerate(x for x in get_leaves(root)))
    nnodes = length(post_order(root))
    sizehint!(height, nnodes)
    findheights!(root)

    depth = Dict{T, Float64}()
    ### Kann man hier node.num stattdessen verwenden, oder Labels hängen später
    ### von dem name Array ab
    names = String[]
    sizehint!(depth, nnodes)
    sizehint!(names, nnodes)
    finddepths!(root)

    depth, height, names
end

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

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

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


text:
# --- PRIVATE METHODS ---------------------------------------------------------- #
function check_orgamism_code(organism_code::String)
end
# ------------------------------------------------------------------------------ #

# --- PUBLIC METHODS ----------------------------------------------------------- #
function get_gene_list_for_organism(organism_code::String)

    # check the organism_id -
    check_orgamism_code(organism_code)

end

function get_pathway_list_for_organism(organism_code::String)::Dict{String,String}

    # check the organism_id -
    check_orgamism_code(organism_code)
end

function build_metabolite_matching_table()::Union{Dict{String,String},Nothing}

    # initialize -
    metabolite_matching_table = Dict{String,String}()

    # build the url -
    url_string = "http://rest.kegg.jp/list/compound"

    # make the call -
    http_body = http_get_call_with_url(url_string)

    # ok, so no body?
    if (length(http_body) == 0)
        return nothing
    end

    # populate the table -
    records_array = split(http_body,"\n")
    for record in records_array

        # records take the form:
        # kegg_compound_id\t{Names;*}

        if (record != "")
            # split around the \t -
            kegg_compound_id = string(split(record,"\t")[1])
            compound_name_array = split(record,"\t")[2]
            if (occursin(";",compound_name_array) == true)

                # ok, so we have multiple possible names -
                synomym_array = split(compound_name_array,";")
                for synonum in synomym_array
                    metabolite_matching_table[lstrip(string(synonum))] = kegg_compound_id
                end
            else
                metabolite_matching_table[lstrip(string(compound_name_array))] = kegg_compound_id
            end
        end
    end

    # return -
    return metabolite_matching_table
end

function get_organism_codes()::Union{Dict{String,Organism},Nothing}

    # initialize -
    organism_dictionary = Dict{String,Organism}()

    # formulate the URL -
    url_string = "http://rest.kegg.jp/list/organism"

    # make the call -
    http_body = http_get_call_with_url(url_string)

    # ok, so no body?
    if (length(http_body) == 0)
        return nothing
    end

    # split along the newline -
    tokenized_body = split(http_body,"\n")
    for record in tokenized_body

        # record -
        # {Txxx}\t{organism code}\t{description}\t{species_taxonomy}

        # ok, we need to cut again, along the \t -
        field_array = split(record,"\t")

        # check, do we have 4 elements -
        if (length(field_array) == 4)

            # create new organism -
            organism_wrapper = Organism()
            organism_wrapper.organism_id = field_array[1]
            organism_wrapper.species_description = field_array[3]
            organism_wrapper.species_taxonomy = field_array[4]

            # # grab the organism_code -
            organism_code = field_array[2]
            organism_wrapper.organism_code = organism_code

            # cache -
            organism_dictionary[organism_code] = organism_wrapper
        end
    end

    # return -
    return organism_dictionary
end
# ------------------------------------------------------------------------------ #


text:
"""
    GaborKernel(; ell::Real=1.0, p::Real=1.0)

Gabor kernel with length scale ell and period p. Given by
```math
    κ(x,y) =  h(x-z), h(t) = exp(-sum(t.^2./(ell.^2)))*cos(pi*sum(t./p))
```

"""
struct GaborKernel{K<:Kernel} <: BaseKernel
    kernel::K
    function GaborKernel(;ell=nothing, p=nothing)
        k = _gabor(ell=ell, p=p)
        new{typeof(k)}(k)
    end
end

(κ::GaborKernel)(x, y) = κ.kernel(x ,y)

function _gabor(; ell = nothing, p = nothing)
    if ell === nothing
        if p === nothing
            return SqExponentialKernel() * CosineKernel()
        else
            return SqExponentialKernel() * transform(CosineKernel(), 1 ./ p)
        end
    elseif p === nothing
        return transform(SqExponentialKernel(), 1 ./ ell) * CosineKernel()
    else
        return transform(SqExponentialKernel(), 1 ./ ell) * transform(CosineKernel(), 1 ./ p)
    end
end

function Base.getproperty(k::GaborKernel, v::Symbol)
    if v == :kernel
        return getfield(k, v)
    elseif v == :ell
        kernel1 = k.kernel.kernels[1]
    	if kernel1 isa TransformedKernel
    	    return 1 ./ kernel1.transform.s[1]
    	else
    	    return 1.0
    	end
    elseif v == :p
        kernel2 = k.kernel.kernels[2]
        if kernel2 isa TransformedKernel
            return 1 ./ kernel2.transform.s[1]
        else
            return 1.0
        end
    else
        error("Invalid Property")
    end
end

Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")

function kernelmatrix(
    κ::GaborKernel,
    X::AbstractMatrix;
    obsdim::Int=defaultobs)
    kernelmatrix(κ.kernel, X; obsdim=obsdim)
end

function kernelmatrix(
    κ::GaborKernel,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int=defaultobs)
    kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
end

function kerneldiagmatrix(
    κ::GaborKernel,
    X::AbstractMatrix;
    obsdim::Int=defaultobs) #TODO Add test
    kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
end