text:
#
# Plot Restricted Two-body Problem trajectories
#

"""
Plot the orbital positions of a R2BP orbit.
"""
function plotpositions(traj::Trajectory{<:RestrictedTwoBodyOrbit}; lengthunit = lengthunit(first(traj)), kwargs...)

    pos = transpose(hcat(position_vector.(traj)...))
    pos = ustrip.(lengthunit, pos)
    LU  = string(lengthunit)

    defaults = (; title     = "R2BP Positions" * (first(traj).system.name == "" ? "" : " w.r.t. $(first(traj).system.name)'s Center"),
                  xlabel    = "X ($LU)",
                  ylabel    = "Y ($LU)",
                  zlabel    = "Z ($LU)",
                  label     = "Orbit Position",
                  formatter = :plain,
                  grid      = :on,
                  linewidth = 2,
                  dpi       = 150)

    options = merge(defaults, kwargs)

    plot(pos[:,1], pos[:,2], pos[:,3]; options...)

end

"""
Plot the orbital velocities of a R2BP orbit.
"""
function plotvelocities(traj::Trajectory{<:RestrictedTwoBodyOrbit}; velocityunit = velocityunit(first(traj)), kwargs...)

    vel = transpose(hcat(velocity_vector.(traj)...))
    vel = ustrip.(velocityunit, vel)
    VU  = string(velocityunit)

    defaults = (; title     = "R2BP Velocities" * (first(traj).system.name == "" ? "" : " w.r.t. $(first(traj).system.name)'s Center"),
                  xlabel    = "X ($VU)",
                  ylabel    = "Y ($VU)",
                  zlabel    = "Z ($VU)",
                  label     = "Orbit Velocity",
                  formatter = :plain,
                  grid      = :on,
                  linewidth = 2,
                  dpi       = 150)

    options = merge(defaults, kwargs)

    plot(vel[:,1], vel[:,2], vel[:,3]; options...)

end

"""
Plot the positions of an orbit.
"""
function plotpositions(pos::AbstractVector{V}; lengthunit = u"km", exclude_z = false, kwargs...) where V <: AbstractVector{<:Real}

    pos = transpose(hcat(pos...))
    LU  = string(lengthunit)

    defaults = (; title     = "Orbit Positions",
                  xlabel    = "X ($LU)",
                  ylabel    = "Y ($LU)",
                  zlabel    = "Z ($LU)",
                  label     = "Orbit Position",
                  formatter = :plain,
                  grid      = :on,
                  linewidth = 2,
                  dpi       = 150)

    options = merge(defaults, kwargs)

    if !exclude_z
        plot(pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        plot(pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the positions of an orbit.
"""
function plotpositions(pos::AbstractVector{<:V}; lengthunit = unit(pos[1][1]), exclude_z = false, kwargs...) where V <: AbstractVector{<:Unitful.Length}
    return plotpositions(map(p->ustrip.(lengthunit, p), pos); lengthunit = lengthunit, exclude_z = exclude_z, kwargs...)
end

"""
Plot the positions of an orbit.
"""
function plotpositions!(fig, pos::AbstractVector{V}; lengthunit = u"km", exclude_z = false, kwargs...) where V <: AbstractVector{<:Real}

    pos = transpose(hcat(pos...))
    LU  = string(lengthunit)

    defaults = (; label=:none)

    options = merge(defaults, kwargs)

    if !exclude_z
        plot!(fig, pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        plot!(fig, pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the positions of an orbit.
"""
function plotpositions!(fig, pos::AbstractVector{V}; lengthunit = unit(pos[1][1]), exclude_z = false, kwargs...) where V <: AbstractVector{<:Unitful.Length}
    return plotpositions!(fig, map(p->ustrip.(lengthunit, p), pos); lengthunit = lengthunit, exclude_z = exclude_z, kwargs...)
end

"""
Plot the positions of an orbit.
"""
function plotpositions(pos::AbstractMatrix{<:Real}; lengthunit = u"km", exclude_z = false, kwargs...)

    size(pos,2) == 3 || (pos = transpose(pos))
    LU  = string(lengthunit)

    defaults = (; title     = "Orbit Positions",
                  xlabel    = "X ($LU)",
                  ylabel    = "Y ($LU)",
                  zlabel    = "Z ($LU)",
                  label     = "Orbit Position",
                  formatter = :plain,
                  grid      = :on,
                  linewidth = 2,
                  dpi       = 150)

    options = merge(defaults, kwargs)

    if !exclude_z
        plot(pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        plot(pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the positions of an orbit.
"""
function plotpositions(pos::AbstractMatrix{<:Unitful.Length}; lengthunit = unit(pos[1]), exclude_z = false, kwargs...)
    return plotpositions(ustrip.(lengthunit, pos); lengthunit = lengthunit, exclude_z = exclude_z, kwargs...)
end

"""
Plot the positions of an orbit.
"""
function plotpositions!(fig, pos::AbstractMatrix{<:Real}; lengthunit = u"km", exclude_z = false, kwargs...)

    size(pos,2) == 3 || (pos = transpose(pos))
    LU  = string(lengthunit)

    defaults = (; label=:none)

    options = merge(defaults, kwargs)

    if !exclude_z
        plot!(fig, pos[:,1], pos[:,2], pos[:,3]; options...)
    else
        plot!(fig, pos[:,1], pos[:,2]; options...)
    end
end

"""
Plot the positions of an orbit.
"""
function plotpositions!(fig, pos::AbstractMatrix{<:Unitful.Length}; lengthunit = unit(pos[1]), exclude_z = false, kwargs...)
    return plotpositions!(fig, ustrip.(lengthunit(pos[1]), pos); lengthunit = lengthunit, exclude_z = exclude_z, kwargs...)
end


text:
module DistanceTransforms

using ImageMorphology
using CUDA
using FLoops
using FoldsCUDA

include("./chamfer.jl")
include("./euclidean.jl")
include("./squared_euclidean.jl")
include("./utils.jl")

export transform,
    transform!,

    # Export chamfer.jl functions
    Chamfer,

    # Export euclidean.jl functions
    euclidean,

    # Export squared_euclidean.jl functions
    SquaredEuclidean,

    # Export utils.jl functions
    boolean_indicator
end

ext: jl


text: # import PGFPlots


text:
using HerbLanguage
using HerbProlog
using Test

include("swipl_test.jl")


text:
## materialize
@inline mtab_materialize(x, ::Integer) = x
@inline mtab_materialize(bc::Broadcasted, TN::Integer) = mtab_copy(instantiate(bc), TN)

@inline mtb_materialize!(dest::Tuple{Vararg{AbstractArray}}, bc::Broadcasted, TN::Integer) =
    (mtb_materialize!(TupleDummy(dest), bc, TN); dest)

## copyto
@inline mtab_copy(bc::Broadcasted{Style{Tuple}}, ::Integer) = tab_copy(bc)
@inline mtab_copy(bc::FilledBC, ::Integer) = copy(bc)
@inline mtab_copy(bc::Broadcasted, TN::Integer) = (mtb_copyto!(toa_similar(bc), bc, TN); dest)


text:
# See Table 4.1 of Pozar, Microwave Engineering

"""
    series_network(Z::CircuitParams)
Promote `CircuitParams` to `ABCDparams` assuming series connection.
```
        ┌────┐
○───────┤ Z  ├───────○
port 1  └────┘  port 2
○────────────────────○
```
"""
series_network(Z::CircuitParams{T}) where {T<:Real} =
    ABCDparams([1 convert(Impedance, Z).data; 0 1])
"""
    series_network(D::CircuitData{T}) where {T<:CircuitParams}
Promote `Vector{CircuitParams}` to `Vector{ABCDparams}` assuming series
connection.
"""
series_network(Z::Vector{T}) where {T<:CircuitParams} =
    [series_network(Z[idx]) for idx in 1:length(Z)]
"""
    series_network(D::CircuitData{T}) where {T<:CircuitParams}
Promote `CircuitData` to `NetworkData{ABCDparams}` assuming series connection.
"""
series_network(D::CircuitData{S, T}; Z0=50.0) where {S<:Real, T<:CircuitParams} =
    NetworkData(D.frequency, series_network(D.params); Z0=Z0)

"""
    parallel_network(Y::CircuitParams)
Promote `CircuitParams` to `ABCDparams` assuming parallel connection.
```
○──────────┬─────────○
         ┌─┴─┐
port 1   | Y |  port 2
         └─┬─┘
○──────────┴─────────○
```
"""
parallel_network(Y::CircuitParams) =
    ABCDparams([1 0; convert(Admittance, Y).data 1])
"""
    parallel_network(D::CircuitData{T}) where {T<:CircuitParams}
Promote `Vector{CircuitParams}` to `Vector{ABCDparams}` assuming parallel
connection.
"""
parallel_network(Y::Vector{T}) where {T<:CircuitParams} =
    [parallel_network(Y[idx]) for idx in 1:length(Y)]
"""
    parallel_network(D::CircuitData{T}) where {T<:CircuitParams}
Promote `CircuitData` to `NetworkData{ABCDparams}` assuming parallel connection.
"""
parallel_network(D::CircuitData{S, T}; Z0=50.0) where {S<:Real, T<:CircuitParams} =
    NetworkData(D.frequency, parallel_network(D.params); Z0=Z0)

"""
    terminated_network(Z::CircuitParams; Z0=50.0)
Promote `CircuitParams` to `Sparams` assuming termination with impedance `Z`.
```
     ┌─────────○
   ┌─┴─┐
   | Z |  port 1
   └─┬─┘
◁───┴─────────○
```
"""
function terminated_network(Z::Impedance{T}; Z0=50.0) where {T<:Number}
    d = [reflection_coefficient(Z0, convert(Impedance, Z).data) for i in 1:1,
        j in 1:1]
    return Sparams(d)
end
"""
    terminated_network(Z::Vector{T}; Z0=50.0) where {T<:CircuitParams}
Promote `Vector{CircuitParams}` to `Vector{Sparams}` assuming termination
with impedance `Z0`.
"""
terminated_network(Z::Vector{T}; Z0=50.0) where {T<:CircuitParams} =
    [terminated_network(Z[idx]; Z0=Z0) for idx in 1:length(Z)]
"""
    terminated_network(D::CircuitData{T}; Z0=50.0) where {T<:CircuitParams}
Promote `CircuitData` to `NetworkData{Sparams}` assuming terminated
connection to ground.
"""
terminated_network(D::CircuitData{S, T}; Z0=50.0) where {S<:Real, T<:CircuitParams} =
    NetworkData(D.frequency, terminated_network(D.params; Z0=Z0))

"""
    π_network(Y1::CircuitParams, Y2::CircuitParams, Y3::CircuitParams)
Promote `CircuitParams` Y1, Y2, and Y3 to `ABCDparams` assuming Π-network
configuration.
```
              ┌────┐
   ○──────┬───┤ Y3 ├───┬──────○
        ┌─┴─┐ └────┘ ┌─┴─┐
 port1  |Y1 |        |Y2 |  port2
        └─┬─┘        └─┬─┘
   ○──────┴────────────┴──────○
```
"""
function π_network(Y1::CircuitParams, Y2::CircuitParams, Y3::CircuitParams)
    _Y1, _Y2, _Y3 = [convert(Admittance, Y).data for Y in [Y1, Y2, Y3]]
    return ABCDparams([(1 + _Y2/_Y3) 1/_Y3;
        (_Y1 + _Y2 + _Y1 * _Y2/_Y3) (1 + _Y1/_Y3)])
end

"""
    π_network(Y1::Array{T, 1}, Y2::Array{S, 1}, Y3::Array{U, 1}) where
        {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams}
Promote arrays of `CircuitParams` Y1, Y2, and Y3 to those of `ABCDparams`
assuming Π-network configuration.
"""
π_network(Y1::Array{T, 1}, Y2::Array{S, 1}, Y3::Array{U, 1}) where
    {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams} =
    [π_network(Y1[idx], Y2[idx], Y3[idx]) for idx in 1:length(Y1)]

"""
    π_network(D1::CircuitData{T}, D2::CircuitData{S}, D3::CircuitData{U}) where
        {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams}
Promote `CircuitData` D1, D2, and D3 to `NetworkData{ABCDparams}` assuming
Π-network configuration. Supported only for data with same frequency range.
"""
π_network(D1::CircuitData, D2::CircuitData, D3::CircuitData) =
    check_frequency_identical(D1, D2, D3) ?
    NetworkData(D1.frequency, π_network(D1.params, D2.params, D3.params)) :
    error("`NetworkData` can only be constructed for `CircuitData` defined in same frequencies")

"""
    t_network(Z1::CircuitParams, Z2::CircuitParams, Z3::CircuitParams)
Promote `CircuitParams` Z1, Z2, and Z3 to `ABCDparams` assuming T-network
configuration.
```
      ┌────┐       ┌────┐
 ○────┤ Z1 ├───┬───┤ Z2 ├───○
      └────┘ ┌─┴─┐ └────┘
 port1       |Z3 |      port2
             └─┬─┘
 ○─────────────┴────────────○
```
"""
function t_network(Z1::CircuitParams, Z2::CircuitParams, Z3::CircuitParams)
    _Z1, _Z2, _Z3 = [convert(Impedance, Z).data for Z in [Z1, Z2, Z3]]
    return ABCDparams([(1 + _Z1/_Z3) (_Z1 + _Z2 + _Z1 * _Z2 / _Z3);
        1/_Z3 (1 + _Z2/_Z3)])
end

"""
    t_network(D1::Array{T, 1}, D2::Array{S, 1}, D3::Array{U, 1}) where
        {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams}
Promote arrays of `CircuitParams` Z1, Z2, and Z3 to those of `ABCDparams`
assuming T-network configuration.

"""
t_network(Z1::Array{T, 1}, Z2::Array{S, 1}, Z3::Array{U, 1}) where
    {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams} =
    [t_network(Z1[idx], Z2[idx], Z3[idx]) for idx in 1:length(Z1)]

"""
    t_network(D1::CircuitData{T}, D2::CircuitData{S}, D3::CircuitData{U}) where
        {T<:CircuitParams, S<:CircuitParams, U<:CircuitParams}
Promote `CircuitData` D1, D2, and D3 to `NetworkData{ABCDparams}` assuming
T-network configuration. Supported only for data with same frequency range.
"""
t_network(D1::CircuitData, D2::CircuitData, D3::CircuitData) =
    check_frequency_identical(D1, D2, D3) ?
    NetworkData(D1.frequency, Tnetwork(D1.params, D2.params, D3.params)) :
    error("`NetworkData` can only be constructed for `CircuitData` defined in same frequencies")


text:
using LinearAlgebra

function plus(a, b)
	return a + b
end

minus(a, b) = a - b

# TODO: Overloading is not implemented in the core yet, this should be reviewed or
# implemented ad-hoc in Julia loader when resolving the calls in the invoke
function plus(a)
	return a+1
end

function mul(a,b)
	return a*b
end

max_stars_repo_name: gargakshit/core

max_issues_repo_name: gargakshit/core

max_forks_repo_name: gargakshit/core


text:
using MolecularGraph:
    fastidentityfilter, fastsubstrfilter, emaptonmap

@testset "structurematch" begin

@testset "prefilter" begin
    BuOH = smilestomol("CCCO")
    butane = smilestomol("CCCC")
    propane = smilestomol("CCC")
    @test fastidentityfilter(butane, BuOH)
    @test !fastidentityfilter(propane, butane)
    @test fastsubstrfilter(butane, propane)
    @test !fastsubstrfilter(propane, butane)
end

@testset "structmatch" begin
    null = smilestomol("")
    @test !hasexactmatch(null, null)
    @test !hassubstructmatch(null, null)
    @test !hasnodeinducedmatch(null, null)
    @test !hasedgeinducedmatch(null, null)

    hexane = smilestomol("CCCCCC")
    iso = smilestomol("CC(C)CCC")
    iso2 = smilestomol("CCCC(C)C")
    cyclohexane = smilestomol("C1CCCCC1")
    @test !hasexactmatch(hexane, iso, prefilter=false)
    @test hasexactmatch(iso, iso2, prefilter=false)
    @test !hasexactmatch(cyclohexane, hexane, prefilter=false)
    @test hassubstructmatch(cyclohexane, hexane, prefilter=false)
    @test !hasnodeinducedmatch(cyclohexane, hexane, prefilter=false)
    @test hasedgeinducedmatch(cyclohexane, hexane, prefilter=false)

    tms = smilestomol("C[Si](C)(C)C")
    tsm = smilestomol("[Si]C([Si])([Si])[Si]")
    @test !hasexactmatch(tms, tsm, prefilter=false)
    @test !hasedgeinducedmatch(tms, tsm, prefilter=false)
    @test hasexactmatch(tms, tsm, prefilter=false, atommatcher=(m,q)->((m,q)->true))
    @test hasedgeinducedmatch(tms, tsm, prefilter=false, atommatcher=(m,q)->((m,q)->true))

    sulfide = smilestomol("CSSC")
    disconn = smilestomol("CS.SC")
    @test !hasexactmatch(sulfide, disconn, prefilter=false)
    @test hassubstructmatch(sulfide, disconn, prefilter=false)
    @test !hasnodeinducedmatch(sulfide, disconn, prefilter=false)
    @test hasedgeinducedmatch(sulfide, disconn, prefilter=false)

    # No edges
    nacl = smilestomol("[Na+].[Cl-]")
    na = smilestomol("[Na]")
    @test hassubstructmatch(nacl, na, prefilter=false)
    @test !hasedgeinducedmatch(nacl, na, prefilter=false)

    tetrahedrane = smilestomol("C12C3C1C23")
    fused = smilestomol("C1C2C1C2")
    spiro = smilestomol("C1CC12CC2")
    @test hassubstructmatch(tetrahedrane, fused, prefilter=false)
    @test !hassubstructmatch(tetrahedrane, spiro, prefilter=false)
    @test hasedgeinducedmatch(tetrahedrane, fused, prefilter=false)
    @test !hasedgeinducedmatch(tetrahedrane, spiro, prefilter=false)
end

@testset "connectedquery" begin
    hexane = smilestomol("CCCCCC")
    anyatom = parse(SMARTS, "*")
    @test hassubstructmatch(hexane, anyatom)

    aniline = smilestomol("C1=CC=CC=C1N")
    dieamine = smilestomol("CCNCC")
    priamine = parse(SMARTS, "[NX3;H2]")
    @test hassubstructmatch(aniline, priamine)
    @test !hassubstructmatch(dieamine, priamine)

    diether = smilestomol("CCOCC")
    phenol = smilestomol("c1ccccc1O")
    glycerol = smilestomol("OCC(O)CO")
    alcohol = parse(SMARTS, "[#6][OD]")
    @test !hassubstructmatch(diether, alcohol)
    @test hassubstructmatch(phenol, alcohol)
    @test hassubstructmatch(glycerol, alcohol)

    cyclopentane = smilestomol("C1CCCC1")
    pyrrole = smilestomol("N1C=CC=C1")
    aliphring = parse(SMARTS, "*@;!:*")
    @test hassubstructmatch(cyclopentane, aliphring)
    @test !hassubstructmatch(pyrrole, aliphring)

    triamine = smilestomol("CCN(CC)CC")
    acetamide = smilestomol("CC(=O)N")
    notamide = parse(SMARTS, raw"[NX3;!$(NC=O)]")
    @test hassubstructmatch(triamine, notamide)
    @test !hassubstructmatch(acetamide, notamide)

    po1 = smilestomol("COOC")
    npo1 = smilestomol("COCOCOC")
    peroxide = parse(SMARTS, "[OX2][OX2]")
    @test hassubstructmatch(po1, peroxide)
    @test !hassubstructmatch(npo1, peroxide)

    pyridine = smilestomol("n1ccccc1")
    pyrrole = smilestomol("[nH]1cccc1")
    sixmem = parse(SMARTS, "[*r6]1[*r6][*r6][*r6][*r6][*r6]1")
    @test hassubstructmatch(pyridine, sixmem)
    @test !hassubstructmatch(pyrrole, sixmem)
end

@testset "disconnectedquery" begin
    hetero1 = smilestomol("CCN(O)[O-]")
    hetero2 = smilestomol("CCN=N")
    hetero3 = smilestomol("BrCCOC#N")
    hetero4 = smilestomol("CCS(=O)(=O)O")
    disconn = parse(SMARTS, "[#7,#8].[!#6].N")
    @test hassubstructmatch(hetero1, disconn)
    @test !hassubstructmatch(hetero2, disconn)
    @test hassubstructmatch(hetero3, disconn)
    @test !hassubstructmatch(hetero4, disconn)
end

@testset "node matching" begin
    function nodematch(mol, query, idx=1, len=1)
        afunc = atommatch(mol, query)
        bfunc = bondmatch(mol, query)
        matches = edgesubgraph_isomorphisms(mol, query, nodematcher=afunc, edgematcher=bfunc)
        emaps = collect(matches)
        @test length(emaps) == len
        emap = emaps[idx]
        return emaptonmap(emap, mol, query)
    end
    function single(v)   # like `only` on Julia 1.4+
        length(v) == 1 || error("expected a single entry")
        return v[1]
    end
    # 6-membered carbon ring fused to a 5-membered carbon ring with an attached oxygen
    query = smartstomol("[#6]~1~[#6]~[#6]~[#6]~2~[#6](~[#6]~1)~[#6]~[#6]~[#6]~2~[#8]")
    # The carbon numbering looks like this:
    #         3     9
    #      2     4     8
    #      1     5     7
    #         6
    # 1-Indanone, PubChem CID6735
    mol = smilestomol("C1CC(=O)C2=CC=CC=C21")
    nmap = nodematch(mol, query)
    @test single.(nmap) == [8,7,6,5,10,9,1,2,3,4]  # indices in mol corresponding to atoms in query
    # 1-Acenaphthenone, PubChem CID75229
    mol = smilestomol("C1C2=CC=CC3=C2C(=CC=C3)C1=O")
    nmap = nodematch(mol, query)
    @test single.(nmap) == [11,10,9,8,7,6,2,1,12,13]
    # Ninhydrin, PubChem CID10236 (2 matches)
    mol = smilestomol("C1=CC=C2C(=C1)C(=O)C(C2=O)(O)O")
    nmaps = ([1,2,3,4,5,6,7,9,10,11], [2,1,6,5,4,3,10,9,7,8])
    nmap = nodematch(mol, query, 1, 2)
    snmap = single.(nmap)
    @test snmap ∈ nmaps
    nmap = nodematch(mol, query, 2, 2)
    snmap = single.(nmap)
    @test snmap ∈ nmaps

    # Ambiguous match
    query = smartstomol("[H][H]")
    mol = smilestomol("[H][H]")
    nmap = nodematch(mol, query)
    @test nmap == [[1,2], [1,2]]
end

@testset "MCS" begin
    null = smilestomol("")
    @test size(disconnectedmcis(null, null)) == 0

    BuOH = smilestomol("CCCO")
    butane = smilestomol("CCCC")
    @test size(disconnectedmcis(butane, BuOH)) == 3
    @test size(disconnectedmces(butane, BuOH)) == 2

    hexane = smilestomol("CCCCCC")
    cyclohexane = smilestomol("C1CCCCC1")
    @test size(disconnectedmcis(hexane, cyclohexane)) == 5
    @test size(disconnectedmces(hexane, cyclohexane)) == 5

    tms = smilestomol("C[Si](C)(C)C")
    tsm = smilestomol("[Si]C([Si])([Si])[Si]")
    @test size(disconnectedmcis(tms, tsm)) == 2
    @test size(disconnectedmces(tms, tsm)) == 1

    sulfide = smilestomol("CSSC")
    disconn = smilestomol("CS.SC")
    @test size(disconnectedmcis(sulfide, disconn)) == 3
    @test size(disconnectedmces(sulfide, disconn)) == 2

    # Delta-Y
    cyclopropane = smilestomol("C1CC1")
    isopropane = smilestomol("CC(C)C")
    @test size(disconnectedmcis(cyclopropane, isopropane)) == 2
    @test size(disconnectedmces(cyclopropane, isopropane)) == 2

    tetrahedrane = smilestomol("C12C3C1C23")
    fused = smilestomol("C1C2C1C2")
    spiro = smilestomol("C1CC12CC2")
    @test size(disconnectedmcis(tetrahedrane, fused)) == 3
    @test size(disconnectedmcis(tetrahedrane, spiro)) == 3
    @test size(disconnectedmces(tetrahedrane, fused)) == 5
    @test size(disconnectedmces(tetrahedrane, spiro)) == 4
end

end # structurematch


text:
module CatArrays
using Base

mutable struct CatArray{T,N,M,TA} <: AbstractArray{T,N} where TA <: AbstractArray
    # dimension over which the sub-arrays are concatenated
    dim::Int
    # tuple of all sub-arrays
    arrays::NTuple{M,TA}
    # size of all subarrays
    asize::Array{Int,2}
    # start indices of every subarrays in the combined array
    start::Array{Int,2}
    # size of the combined array
    sz::NTuple{N,Int}
end

"""
C = CatArray(dim,array1,array2,...)

Create a concatenated array view from a list of arrays. Individual
elements can be accessed by subscribs
"""
function CatArray(dim::Int,arrays...)
    na = length(arrays)

    # number of dimensions
    nd = ndims(arrays[1])
    if dim > nd
        nd = dim
    end

    asize = ones(Int,na,nd)
    start = ones(Int,na,nd)

    # get size of all arrays
    for i=1:na
        tmp = size(arrays[i])
        asize[i,1:length(tmp)] = [tmp...]
    end

    # check if dimensions are consistent
    ncd = 1:nd .!= dim

    for i=2:na
        if any(asize[i,ncd] .!= asize[1,ncd])
            error("Array number $i has inconsistent size")
        end
    end

    # start index of each sub-array

    for i=2:na
        start[i,:] = start[i-1,:]
        start[i,dim] = start[i,dim] + asize[i-1,dim]
    end

    sz = asize[1,:]
    sz[dim] = sum(asize[:,dim])

    return CatArray{eltype(arrays[1]),nd,na,typeof(arrays[1])}(
                    dim,
                    arrays,
                    asize,
                    start,
                    (sz...,))

end

function normalizeindexes(sz,indexes)
    ndims = length(sz)
    ind = Vector{StepRange}(undef,ndims)
    squeezedim = falses(ndims)

    # normalize indexes
    for i = 1:ndims
        indT = typeof(indexes[i])
        # :
        if indT == Colon
            ind[i] = 1:1:sz[i]
            # just a number
        elseif indT == Int
            ind[i] = indexes[i]:1:indexes[i]
            squeezedim[i] = true
            # range with a step equal to 1
        elseif indT == UnitRange{Int}
            ind[i] = first(indexes[i]):1:last(indexes[i])
        elseif indT == StepRange{Int,Int}
            ind[i] = indexes[i]
        else
            error("$indT: unsupported index")
        end
    end

    return (ind...,),squeezedim
end

#function Base.getindex(CA::CatArray{T,N},idx::Numbers...) where {T,N}

function Base.getindex(CA::CatArray{T,N},idx...) where {T,N}
    checkbounds(CA,idx...)

    ind,squeezedim = normalizeindexes(size(CA),idx)

    idx_global,idx_local,sz = idx_global_local_(CA,ind)

    B = Array{T,length(sz)}(undef,sz...)

    for j=1:length(CA.arrays)
        if prod(length.(idx_local[j])) > 0
            # get subset from j-th array
            subset = CA.arrays[j][idx_local[j]...]

            B[idx_global[j]...] = subset
        end
    end

    if all(squeezedim)
        # return just a scalar
        return B[1]
    elseif any(squeezedim)
        return dropdims(B,dims = (findall(squeezedim)...,))
    else
        return B
    end
end

Base.size(CA::CatArray) = CA.sz

function idx_global_local_(CA::CatArray,idx::NTuple{N}) where N

    n = length(size(CA))

    # number of indices must be equal to dimension
    @assert(length(idx) == n)

    lb = ones(Int,1,n) # lower bound
    ub = ones(Int,1,n) # upper bound

    # transform all colons to index range
    for i=1:n
        if (idx[i] == :)
            idx[i] = 1:CA.sz[i]
        end

        lb[i] = minimum(idx[i])
        ub[i] = maximum(idx[i])
    end

    sz = ntuple(i -> length(idx[i]),Val(N))

    idx_local  = Vector{NTuple{N,StepRange{Int,Int}}}(undef,length(CA.arrays))
    idx_global = Vector{NTuple{N,StepRange{Int,Int}}}(undef,length(CA.arrays))

    # loop over all arrays
    for j = 1:length(CA.arrays)
        idx_local_tmp = Vector{StepRange{Int,Int}}(undef,n)
        idx_global_tmp = Vector{StepRange{Int,Int}}(undef,n)

        # loop over all dimensions
        for i=1:n
            # rebase subscribt at CA.start[j,i]
            tmp = idx[i] .- (CA.start[j,i] - 1)

            # only indeces within bounds of the j-th array
            sel = (1 .<= tmp) .& (tmp .<= CA.asize[j,i])

            if sum(sel) == 0
                idx_local_tmp[i] = 1:0
                idx_global_tmp[i] = 1:0
            else
                # index for getting the data from the local j-th array
                idx_local_tmp[i] = tmp[findfirst(sel):findlast(sel)]
                idx_global_tmp[i] = (1:sum(sel)) .+ (findfirst(sel) - 1)
            end
        end

        idx_local[j] = (idx_local_tmp...,)
        idx_global[j] = (idx_global_tmp...,)
    end

    return idx_global,idx_local,sz

end

function Base.setindex!(CA::CatArray{T,N},data,idx...) where {T,N}
    ind,squeezedim = normalizeindexes(size(CA),idx)
    idx_global,idx_local,sz = idx_global_local_(CA,ind);
    @debug ind,idx_global,idx_local,sz

    data2 = reshape(data,sz)
    for j = 1:length(CA.arrays)
        # get subset from global array x
        subset = @view data2[idx_global[j]...]

        @debug idx_local[j]

        # set subset in j-th array
        CA.arrays[j][idx_local[j]...] = subset;
    end

    return data
end

export CatArray
end


text:
@testset "UnitRange" begin
    for R in (UnitMRange, UnitSRange)
        @testset "$R" begin
            r = R(1, 10)
            b = UnitRange(1, 10)
            rfloat = AbstractUnitRange{Float64}(r)
            @test eltype(rfloat) == Float64
            @test isa(rfloat, R)
            @test R{Int}(r) === r
            @test R{Float64}(r) == R(1., 10.)
            @test eltype(R{Int}(UnitRange(UInt(1), UInt(10)))) == Int
            @test R(UnitRange(UInt(1), UInt(10))) == R(UInt(1), UInt(10))
            @test first(r) == r.start
            @test last(r) == r.stop

            @test intersect(r, r[2]) == intersect(b, b[2])
            @test intersect(r, r[2]) == intersect(b, b[2])

            @test_throws ErrorException r.notfield
            if R == UnitMRange
                setproperty!(r, :start, 2)
                @test r == UnitMRange(2, 10)

                setproperty!(r, :stop, 8)
                @test r == UnitMRange(2, 8)

                @test_throws ErrorException setproperty!(r, :anything, 3)
            end
        end
    end
end


text:
# normalizecolumns(m::CuArray{T, 2}) where {T} = m ./ sqrt.(sum(m .^ 2, dims = 1)) .+ eps(T)
# normalizecolumns(m::TrackedArray{T, 2, CuArray{T, 2}}) where {T} = m ./ sqrt.(sum(m .^ 2, dims = 1)) .+ eps(T)

GPUArray{T, N} = Union{CuArray{T, N}, TrackedArray{T, N, CuArray{T, N}}}

samplehsuniform_gpu(size::Int...) = samplehsuniform_gpu(Float32, size...)
function samplehsuniform_gpu(T::Type, size::Int...)
    v = cuzeros(T, size...)
    randn!(v)
	normalizecolumns!(v)
end

function wloss(m::SVAEbase, x::GPUArray{T, 2}, β::T, d) where {T}
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	zp = samplehsuniform_gpu(T, size(z)...)
	Ω = d(z, zp)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * Ω
end

function samplez(m::SVAE, μz::GPUArray{T, 2}, κz::GPUArray{T, 2}) where {T}
	ω = sampleω(m, κz)
    v = samplehsuniform_gpu(T, size(μz, 1) - 1, size(μz, 2))
    ω2 = 1 .- CUDAnative.pow.(ω, 2f0) .+ eps(T)
	z = householderrotation(vcat(ω, CUDAnative.sqrt.(ω2) .* v), μz)
	return z
end

rand_gamma_cpu(k::Number, θ::Number, size::Int...) = rand_gamma_cpu(Float32, k, θ, size...)
rand_gamma_cpu(T::Type, k::Number, θ::Number, size::Int...) = rand_gamma(T::Type, k, θ, fill, zeros, size...)
function rand_gamma(T::Type, k, θ, fillfun, zerosfun, size...)
    k = T(k)
    θ = T(θ)
    if k < 1
        u = zerosfun(T, size...)
        rand!(u)
        return rand_gamma(T, 1 + k, θ, fillfun, zerosfun, size...) .* (u .^ (1 / k))
    end

    x = zerosfun(T, size...)
    v = zerosfun(T, size...)
    u = zerosfun(T, size...)

    d = k - 1 / T(3)
    c = 1 / (3 * sqrt(d))

    masku = fillfun(true, size...)

    while true
        maskv = copy(masku)
        while any(maskv)
            x[maskv] .= randn!(x[maskv])
            v[maskv] .= 1 .+ c .* x[maskv]
            maskv[maskv] .= v[maskv] .<= 0
        end

        @. v[masku] = v[masku] * v[masku] * v[masku]
        u[masku] .= rand!(u[masku])
        @. masku[masku] = masku[masku] & !(u[masku] < 1 - 0.0331 * x[masku] * x[masku] * x[masku] * x[masku])
        @. masku[masku] = masku[masku] & !(log(u[masku]) < (0.5 * x[masku] * x[masku] + d * (1 - v[masku] + log(v[masku]))))
        if !any(masku)
            break
        end
    end
    return θ .* d .* v
end

rand_gamma_gpu(k::Number, θ::Number, size::Int...) = rand_gamma_gpu(Float32, k, θ, size...)
function rand_gamma_gpu(T::Type, k, θ, size...)
    k = T(k)
    θ = T(θ)
    if k < 1
        u = cuzeros(T, size...)
        rand!(u)
        return rand_gamma_gpu(T, 1 + k, θ, size...) .* CUDAnative.pow.(u, (1 / k))
    end

    x = cuzeros(T, size...)
    v = cuzeros(T, size...)
    u = cuzeros(T, size...)

    d = k - 1 / T(3)
    c = 1 / (3 * sqrt(d))

    masku = cufill(true, size...)
    while any(masku)
        maskv = copy(masku)
        while any(maskv)
            v .= map((m, v, nv) -> m ? nv : v, maskv, v, next_gamma_v_sample(x, c))
            maskv .= maskv .& (v .<= 0)
        end
        x .= (v .- 1) ./ c

        v .= map((m, v, nv) -> m ? nv : v, masku, v, v .* v .*v)
        rand!(u)
        @. masku = masku & !(u < 1 - 0.0331 * x * x * x * x)
        @. masku = masku & !(CUDAnative.log(u) < (0.5 * x * x + d * (1 - v + CUDAnative.log(v))))
    end
    return θ .* d .* v
end

function next_gamma_v_sample(x, c)
    randn!(x)
    return 1 .+ c .* x
end

rand_beta_gpu(α::Number, β::Number, size::Int...) = rand_beta_gpu(Float32, α, β, size...)
function rand_beta_gpu(T::Type, α, β, size...)
    α = T(α)
    β = T(β)

    if (α > 1) || (β > 1)
        g1 = rand_gamma_gpu(T, α, 1, size...)
        g2 = rand_gamma_gpu(T, β, 1, size...)
        return g1 ./ (g1 .+ g2)
    end

    u = cuzeros(T, size...)
    v = cuzeros(T, size...)
    x = cuzeros(T, size...)
    y = cuzeros(T, size...)

    mask = cufill(true, size...)
    while any(mask)
        rand!(u)
        rand!(v)
        x .= map((m, x, nx) -> m ? nx : x, mask, x, CUDAnative.pow.(u, (1 / α)))
        y .= map((m, y, ny) -> m ? ny : y, mask, y, CUDAnative.pow.(v, (1 / β)))
        @. mask = mask & ((x + y) > 1)
    end

    return map((x, y) -> (x + y) > 0 ? x / (x + y) : log_beta_expression(CUDAnative.log(x), CUDAnative.log(y)), x, y)
end

function log_beta_expression(logX, logY)
    logM = logX > logY ? logX : logY
    logX -= logM
    logY -= logM
    return CUDAnative.exp(logX - CUDAnative.log(CUDAnative.exp(logX) + CUDAnative.exp(logY)))
end

function sampleω(model::SVAE, κ::GPUArray{T, N}) where {T, N}
    m = CuArray([T(model.zdim)])
    c2 = @. CUDAnative.pow(4κ, 2f0) + CUDAnative.pow((m - 1), 2f0)
    c = CUDAnative.sqrt.(c2)
	b = @. (-2κ + c) / (m - 1)
    a = @. (m - 1 + 2κ + c) / 4
    logm1 = CUDAnative.log.(m .- 1)
    d = @. (4 * a * b) / (1 + b) - (m - 1) * logm1
	ω = rejectionsampling(m[1], a, b, d, κ)
	return ω
end

function rejectionsampling(m, a, b, d, κ::GPUArray{T, N}) where {T, N}
    nϵ = cuzeros(T, size(a)...)
    ϵ = cuzeros(T, size(a)...)
    u = cuzeros(T, size(a)...)
    # ω = cuzeros(T, size(a)...)
    t = cuzeros(T, size(a)...)
    mask = cufill(true, size(a)...)

    while any(mask)
        nϵ .= rand_beta_gpu(T, (m - 1) / 2, (m - 1) / 2, size(a)...)
        ϵ .= map((m, x, nx) -> m ? nx : x, mask, ϵ, nϵ)
        rand!(u)
        t = @. 2 * a * b / (1 - (1 - b) * ϵ)
        mask = @. mask & ((m - 1) * CUDAnative.log(t) - t + d >= CUDAnative.log(u))
	end
	return @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

df = DataFrame!(CSV.File(scm_path("..", "data", "WaffleDivorce.csv"), delim=';'));
df = DataFrame(
  :s => df[:, :South],
  :a => df[:, :MedianAgeMarriage],
  :m => df[:, :Marriage],
  :w => df[:, :WaffleHouses],
  :d => df[:, :Divorce]
);

d = OrderedDict(
  :s => [:a, :m, :w],
  :a => [:m, :d],
  :m => [:d],
  :w => [:d]
);
u = []

dag = DAG("sr6.4.3", d; df=df);
fname = scm_path("..", "examples", "SR", "SR6.4.3", "sr6.4.3.dot")
#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

allpaths  = all_paths(dag, :w, :d)
backdoorpaths = backdoor_paths(dag, allpaths, :w)
adjustmentsets = adjustment_sets(dag, :w, :d)

@testset "sr6_4_3" begin

  @test length(allpaths) == 4
  @test allpaths == [
    [:w, :s, :a, :d],
    [:w, :s, :a, :m, :d],
    [:w, :s, :m, :d],
    [:w, :s, :m, :a, :d]
  ]
  @test length(backdoorpaths) == 0
  @test length(adjustmentsets) == 0

end


text:
# ---
# title: 542. 01 Matrix
# id: problem542
# author: Indigo
# date: 2021-06-23
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/01-matrix/description/>
# hidden: true
# ---
#
# Given a matrix consists of 0 and 1, find the distance of the nearest 0 for
# each cell.
#
# The distance between two adjacent cells is 1.
#
#
#
# **Example 1:**
#
#
#
#     Input:
#     [[0,0,0],
#      [0,1,0],
#      [0,0,0]]
#
#     Output:
#     [[0,0,0],
#      [0,1,0],
#      [0,0,0]]
#
#
# **Example 2:**
#
#
#
#     Input:
#     [[0,0,0],
#      [0,1,0],
#      [1,1,1]]
#
#     Output:
#     [[0,0,0],
#      [0,1,0],
#      [1,2,1]]
#
#
#
#
# **Note:**
#
#   1. The number of elements of the given matrix will not exceed 10,000.
#   2. There are at least one 0 in the given matrix.
#   3. The cells are adjacent in only four directions: up, down, left and right.
#
#
## @lc code=start
using LeetCode

function update_matrix(mat::Matrix{Int})
    m, n = size(mat)
    dist = fill(typemax(Int) >> 1, m, n)
    for I in CartesianIndices(mat)
        (mat[I] == 0) && (dist[I] = 0)
    end
    for i in 1:m, j in 1:n
        (i - 1 >= 1) && (dist[i, j] = min(dist[i, j], dist[i - 1, j] + 1))
        (j - 1 >= 1) && (dist[i, j] = min(dist[i, j], dist[i, j - 1] + 1))
    end
    for i in m:-1:1, j in n:-1:1
        (i + 1 <= m) && (dist[i, j] = min(dist[i, j], dist[i + 1, j] + 1))
        (j + 1 <= n) && (dist[i, j] = min(dist[i, j], dist[i, j + 1] + 1))
    end
    return dist
end
## @lc code=end


text:
############################################################# map to julia
m, d  = mj.mapmujoco(pm, pd)

@test mj.get(m, :nq) == 28
@test mj.get(m, :nv) == 27
@test mj.get(m, :nbody) == 14
@test mj.get(m, :opt, :timestep) == 0.002

mj.set(m, :nq, -28)
mj.set(m, :nv, -27)
mj.set(m, :nbody, -14)
mj.set(m, :opt, :timestep, -0.002)

@test mj.get(m, :nq) == -28
@test mj.get(m, :nv) == -27
@test mj.get(m, :nbody) == -14
@test mj.get(m, :opt, :timestep) == -0.002

mj.set(m, :nq, 28)
mj.set(m, :nv, 27)
mj.set(m, :nbody, 14)
mj.set(m, :opt, :timestep, 0.002)

@test mj.get(m, :nq) == 28
@test mj.get(m, :nv) == 27
@test mj.get(m, :nbody) == 14
@test mj.get(m, :opt, :timestep) == 0.002

############################################################# mujoco basics

qpos0 = copy(d.qpos)

mj_step(m, d)

@test d.qpos != qpos0
@test mj.get(d, :time) == 0.002

mj.set(m, :opt, :timestep, -0.002) # backwards step
mj_step(m, d)

@test isapprox(qpos0, d.qpos; atol = 1e-3)
@test mj.get(d, :time) == 0.0

mj.set(d, :time, -100.0)
@test mj.get(d, :time) == -100.0

mj_resetData(m, d)

@test qpos0 == d.qpos
@test mj.get(d, :time) == 0.00

mj_deleteModel(m)
mj_deleteData(d)


text:
#############################################################################
#############################################################################
#
# This file implements polynomial addition
#
#############################################################################
#############################################################################

"""
Add two polynomials.
"""
function +(p1::Polynomial, p2::Polynomial)::Polynomial
    p1, p2 = deepcopy(p1), deepcopy(p2)
    p3 = Polynomial()
    while !iszero(p1) && !iszero(p2)
        t1, t2 = leading(p1), leading(p2)
        if t1.degree == t2.degree
            push!(p3, pop!(p1)+pop!(p2))
        elseif t1.degree < t2.degree
            push!(p3,pop!(p2))
        else
            push!(p3,pop!(p1))
        end
    end
    while !iszero(p1)
        push!(p3,pop!(p1))
    end
    while !iszero(p2)
        push!(p3,pop!(p2))
    end
    return p3
end

"""
Add a polynomial and a term.
"""
+(p::Polynomial, t::Term) = p + Polynomial(t)
+(t::Term, p::Polynomial) = p + t

"""
Add a polynomial and an integer.
"""
+(p::Polynomial, n::Int) = p + Term(n,0)
+(n::Int, p::Polynomial) = p + Term(n,0)


text:
#
#  7 ---9--- 8 ---12-- 9
#
#  |         |         |
#  |         |         |
#  |         |         |
#  10   4    11   5    13
#  |         |         |
#  |         |         |
#  |         |         |
#
#  4 ---2--- 5 ---8--- 6
#
#  |         | \   3   |
#  |         |  \      |
#  |         |   \     |
#  3    1    4    6    7
#  |         |     \   |
#  |         |  2   \  |
#  |         |       \ |
#
#  1 ---1--- 2 ---5--- 3
#

struct GridTopologyMock <: GridTopology{2,2} end

function get_faces(g::GridTopologyMock,dimfrom::Integer,dimto::Integer)
  if dimfrom == 0
    if dimto == 0
      nvertex = 9
      return identity_table(Int,Int32,nvertex)
    elseif dimto == 1
      vertex_to_faces_data = [1,3,1,4,5,5,6,7,2,3,10,2,4,6,8,11,7,8,13,9,10,9,11,12,12,13]
      vertex_to_faces_ptrs = Int32[1,3,6,9,12,17,20,22,25,27]
      return Table(vertex_to_faces_data,vertex_to_faces_ptrs)
    elseif dimto == 2
      vertex_to_cells_data = [1,1,2,2,3,1,4,1,2,3,4,5,3,5,4,4,5,5]
      vertex_to_cells_ptrs = Int32[1,2,4,6,8,13,15,16,18,19]
      return Table(vertex_to_cells_data,vertex_to_cells_ptrs)
    else
      @unreachable
    end

  elseif dimfrom == 1
    if dimto == 0
      face_to_vertices_data = [1,2,4,5,1,4,2,5,2,3,3,5,3,6,6,5,7,8,4,7,5,8,8,9,6,9]
      face_to_vertices_ptrs = Int32[1,3,5,7,9,11,13,15,17,19,21,23,25,27]
      return Table(face_to_vertices_data, face_to_vertices_ptrs)
    elseif dimto == 1
      nedges = 13
      return identity_table(Int,Int32,nedges)
    elseif dimto == 2
      face_to_cells_data = [1,1,4,1,1,2,2,2,3,3,3,5,4,4,4,5,5,5]
      face_to_cells_ptrs = Int32[1,2,4,5,7,8,10,11,13,14,15,17,18,19]
      return Table(face_to_cells_data,face_to_cells_ptrs)
    else
      @unreachable
    end

  elseif dimfrom == 2
    if dimto == 0
      cell_to_vertices_data = [1,2,4,5,2,3,5,3,6,5,4,5,7,8,5,6,8,9]
      cell_to_vertices_ptrs = Int32[1,5,8,11,15,19]
      return Table(cell_to_vertices_data,cell_to_vertices_ptrs)
    elseif dimto == 1
      cell_to_faces_data = [1,2,3,4,5,4,6,7,6,8,2,9,10,11,8,12,11,13]
      cell_to_faces_ptrs = Int32[1,5,8,11,15,19]
      return Table(cell_to_faces_data, cell_to_faces_ptrs)
    elseif dimto == 2
      ncells = 5
      return identity_table(Int,Int32,ncells)
    else
      @unreachable
    end

  else
    @unreachable
  end
end

function get_polytopes(g::GridTopologyMock) where d
  [QUAD, TRI]
end

function get_cell_type(g::GridTopologyMock)
  quad = Int8(1)
  tri = Int8(2)
  [quad, tri, tri, quad, quad]
end

function get_vertex_coordinates(g::GridTopologyMock)
  Point{2,Float64}[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)]
end


text:
import TowardsNotebook

function julia_main()::Cint
  TowardsNotebook.notebook()
  return 0 # if things finished successfully
end

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Registers

using QuantumCircuits.QCircuits.QBase

export Bit, Qubit, Cbit, Register, ClassicalRegister, QuantumRegister,
       getid, setid

"The abstract bit"
abstract type Bit end

"The abstract register"
abstract type Register <: AbstractVector{Int} end

"Classical bit"
mutable struct Cbit <: Bit
    index::Union{Int, Nothing}
    regIndex::Int

    Cbit(idx::Integer) = new(nothing, idx)
end

"Quantum bit"
mutable struct Qubit <: Bit
    index::Union{Int, Nothing}
    regIndex::Int

    Qubit(idx::Integer) = new(nothing, idx)
end

function Base.show(io::IO, q::Qubit)
    if isnothing(q.index)
        print(io, "$(q.regIndex)")
    else
        print(io, "$(q.index)")
    end
end

function Base.:(==)(q1::Qubit, q2::Qubit)
    @assert !isnothing(q1.index) "Unable compare un assigned register."
    @assert !isnothing(q2.index) "Unable compare un assigned register."

    return q1.index == q2.index
end

function Base.hash(q::Qubit, h::UInt)
    @assert !isnothing(q.index) "Unable compare un assigned register."

    return hash(q.index, h)
end

"Get bit index"
function getid(b::Bit)
    @assert !isnothing(b.index) "The bit is not assigned to circuit."

    return b.index
end
"Get bit index"
function setid(b::Bit, id::Integer)
    @assert isnothing(b.index) "Unbale to assig bit to circuit twice."

    b.index = id
end

"Register macro"
macro register(name, bit)
    eval(quote
        struct $name <: Register
            name::Union{String, Nothing}
            bits::Vector{$bit}

            $name(n::Integer, name::Union{String, Nothing}) = new(name, [$bit(i) for i in 0:(n-1)])
        end
        $name(n::Integer) = $name(n, nothing)
    end)
end

# Classical Register
@register(ClassicalRegister, Cbit)
# Quantum Register
@register(QuantumRegister, Qubit)

"The size of register."
Base.length(reg::Register) = length(reg.bits)

"Size"
Base.size(reg::Register) = (length(reg),)

"Each index method"
Base.eachindex(reg::Register) = 0:(length(reg)-1)

"Get bit from register on given index."
function Base.getindex(reg::Register, idx::Integer)
    @assert idx < length(reg) "The index $idx is out of bound for register $(reg.name)."

    return reg.bits[idx + 1]
end

"Iterable method"
Base.iterate(reg::Register) = Base.iterate(reg.bits)
Base.iterate(reg::Register, state) = Base.iterate(reg.bits, state)

end  # module Registers


text:
#
# Protein - Water
#

using ComplexMixtures ; const CM = ComplexMixtures
using PDBTools
using Plots
ENV["GKSwstype"] = "nul"

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = readPDB("../../test/data/NAMD/structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = select(atoms,"protein")
solute = CM.Selection( solute_indexes, nmols=1 )

# The solvent is Water, which has 3 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent_indexes = select(atoms,"resname TIP3")
solvent = CM.Selection( solvent_indexes, natomspermol=3 )

# Initialize trajectroy data structure and open input stream
trajectory = CM.Trajectory("../../test/data/NAMD/trajectory.dcd",solute,solvent)

# Input options for the calcualtion
options = CM.Options(binstep=0.2)

# Run CM calculation, and get the resutls in the R structure
@time R = CM.mddf(trajectory,options)

CM.save(R,"example1.json")
CM.write(R,"example1.dat",solute,solvent)

plot(layout=(4,1))

x = [0,10]
y = [1, 1]

sp=1
plot!(ylabel="MDDF or RDF",subplot=sp)
plot!(R.d,R.mddf,subplot=sp,label="mddf")
plot!(R.d,R.rdf,subplot=sp,label="rdf")
plot!(x,y)
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(R.d,R.kb,subplot=sp,label="mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="atom contrib",subplot=sp)
for i in 1:solvent.natomspermol
  plot!(R.d,R.solvent_atom[:,i],subplot=sp,label="",linewidth=2)
end
plot!(legend=:topright,subplot=sp)

sp=4
plot!(ylabel="count",subplot=sp)
plot!(R.d,R.md_count,subplot=sp,label="R md")
plot!(R.d,R.md_count_random,subplot=sp,label="R rand")
plot!(legend=:topleft,subplot=sp)

plot!(size=(600,800))
savefig("./example2.pdf")


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here
set_line_width(cr, 30.0);
set_line_cap(cr, Cairo.CAIRO_LINE_CAP_BUTT); # default
move_to(cr, 64.0, 50.0); line_to(cr, 64.0, 200.0);
stroke(cr);
set_line_cap(cr, Cairo.CAIRO_LINE_CAP_ROUND);
move_to(cr, 128.0, 50.0); line_to(cr, 128.0, 200.0);
stroke(cr);
set_line_cap(cr, Cairo.CAIRO_LINE_CAP_SQUARE);
move_to(cr, 192.0, 50.0); line_to(cr, 192.0, 200.0);
stroke(cr);

# draw helping lines
set_source_rgb(cr, 1, 0.2, 0.2);
set_line_width(cr, 2.56);
move_to(cr, 64.0, 50.0);  line_to(cr, 64.0, 200.0);
move_to(cr, 128.0, 50.0); line_to(cr, 128.0, 200.0);
move_to(cr, 192.0, 50.0); line_to(cr, 192.0, 200.0);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_set_line_cap.png");


text:
@testset "model" begin
    N = 10
    sm = SellmeierModel(rand(N), rand(N))

    @test length(sm) == N
    @test sm(0.0) == 1.0
    @test sm(Inf) == 1.0 + sum(sm.str)
end  # @testset "model"


text:
using HPAT

function get_view_items(A::Vector{Tuple{Int64,Int64,Int64}}, days_in_sec_before_purchase, views_before_purchase, purchased_item)
    sort!(A, by=a->a[1], rev=true)
    out = Int64[]
    sale = false
    sale_time = 0
    count = 0
    for (t,i,s) in A
        # new sale happened, reset
        if s!=0
            count = 0
            if i==purchased_item
                sale = true
                sale_time = t
            else
                sale = false
            end
        # just a click by user
        else
            if sale==true && count<views_before_purchase && t>=sale_time-days_in_sec_before_purchase
                push!(out, i)
                count += 1
            end
        end
    end
    return out
end

@acc hpat function q03(days_in_sec_before_purchase, views_before_purchase, purchased_item, my_categories, file_name)

    web_clickstreams = DataSource(DataTable{:wcs_user_sk=Int64,:wcs_sales_sk=Int64,:wcs_item_sk=Int64,wcs_click_time_sk=Int64,wcs_click_date_sk=Int64}, HDF5, file_name)
    web_clickstreams[:tstamp_inSec] = web_clickstreams[:wcs_click_date_sk]*24*60*60 .+ web_clickstreams[:wcs_click_time_sk]
    item = DataSource(DataTable{:item_sk=Int64,:i_category_id=Int64}, HDF5, file_name)

    user_clicks = aggregate(web_clickstreams, :wcs_user_sk, collect(:tstamp_inSec),collect(:wcs_item_sk),collect(:wcs_sales_sk))
    view_items = map(a->get_view_items(a, days_in_sec_before_purchase, views_before_purchase, purchased_item), user_clicks)

    view_items = flatten(view_items)
    view_items_category = join(view_items, item, :view_items==:i_item_sk, :item_sk)
    view_items_category = view_items_category[:i_category_id in my_categories]
    item_counts = collect(countmap(view_items_category[:item_sk]))
    sort!(item_counts, rev=true)
    return item_counts[1:100]

end

println(q02(10*24*60*60, 5, 11, [3,4,5], "data.hdf5"))


text:
using SMC
using Base.Test

@testset "particles" begin include("particleset_test.jl")  end
@testset "hmm"       begin include("hmm_test.jl")          end
@testset "kalman"    begin include("kalman_test.jl")       end
@testset "pf,ps"     begin include("particlef+s_test.jl")  end
@testset "ps,subq"   begin include("particlesubq_test.jl") end


text:
# Feynman-Kleinert-1986-PRA
# Replot fits from Feynman-Kleinert Table II
# - basically checking whether the Wtilde equation is reporting correctly,
# and whether (as described in the caption) only Wtilde for g=0.1976 has double well structure

macro doublewell(g)
    return :(x -> -0.5*x^2 + 0.25*g*x^4 + 0.25/g)
end

@printf("\n\nFeynman-Kleinert Table II reads:\n g      E0=Wtilde\n 0.1976 0.650 \n 0.4    0.549 \n 4.0    0.598 \n 40.0   1.409\n")
β=3000
# Much higher than this, and results collapse to Inf. But essentially agrees with Table II now!

@printf("\nI calculate, β=%g: \n",β)
# Wtilde(x0,a2,Ω)
xrange=-4:0.1:4
# These values, from Feynman and Kleinhert, table II
g=0.1976
myW=Wtilde(1.943,0.397,1.255,@doublewell(g),β)
@printf("%f %f\n",g,myW)
@test myW ≈ 0.650 atol=0.02

g=0.4
myW=Wtilde(0.0,1.030,0.486,@doublewell(g),β)
@printf("%f %f\n",g,myW)
@test  myW ≈ 0.549 atol=0.01

g=4.0
myW=Wtilde(0.0,0.3059,1.634,@doublewell(g),β)
@printf("%f %f\n",g,myW)
@test myW ≈ 0.598 atol=0.01

g=40.0
myW=Wtilde(0.0,0.1306,3.829,@doublewell(g),β)
@printf("%f %f\n",g,myW)
@test myW ≈ 1.409 atol=0.01


text:
#######
# Here we store the parameters corresponding to the SQUID-like geometry (`mask = true`)
######

####### Flux to B conversion
Φ0 = 2067.8338 # T*nm^2
function Bprange(p, Φrange)
    Bperp = collect(Φrange./(p.Ln * p.W)) #conversion from flux to Bfield (size dependent)
    Bcriteria(Bperp, p.scale)
    return Bperp
end

"We need to make sure that we are not in the Hofstadter regime, for the scaled and
non-scaled system"
function Bcriteria(B::Array, scale; threshold = 50.)
    ħc = ustrip(u"eV*nm",ħ*c_0)
    a0 = 0.246    #nm
    Bmax = findmax(abs.(B))[1]
    l_B = ħc/Bmax
    l_scale = scale * a0
    println(l_B/ l_scale)
    println(l_B/a0)
    l_B/ l_scale > threshold ? nothing : @warn("lB ≈ a0*scale, decrease scaling")
    l_B/ a0 > threshold ? nothing : @warn("lB ≈ a0, increase BLG area (we are in t
        he Hofstadter regime)")
end

#################
# Params
#################
p = Params(scale = 10, Δ=.0, Ln=200, Ls=0, W = 320, U = 0.00005, μSL = .8, μSR=.8,
    μN = 0*-0.0034, V_perp =  0., V_par = 1e-6, B_par = 0, B_perp = 0., Vu = 0.,
    edges = :armchair_zigzag, smooth = false, dir = :x, U_dimer = 0.0, α = 0.,
    smoothμ = false, localselfenergy = false, analyticalselfenergy = false, normal=true);

p_inv = reconstruct(p, λ = 0.02,
    scale = 10, Ln = 70, W = 74,
    smoothμ = false, mask = true, Δx_mask = 7, Δy_mask = 5,
    analyticalselfenergy = false, normal = true)

p_inv_sc = reconstruct(p_inv,
    localselfenergy = true, normal = false, Δ = 0.1)

p_inv_sc_smooth = reconstruct(p_inv_sc,
    smoothμ = true, Δxμ = 7, λμ = 1, λμbarrier = 1, Δμbarrier = -3)

    ############################################################
p_norm_sc_smooth = reconstruct(p_inv_sc_smooth,
    Δx_mask = 14,  μN = -0.34)

##############################
p_inv_2 = reconstruct(p, λ = 0.02,
    scale = 10, Ln = 50, W = 50,
    smoothμ = false, mask = true, Δx_mask = 7, Δy_mask = 5,
    analyticalselfenergy = false, normal = true)

p_inv_sc_2 = reconstruct(p_inv_2,
    localselfenergy = true, normal = false, Δ = 0.1)

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

p_inv_3 = reconstruct(p, λ = 0.02,
    scale = 20, Ln = 50, W = 80,
    smoothμ = false, mask = true, Δx_mask = 5, Δy_mask = 8,
    analyticalselfenergy = false, normal = true)

p_inv_sc_3 = reconstruct(p_inv_3,
    localselfenergy = true, normal = false, Δ = 0.1)

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

p_nodeg = reconstruct(p_inv_sc_smooth, V_par = 1e-4, α = 1e-4, Ex = 1e-4, Ey = 1e-4,
    Ez = 1e-4)


text:
module Data
    function parse_data_list(source::Array)
        result = Dict(
            "type"=> nothing,
            "data"=> []
        )
        list_type = Dict(
            "Float64"=> 0,
            "Int"=> 0,
            "String"=> 0
        )

        for item in source
            item = string(item)
            if !isnull(tryparse(Float64, item))
                if contains(item, ".")
                    list_type["Float64"] += 1
                else
                    list_type["Int"] += 1
                end
            else
                list_type["String"] += 1
            end
        end
        if list_type["String"] > 0
            result["type"] = String
            result["data"] = [string(item) for item in source]
        else
            if list_type["Float64"] > 0
                result["type"] = Float64
            else
                result["type"] = Int
            end

            result["data"] = [parse(result["type"], string(item)) for item in source]
        end

        return result
    end
end


text:
#!/bin/julia

#=
NAME: download_julia_packages - download all registered Julia packages in one go.

SYNOPSYS
    download_julia_packages [--dir <path>] [--archive <option>] [--julia-dir <path>]
                            [--registry <path to 'Registry.toml'>]

DESCRIPTION
	The script downloads the latest Julia nightly from http://julialang.org and installs it in
	/opt if --dir is not used. It also generates a sylink named 'julia_nightly' in /usr/local/bin.
OPTIONS
	--dir <install path>
        Downloads all Julia packages to the path specified. If the path does not exist,
        it is created.
	--archive
		Additionally creates a '.tar.gz' archive with all packages.
    --julia-dir <path>
        The path to the '.julia' folder; this is used for accessing package information.
    --registry <path>
        Custom path for the 'Registry.toml' file.

AUTHOR

Written by Corneliu Cofaru, (c) 2018
=#

# TODO(Corneliu): Add update option (i.e. incremental download + update existing)
# TODO(Corneliu): Integrate with LibGit2 to perform operations.
# TODO(Corneliu): Get/update Registry.toml by default without needing a Julia
# installation path or path to registry toml; the options can be used as
# overrides only.

using Pkg
installed = Pkg.installed()
!("ArgParse" in keys(installed)) &&  Pkg.add("ArgParse")
using ArgParse

# Function that parses input arguments
function get_commandline_arguments(args::Vector{String})
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dir"
            help = "download path"
            default = "./"
        "--archive"
            help = "archive the downloaded repositories into a '.tar.gz' archive"
            action = :store_true
        "--julia-dir"
            help = "path to the '.julia' folder"
            default = abspath(expanduser("~/.julia"))
        "--registry"
            help = "path to the 'Registry.toml' file"
            default = abspath(expanduser("~/.julia/registries/General.toml"))
    end
    return parse_args(args,s)
end

function main(ARGS::Vector{String})

	# Parse command line arguments
	args = get_commandline_arguments(ARGS)

	DIR = args["dir"]
	ARCHIVE = args["archive"]
    JULIA_DIR = args["julia-dir"]
    REGISTRY = args["registry"]
    ARFILE = "General.tar.gz"

    # Create installation path if it does not exist
    directory = abspath(DIR)
    !isdir(directory) && mkpath(directory)

    # Check the Julia directory for registries
    pkg_reg_path = abspath(joinpath(JULIA_DIR, "registries"))
    !isdir(pkg_reg_path) && @error "$pkg_reg_path does not exist."

    # Get the Registry.toml path
    registry_path = joinpath(pkg_reg_path, "_ex")
    !isdir(registry_path) && mkpath(registry_path)
    try
        run(`tar zxvf $(joinpath(pkg_reg_path, ARFILE)) -C $registry_path`)
    catch e
        error("Could not extract Registry.toml: $e")
    end

    # Parse file and extract information in a Vector of tuples (package_name, path/to/Package.toml)
    @info "Parsing $registry_path..."

    # Function that parses the line in the Registry.toml file
    function _parse_registry_line(line)
        _, _tmp1 = split(line, "{");
        _tmp1=replace(_tmp1, r"(\"|})"=>"");
        _name, _path = split(_tmp1, ",")
        _, name = strip.(split(_name, "="))
        _, path = strip.(split(_path, "="))
        return name, path
    end

    function _git_clone(repo, path)
        try
            run(pipeline(`git clone $repo $path`, stdout=devnull))
        catch
        end
    end

    PkgInfos = Vector{Tuple{String, String}}()
    open(joinpath(registry_path, "Registry.toml"), "r") do fid_registry
        # Inside Registry.toml
        for line in eachline(fid_registry)
            if all(occursin(symbol,line) for symbol in ["{", "}", "="])
                #00701ae9-d1dc-5365-b64a-a3a3ebf5695e = { name = "BioAlignments", path = "B/BioAlignments" }
                name, regdir = _parse_registry_line(line)
                open(joinpath(registry_path, regdir, "Package.toml"), "r") do fid_package
                    # Inside Package.toml
                    for _line in eachline(fid_package)
                        key, val = strip.(split(_line, "="))
                        if key == "repo"
                            push!(PkgInfos, (name, replace(val, "\""=>"")))
                        end
                    end
                end
            end
        end
    end

    @info "Cleaning up $registry_path"
    run(`rm -rf $registry_path`)

    @info "Cleaning up $directory ..."
    directories = readdir(directory)
    if !isempty(directories)
        for _dir in directories
            @info "removing $_dir"
            Base.Filesystem.rm(joinpath(directory,_dir), recursive=true, force=true)
        end
    end

    @info "Cloning repositories..."
    CONCURRENT = 16
    asyncmap(PkgInfos, ntasks=CONCURRENT) do x
        _git_clone(x[2], joinpath(DIR, x[1]))
    end

    if ARCHIVE
        @info "Compressing repositories"
        run(pipeline(`tar -I pigz -cf ./all_julia_repositories.tgz $(joinpath(directory, "/"))`,
                     stdout=devnull))
    end
end

#########################
# Run the main function #
#########################
main(ARGS)


text:
latexinline(args...;kwargs...) = latexify(args...;kwargs...,env=:inline)

function _latexinline(x; kwargs...)
    latexstr = latexstring( latexify(x; kwargs...,env=:raw) )
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end


text:
mutable struct QGateBlock
	gates::Vector{Union{QGate, QGateBlock}}
	qubits::Vector{Vector{Int}}
	QGateBlock() = new(Vector{T where {T<:Union{QGate,QGateBlock}}}(undef, 0), Vector{Vector{Int}}(undef, 0))
	QGateBlock(gates::Vector{T} where {T<:Union{QGate,QGateBlock}}, qubits::Vector{Vector{Int}}) = new(gates, qubits)
end # struct

const Operator = Union{QGate, QGateBlock}

gates(block::QGateBlock) = block.gates
qubits(block::QGateBlock) = block.qubits
size(block::QGateBlock) = size(block.gates)
iterate(block::QGateBlock, state::Int=1) = iterate(block.gates, state)
copy(block::QGateBlock) = QGateBlock(copy.(gates(block)), copy(qubits(block)))
length(block::QGateBlock) = length(gates(block))
function ==(block1::QGateBlock, block2::QGateBlock)
	length(block1)!=length(block2) && (return false)
	for i =1:length(block1)
		if (block1.gates[i]!= block2.gates[i]) || (block1.qubits[i]!= block2.qubits[i])
			return false
		end
	end
	return true
end
function add!(block::QGateBlock, gate::Operator, pos::Vector{Int})
	if isa(gate, VarGate)
		gate = copy(gate)
	end
	push!(gates(block), gate)
	push!(qubits(block), pos)
	return block
end

add!(block::QGateBlock, gate::Operator, qubits::Int...) = add!(block, gate, [qubits...])

const GatePosTuple = Tuple{T, Vector{Int}} where {T <: Operator}

function add!(block::QGateBlock, tuples::GatePosTuple...)
	for tuple ∈ tuples
		add!(block, tuple[1], tuple[2])
	end
	return block
end

addCopy!(block::QGateBlock, gate::Operator, pos::Vector{Int}) = add!(block, copy(gate), pos)
function addCopy!(block::QGateBlock, tuples::GatePosTuple...)
	for tuple ∈ tuples
		add!(block, copy(tuple[1]), copy(tuple[2]))
	end
end

function flatten(block::QGateBlock, pos::Vector{Int}, flatttened = nothing)
	flatttened===nothing && (flatttened = QGateBlock())
	for i = 1:length(block)
		gateOrBlock = gates(block)[i]
		gateOrBlockPos = qubits(block)[i]
		if isa(gateOrBlock, QGateBlock)
			flatten(gateOrBlock, [pos[i] for i ∈ gateOrBlockPos],flatttened)
		else
			add!(flatttened, gateOrBlock, [pos[i] for i ∈ gateOrBlockPos])
		end
	end
	return flatttened
end

function extractParams(block::QGateBlock, params= nothing)
	params===nothing && (params= Real[])
	for gateOrBlock ∈ block
		if isa(gateOrBlock, QGateBlock)
			extractParams(gateOrBlock, params)
		elseif isa(gateOrBlock, VarGate)
			push!(params, param(gateOrBlock)...)
		end
	end
	return params
end

function insertParams!(block::QGateBlock, params::Vector{T} where {T<:Real}, i::Int = 1)
	for gateOrBlock ∈ block
		if isa(gateOrBlock, QGateBlock)
			i = insertParams!(gateOrBlock, params, i)
		elseif isa(gateOrBlock, VarGate)
			n = paramCount(gateOrBlock)
			setParam!(gateOrBlock, params[i:i+n-1])
			i+=n
		end
	end
	return i
end

function inverse(block::QGateBlock)
	inversed = QGateBlock()
	reverseOrderGates = reverse(block.gates)
	reverseOrderPos = reverse(block.qubits)
	for i = 1:length(block)
		add!(inversed, inverse(reverseOrderGates[i]), reverseOrderPos[i])
	end
	return inversed
end

function randomQGateBlock(size::Int, depth::Int=rand(1:10))
	block = QGateBlock()
	for i=1:depth
		type = rand(1:16)
		if type <=12
			subtype = rand(1:12)
			if subtype<=8
				numqubits=1
			elseif subtype<=10
				numqubits=2
			else
				numqubits = 3
			end
			add!(block,randomConstGate(numqubits), Random.randperm(size)[1:numqubits])
		else
			add!(block, randomVarGate(), rand(1:size))
		end
	end
	return block
end

function show(io::IO, block::QGateBlock)
	for i = 1:length(block)
		if isa(block.gates[i], QGateBlock)
			printstyled(io, "BLOCK\n"; color = :red)
			print(io, block.qubits[i],"\n")
			print(io, block.gates[i])
		else
			print(io, block.gates[i]," , ", block.qubits[i], "\n")
		end
	end
end


text:
"""
This type is a container type that we hold a reference to in order to make sure
that the garbage collector doesn't free up the memory that backs up our ring
buffers.
"""
type LockFreeRingBufferData{T}
    buf::Vector{T}
    nread::UInt
    nwritten::UInt
    cond::AsyncCondition
end

# keep all our references here
const ringbufs = Set{LockFreeRingBufferData}()

"""
A LockFreeRingBuffer implements a single-producer single-consumer lock-free ring
buffer. It doesn't provide any of the queuing behavior.

This also provides a mechanism to read and write from a ring buffer without
interacting with Julia's memory management, which makes this type useful for
cases where data is produced and/or consumed from a separate thread's context,
such as when interfacing with an audio library.

Note that this is an immutable object with pointers to heap-allocated data.
This way we can pass it around as a stack-allocated object without hitting
the garbage collector infrastructure.
"""
immutable LockFreeRingBuffer{T}
    size::UInt
    buf::Ptr{T}
    nread::Ptr{UInt}
    nwritten::Ptr{UInt}
    container::Ptr{LockFreeRingBufferData{T}}
    cond::Ptr{AsyncCondition}
    condhandle::Ptr{Void}
end

function LockFreeRingBuffer(T, size)
    size = nextpow2(size)
    cond = AsyncCondition()
    container = LockFreeRingBufferData(Array(T, size), UInt(0), UInt(0), cond)
    push!(ringbufs, container)
    rootaddr = pointer_from_objref(container)
    bufptr = Ptr{T}(pointer(container.buf))
    nreadptr = Ptr{UInt}(rootaddr + fieldnameoffset(LockFreeRingBufferData{T}, :nread))
    nwrittenptr = Ptr{UInt}(rootaddr + fieldnameoffset(LockFreeRingBufferData{T}, :nwritten))
    containerptr = Ptr{LockFreeRingBufferData{T}}(rootaddr)

    LockFreeRingBuffer(UInt(size), bufptr, nreadptr, nwrittenptr, containerptr,
            Ptr{AsyncCondition}(pointer_from_objref(cond)), cond.handle)
end

function nreadable(buf::LockFreeRingBuffer)
    nread = unsafe_load(buf.nread)
    nwritten = unsafe_load(buf.nwritten)

    nwritten - nread
end

nwritable(buf::LockFreeRingBuffer) = buf.size - nreadable(buf)

function write{T}(buf::LockFreeRingBuffer{T}, data::Ptr{T}, n::UInt)
    n = min(n, nwritable(buf))
    writepos = unsafe_load(buf.nwritten) + 1
    sizemask = buf.size - 1
    for i in 1:n
        # note unsafe_store! wants 1-indexed data, but the masking assumes 0-indexed
        unsafe_store!(buf.buf, unsafe_load(data, i), (writepos + i - 1) & sizemask + 1)
    end

    unsafe_store!(buf.nwritten, unsafe_load(buf.nwritten) + n)
    wakewaiters(buf)

    n
end

# convert the sample count to a UInt if necessary
write{T}(buf::LockFreeRingBuffer{T}, data::Ptr{T}, n::Real) = write(buf, data, UInt(n))

# also handle Vectors, possibly detecting the length
write{T}(buf::LockFreeRingBuffer{T}, data::Array{T}, n=length(data)) = write(buf, pointer(data), UInt(n))

function read!{T}(buf::LockFreeRingBuffer{T}, data::Ptr{T}, n::UInt)
    n = min(n, nreadable(buf))
    readpos = unsafe_load(buf.nread) + 1
    sizemask = buf.size - 1
    for i in 1:n
        # note unsafe_load wants 1-indexed data, but the masking assumes 0-indexed
        unsafe_store!(data, unsafe_load(buf.buf, (readpos + i - 1) & sizemask + 1), i)
    end

    unsafe_store!(buf.nread, unsafe_load(buf.nread) + n)
    wakewaiters(buf)

    n
end

# convert the sample count to a UInt if necessary
read!{T}(buf::LockFreeRingBuffer{T}, data::Ptr{T}, n::Real) = read!(buf, data, UInt(n))

# also handle Vectors, possibly detecting the length
read!{T}(buf::LockFreeRingBuffer{T}, data::Array{T}, n=length(data)) = read!(buf, pointer(data), UInt(n))

# this just removes the container from our ringbufs list, so the GC can do its
# thing. After this all the pointers should be considered invalid
function Base.close(buf::LockFreeRingBuffer)
    pop!(ringbufs, unsafe_pointer_to_objref(buf.container))
end

Base.length(buf::LockFreeRingBuffer) = buf.size

wait(buf::LockFreeRingBuffer) = wait(unsafe_load(buf.cond))
wakewaiters(buf::LockFreeRingBuffer) = ccall(:uv_async_send, Cint, (Ptr{Void}, ), buf.condhandle)

"Gives the memory offset of the given field, given as a symbol"
fieldnameoffset(T, fname) = fieldoffset(T, findfirst(fieldnames(T), fname))


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

const d_lng = SVector(0, -2, 0, 0, 0, 0, -2, 0, 0, -2, -2, -2, 0, 2, 0, 2, 0, 0, -2, 0,
                      2, 0, 0, -2, 0, -2, 0, 0, 2,  -2, 0, -2, 0, 0, 2, 2, 0, -2, 0, 2,
                      2, -2, -2, 2, 2, 0, -2, -2, 0, -2, -2, 0, -1, -2, 1, 0, 0, -1, 0,
                      0,  2, 0, 2)
const M_lng = SVector(0, 0, 0, 0, 1, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                      0, 0, 0, 0, 0, 0, 2, 0, 2, 1, 0, -1, 0, 0, 0, 1, 1, -1, 0, 0, 0,
                      0, 0, 0, -1, -1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, -1, 1, -1, -1, 0,
                      -1)
const Mprime_lng = SVector(0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, -1, 0, 1, -1, -1, 1, 2,
                           -2, 0, 2, 2, 1, 0, 0, -1, 0, -1,  0, 0, 1, 0, 2, -1, 1, 0,
                           1, 0, 0, 1, 2, 1, -2, 0, 1, 0, 0, 2, 2, 0, 1, 1, 0, 0, 1,
                           -2, 1, 1, 1, -1, 3, 0)
const F_lng = SVector(0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 2, 2, 0, 0, 2, 0, 2, 0, 2, 2, 2,
                      0, 2, 2, 2, 2, 0, 0, 2, 0, 0,  0, -2, 2, 2, 2, 0, 2, 2, 0, 2, 2,
                      0, 0, 0, 2, 0, 2, 0, 2, -2, 0, 0, 0, 2, 2, 0, 0, 2, 2, 2, 2)
const ω_lng = SVector(1, 2, 2, 2, 0, 0, 2, 1, 2, 2, 0, 1, 2, 0, 1, 2, 1, 1, 0, 1, 2, 2,
                      0, 2, 0, 0, 1, 0, 1, 2, 1, 1, 1, 0, 1, 2, 2, 0, 2, 1, 0, 2, 1, 1,
                      1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 2, 0, 0, 2, 2, 2, 2)
const sin_lng = SVector(-171996, -13187, -2274, 2062, 1426, 712, -517, -386, -301, 217,
                        -158, 129, 123, 63, 63, -59, -58, -51, 48, 46, -38, -31, 29,
                        29, 26, -22, 21, 17, 16, -16, -15, -13, -12, 11, -10, -8, 7,
                        -7, -7, -7, 6, 6, 6, -6, -6, 5, -5, -5, -5, 4, 4, 4, -4, -4,
                        -4, 3, -3, -3, -3, -3, -3, -3, -3)
const sdelt = SVector(-174.2, -1.6, -0.2, 0.2, -3.4, 0.1, 1.2, -0.4, 0, -0.5, 0, 0.1,
                      0.0, 0.0, 0.1, 0.0, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, -0.1, 0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
const cos_lng = SVector(92025, 5736, 977, -895, 54, -7, 224, 200, 129, -95, 0, -70,
                        -53, 0, -33, 26, 32, 27, 0, -24, 16, 13, 0, -12, 0, 0, -10, 0,
                        -8, 7, 9, 7, 6, 0, 5, 3, -3, 0, 3, 3, 0, -3, -3, 3, 3, 0, 3, 3,
                        3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
const cdelt = SVector(8.9, -3.1, -0.5, 0.5, -0.1, 0.0, -0.6, 0.0, -0.1, 0.3, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

function nutate(jd::AbstractFloat)
    # Number of Julian centuries since 2000-01-01T12:00:00
    t = (jd - J2000) / JULIANCENTURY
    # Mean elongation of the Moon
    d = deg2rad(mod(@evalpoly(t, 297.85036, 445267.111480, -0.0019142, inv(189474)), 360))
    # Sun's mean anomaly
    M = deg2rad(mod(@evalpoly(t, 357.52772, 35999.050340, -0.0001603, -inv(3e5)), 360))
    # Moon's mean anomaly
    Mprime = deg2rad(mod(@evalpoly(t, 134.96298, 477198.867398, 0.0086972, inv(5.625e4)), 360))
    # Moon's argument of latitude
    F = deg2rad(mod(@evalpoly(t, 93.27191, 483202.017538, -0.0036825, -inv(3.27270e5)), 360))
    # Longitude of the ascending node of the Moon's mean orbit on the ecliptic,
    # measured from the mean equinox of the date
    ω = deg2rad(mod(@evalpoly(t, 125.04452, -1934.136261, 0.0020708, inv(4.5e5)), 360))
    arg = d_lng * d + M_lng * M + Mprime_lng * Mprime + F_lng * F + ω_lng * ω
    sc_arg = sincos.(arg)
    s = getindex.(sc_arg, 1)
    c = getindex.(sc_arg, 2)
    long  = dot((sdelt .* t .+ sin_lng), s) / 10000
    obliq = dot((cdelt .* t .+ cos_lng), c) / 10000
    return long, obliq
end

"""
    nutate(jd) -> long, obliq

### Purpose ###

Return the nutation in longitude and obliquity for a given Julian date.

### Arguments ###

* `jd`: Julian ephemeris date, it can be either a scalar or a vector

### Output ###

The 2-tuple `(long, obliq)`, where

* `long`: the nutation in longitude
* `obl`: the nutation in latitude

If `jd` is an array, `long` and `obl` are arrays of the same length.

### Method ###

Uses the formula in Chapter 22 of ``Astronomical Algorithms'' by Jean Meeus
(1998, 2nd ed.) which is based on the 1980 IAU Theory of Nutation and includes
all terms larger than 0.0003".

### Example ###

(1) Find the nutation in longitude and obliquity 1987 on Apr 10 at Oh.  This is
example 22.a from Meeus

```jldoctest
julia> using AstroLib

julia> jd = jdcnv(1987, 4, 10);

julia> nutate(jd)
(-3.7879310771104917, 9.442520698644492)
```

(2) Plot the daily nutation in longitude and obliquity during the 21st century.
Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.

```julia
using PyPlot
years = DateTime(2000):DateTime(2100);
long, obl = nutate(jdcnv.(years));
plot(years, long); plot(years, obl)
```

You can see both the dominant large scale period of nutation, of 18.6 years, and
smaller oscillations with shorter periods.

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
nutate(jd::Real) = nutate(float(jd))

function nutate(jd::AbstractArray{J}) where {J<:Real}
    typej = float(J)
    long = similar(jd, typej)
    obliq = similar(jd, typej)
    for i in eachindex(jd)
        long[i], obliq[i] = nutate(jd[i])
    end
    return long, obliq
end


text:
module Plotter

include("pitch.jl")
include("passmap.jl")

using .StatsbombPitch
using .Passmap

end ##end of module


text:
using CMPUtils
recodeaudio("wav/W")

#recodeaudio("wav/W")
#recodeaudio("wav/1")
#recodeaudio("wav/2")
#recodeaudio("wav/34")
#recodeaudio("wav/R")

# Do the above one at a time
# and restart VSCode in between
# And don't put anything else in the terminal, including the comments!
# Only run the two lines "using CMPUtils" and "recodeaudio(...)"
# Otherwise the terminal input gets mangled

# Features:  5 = Recommended
#           13 = Maximum (see CMPUtils)

# TODO: try this with 13 and see how the classifiers like it


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# nominal model
include(joinpath(module_dir(), "src/dynamics/hopper_2D/model.jl"))
model_nom = get_model("hopper_2D")

model_sim = Hopper2D(Dimensions(nq, nu, nw, nc, nb),
			   1.5 * mb, ml, 1.5 * Jb, Jl,
			   μ_world, μ_joint, g,
			   BaseMethods(), DynamicsMethods(), ContactMethods(),
			   ResidualMethods(), ResidualMethods(),
			   SparseStructure(spzeros(0, 0), spzeros(0, 0)),
			   SVector{4}(zeros(4)),
			   environment_2D_flat())

expr_base = generate_base_expressions(model, M_analytical=true)
instantiate_base!(model_sim.base, expr_base)
expr_dyn = generate_dynamics_expressions(model_sim)
instantiate_dynamics!(model_sim.dyn, expr_dyn)
expr_res, rz_sp, rθ_sp = generate_residual_expressions(model_sim)
instantiate_residual!(model_sim.res, expr_res)
model_sim.spa.rz_sp = rz_sp
model_sim.spa.rθ_sp = rθ_sp

nq = model_nom.dim.q
nu = model_nom.dim.u
nc = model_nom.dim.c
nb = model_nom.dim.b

# time
H = 92#traj.H
h = 0.01#traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 5000

# Select the initial speed
v0 = 0.2
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
include(joinpath(module_dir(), "src/controller/raibert_policy.jl"))
p = raibert_policy(model_nom, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h)

off0 = SVector{nq,T}([0.0, 0.5, 0.0, 0.0])
off1 = SVector{nq,T}([v0*h_sim, 0.5, 0.0, 0.0])
q_ref = SVector{nq,T}([0.0, 0.5, 0.0, 0.5])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

sim = simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	uL = [-0.25, -0.3], uU = [0.25, 0.3],
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = SimulatorOptions(warmstart = true)
    )

@time status = simulate!(sim)

anim = visualize_robot!(vis, sim.model, sim.traj, sample=20)

plot(hcat(sim.traj.u...)[1:2, 1:10:end]',
	label = ["torque" "force"], xlabel = "time step", ylabel = "control")


text:
import MittagLefflerFunctions: rΓ, MLQuad
using Printf
import SpecialFunctions: gamma

Γ(x) = gamma(x)

α = 0.7
β = 1.0
tol = 1e-12
Eref = MLQuad(α, β, 14, :hyperbola)
xvals = collect(5.0:10.0:55.0)

function asymp_coeffs(α, β, N)
    s = 0.0
    σ = zeros(N)
    τ = zeros(N)
    sgn = 1.0
    s = 0.0
    for n = 1:N
        if n*α < β
            σ[n] = 1.0
            τ[n] = 1 / Γ(β-n*α)
        else
            σ[n] = -sin(π*(n*α-β))
            τ[n] = Γ(1+n*α-β) / π
        end
    end
    return σ, τ
end

function table_row(α, β, x, tol)
    Nopt = x^(1/α) / α
    N = ceil(Int64, Nopt)
    σ, τ = asymp_coeffs(α, β, N)
    z = -x
    powz = 1.0
    m = 0
    next = NaN
    last = NaN
    s = 0.0
    for n = 1:N-1
        powz *= z
        s -= σ[n] * τ[n] / powz
        last = τ[n] / abs(powz)
        m = n + 1
        if last < tol
            break
        end
    end
    powz *= z
    next = τ[m] / abs(powz)
    return s, m, Nopt, last, next
end

@printf("α = %0.2f, β = %0.2f, tol = %0.2g\n\n", α, β, tol)
@printf("%6s  %5s  %8s  %10s  %10s  %10s\n\n",
        "x", "m", "Nx", "error", "m-1", "m")
for x in xvals
    s, m, Nopt, last, next = table_row(α, β, x, tol)
    err = Eref(-x) - s
    @printf("%6.2g& %5d& %8.1f& %10.2e& %10.2e& %10.2e\\\\\n",
            x, m, Nopt, err, last, next)
end


text:
"""
    FENodeSetModule

Module for the finite element node set.
"""
module FENodeSetModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.count

"""
    Finite element node set type.

xyz =
    Array of node locations.
    Array of coordinates, the number of rows corresponds to the number of
    nodes in the set and the columns corresponds to the space dimensions.
    The location of node j is given by x[j,:].

"""
mutable struct FENodeSet
  xyz::Array{FFlt, 2}

  function FENodeSet(xyz::FFltMat)
    self = new(deepcopy(xyz)); # Need to make a COPY of the input array!
    return self
  end
end

"""
    spacedim(self::FENodeSet)

Number of dimensions of the space in which the node lives, 1, 2, or 3.
"""
spacedim(self::FENodeSet) =size(self.xyz,2)

"""
    xyz3(self::FENodeSet)

Get the  3-D coordinate that define the location  of the node.
Even if the nodes  were specified in  lower dimension (1-D, 2-D)
this function returns  a 3-D coordinate  by padding with zeros.
"""
function xyz3(self::FENodeSet)
    if (size(self.xyz,2)==1)
        val = [self.xyz zeros(size(self.xyz,1),2)];
    elseif (size(self.xyz,2)==2)
        val = [self.xyz zeros(size(self.xyz,1),1)];
    else
        val = deepcopy(self.xyz);
    end
end

"""
    count(self::FENodeSet)

Get the number of finite element nodes in the node set.
"""
function count(self::FENodeSet)
    return size(self.xyz,1)::FInt
end

end


text:
using JLD
using DiffPrecTest
using LaTeXStrings

include("../qq_plot.jl")

pArr = [100, 200]
elemArr = [(5,5), (8, 7), (50, 25)]
methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"]

## qq-plots

fig_num = 0

fig = figure(figsize=(7.08661, 4.72441), dpi=1000)
for ip=1:2
  for iElem=1:3
    global fig_num

    res = load("results/sim1_res_$(ip)_$(iElem).jld", "results")
    fig_num += 1

    qq_my     = [res[1, i].p / res[1, i].std for i=1:1000]
    qq_oracle = [res[6, i].p / res[6, i].std for i=1:1000]

    ax = subplot(2, 3, fig_num)
    qqplot(qq_oracle, color="grey")
    qqplot(qq_my)
    title(
     latexstring("p = $(pArr[ip]), \\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"),
     size="small"
    )
    ax[:tick_params]("both", labelsize="xx-small", length=2, pad=2)
  end
end
tight_layout()
savefig("qq_plot_sim1.pdf")
close(fig)

## hist-plots

fig_num = 0

x = range(-3.290, stop=3.290, length=100)
y = pdf.(Normal(), x)

fig = figure(figsize=(7.08661, 4.72441), dpi=1000)
for ip=1:2
  for iElem=1:3
    global fig_num

    res = load("results/sim1_res_$(ip)_$(iElem).jld", "results")
    fig_num += 1

    ax = subplot(2, 3, fig_num)
    plt[:hist]([res[1, i].p / res[1, i].std for i=1:1000], 100, density=true)
    plot(x,y)
    xlim(-3.290,3.290)
    title(
     latexstring("p = $(pArr[ip]), \\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0"),
     size="small"
    )
    ax[:tick_params]("both", labelsize="xx-small", length=2, pad=2)
  end
end
tight_layout()
savefig("hist_plot_sim1.pdf")
close(fig)

## irina plots


text:
function interpolate(c::EllipticalArc, p0, t)
    ϕ, rx, ry = c.angle, c.r1, c.r2
    θ = c.a1 + (c.a2 - c.a1) * t
    m = Mat2(cos(ϕ), sin(ϕ), -sin(ϕ), cos(ϕ))
    m * Point(rx * cos(θ), ry * sin(θ)) + c.c
end

function interpolate(c::CurveTo, p0, t)
    p1, p2, p3 = c.c1, c.c2, c.p
    (1 - t)^3 * p0 + 3(t - 2t^2 + t^3) * p1 + 3(t^2 -t^3) * p2 + t^3 * p3
end

discretize!(v::Vector{<:AbstractPoint}, c::Union{MoveTo, LineTo}) = push!(v, c.p)
function discretize!(v::Vector{<:AbstractPoint}, c::ClosePath)
    if v[end] !== v[begin]
        push!(v, v[begin])
    end
end

function discretize!(v::Vector{<:AbstractPoint}, c::Union{EllipticalArc, CurveTo})
    N0 = length(v)
    p0 = v[end]
    N = 100
    resize!(v, N0 + N)
    dt = 1.0/N
    for (i, t) in enumerate(dt:dt:1.0)
        v[N0 + i] = interpolate(c, p0, t)
    end
end

export discretize
function discretize(path::BezierPath)
    T = typeof(path.commands[1].p) #XXX BezierPath{T} ?
    v = Vector{T}()
    for c in path.commands
        discretize!(v, c)
    end
    return v
end

#=
using LinearAlgebra: norm, cross, ⋅, normalize

"""continous parametric from of bezier path"""
function (path::BezierPath)(t)
    N = length(path.commands) - 1
    t = t * N
    seg = round(Int, t % N) + 1
    tseg = rem(t, N)
    p0 = endpos(path.commands[seg]) # get last point from prev
    c = path.commands[seg+1]
    if c isa ClosePath
        return interpolate(LineTo(endpos(path.commands[1])), p0, tseg)
    else
        return interpolate(path.commands[seg+1], p0, tseg)
    end
end

endpos(c::MoveTo) = c.p
endpos(c::LineTo) = c.p
endpos(c::CurveTo) = c.p
endpos(c::EllipticalArc) = c(zero(Point2), 1)

interpolate(c::MoveTo, p0, t) = endpos(c)
interpolate(c::LineTo, p0, t) = p0 + (c.p - p0) * t

# experiments on adaptive discretization
function tangent(c::CurveTo, p0, t)
    p1, p2, p3 = c.c1, c.c2, c.p
    normalize(-3(1 - t)^2 * p0 + 3(1 - 4t + 3t^2) * p1 + 3(2t -3t^2) * p2 + 3t^2 * p3)
end

function discretize_adaptive!(v::Vector{<:AbstractPoint}, c::Union{EllipticalArc, CurveTo})
    α = 2 / 360 * 2π
    dtstep = 0.01

    p0 = v[end]

    t = 0.0
    while true
        oldT = tangent(c, p0, t)
        T = oldT
        steps = 0
        while abs(oldT ⋅ T) > cos(α)
            t += dtstep
            steps += 1
            T = tangent(c, p0, t)
        end
        if steps == 1 # if a single step was to much than adjust stepsize
            dtstep = dtstep/2
            continue
        end
        t >= 1.0 && break

        t -= dtstep
        push!(v, interpolate(c, p0, t))
    end
    push!(v, interpolate(c, p0, 1.0))
end
=#


text:
"Parse a tree grid from a string."
function parsetreegrid(s)
    lines = split(s, "\n", keepempty=false)
    width = length(lines[1])
    height = length(lines)
    grid = zeros(Int8, height, width)
    for (i, line) in enumerate(lines)
        for (j, c) in enumerate(line)
            if c == '#'
                grid[i, j] = 1
            end
        end
    end
    grid
end

function walktreegrid(grid, right, down)
    count = 0
    pos = (1,1)
    sz = size(grid)
    while pos[1] <= sz[1]
        if grid[pos...] == 1
            count += 1
        end
        pos = (pos[1]+down, pos[2]+right)
        if pos[2] > sz[2]
            pos = (pos[1], pos[2] % sz[2])
        end
    end
    count
end

day03part1(s) = s |> parsetreegrid |> x->walktreegrid(x, 3, 1)

function day03part2(s)
    grid = parsetreegrid(s)
    dirs = [(1,1), (3,1), (5,1), (7,1), (1,2)]
    mapreduce(x->walktreegrid(grid, x...), *, dirs)
end

function day03(input)
    s = read(input, String)
    p1 = day03part1(s)
    p2 = day03part2(s)
    (p1, p2)
end

#=

s = """
..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#
"""
grid = parsetreegrid(s)

walktreegrid(grid, (1,1)...)

input = read("day03/input.txt", String)
parsetreegrid(input)

day03part1(input)
day03part2(input)

=#


text:
using CorpusLoaders
using CorpusLoaders: @NestedVector
using Test
using Base.Iterators
using MultiResolutionIterators
using DataDeps

@testset "basic use $category" for category in ["train_pos", "train_neg", "test_pos", "test_neg"]
    twitter_gen = load(Twitter(category))
    docs = collect(take(twitter_gen, 10))

    @test typeof(docs) == @NestedVector(String, 3)

    words = flatten_levels(docs, (!lvls)(Twitter, :words))|>full_consolidate
    @test typeof(words) == Vector{String}

    sentences = flatten_levels(docs, (lvls)(Twitter, :documents))|>full_consolidate
    @test length([sent for sent in sentences if length(sent) != 0]) == length(sentences)
    @test sum(length.(sentences)) == length(words)

    @test typeof(full_consolidate(flatten_levels(docs, (!lvls)(Twitter, :documents)))) == Vector{String}
    @test typeof(flatten_levels(docs, (!lvls)(Twitter, :sentences))|>full_consolidate) == Vector{String}

end


text:
module mna

using Reduce
@force using Reduce.Algebra

export nodes, volts, analyze, transferFunction

"""
    nodes(netlist::Dict)::Array{Int64,1}

Takes a netlist dictionary and returns a list of all
nodes in the netlist.
"""
function nodes(netlist::Dict)::Array{Int64,1}
    return ( values(netlist)
                 |> x -> reduce((n,p) ->  [n...,p[2:(end-1)]...], x, init=[])
                 |> unique
                 |> x -> filter((n) -> n != 0, x))
end

"""
    volts(netlist::Dict)::Array{Any,1}

Takes a netlist dictionary and returns all the voltages
in the netlist.
"""
function volts(netlist::Dict)::Array{Any,1}
    return ( keys(netlist)
           |> x -> filter( c -> in( first(netlist[c])
                                  , [:V :VCVS :VCCS :CCVS :CCCS] )
                         , x)
           |> x -> identity.(x) )
end

"""
    systemMatrix(netlist::Dict, ω::Any)::Tuple

Takes a netlist and a angle freuqency and returns the
system matrix, as well as the y vector.
"""
function systemMatrix(netlist::Dict, ω::Any)::Tuple
    m = length(nodes(netlist));
    n = length(volts(netlist));
    M = RExpr(zeros(m+n , m+n)) |> parse |> Reduce.mat;
    y = RExpr(zeros(m+n , 1)) |> parse |> Reduce.mat;

    v = volts(netlist);

    expr = function(type, value)
        if type in [:R :Z :VCCS]
            return (1 / value)
        elseif type == :C
            return (im * (ω == nothing ? :s : ω) * value)
        elseif type == :L
            return (1 / (im * (ω == nothing ? :s : ω) * value))
        else
            return value
        end
    end

    for (id, component) in netlist
        type = first(component);
        nodes = component[2:(end-1)];
        value = expr(first(component), last(component));

        if in(type, [:R :G :L :C :Z :Y])
            idx = ( repeat([nodes], 2)
                  |> x -> Iterators.product(x...)
                  |> collect
                  |> x -> reshape(x, (sum(size(x)),1))
                  |> x -> filter(y -> !in(0,y), x) );
            for (i,j) in idx
                M[i,j] = M[i,j] + (i == j ? value : -value);
            end
        elseif type == :I
            iv = ( identity.(zip(nodes, [-value,value]))
                  |> x -> filter(y -> y[1] != 0, x) );
            for i in iv
                y[i[1]] = y[i[1]] + i[2];
            end
        elseif type == :V
            i_idx = m + findfirst(x -> x == id, v);
            y[i_idx] = y[i_idx] + value;
            m_idx = ( identity.(zip(nodes, [1,-1]))
                    |> x -> filter(y -> y[1] != 0, x) );
            for m in m_idx
                M[m[1], i_idx] = M[m[1], i_idx] + m[2];
                M[i_idx, m[1]] = M[i_idx, m[1]] + m[2];
            end
        elseif type == :VCCS
            i_idx = m + findfirst(x -> x == id, v);
            b_idx = ( identity.(zip(nodes[1:2], [-1 1]))
                    |> x -> filter(y -> y[1] != 0, x) );
            c_idx = ( identity.(zip(nodes[3:4], [1 -1]))
                    |> x -> filter(y -> y[1] != 0, x) );

            M[i_idx, i_idx] = M[i_idx, i_idx] + value;
            for b in b_idx
                M[b[1], i_idx] = M[b[1], i_idx] + b[2];
            end
            for c in c_idx
                M[i_idx, c[1]] = M[i_idx, c[1]] + c[2];
            end
        elseif type == :VCVS
            i_idx = m + findfirst(x -> x == id, v);
            b_idx = ( identity.(zip(nodes[3:4], [1 -1]))
                    |> x -> filter(y -> y[1] != 0, x) );
            c_idx = ( identity.(zip(nodes, [-value value 1 -1]))
                    |> x -> filter(y -> y[1] != 0, x) );
            for b in b_idx
                M[b[1], i_idx] = m[b[1], i_idx] + b[2];
            end
            for c in c_idx
                M[i_idx, c[1]] = m[i_idx, c[1]] + c[2];
            end
        elseif type == :CCVS
            nothing
        elseif type == :CCCS
            nothing
        else
            nothing
        end
    end

    return (M, y)
end

"""
    analyze(netlist::Dict)::Dict

Performs modified nodal analysis (mna) on the given
`netlist` and returns the vector containing node voltages
and branch currents, as well as the transfer function H(ω).
"""
function analyze(netlist::Dict, ω = nothing)

    Reduce.rounded(true);
    Algebra.scientific_notation(2,2);
    M, y = systemMatrix(netlist, ω);

    if ω == nothing
        x = RExpr(M \ y);
    else
        x = convert(Array{Complex, 2}, M) \ convert(Array{Complex, 1}, y);
    end

    return x
end

"""
    transferFunciton(::RExpr) => ::Function

Takes an RExpr with only one symbolic variable `s` and
returns a function that evaluates the small signal model
at the given angle frequency ω.
"""
function transferFunction(x::RExpr)
    return (ω) -> Algebra.sub(:(s = $ω), x) |> parse |> eval
end

end # module


text:
# Pseudosampling experiment
#
# Compares the Stein discrepancy convergence behavior of herding and
# quasi Monte Carlo pseudosampling with that of independent random sampling.
#
# Before running this script, run src/experiments/pseudosample/pseudosample.m from the
# repository base directory to generate pseudosamples.

include("experiment_utils.jl")

using MAT
using SteinDistributions: SteinUniform, gradlogdensity, supportlowerbound, supportupperbound
using SteinDiscrepancy: gsd

# Load generated pseudosamples from experiment results directory
experdir = joinpath("src", "experiments", "pseudosample", "results", "pseudosample")
samplefile = joinpath(experdir,"samples.mat")
samples = matread(samplefile)

# Select an optimization problem solver
solver = "clp"
#solver = "gurobi"
# Specify number of samples and dimension of sampled vectors
n = length(samples["samples_herding"])
# dimension of x
d = 1
# distribution sampled from
distname = "uniform"
# Independent Uniform([0.0,1.0]) with best known Stein constants
(c1,c2,c3) = (0.5,0.5,1.0)
target = SteinUniform(d)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# get the support
lowerbounds = [supportlowerbound(target, 1)]
upperbounds = [supportupperbound(target, 1)]
# Sample sizes at which optimization problem will be solved
ns = 1:200
# Solve optimization problem for each sampler at each sample size
samplers = ["herding","sobol","independent"]

# Collect all objective values for interactive plotting
for sampler in samplers, i = ns
    @printf("Beginning optimization for %s, n=%d, d=%d\n", sampler, i, d)
    # Create program with data subset
    # each column represents a separate trial
    simulateddata = samples["samples_$sampler"]
    numtrials = size(simulateddata, 2)

    for trial in 1:numtrials
        # Compute Stein discrepancy for first i points in this trial
        res = gsd(points=simulateddata[1:i,trial]'',
                  gradlogdensity=gradlogp,
                  c1=c1,
                  c2=c2,
                  c3=c3,
                  supportlowerbounds=lowerbounds,
                  supportupperbounds=upperbounds,
                  solver=solver)
        println("\tn = $(i), objective = $(res.objectivevalue)")

        # Package and save results
        instance_data = Dict{Any, Any}(
            "sampler" => sampler,
            "distname" => distname,
            "n" => i,
            "trial" => trial,
            "d" => d,
            "seed" => samples["seed"],
            "X" => res.points,
            "q" => res.weights,
            "g" => res.g,
            "gprime" => res.gradg,
            "objectivevalue" => res.objectivevalue,
            "edgetime" => res.edgetime,
            "solvetime" => res.solvetime
        )

        save_json(
            instance_data;
            dir="convergence_rates-pseudosamplers-uniform",
            distname=distname,
            sampler=sampler,
            n=i,
            trial=trial,
            d=d,
            seed=samples["seed"]
        )
    end
end

@printf("COMPLETE!")


text:
export iauDat
"""
For a given UTC date, calculate Delta(AT) = TAI-UTC.

   :------------------------------------------:
   :                                          :
   :                 IMPORTANT                :
   :                                          :
   :  A new version of this function must be  :
   :  produced whenever a new leap second is  :
   :  announced.  There are four items to     :
   :  change on each such occasion:           :
   :                                          :
   :  1) A new line must be added to the set  :
   :     of statements that initialize the    :
   :     array "changes".                     :
   :                                          :
   :  2) The constant IYV must be set to the  :
   :     current year.                        :
   :                                          :
   :  3) The "Latest leap second" comment     :
   :     below must be set to the new leap    :
   :     second date.                         :
   :                                          :
   :  4) The "This revision" comment, later,  :
   :     must be set to the current date.     :
   :                                          :
   :  Change (2) must also be carried out     :
   :  whenever the function is re-issued,     :
   :  even if no leap seconds have been       :
   :  added.                                  :
   :                                          :
   :  Latest leap second:  2016 December 31   :
   :                                          :
   :__________________________________________:

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  user-replaceable support function.

Given:
   iy     int      UTC:  year (Notes 1 and 2)
   im     int            month (Note 2)
   id     int            day (Notes 2 and 3)
   fd     double         fraction of day (Note 4)

Returned:
   deltat double   TAI minus UTC, seconds

Returned (function value):
         int      status (Note 5):
                     1 = dubious year (Note 1)
                     0 = OK
                     -1 = bad year
                     -2 = bad month
                     -3 = bad day (Note 3)
                     -4 = bad fraction (Note 4)
                     -5 = internal error (Note 5)

Notes:

   1. UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper
   to call the function with an earlier date.  If this is attempted,
   zero is returned together with a warning status.

   Because leap seconds cannot, in principle, be predicted in
   advance, a reliable check for dates beyond the valid range is
   impossible.  To guard against gross errors, a year five or more
   after the release year of the present function (see the constant
   IYV) is considered dubious.  In this case a warning status is
   returned but the result is computed in the normal way.

   For both too-early and too-late years, the warning status is +1.
   This is distinct from the error status -1, which signifies a year
   so early that JD could not be computed.

   2. If the specified date is for a day which ends with a leap second,
   the TAI-UTC value returned is for the period leading up to the
   leap second.  If the date is for a day which begins as a leap
   second ends, the TAI-UTC returned is for the period following the
   leap second.

   3. The day number must be in the normal calendar range, for example
   1 through 30 for April.  The "almanac" convention of allowing
   such dates as January 0 and December 32 is not supported in this
   function, in order to avoid confusion near leap seconds.

   4. The fraction of day is used only for dates before the
   introduction of leap seconds, the first of which occurred at the
   end of 1971.  It is tested for validity (0 to 1 is the valid
   range) even if not used;  if invalid, zero is used and status -4
   is returned.  For many applications, setting fd to zero is
   acceptable;  the resulting error is always less than 3 ms (and
   occurs only pre-1972).

   5. The status value returned in the case where there are multiple
   errors refers to the first error detected.  For example, if the
   month and day are 13 and 32 respectively, status -2 (bad month)
   will be returned.  The "internal error" status refers to a
   case that is impossible but causes some compilers to issue a
   warning.

   6. In cases where a valid result is not available, zero is returned.

References:

   1. For dates from 1961 January 1 onwards, the expressions from the
   file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used.

   2. The 5ms timestep at 1961 January 1 is taken from 2.58.1 (p87) of
   the 1992 Explanatory Supplement.

Called:
   iauCal2jd    Gregorian calendar to JD

This revision:  2017 October 7

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauDat(int iy, int im, int id, double fd, double *deltat )

function iauDat(iy::Real, im::Real, id::Real, fd::Real)
   # Allocate return values
   ref_dat = Ref{Float64}(0.0)

   status = ccall((:iauDat, libsofa_c), Cint,
                  (Cint, Cint, Cint, Cdouble, Ref{Cdouble}),
                  convert(Int32, iy), convert(Int32, im),
                  convert(Int32, id), convert(Float64, fd),
                  ref_dat)

   if status != 0
       @warn "Non-zero return code form iauD2dtf: $status"
   end

   return status, ref_dat[]
end


text:
using Printf

a = readlines()

ROWS = lastindex(a)
COLS = lastindex(a[begin])

digits = []
for row in 1:ROWS
    for col in 1:COLS
        append!(digits, parse(Int, a[row][col]))
    end
end

basin = []

function contains(what, lst)
    # @printf("contains %d in %s\n", what, lst)
    what in Set(lst)
end

function candidates(i, previous)
    out = []

    col = (i - 1) % COLS
    row = div((i - 1), COLS)

    if col > 0
        append!(out, i - 1)
    end
    if col < (COLS-1)
        append!(out, i + 1)
    end
    if row > 0
        append!(out, i - COLS)
    end
    if row < (ROWS-1)
        append!(out, i + COLS)
    end
    ret = []
    for x in out
        if ! (x in previous)
            append!(ret, x)
        end
    end
    ret
end

for i = 1:lastindex(digits)
    digit = digits[i]
    cands = candidates(i, [])
    low = true
    for c in cands
        if digit >= digits[c]
            low = false
        end
    end
    if low
        append!(basin, i)
    end
end

total = 0
for x in basin
    global total += digits[x] + 1
end
println(total)

basin_sizes = []
additional = []
basin = Set(basin)
seen = Set()

lengths = []

for start in basin
    b = Set([start])
    q = []
    for c in candidates(start, [])
        append!(q, c)
    end
    while 0 < length(q)
        c = popfirst!(q)
        if digits[c] < 9
            ns = candidates(c, b)
            good = true
            for n in ns
                if digits[n] < digits[c]
                    good = false
                end
            end
            if good
                push!(b, c)
                for n in ns
                    append!(q, n)
                end
            end
        end
    end
    append!(lengths, length(b))
end

sorted_lengths = sort(lengths)
# println(sorted_lengths)

println(sorted_lengths[end] * sorted_lengths[end-1] * sorted_lengths[end-2])