text:
module BSON

export BSONObject, ObjectID,
       BSON_OK, BSON_ERROR,
       dict

import Base.show,
       Base.start,
       Base.next,
       Base.done

const BSON_LIB = "libmongoc"
const BSON_OK = 0
const BSON_ERROR = -1

## BSON types ##
const BSON_EOO        = 0
const BSON_DOUBLE     = 1
const BSON_STRING     = 2
const BSON_OBJECT     = 3
const BSON_ARRAY      = 4
const BSON_BINDATA    = 5
const BSON_UNDEFINED  = 6
const BSON_OID        = 7
const BSON_BOOL       = 8
const BSON_DATE       = 9
const BSON_NULL       = 10
const BSON_REGEX      = 11
const BSON_DBREF      = 12
const BSON_CODE       = 13
const BSON_SYMBOL     = 14
const BSON_CODEWSCOPE = 15
const BSON_INT        = 16
const BSON_TIMESTAMP  = 17
const BSON_LONG       = 18

type BSONObject
    _bson::Ptr{Void}

    function BSONObject(_bson::Ptr{Void})
        bson = new(_bson)
        finalizer(bson, destroy)
        return bson
    end
end

type ObjectID
    id::String
    _oid::Ptr{Void}

    function ObjectID(id::String, _oid::Ptr{Void})
        new(id, _oid)
    end
end

ObjectID(id::String) = begin
    _oid = Array(Uint8, 12)
    ccall((:bson_oid_from_string, BSON_LIB), Void, (Ptr{Void}, Ptr{Uint8}), _oid, bytestring(id))
    ObjectID(id, convert(Ptr{Void}, _oid))
end

ObjectID(_oid::Ptr{Void}) = begin
    _str = Array(Uint8, 25)
    ccall((:bson_oid_to_string, BSON_LIB), Void, (Ptr{Void}, Ptr{Uint8}), _oid, _str)
    id = bytestring(_str[1:24])
    ObjectID(id, _oid)
end

function show(io::IO, oid::ObjectID)
    print(io, "ObjectID(\"$(oid.id)\")")
end

BSONObject(bson::BSONObject) = bson
BSONObject(dict::Associative) = begin
    _bson = ccall((:bson_alloc, BSON_LIB), Ptr{Void}, ())
    if _bson == C_NULL
        error("Unable to create BSON object")
    end
    ccall((:bson_init, BSON_LIB), Void, (Ptr{Void},), _bson)

    build(_bson, dict)

    if ccall((:bson_finish, BSON_LIB), Int32, (Ptr{Void},), _bson) == BSON_ERROR
        error("Unable to build BSON object")
    end

    BSONObject(_bson)
end

BSONObject() = BSONObject(Dict{Any,Any}())

function show(io::IO, bson::BSONObject)
    show(io, dict(bson))
end

function dict(bson::BSONObject)
    d = Dict{Any, Any}()
    for (k,v) in bson
        if v == BSONObject || v == Array
            d[k] = dict(v)
        else
            d[k] = v
        end
    end
    return d
end

immutable Sentinel end
function Base.get(bson::BSONObject, key::String, default)
    for (k,v) in bson
        if k == key
            return v
        end
    end
    default
end
function Base.get(bson::BSONObject, key::String)
    val = get(bson, key, Sentinel())
    if is(val, Sentinel())
        error("key not found: $(repr(key))")
    end
    val
end
Base.getindex(bson::BSONObject, key) = get(bson, key)

## Iterator ##

start(b::BSONObject) = begin
    _iterator = ccall((:bson_iterator_alloc, BSON_LIB), Ptr{Void}, ())
    ccall((:bson_iterator_init, BSON_LIB), Void, (Ptr{Void}, Ptr{Void}), _iterator, b._bson)
    _iterator
end

next(b::BSONObject, _iterator) = begin
    ((key(_iterator), value(_iterator)), _iterator)
end

done(b::BSONObject, _iterator) = begin
    bson_type = ccall((:bson_iterator_next, BSON_LIB), Int32, (Ptr{Void},), _iterator)
    if bson_type == BSON_EOO
        ccall((:bson_iterator_dealloc, BSON_LIB), Void, (Ptr{Void},), _iterator)
        return true
    end
    false
end

## Private methods ##

function key(_iterator::Ptr{Void})
    str = ccall((:bson_iterator_key, BSON_LIB), Ptr{Uint8}, (Ptr{Void},), _iterator)
    bytestring(str)
end

function value(_iterator::Ptr{Void})
    bson_type = ccall((:bson_iterator_type, BSON_LIB), Int32, (Ptr{Void},), _iterator)

    if BSON_DOUBLE == bson_type
        ccall((:bson_iterator_double, BSON_LIB), Float64, (Ptr{Void},), _iterator)
    elseif BSON_INT == bson_type
        ccall((:bson_iterator_int, BSON_LIB), Int32, (Ptr{Void},), _iterator)
    elseif BSON_LONG == bson_type
        ccall((:bson_iterator_long, BSON_LIB), Int64, (Ptr{Void},), _iterator)
    elseif BSON_STRING == bson_type || BSON_SYMBOL == bson_type
        _str = ccall((:bson_iterator_string, BSON_LIB), Ptr{Uint8}, (Ptr{Void},), _iterator)
        bytestring(_str)
    elseif BSON_OID == bson_type
        _oid = ccall((:bson_iterator_oid, BSON_LIB), Ptr{Void}, (Ptr{Void},), _iterator)
        ObjectID(_oid)
    elseif BSON_BOOL == bson_type
        _bool = ccall((:bson_iterator_bool, BSON_LIB), Int32, (Ptr{Void},), _iterator)
        _bool == 0 ? false : true
    elseif BSON_DATE == bson_type
        ccall((:bson_iterator_date, BSON_LIB), Int64, (Ptr{Void},), _iterator)
    elseif BSON_REGEX == bson_type
        _str = ccall((:bson_iterator_regex, BSON_LIB), Ptr{Uint8}, (Ptr{Void},), _iterator)
        Regex(bytestring(_str))
    elseif BSON_NULL == bson_type || BSON_UNDEFINED == bson_type
        nothing
    elseif BSON_OBJECT == bson_type
        _bson = ccall((:bson_alloc, BSON_LIB), Ptr{Void}, ())
        ccall((:bson_iterator_subobject, BSON_LIB), Void, (Ptr{Void}, Ptr{Void}), _iterator, _bson)
        BSONObject(_bson)
    elseif BSON_ARRAY == bson_type
        a = {}
        _subiterator = ccall((:bson_iterator_alloc, BSON_LIB), Ptr{Void}, ())
        ccall((:bson_iterator_subiterator, BSON_LIB), Void, (Ptr{Void}, Ptr{Void}), _iterator, _subiterator)
        while(ccall((:bson_iterator_next, BSON_LIB), Int32, (Ptr{Void},), _subiterator) != BSON_EOO)
            push!(a, value(_subiterator))
        end
        ccall((:bson_iterator_dealloc, BSON_LIB), Void, (Ptr{Void},), _subiterator)
        a
    else
        # Not supported:
        #   BSON_TIMESTAMP
        #   BSON_BINDATA
        #   BSON_CODE
        #   BSON_CODEWSCOPE
        #
        error("Unsupported type: $bson_type")
    end
end

function build(_bson::Ptr{Void}, dict::Dict{Any,Any})
    for (k,v) in dict
      append(_bson, k, v)
    end
end

append(_bson::Ptr{Void}, k::String, v::Union(ASCIIString, UTF8String)) =
    ccall((:bson_append_string, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8}), _bson, bytestring(k), bytestring(v))
append(_bson::Ptr{Void}, k::String, v::Int32) =
    ccall((:bson_append_int, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Int32), _bson, bytestring(k), v)
append(_bson::Ptr{Void}, k::String, v::Int64) =
    ccall((:bson_append_long, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Int64), _bson, bytestring(k), v)
append(_bson::Ptr{Void}, k::String, v::Union(Float32, Float64)) =
    ccall((:bson_append_double, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Float64), _bson, bytestring(k), v)
append(_bson::Ptr{Void}, k::String, v::ObjectID) =
    ccall((:bson_append_oid, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}), _bson, bytestring(k), v._oid)
append(_bson::Ptr{Void}, k::String, v::Bool) =
    ccall((:bson_append_bool, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Int32), _bson, bytestring(k), v == true ? 1 : 0)
function append(_bson::Ptr{Void}, k::String, v::Dict)
    ccall((:bson_append_start_object, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}), _bson, bytestring(k))
    build(_bson, v)
    ccall((:bson_append_finish_object, BSON_LIB), Int32, (Ptr{Void},), _bson)
end
function append(_bson::Ptr{Void}, k::String, v::Vector)
    ccall((:bson_append_start_array, BSON_LIB), Int32, (Ptr{Void}, Ptr{Uint8}), _bson, bytestring(k))
    for i in 1:length(v)
        append(_bson, string(i-1), v[i])
    end
    ccall((:bson_append_finish_array, BSON_LIB), Int32, (Ptr{Void},), _bson)
end
function append(_bson::Ptr{Void}, k::String, v::Any)
    # Not supported:
    #   BSON_REGEX
    #   BSON_DATE
    #   BSON_TIMESTAMP
    #   BSON_BINDATA
    #   BSON_CODE
    #   BSON_CODEWSCOPE
    #
    error("Unsupported type: $t")
end

function destroy(bson::BSONObject)
    ccall((:bson_destroy, BSON_LIB), Void, (Ptr{Void},), bson._bson)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5720d27a-f176-11ea-35ec-efcc84b4e2a9
using NiLang, Plots

# ╔═╡ d87a3c74-f17a-11ea-1ce7-31659781efee
using ReversibleSeismic: i_one_step!, AcousticPropagatorParams, Ricker

# ╔═╡ b0d4d10c-f191-11ea-2f84-1b19dcb85d6a
using Optim

# ╔═╡ 701eab34-f17f-11ea-17e9-fbfa0b20ab76
md"# To define a loss"

# ╔═╡ 937e8aa6-f17a-11ea-1d5d-6f495a29d584
@i function i_solve!(param::AcousticPropagatorParams, src,
            srcv::Array{T, 1}, c::Array{T, 2},
            tu::Array{T,3}, tφ::Array{T,3}, tψ::Array{T,3}) where T

    @routine begin
        d2 ← zero(param.DELTAT)
        d2 += param.DELTAT^2
    end
    for i = 3:param.NSTEP+1
        i_one_step!(param, view(tu,:,:,i), view(tu,:,:,i-1), view(tu,:,:,i-2),
            view(tφ,:,:,i), view(tφ,:,:,i-1), view(tψ,:,:,i), view(tψ,:,:,i-1), param.Σx, param.Σy, c)
        tu[src[1], src[2], i] += srcv[i-2] * d2
    end
    ~@routine
end

# ╔═╡ 453bb1f2-f18e-11ea-086f-d18ed1055e0f
nstep = 1500

# ╔═╡ 0fa21e04-f192-11ea-03db-7d8b07acfa91
c = 700 * (1 .+ sin.(LinRange(0, 5π, 82))' .* cos.(LinRange(0, 3π, 82)));

# ╔═╡ b4a4159a-f178-11ea-07fc-059d9c6119cb
"""
the reversible loss
"""
@i function i_loss!(out::T, param, src, srcv::AbstractVector{T}, c0::AbstractMatrix{T}, tu::AbstractArray{T,3}, tφ::AbstractArray{T,3}, tψ::AbstractArray{T,3}) where T
	@routine begin
		c ← zero(c0)
		for i=1:length(c0)
			c[i] += abs(c0[i])
		end
	end
    i_solve!(param, src, srcv, c, tu, tφ, tψ)
	out -= tu[size(tu,1)÷2,size(tu,2)÷2+20,end]
	~@routine
end

# ╔═╡ 66d4e274-f180-11ea-0831-a946570eca63
function loss(c; nstep)
	nx, ny = size(c) .- 2
	param = AcousticPropagatorParams(nx=nx, ny=ny,
	Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)

 	tu = zeros(nx+2, ny+2, nstep+1)
 	tφ = zeros(nx+2, ny+2, nstep+1)
 	tψ = zeros(nx+2, ny+2, nstep+1)

 	src = (nx ÷ 2, ny ÷ 2)
 	srcv = Ricker(param, 100.0, 500.0)

 	res = i_loss!(0.0, param, src, srcv, copy(c), tu, tφ, tψ)
	res[1], res[7]
end

# ╔═╡ a8f5a20a-f195-11ea-167f-7d408ebe0366
heatmap(c)

# ╔═╡ 3b10da96-f18d-11ea-16fa-bf8437f54903
tu_seq = loss(c, nstep=nstep)[2];

# ╔═╡ cd4c73b6-f17c-11ea-08be-5bf044147749
ani = @animate for i=1:100:nstep
	heatmap(tu_seq[2:end-1,2:end-1,i], clim=(-0.005, 0.005))
end

# ╔═╡ cea59342-f17d-11ea-11c3-5ddc3062b120
gif(ani, fps=5)

# ╔═╡ 03c0a192-f198-11ea-3a75-f7b26638fc6a
heatmap(tu_seq[:,:,end])

# ╔═╡ 5b9381f2-f17f-11ea-3fb0-c5f5314c3c69
md"# Obtaining gradients"

# ╔═╡ 0681da22-f17a-11ea-1fca-a7d4004af233
"""
obtain gradients with NiLang.AD
"""
function getgrad(c::AbstractMatrix{T}; nstep::Int) where T
     param = AcousticPropagatorParams(nx=size(c,1)-2, ny=size(c,2)-2,
          Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=nstep)

     c = copy(c)
     tu = zeros(T, size(c)..., nstep+1)
     tφ = zeros(T, size(c)..., nstep+1)
     tψ = zeros(T, size(c)..., nstep+1)

     src = size(c) .÷ 2 .- 1
     srcv = Ricker(param, 100.0, 500.0)
     NiLang.AD.gradient(Val(1), i_loss!, (0.0, param, src, srcv, c, tu, tφ, tψ))[end-3]
end;

# ╔═╡ 5f8782ca-f17a-11ea-082a-a98ea089b700
gradc = getgrad(c; nstep=nstep);

# ╔═╡ 7bea3116-f18d-11ea-19d8-e5ef06dd114a
heatmap(gradc)

# ╔═╡ 16504abe-f180-11ea-33b6-bfb2937f0846
res = optimize(c->(@show loss(c, nstep=nstep)[1]), (g, c)->(g.=getgrad(c; nstep=nstep)), c, BFGS(), Optim.Options(iterations=5, g_tol=1e-30, f_tol=1e-30))

# ╔═╡ 07a889aa-f193-11ea-3a0d-1bacdd6cc2cb
heatmap(abs.(res.minimizer) - abs.(c))

# ╔═╡ 5a869532-f18d-11ea-2f3e-2b5886e7fa70
tu_seq2 = loss(res.minimizer, nstep=nstep)[2];

# ╔═╡ 5cd02706-f181-11ea-217e-c167fc0a4164
ani2 = @animate for i=1:nstep
	heatmap(tu_seq2[2:end-1,2:end-1,i], clim=(-0.005, 0.005))
end every 100

# ╔═╡ af3428da-f181-11ea-0ea8-f9bfabc320d3
gif(ani2, fps=5)

# ╔═╡ f36eea06-f192-11ea-003e-b91e794f3a6f
heatmap(tu_seq2[2:end-1,2:end-1,end], clim=(-0.005, 0.005))

# ╔═╡ Cell order:
# ╠═5720d27a-f176-11ea-35ec-efcc84b4e2a9
# ╠═d87a3c74-f17a-11ea-1ce7-31659781efee
# ╟─701eab34-f17f-11ea-17e9-fbfa0b20ab76
# ╠═937e8aa6-f17a-11ea-1d5d-6f495a29d584
# ╠═b4a4159a-f178-11ea-07fc-059d9c6119cb
# ╠═66d4e274-f180-11ea-0831-a946570eca63
# ╠═453bb1f2-f18e-11ea-086f-d18ed1055e0f
# ╠═0fa21e04-f192-11ea-03db-7d8b07acfa91
# ╠═a8f5a20a-f195-11ea-167f-7d408ebe0366
# ╠═3b10da96-f18d-11ea-16fa-bf8437f54903
# ╠═cd4c73b6-f17c-11ea-08be-5bf044147749
# ╠═cea59342-f17d-11ea-11c3-5ddc3062b120
# ╠═03c0a192-f198-11ea-3a75-f7b26638fc6a
# ╟─5b9381f2-f17f-11ea-3fb0-c5f5314c3c69
# ╠═0681da22-f17a-11ea-1fca-a7d4004af233
# ╠═5f8782ca-f17a-11ea-082a-a98ea089b700
# ╠═7bea3116-f18d-11ea-19d8-e5ef06dd114a
# ╠═b0d4d10c-f191-11ea-2f84-1b19dcb85d6a
# ╠═16504abe-f180-11ea-33b6-bfb2937f0846
# ╠═07a889aa-f193-11ea-3a0d-1bacdd6cc2cb
# ╠═5a869532-f18d-11ea-2f3e-2b5886e7fa70
# ╠═5cd02706-f181-11ea-217e-c167fc0a4164
# ╠═af3428da-f181-11ea-0ea8-f9bfabc320d3
# ╠═f36eea06-f192-11ea-003e-b91e794f3a6f

ext: jl


text:
"""
    hopper 3D
        orientation representation: modified rodrigues parameters
		similar to Raibert hopper, all mass is located at the body
		s = (px, py, pz, tx, ty, tz, r) = (3d_position, MRP(m1, m2, 0.0), leg_length)
"""
struct Hopper3D{T} <: Model{T}
    nq::Int
	nu::Int
	nw::Int
	nc::Int

	mb::T # mass of body
    ml::T # mass of leg
    Jb::T # inertia of body
    Jl::T # inertia of leg

    μ_world::T  # coefficient of friction
	μ_joint::T
    g::T # gravity

	orientation::Symbol

	base::BaseMethods
	dyn::DynamicsMethods

	joint_friction::SVector
end

lagrangian(model::Hopper3D, q, q̇) = 0.0

# Kinematics
function kinematics(model::Hopper3D, q)
	p = view(q, 1:3)
	R = eval(model.orientation)(view(q, 4:6)...)
	p + R * [0.0; 0.0; -1.0 * q[7]]
end

# Methods
function M_func(model::Hopper3D, q)
	Diagonal(@SVector [model.mb + model.ml, model.mb + model.ml, model.mb + model.ml,
					   model.Jb + model.Jl, model.Jb + model.Jl, model.Jb + model.Jl,
					   model.ml])
end

function C_func(model::Hopper3D, q, q̇)
	@SVector [0.0, 0.0, (model.mb + model.ml) * model.g, 0.0, 0.0, 0.0, 0.0]
end

function ϕ_func(model::Hopper3D, env::Environment, q)
	SVector{1}(kinematics(model, q)[3] - env.surf(kinematics(model, q)[1:2]))
end

function B_func(::Hopper3D, q)
    rot = view(q, 4:6)
    R = eval(model.orientation)(rot...)
    @SMatrix [0.0 0.0 0.0 R[1,1] R[2,1] R[3,1] 0.0;
              0.0 0.0 0.0 R[1,2] R[2,2] R[3,2] 0.0;
			  R[1,3] R[2,3] R[3,3] 0.0 0.0 0.0 1.0]
end

function A_func(::Hopper3D, q)
    rot = view(q, 4:6)
    R = eval(model.orientation)(rot...)
    @SMatrix [1.0 0.0 0.0 0.0 0.0 0.0 0.0;
              0.0 1.0 0.0 0.0 0.0 0.0 0.0;
			  0.0 0.0 1.0 0.0 0.0 0.0 0.0]
end

function J_func(model::Hopper3D, env::Environment, q)
    k(z) = kinematics(model, z)
    ForwardDiff.jacobian(k, q)
end

function contact_forces(model::Hopper3D, env::Environment{<:World, LinearizedCone}, γ1, b1, q2, k)
	m = friction_mapping(env)

	SVector{3}(transpose(rotation(env, k)) * [m * b1; γ1])
end

function velocity_stack(model::Hopper3D, env::Environment{<:World, LinearizedCone}, q1, q2, k, h)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]

	v1_surf = rotation(env, k) * v

	SVector{4}(friction_mapping(env)' * v1_surf[1:2])
end

function get_stride(model::Hopper3D, traj)
    stride = zeros(model.nq)
    stride[1:2] = traj.q[end-1][1:2] - traj.q[1][1:2]
    return stride
end

# Dimensions
nq = 7 # configuration dimension
nu = 3 # control dimension
nw = 3 # disturbance dimension
nc = 1 # number of contact points
nquat = 0 # number of quaternions

# Parameters
g = 9.81 # gravity
μ_world = 1.5 # coefficient of friction
μ_joint = 0.0

# TODO: change to Raibert parameters
mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

hopper_3D = Hopper3D(nq,nu,nw,nc,
			mb, ml, Jb, Jl,
			μ_world, μ_joint, g,
			:MRP,
			BaseMethods(), DynamicsMethods(),
			SVector{7}(zeros(7)))

function friction_coefficients(model::Hopper3D)
	return [model.μ_world]
end

function initialize_z!(z, model::Hopper3D, idx::RoboDojo.IndicesZ, q)
    z .= 1.0
    z[idx.q] .= q
end


text:
# This is for the PriorityQueue
using DataStructures

# creating abstract type for binary trees
abstract type BT end
type Empty <: BT end

type Node <: BT
    right::BT
    left::BT
    weight::Int64
    key::String
    bitstring::String
end

# This will search through the tree to create bitstrings for each character
# This outputs encoding Dict to be used for encoding
function create_codebook(n::Node)
    s = Stack(Node)
    push!(s, n)

    codebook = Dict{Char, String}()

    while (length(s) > 0)
        temp = pop!(s)
        if (typeof(temp.left) != Empty)
            temp.left.bitstring = temp.bitstring * "0"
            push!(s, temp.left)
        end
        if (typeof(temp.right) != Empty)
            temp.right.bitstring = temp.bitstring * "1"
            push!(s, temp.right)
        end
        if (typeof(temp.right) == Empty &&
            typeof(temp.left) == Empty)
            codebook[temp.key[1]] = temp.bitstring
        end
    end

    return codebook
end

# This outputs huffman tree to generate dictionary for encoding
function create_tree(phrase::String)

    # creating weights
    weights = PriorityQueue()
    for i in phrase
        temp_string = string(i)
        if (haskey(weights, temp_string))
            weights[temp_string] += 1
        else
            weights[temp_string] = 1
        end
    end

    # Creating all nodes to iterate through
    nodes = PriorityQueue{Node, Int64}()
    while(length(weights) > 0)
        weight = peek(weights)[2]
        key = dequeue!(weights)
        temp_node = Node(Empty(), Empty(), weight, key, "")
        enqueue!(nodes, temp_node, weight)
    end

    while(length(nodes) > 1)
        node1 = dequeue!(nodes)
        node2 = dequeue!(nodes)
        temp_node = Node(node1, node2, node1.weight + node2.weight, "", "")
        enqueue!(nodes, temp_node, temp_node.weight)
    end

    huffman_tree = dequeue!(nodes)
    return huffman_tree

end

function encode(codebook::Dict{Char, String}, phrase::String)
    final_bitstring = ""
    for i in phrase
        final_bitstring = final_bitstring * codebook[i]
    end

    return final_bitstring
end

function decode(huffman_tree::Node, bitstring::String)
    current = huffman_tree
    final_string = ""
    for i in bitstring
        if (i == '0')
            current = current.left
        else
            current = current.right
        end
        if(typeof(current.left) == Empty && typeof(current.right) == Empty)
            final_string = final_string * current.key
            current = huffman_tree
        end
    end

    return final_string
end

function two_pass_huffman(phrase::String)
    huffman_tree = create_tree(phrase)
    codebook = create_codebook(huffman_tree)
    println(codebook)
    bitstring = encode(codebook, phrase)
    final_string = decode(huffman_tree, bitstring)
    println(bitstring)
    println(final_string)
end

two_pass_huffman("aaaabbbccd")

max_issues_repo_name: shiffman/simuleios

max_forks_repo_name: shiffman/simuleios


text:
using OnlineAI

function runtest(net, x, y, niter)
  for i in 1:niter
    fit!(net, x, y)
  end
end

nin, nout = 2, 1
hidden = [100,100]
params = NetParams(
    # updater = SGDModel()
    updater = AdaMaxModel()
  )
net = buildTanhClassificationNet(nin, nout, hidden; params = params)
x = randn(nin)
y = randn(nout)

runtest(net, x, y, 1)

Profile.clear()

# note: last pass this ran in ~8 secs
@time @profile runtest(net, x, y, 30000)


text:
include("ECapeExceptions.jl")
include("CapeOpen.jl")
module ThermodynamicsTable
  using ECapeExceptions
  import CapeOpen.perryanalytic
  thisfiledirname=dirname(@__FILE__())

  pabf=open(joinpath(thisfiledirname,"..","lib",perryanalytic.datafilename*".binary"),"r");
  pata=perryanalytic.tableaddreses
  adrs=readdlm(joinpath(thisfiledirname,"..","lib",perryanalytic.datafilename*".addresses"),';')
  for i=1:size(adrs)[1]
    pata[adrs[i]][1][1:3]=adrs[i,2:4]
  end
  function getdatamatrix(path::AbstractString)
    global floattables,compondtable
    file=open(thisfiledirname * "/Tables/"*path);
    matrix=readdlm(file,';',header=false);
    close(file)
    return matrix
  end

  #for LiquidsDensities, Cp, LiquidsCp props Max of two data row per prop may exists
  function readbinarydatabase(id::UInt16, table::ASCIIString, skipdata::UInt8=0%UInt8)
    if haskey(pata,table)
      v=pata[table]
    else
      throw(ECapeInvalidArgument())
    end
    seek(pabf,v[1][1])
    skipsize=v[1][3]-sizeof(UInt16)
    j=1
    while(id!=(read(pabf,UInt16)) && j<=v[1][2])
      skip(pabf,skipsize)
      j+=1
    end
    if(skipdata>0 && j+skipdata<=v[1][2])
      skip(pabf,skipdata*skipsize)
      id!=(read(pabf,UInt16)) && throw(ECapeInvalidArgument())
    end
    if(j<=v[1][2])
      for jj in 2:length(v)
        if isa(v[jj],Array)
          read!(pabf,v[jj])
        else
          v[jj]=read(pabf,typeof(v[jj]))
        end
      end
      return v[2:end]
    end
    throw(ECapeLimitedImpl())
  end
end

include("PhysicalPropertyCalculator.jl")
include("ICapeThermoCompounds.jl")
include("ICapeThermoMaterial.jl")
include("ICapeThermoPhases.jl")
include("ICapeThermoPropertyRoutine.jl")
include("ICapeThermoUniversalConstants.jl")
using ICapeThermoUniversalConstants,ICapeThermoCompounds,ECapeExceptions
perryanalytic = CapeOpen.perryanalytic
println("`perryanalytic` Property Package was created and is ready to use.")
include("CoolProp.jl")
include("FreeSteam.jl")


text:
using DrWatson
@quickactivate "OneVision"

using PackageCompiler
include("pkg_list.jl")
println("Creating Sysimage...")
time_taken = @elapsed create_sysimage(pkg_list,
    sysimage_path="JuliaSysimage.dylib",
    precompile_execution_file="scripts/precompile/run_tests.jl"
)
println(time_taken)


text:
using LinearAlgebra
using GLMakie

# LSP indexing solution
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/SpinDoctor.jl")
    using .SpinDoctor
else
    using SpinDoctor
end

set_theme!(theme_dark())

## Choose setup recipe
# include("setups/axon.jl")
# include("setups/sphere.jl")
include("setups/plates.jl")
# include("setups/cylinders.jl")
# include("setups/spheres.jl")
# include("setups/neuron.jl")

## Create model from setup recipe
femesh, = @time create_geometry(setup; recreate = true)
model = Model(; mesh = femesh, coeffs...)
volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)
ncompartment = length(model.mesh.points)
matrices = @time assemble_matrices(model);
@info "Number of nodes per compartment:" length.(model.mesh.points)

## Magnetic field gradient
dir = [1.0, 1.0, 1.0]
profile = PGSE(2500.0, 4000.0)
b = 1000
g = √(b / int_F²(profile)) / model.γ
gradient = ScalarGradient(dir, profile, g)

## Short term approximation (STA)
adc_sta_cmpts = compute_adc_sta(model, gradient)

## Solve BTPDE and fit ADC for different b-values
bvalues = 0:100:1000
gvalues = map(b -> √(b / int_F²(profile)) / coeffs.γ, bvalues)
gradients = [ScalarGradient(gradient.dir, gradient.profile, g) for g ∈ gvalues]
# btpde = GeneralBTPDE(;
#     model, matrices, reltol = 1e-4, abstol = 1e-6, odesolver = QNDF(autodiff = false),
# )
btpde = IntervalConstanBTPDE{T}(; model, matrices, θ = 0.5, timestep = 5)
ξ, = solve_multigrad(btpde, gradients)
signals = [compute_signal(matrices.M, ξ) for ξ ∈ ξ]
signals_cmpts = [compute_signal.(matrices.M_cmpts, split_field(model.mesh, ξ)) for ξ ∈ ξ]
adc_fit = fit_adc(bvalues, signals)
adc_fit_cmpts =
    [fit_adc(bvalues, [s[icmpt] for s ∈ signals_cmpts]) for icmpt = 1:ncompartment]

## Solve HADC
hadc = HADC(; model, matrices, odesolver = QNDF(), reltol = 1e-4, abstol = 1e-6)
adc_homogenized_cmpts = @time solve(hadc, gradient)

## Matrix Formalism ADC
# Perform Laplace eigendecomposition
laplace = Laplace{T}(; model, matrices, neig_max = 400)
lap_eig = @time solve(laplace)

D_mf = compute_mf_diffusion_tensor(model.mesh, matrices.M, lap_eig, gradient)
adc_mf = dir'D_mf * dir / dir'dir


text:
module CodeGen

using MacroTools
using JSON
using CxxWrap

include("stub.jl")

const JC_TYPEMAP = Dict(
    :Int8          => "int8_t",
    :UInt8         => "uint8_t",
    :Int16         => "int16_t",
    :UInt16        => "uint16_t",
    :Int32         => "int32_t",
    :UInt32        => "uint32_t",
    :Int64         => "int64_t",
    :UInt64        => "uint64_t",
    :CxxBool       => "bool",
    :CxxChar       => "char",
    :CxxWchar      => "wchar_t",
    :CxxSignedChar => "signed char",
    :CxxUChar      => "unsigned char",
    :CxxShort      => "short",
    :CxxUShort     => "unsigned short",
    :CxxInt        => "int",
    :CxxULongLong  => "unsigned long long",
    :CxxLong       => "long",
    :CxxULong      => "unsigned long",
    :CxxLongLong   => "long long",
    :CxxULongLong  => "unsigned long long",
    :Cstring       => "char*",
    :Cvoid         => "void",
    :Cdouble       => "double",
    :StdString => "std::string",
    :StdVector => "std::vector",
    :SharedPtr => "std::shared_ptr",
)

const DEREFERENCE_BLACKLIST = []

const REFERENCE_BLACKLIST = [
    :SharedPtr
    :UniquePtr
    :WeakPtr
]

const LTYPE_TO_RTYPE = Dict(:StdString => :Cstring)

function gen_call(f, ex...)
    @info typeof(f)
    callex = Expr(:call, f)
    for e in ex
        if @capture(e, a_ = b_)
            push!(callex.args, Expr(:kw, a, b))
        else
            push!(callex.args, e isa Symbol ? QuoteNode(e) : e)
        end
    end
    return callex
end

struct PType{N}
    name::Symbol
    parameters::NTuple{N,<:Union{Symbol,<:PType}}
end
PType(name::Symbol, parameters::Union{Symbol,<:PType}...) = PType(name, parameters)
const CType = Union{Symbol,<:PType}

struct Field
    name::Symbol
    type::CType
    dereference::Bool
    mutable::Bool
end
function Field(name, type; dereference=should_dereference(type), mutable=true)
    return Field(name, type, dereference, mutable)
end
macro field(ex...)
    return gen_call(Field, ex...)
end

struct Method
    name::Symbol
    cppmethod::String
end
Method(name, cppmethod) = Method(name, cppmethod)

struct Struct
    name::Symbol
    fields::Vector{Field}
    methods::Vector{Method}
end
Struct(name; fields=Field[], methods=Method[]) = Struct(name, fields, methods)
macro type(ex...)
    return gen_call(Struct, ex...)
end

struct GlobalMethod
    name::Symbol
    cppmethod::String
    override_module::Union{Nothing,Symbol}
    GlobalMethod(name, cppmethod, override_module=nothing) = new(name, cppmethod, override_module)
end

function operator(m::Symbol, A::Struct, B::Struct=A)
    @assert m in (:(==), :(>=), :(<=), :(>), :(<), :(-), :(+))
    return GlobalMethod(m, "[]($(A.name)& a, $(B.name)& b) { return a $m b; }", :Base)
end

cname(type::Symbol) = get(JC_TYPEMAP, type, type)
function cname(type::PType)
    return "$(cname(type.name))<$(join(map(cname, type.parameters), ", "))>"
end

write_jl_typestub(path::AbstractString, args...; kwargs...) =
    _write(write_jl_typestub, path, args, kwargs)
function write_jl_typestub(io::IO, type::Struct)
    function gen(i)
        field = type.fields[i]
        return (
            field,
            :(name === $(QuoteNode(field.name))),
            if field.dereference
                :(CxxWrap.CxxWrapCore.dereference_argument(
                    $(Symbol(type.name, :_get_, field.name))(x)
                ))
            else
                :($(Symbol(type.name, :_get_, field.name))(x))
            end,
            :($(Symbol(type.name, :_set_, field.name))(x, value)),
        )
    end

    get_ex = Expr(:block, :($Base.getfield(x, name)))
    set_ex = Expr(:block, :($Base.setfield!(x, name, value)))
    for i in 2:length(type.fields)
        field, case, getfield, setfield! = gen(i)
        get_ex = Expr(:elseif, case, getfield, get_ex)
        if field.mutable
            set_ex = Expr(:elseif, case, setfield!, set_ex)
        end
    end
    field, case, getfield, setfield! = gen(1)
    get_ex = Expr(:if, case, getfield, get_ex)
    if field.mutable
        set_ex = Expr(:if, case, setfield!, set_ex)
    end

    T = :(Union{$(type.name),SmartPointer{$(type.name)}})
    ex = quote
        function Base.getproperty(x::$T, name::Symbol)
            Base.@_inline_meta
            return $get_ex
        end
        function Base.setproperty!(x::$T, name::$Symbol, value)
            Base.@_inline_meta
            $set_ex
            return x
        end
        function Base.propertynames(x::$T)
            Base.@_inline_meta
            # return $(Expr(:tuple, map(QuoteNode, fields.args)...))
            return $(Tuple(f.name for f in type.fields))
        end
    end

    write(io, read(joinpath(@__DIR__, "stub.jl")))
    for line in prettify(ex).args
        write(io, string(line), '\n')
    end

    return nothing
end

write_cpp_typestub(path::AbstractString, args...; kwargs...) =
    _write(write_cpp_typestub, path, args, kwargs)
function write_cpp_typestub(io::IO, mod::Symbol, type::Struct; extra=[])
    name = type.name
    buf = IOBuffer()
    write(buf, "$(mod).add_type<$name>(\"$name\")")

    for field in type.fields
        fname = field.name
        T = cname(field.type)
        R = isreferencetype(field.type) ? "auto&" : "auto"
        write(
            buf,
            ".method(\"$(name)_get_$(fname)\", [](const $(name)& x) -> $R { return x.$(fname); })",
        )
        if field.mutable
            write(
                buf,
                ".method(\"$(name)_set_$(fname)\", []($(name)& x, const $T y) { x.$(fname) = y; })",
            )
        end
    end

    for method in type.methods
        write(buf, ".method(\"$(method.name)\", $(method.cppmethod))")
    end

    extra = extra isa AbstractString ? [extra] : extra
    for l in extra
        !startswith(l, '.') && write(buf, '.')
        write(buf, l)
    end

    format_cpp(buf)
    write(io, seekstart(buf), ';', '\n')
    return nothing
end

write_cpp_methodstub(path::AbstractString, args...; kwargs...) =
    _write(write_cpp_methodstub, path, args, kwargs)
function write_cpp_methodstub(io::IO, mod::Symbol, methods::Vector{GlobalMethod}; extra=[])
    buf = IOBuffer()
    for method in methods
        if method.override_module === mod
            write(buf, "$mod.set_override_module($mod.julia_module());")
        elseif method.override_module === :Base
            write(buf, "$mod.set_override_module(jl_base_module);")
        elseif method.override_module !== nothing
            error("Unknown override_module $(method.override_module)")
        end
        write(buf, "$mod.method(\"$(method.name)\", $(method.cppmethod));")
        if method.override_module !== nothing
            write(buf, "$mod.unset_override_module();")
        end
    end

    extra = extra isa AbstractString ? [extra] : extra
    for l in extra
        write(buf, l)
        !endswith(l, ';') && write(buf, ';')
    end

    format_cpp(buf)
    write(io, seekstart(buf), '\n')
    return nothing
end

function _write(f::Function, path::AbstractString, args, kwargs)
    @info "Generating $path"
    return open(io -> f(io, args...; kwargs...), path, "w")
end

writeln(io::IO, xs...) = write(io, xs..., '\n')

function format_cpp(stdin::IO)
    style = JSON.json(Dict("BasedOnStyle" => "google", "ColumnLimit" => 120))
    clang_cmd = `clang-format --style=$style`
    if success(`nix --help`)
        cmd = `nix shell nixpkgs'#'libclang -c $clang_cmd`
    elseif success(`clang-format --help`)
        cmd = `clang-format`
    else
        return nothing
    end

    stdout = IOBuffer()
    seekstart(stdin)
    x = run(pipeline(seekstart(stdin), `nix shell nixpkgs'#'libclang -c $clang_cmd`, stdout))
    write(seekstart(stdin), seekstart(stdout))

    return nothing
end

function isreferencetype(type::Symbol)
    return !isprimitivetype(getfield(@__MODULE__, type)) && !(type in REFERENCE_BLACKLIST)
end
isreferencetype(type::PType) = isreferencetype(type.name)

should_dereference(type::Symbol) = isreferencetype(type) && !(type in DEREFERENCE_BLACKLIST)
should_dereference(type::PType) = should_dereference(type.name)

end


text:
include("AlgAssAbsOrd/Order.jl")
include("AlgAssAbsOrd/Ideal.jl")
include("AlgAssAbsOrd/Elem.jl")
include("AlgAssAbsOrd/PicardGroup.jl")
include("AlgAssAbsOrd/ResidueRingMultGrp.jl")
include("AlgAssAbsOrd/LocallyFreeClassGroup.jl")
include("AlgAssAbsOrd/ICM.jl")
include("AlgAssAbsOrd/NEQ.jl")
include("AlgAssAbsOrd/Eichler.jl")


text:
include("src/SpikingNeuralNetworks.jl")
include("src/units.jl")

using WaspNet

include("src/plot.jl")

unicodeplots()

SNN = SpikingNeuralNetworks
adparam = SNN.ADEXParameter(;a = 6.050246708405076, b = 7.308480222357973,
    cm = 803.1019662706587,
    v_rest= -63.22881649139353,
    tau_m=19.73777028610565,
    tau_w=351.0551915202058,
    v_thresh=-39.232165554444265,
    delta_T=6.37124632135508,
    v_spike = -36.58205819488362,
    v_reset = -59.18792270568965,
    spike_delta = 16.33506432689027)
E = SNN.AD(;N = 1, param)
E.I = [795.57128906]
#{'value': array(460.57128906) * pA}

SNN.monitor(Etest, [:v,:I])
#SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.vecplot(Etest, :v) |> display

Ne = 1000;
Ni = 600#int(200.0*(10.0/8.0))
#E = SNN.IZ(;N = Ne, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8));#, C=1))
#I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2));#, C=1))

param = SNN.IZParameter_more(;a = 0.02, b = 0.2, c = -65.0, d = 8.0, C=100.0, vr=-65.0, k=4.0, vt=-50.0, vPeak=20.0, cellnumber=7)
E = SNN.IZ_more(;N = 1, param)
#E = SNN.IZ_more(;N = 1, param)
E.I = [190.0015]
SNN.monitor(E, [:v,:I])

SNN.sim!([E],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.vecplot(E, :v) |> display

param = SNN.IZParameter_more(;a = 0.02, b = 0.2, c = -65.0, d = 8.0, C=90.0, vr=-65.0, k=1.6, vt=-62.0, vPeak=20.0, cellnumber=6)
E = SNN.IZ_more(;N = Ne, param)

param = SNN.IZParameter_more(;a = 0.02, b = 0.2, c = -65.0, d = 8.0, C=190.0, vr=-75.0, k=1.6, vt=-62.0, vPeak=20.0, cellnumber=5)

I = SNN.IZ_more(;N = Ni, param)
I.I = [160.0015]

SNN.monitor(I, [:v,:I])
SNN.sim!([I],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.vecplot(I, :v) |> display

SNN.vecplot(I, :I) |> display

Etest =SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))#, C=1))
Etest.I = [3.1]
SNN.monitor(Etest, [:v])
SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms)
SNN.vecplot(Etest, :v) |> display

#Etest = SNN.IZ_more(;N = 1, param = SNN.IZParameter_more(;a = 0.02, b = 0.2, c = -65, d = 8, C=1, vr=-70))

EE = SNN.SpikingSynapse(E, E, :v; σ = 0.25,  p = 0.25);
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.25,  p = 0.25);
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0);
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.25);

#EE = SNN.SpikingSynapse(E, E, :v; σ = 1,  p = 1);
#EI = SNN.SpikingSynapse(E, I, :v; σ = 1,  p = 1);
#IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1);
#II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 1);
println(size(E.N),size(E.N))

P = [E, I];
C = [EE, EI, IE, II];
println("passed !!!")

SNN.monitor([E, I], [:fire])
for t = 1:2000
    E.I .= 5randn(Ne)
    I.I .= 2randn(Ni)
    SNN.sim!(P, C, 1ms)
end
SNN.raster(P) #|> display
#SNN.vecplot(E, :fire) |> display

# println(sum(SNN.raster(P)))
#println(SNN.raster(P))

#{'C': 114.25027724782701, 'k': 0.714430107979936, 'vr': -68.3530583578185, 'vt': -49.88514142246343, 'vPeak': 41.99002172075032, 'a': 0.01223913281254824, 'b': -1.9721425473513066, 'c': -49.74145056836865, 'd': 94.78925041952787}
Etest = SNN.IZ(;N = 1, param = SNN.IZParameter(;a =  0.01223913281254824, b = -1.9721425473513066, c = -49.74145056836865, d = 94.78925041952787))#,C=1.1425,vr=-68.3530583578185));
Etest.I = [168.275];
SNN.monitor(Etest, [:v]);
SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms);

SNN.vecplot(Etest, :v)# |> display

#EE = SNN.SpikingSynapse(E, E, :v; σ = 0.5,  p = 0.8);
#EI = SNN.SpikingSynapse(E, I, :v; σ = 0.5,  p = 0.8);

Ne = 1000;
Ni = 600;
Nb = 50
E = SNN.IZ(;N = Ne, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8));#, C=1))
#E = SNN.IZ_more(;N = Ne, param = SNN.IZ_moreParameter(;a =  0.01223913281254824, b = -1.9721425473513066, c = -49.74145056836865, d = 94.78925041952787))#,C=1.1425,vr=-68.3530583578185));

I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2));
EE = SNN.SpikingSynapse(E, E, :v; σ = 0.25,  p = 0.25);
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.25,  p = 0.25);
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0);
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.25);

E_new = SNN.IZ(;N = Nb, param = SNN.IZParameter(;a =  0.01223913281254824, b = -1.9721425473513066, c = -49.74145056836865, d = 94.78925041952787))#,C=1.1425,vr=-68.3530583578185));
#EE_new = SNN.SpikingSynapse(E_new,E_new, :v; σ = 0.25,  p = 0.25);
#P = [E_new,E, I];
#C = [EE_new,EE, EI, IE, II];

P = [E, I];
C = [EE, EI, IE, II]#,EE_new];

SNN.monitor([E, I], [:fire])
for t = 1:1000
    #E_new.I .= 2randn(Nb)
    E.I .= 0.0112482752(Ne)
    I.I .= 1.68275(Ni)
    SNN.sim!(P, C, 1ms)
end
SNN.raster(P) |> display

typeof(E)

#fieldnames(Plot)

Ne = 1000;      Ni = 200#int(200.0*(10.0/8.0))
Etest = SNN.IZ(;N = 1, param = SNN.IZParameter(;a =  0.07976319652352468, b = -1.545038853189859, c = -52.726764975633166, d = 144.76185252557173))#, C=0.9336,vr=-65.55690901865944))
Etest.I = [1100]
SNN.monitor(Etest, [:v])
SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms)
#SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms)

SNN.vecplot(Etest, :v) |> display

#{'C': 93.36810828872358, 'k': 0.7094612672234731, 'vr': -65.55690901865944, 'vt': -49.93900899948406, 'vPeak': 46.969034033727034, 'a': 0.07976319652352468, 'b': -1.545038853189859, 'c': -52.726764975633166, 'd': 144.76185252557173}
E = SNN.IZ(;N = Ne, param = SNN.IZParameter(;a =  0.07976319652352468, b = -1.545038853189859, c = -52.726764975633166, d = 144.76185252557173))#, C=0.9336,vr=-65.55690901865944))
I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
#I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a =  0.07976319652352468, b = -1.545038853189859, c = -52.726764975633166, d = 144.76185252557173))#, C=0.9336,vr=-65.55690901865944))

EE = SNN.SpikingSynapse(E, E, :v; σ = 0.25,  p = 0.5);
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.25,  p = 0.5);
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0);
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.25);
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor([E, I], [:fire])
for t = 1:2000
    E.I .= 5randn(Ne)
    I.I .= 2randn(Ni)
    SNN.sim!(P, C, 1ms)
end
SNN.raster(P) |> display

println("\n olf mit \n")

Etest = SNN.IZ(;N = 1, param = SNN.IZParameter(;a =  0.015190425633023837, b = -1.989511454925138, c = -54.859511795883556, d = 129.7659026061511))#, C=1.99615880961252,vr=-59.78177430465574))
Etest.I = [400]
SNN.monitor(Etest, [:v])
#SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms)
SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=1000ms)

SNN.vecplot(Etest, :v) #|> display
#E = SNN.HH(;N = 1)
#E.I = [0.0752nA]
#SNN.monitor(E, [:v])
#SNN.sim!([E], []; dt = 0.015ms, delay=100ms, stimulus_duration=1000ms, duration = 1300ms)
#SNN.vecplot(E, :v) |> display

Ne = 1000;      Ni = 200#int(200.0*(10.0/8.0))
#{'C': 93.36810828872358, 'k': 0.7094612672234731, 'vr': -65.55690901865944, 'vt': -49.93900899948406, 'vPeak': 46.969034033727034, 'a': 0.07976319652352468, 'b': -1.545038853189859, 'c': -52.726764975633166, 'd': 144.76185252557173}
E = SNN.IZ(;N = Ne, param = SNN.IZParameter(;a =  0.015190425633023837, b = -1.989511454925138, c = -54.859511795883556, d = 129.7659026061511))#, C=1.99615880961252,vr=-59.78177430465574))
I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
#I = SNN.IZ(;N = Ni, param = SNN.IZParameter(;a =  0.015190425633023837, b = -1.989511454925138, c = -54.859511795883556, d = 129.7659026061511))#, C=1.99615880961252,vr=-59.78177430465574))

EE = SNN.SpikingSynapse(E, E, :v; σ = 0.25,  p = 0.5)
EI = SNN.SpikingSynapse(E, I, :v; σ = 0.25,  p = 0.5)
IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0)
II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.25)
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor([E, I], [:fire])
for t = 1:2000
    E.I .= 5randn(Ne)
    I.I .= 2randn(Ni)
    SNN.sim!(P, C, 1ms)
end
SNN.raster(P) |> display

#using Plots
#=
RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

SNN.monitor(P, [:v])
T = 2second
for t = 1:20000
    for p in [RS, IB, CH, FS, LTS]
        p.I = [10]
    end
    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
    RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
    SNN.sim!(P, [], 0.1ms)
end
=#

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


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

JLLWrappers.@generate_main_file_header("cddlib")
JLLWrappers.@generate_main_file("cddlib", UUID("f07e07eb-5685-515a-97c8-3014f6152feb"))
end  # module cddlib_jll


text:
"""
    slerp(p1, p2, t)

The slerp interpolation at `t` between two vectors.
"""
function slerp(p1,p2,t)
    0 <= t <= 1 || throw(DomainError("Interpolation parameter should be between 0 and 1"))
    m = (1-t)*norm(p1) + t*norm(p2)
    p1 = normalize(p1)
    p2 = normalize(p2)
    d = p1'p2
    if d > 0.99
        return m*((1-t)*p1 + t*p2)
    end
    Ω  = acos(d)
    m*((sin((1-t)*Ω)/sin(Ω))*p1 + (sin(t*Ω)/sin(Ω))*p2)
end

# function sinkhorn_interpolate(a1,c1,a2,c2,ai,ci,α)
#     error("Update to new interface")
#     α == 0 && return a1,c1,0
#     α == 1 && return a2,c2,0
#     # ai = (1-α)*a1 + α*a2
#     # ci  = (1-α)*c1 + α*c2
#     distmat = loss_spectral_ot(a1,c1,a2,c2)[2]
#     function loss(r)
#         ar,cr = r[1:length(ai)], r[length(ai)+1:end]
#         ar[1] = 1
#         (1-α)*loss_spectral_ot!(distmat,a1,c1,ar,cr) +  α*loss_spectral_ot!(distmat,a2,c2,ar,cr)
#     end
#     res = Optim.optimize(loss, [ai;ci], BFGS(), Optim.Options(store_trace=false, show_trace=true, iterations=500, allow_f_increases=false, time_limit=20, x_tol=1e-6, f_tol=1e-6), inplace=false, autodiff=:forward)
#
#     ai,ci = res.minimizer[1:length(ai)], res.minimizer[length(ai)+1:end]
#     ai[1] = 1
#     ai, ci, res
# end

function centraldiff(v::AbstractVector)
    dv = diff(v)./2
    a1 = [dv[1];dv]
    a2 = [dv;dv[end]]
    a = a1+a2
end

"""
    $(SIGNATURES)

Perform displacement interpolation between two models.
"""
function interpolator(d::RationalOptimalTransportDistance,A1,A2)
    @assert d.p == 2 "Interpolation only supported for p=2, you have p=$(d.p)"
    interval   = (0., d.interval[2])
    @assert (interval[1] == 0 || interval[1] == -interval[2])
    e1   = sqrt(2)/sqrt(spectralenergy(domain(d), A1))
    e2   = sqrt(2)/sqrt(spectralenergy(domain(d), A2))
    f1    = w -> evalfr(domain(d), magnitude(d), w, A1, e1)
    f2    = w -> evalfr(domain(d), magnitude(d), w, A2, e2)
    sol1  = c∫(f1,interval...)
    sol2  = c∫(f2,interval...)
    σ1    = sol1(interval[2]) # The total energy in the spectrum
    σ2    = sol2(interval[2]) # The total energy in the spectrum
    F1(w) = sol1(w)/σ1
    F2(w) = sol2(w)/σ2
    iF1   = inv(F1, interval)
    iF2   = inv(F2, interval)
    function (w,t)
        t == 0 && return map(f1,w)
        t == 1 && return map(f2,w)
        r = tmap(w) do w
            finv(e->(1-t)*iF1(e) + t*iF2(e), w, (0,max(σ1, σ2)))
        end
        centraldiff(r) ./ centraldiff(w)
    end
end

function interpolator(d::EuclideanRootDistance,A1,A2; normalize=false)
    p = d.p
    @assert p == 2 "Interpolation only supported for p=2, you have p=$p"
    RT = domain(d) isa Continuous ? ContinuousRoots : DiscreteRoots
    e1,e2 = roots(domain(d), A1), roots(domain(d), A2)
    I1,I2 = d.assignment(e1, e2)
    e1,e2 = RT(e1[I1]),RT(e2[I2])
    w1,w2 = d.weight(e1), d.weight(e2)
    function (w,t)
        A = AR(RT(((1-t)*w1.*e1 + t*w2.*e2)./((1-t).*w1.+t.*w2)))
        e = normalize ? 1/sqrt(spectralenergy(domain(d), A)) : 1
        tmap(w) do w
            evalfr(domain(d), magnitude(d), w, A, e)
        end
    end
end


text:
## Data preparation functions

# eBird data preparation (from DataFrame)

"""
See `prepare_ebd_data!` for details.
"""
function prepare_ebd_data(df::DataFrame)
    # Fix names case & spacing
    newnames =
        names(df) .|>
        string .|>
        titlecase .|>
        lowercasefirst .|>
        x -> replace(x, " " => "") .|> Symbol
    df = rename(df, newnames)
    # Rename species column
    df = rename(df, :scientificName => :species)
    # Separate year-month-day
    df.year = year.(df.observationDate)
    df.month = month.(df.observationDate)
    # Remove entries with missing data
    df = dropmissing(df, :species)
    # Replace spaces by underscores in species names
    df.species .= replace.(df.species, " " .=> "_")
    # Remove not approved observations (exotic species, unvetted data)
    df = filter(:approved => ==(1), df)
    return df
end

"""
    prepare_ebd_data!(df::DataFrame)
    prepare_ebd_data(df::DataFrame)

Prepares the eBird data for the analyses from this project. Currently, this
fixes the column and species names, adds a year and a month column, removes
entries with missing data, and removed unapproved observations (exotic species,
unvetted data).

`prepare_ebd_data(df::DataFrame)` copies the DataFrame, while
`prepare_ebd_data!(df::DataFrame)` operates in-place. The in-place version is
preferred when working on the very large eBird DataFrame.
"""
function prepare_ebd_data!(df::DataFrame)
    # Fix names case & spacing
    newnames =
        names(df) .|>
        string .|>
        titlecase .|>
        lowercasefirst .|>
        x -> replace(x, " " => "") .|> Symbol
    rename!(df, newnames)
    # Rename species column
    rename!(df, :scientificName => :species)
    # Separate year-month-day
    df.year = year.(df.observationDate)
    df.month = month.(df.observationDate)
    # Remove entries with missing data
    dropmissing!(df, :species)
    # Replace spaces by underscores in species names
    df.species .= replace.(df.species, " " .=> "_")
    # Remove not approved observations (exotic species, unvetted data)
    filter!(:approved => ==(1), df)
    return df
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
struct FuncWrap <: AbstractComponent
    func::Function
    p::Vector{Parameter}

    function FuncWrap(func::Function, args...)
        params = Vector{Parameter}()
        for i in 1:length(args)
            push!(params, Parameter(args[i]))
        end
        return new(func, params)
    end
end

function evaluate!(buffer, comp::FuncWrap, domain::AbstractDomain,
                   params...)
    buffer .= comp.func(coords(domain)..., params...)
end


text:
module UpdatableQRFactorizations

using LinearAlgebra
using LinearAlgebra: Givens

export AbstractQR, GivensQ, GivensQR, UpdatableGivensQR, UpdatableQR, UQR,
        add_column!, remove_column!

const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}

# IDEA: could have efficient version for sparse matrices
include("abstract.jl")
include("givensQ.jl")
include("givensQR.jl")
include("updatableQR.jl")

end # module


text:
"""
Bridges `CP.DoublyLexicographicallyLessThan` to `CP.LexicographicallyLessThan`.
"""
struct DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T} <: MOIBC.AbstractBridge
    con::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.LexicographicallyLessThan} # Columns are lexicographically sorted.
    con_transposed::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.LexicographicallyLessThan} # Rows are lexicographically sorted.
end

function MOIBC.bridge_constraint(
    ::Type{DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.DoublyLexicographicallyLessThan,
) where {T}
    return MOIBC.bridge_constraint(
        DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.DoublyLexicographicallyLessThan,
) where {T}
    # Constraint on columns.
    con = MOI.add_constraint(
        model,
        f,
        CP.LexicographicallyLessThan(s.row_dim, s.column_dim)
    )

    # Constraint on rows.
    f_scalars = MOIU.scalarize(f)
    f_matrix = reshape(f_scalars, s.row_dim, s.column_dim)
    f_transposed = MOIU.vectorize(vec(f_matrix)) # vec() does the transposition.

    con_transposed = MOI.add_constraint(
        model,
        f_transposed,
        CP.LexicographicallyLessThan(s.column_dim, s.row_dim)
    )

    return DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge(con, con_transposed)
end

function MOI.supports_constraint(
    ::Type{DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.DoublyLexicographicallyLessThan},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.LexicographicallyLessThan),
    ]
end

function MOI.get(
    ::DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.LexicographicallyLessThan,
    },
) where {T}
    return 2
end

function MOI.get(
    b::DoublyLexicographicallyLessThan2LexicographicallyLessThanBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.LexicographicallyLessThan,
    },
) where {T}
    return [b.con, b.con_transposed]
end


text:
using ExprRules
using Debugger

grammar = @grammar begin
	Real = 1 | 2
	Real = Real + Real
	end

node = RuleNode(3, [RuleNode(2), RuleNode(3, [RuleNode(1), RuleNode(2)])])

_next_state!(node, grammar, 3)


text:
using PyCall
@pyimport matplotlib.pyplot as plt
include("Track.jl")

# Create a Track
track = Track(100)
calculate_C!(track)
generate_y!(track)

# Explore the prior
plt.ion()
steps = 1000
for(i in 1:steps)
	# Do a Metropolis proposal (prior is implicit)
	perturb!(track, theta=0.2)

	plt.hold(false)
	plt.plot(track.y, "bo-")
	plt.xlim([-0.5, track.num_points - 0.5])
	plt.ylim([-4.0, 4.0])
	plt.title("Step $i/$steps")
	plt.draw()
end

plt.ioff()
plt.show()

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
#=
Utilities for Chebyshev polynomial series:
multiplication, dilation, and truncation.
=#
module Cheb

using ..Fct

abstract type ChebSeries{T<:Number} end
struct FullChebSeries{T<:Number} <: ChebSeries{T}
    coeff::AbstractArray{T}
end
struct EvenChebSeries{T<:Number} <: ChebSeries{T}
    coeff::AbstractArray{T}
end
struct OddChebSeries{T<:Number} <: ChebSeries{T}
    coeff::AbstractArray{T}
end

FullChebSeries(a::AbstractArray{T}) where T = FullChebSeries{T}(a)
EvenChebSeries(a::AbstractArray{T}) where T = EvenChebSeries{T}(a)
OddChebSeries(a::AbstractArray{T}) where T = OddChebSeries{T}(a)
ChebSeries(a::AbstractArray{T}) where T = FullChebSeries{T}(a)
ChebSeries{T}(a::AbstractArray{T}) where T = FullChebSeries{T}(a)

FullChebSeries(T::Type,deg::Integer) = FullChebSeries{T}(Vector{T}(undef,deg+1))
EvenChebSeries(T::Type,deg::Integer) = EvenChebSeries{T}(Vector{T}(undef,div(deg,2)+1))
OddChebSeries(T::Type,deg::Integer) = OddChebSeries{T}(Vector{T}(undef,div(deg,2)+1))
ChebSeries(T::Type,deg::Integer) = FullChebSeries(T,deg)

EvenChebSeries(c::FullChebSeries{T}) where T = EvenChebSeries{T}(view(c.coeff,1:2:length(c.coeff)))
OddChebSeries(c::FullChebSeries{T}) where T = OddChebSeries{T}(view(c.coeff,2:2:length(c.coeff)))

#external 0-based indexing; internal 1-based indexing
Base.getindex(c::FullChebSeries{T}, n::Integer) where T = c.coeff[1+n]
Base.getindex(c::EvenChebSeries{T}, n::Integer) where T = iseven(n) ? c.coeff[1+div(n,2)] : zero(T)
Base.getindex(c::OddChebSeries{T}, n::Integer) where T = isodd(n) ? c.coeff[1+div(n,2)] : zero(T)
Base.setindex!(c::FullChebSeries{T}, x, n::Integer) where T = (c.coeff[1+n] = T(x))
Base.setindex!(c::EvenChebSeries{T}, x, n::Integer) where T = (c.coeff[1+div(n,2)] = T(x))
Base.setindex!(c::OddChebSeries{T}, x, n::Integer) where T = (c.coeff[1+div(n,2)] = T(x))

#assumes length(c.coeff) >= 1
degree(c::FullChebSeries) = length(c.coeff) - 1
degree(c::EvenChebSeries) = 2*(length(c.coeff) - 1)
degree(c::OddChebSeries) = 2*length(c.coeff) - 1

support(c::FullChebSeries) = 0:degree(c)
support(c::EvenChebSeries) = 0:2:degree(c)
support(c::OddChebSeries) = 1:2:degree(c)

#assumes deg >= degree(c)
truncate(c::FullChebSeries, deg::Integer) = FullChebSeries(view(c.coeff,1:deg+1))
truncate(c::EvenChebSeries, deg::Integer) = EvenChebSeries(view(c.coeff,1:div(deg,2)+1))

#assumes deg >= degree(c)
truncate(c::OddChebSeries, deg::Integer) = OddChebSeries(view(c.coeff,1:div(deg-1,2)+1))

Base.eltype(c::ChebSeries) = eltype(c.coeff)

FullChebSeries(f::Function) = autotruncate(FullChebSeries(autocheb(f)))
ChebSeries(f::Function) = FullChebSeries(f)
EvenChebSeries(f::Function) = EvenChebSeries(FullChebSeries(f))
OddChebSeries(f::Function) = OddChebSeries(FullChebSeries(f))

import Base: *
function (*)(a::ChebSeries, b::ChebSeries)::FullChebSeries
    L,M=degree(a),degree(b)
    T=promote_type(eltype(a),eltype(b))
    c=FullChebSeries(T,L+M)
    mult!(a,L,b,M,c)
end
function (*)(a::EvenChebSeries, b::EvenChebSeries)::EvenChebSeries
    L,M=degree(a),degree(b)
    T=promote_type(eltype(a),eltype(b))
    c=EvenChebSeries(T,L+M)
    mult!(a,L,b,M,c)
end
function (*)(a::OddChebSeries, b::OddChebSeries)::EvenChebSeries
    L,M=degree(a),degree(b)
    T=promote_type(eltype(a),eltype(b))
    c=EvenChebSeries(T,L+M)
    mult!(a,L,b,M,c)
end
function (*)(a::OddChebSeries, b::EvenChebSeries)::OddChebSeries
    L,M=degree(a),degree(b)
    T=promote_type(eltype(a),eltype(b))
    c=OddChebSeries(T,L+M)
    mult!(a,L,b,M,c)
end
function mult!(a::ChebSeries, L, b::ChebSeries, M, c::ChebSeries)
    T = eltype(c)
    half, init = T(0.5), zero(T)

    for n in support(c)
        #2T_i*T_j==T_{i+j}+T_{|i-j|}
        c_n = init

        itr = intersect(max(0,n-M):min(L,n), support(a))
        length(itr)>0 && (c_n += sum(a[l]*b[n-l] for l in itr))

        itr = intersect(min(L,M+n):-1:n, support(a))
        length(itr)>0 && (c_n += sum(a[l]*b[l-n] for l in itr))

        if n > 0
            itr = intersect(min(L,M-n):-1:0, support(a))
            length(itr)>0 && (c_n += sum(a[l]*b[l+n] for l in itr))
        end

        c[n] = half * c_n
    end
    c
end
(*)(a::EvenChebSeries, b::OddChebSeries)::OddChebSeries = b * a
(*)(a::FullChebSeries, b::EvenChebSeries)::FullChebSeries = b * a
(*)(a::FullChebSeries, b::OddChebSeries)::FullChebSeries = b * a

Base.conj(c::ChebSeries) = typeof(c)(conj(c.coeff))
Base.conj!(c::ChebSeries) = typeof(c)(conj!(c.coeff))

#Outputs o such that sum_n o[n]T_n(x)=sum_n c[n]T_n(kappa*x)
function dilate(c::ChebSeries,kappa::Number)::ChebSeries
    T=eltype(c)
    o=typeof(c)(zeros(eltype(c),size(c.coeff)))
    d=zeros(T,3,div(degree(c),2)+1)
    n0,n1,n2=3,2,1
    d[n0,1]=1
    0 in support(c) && (o[0]=c[0])
    n0,n1,n2 = n2,n0,n1
    d[n0,1]=kappa
    1 in support(c) && (o[1]=c[1]*kappa)
    for n in 2:degree(c)
        n0,n1,n2 = n2,n0,n1
        d[n0,1] = d[n1,1]*kappa
        m0 = div(n,2)
        for m in 1:m0-1
            d[n0,1+m] = (d[n1,1+m] + d[n1,m])*kappa - d[n2,m]
        end
        if isodd(n)
            d[n0,1+m0] = (2d[n1,1+m0] + d[n1,m0])*kappa - d[n2,m0]
        else
            d[n0,1+m0] = d[n1,m0]*kappa - d[n2,m0]
        end
        if n in support(c)
            cn = c[n]
            for m in 0:m0
                o[n-2m] += cn*d[n0,1+m]
            end
        end
    end
    o
end

#evaluate sum_n c[n+1]*T_n(x)
function chebinterp(c::ChebSeries,x::Number)
    old,new=1,x
    out=c[0]+c[1]*x
    for n in 2 : degree(c)
        old, new = new, 2x*new-old
        out+=c[n]*new
    end
    out
end

#=
Outputs coefficient vector c such that f(x) is interpolated
and well approximated by sum_n c[n+1] T_n(x).
=#
function autocheb(f;log2N0::Integer=1,T=Complex{Float64})
    N = 2^max(log2N0,1)
    while true
        samples = Fct.chebsample(f,2N;T=T) #length 2N+1
        coeffs = Fct.chebcoeffs(samples,T)
        head = sum(abs2(coeffs[k]) for k in N+1:-1:1)
        tail = sum(abs2(coeffs[k]) for k in 2N+1:-1:N+2)
        rtol = eps()*(N<=64 ? N<=32 ? 1 : 4 : N<=128 ? 16 : N)
        if tail > head * rtol^2 && N <= 2^19
            samples = Fct.doublesample(f,samples)
            N<<=1
        else
            return view(coeffs,1:min(N+1,2N+1))
        end
    end
end

function autotruncate(c::ChebSeries)
    base = maximum((abs2(c[n]) for n in 0:degree(c)))
    base *= eps()^2
    n = degree(c)
    while abs2(c[n]) < base * n
        n -= 1
    end
    truncate(c,n)
end

end #module


text:
"""
    label = label_components(A; [dims=coords_spatial(A)], [r=1], [bkg])
    label = label_components(A, se; [bkg])

Find and label the connected components of array `A` where the connectivity is defined by
structuring element `se`. Each component is assigned a unique integer value as its label
with `0` representing the background defined by `bkg`.

$(_docstring_se)

# Examples

```jldoctest; setup=:(using ImageMorphology)
julia> A = [false true false true  false;
            true false false  true  true]
2×5 Matrix{Bool}:
 0  1  0  1  0
 1  0  0  1  1

julia> label_components(A) # default diamond shape C4 connectivity
2×5 Matrix{$Int}:
 0  2  0  3  0
 1  0  0  3  3

julia> label_components(A; dims=2) # only the rows are considered
2×5 Matrix{$Int}:
 0  2  0  3  0
 1  0  0  4  4

julia> label_components(A, strel_box((3, 3))) # box shape C8 connectivity
2×5 Matrix{$Int}:
 0  1  0  2  0
 1  0  0  2  2
```

The in-place version is [`label_components!`](@ref). See also [`component_boxes`](@ref),
[`component_lengths`](@ref), [`component_indices`](@ref), [`component_centroids`](@ref) for
basic properties of the labeled components.
"""
label_components(A; dims=coords_spatial(A), r=1, kwargs...) = label_components(A, strel_diamond(A, dims; r); kwargs...)
label_components(A, se; kwargs...) = label_components!(similar(A, Int), A, se; kwargs...)

"""
    label_components!(out, A; [dims], [r] [bkg])
    label_components!(out, A, se; [bkg])

The in-place version of [`label_components`](@ref).
"""
function label_components!(out, A; dims=coords_spatial(A), r=1, kwargs...)
    if !(axes(out) == axes(A))
        msg = "axes of input and output must match, got $(axes(out)) and $(axes(A))"
        hint = if eltype(A) <: Bool || eltype(A) <: CartesianIndex
            "The second argument seems to be a structuring element, it is expected to be the input array."
        else
            ""
        end
        throw(DimensionMismatch("$msg. $hint"))
    end
    return label_components!(out, A, strel_diamond(A, dims; r); kwargs...)
end

function label_components!(out::AbstractArray{T}, A::AbstractArray, se; bkg=zero(eltype(A))) where {T<:Integer}
    axes(out) == axes(A) || throw_dmm(axes(out), axes(A))
    se = _maybe_build_symmetric_strel(se) # compat patch
    is_symmetric(se) || throw(ArgumentError("Non-symmetric structuring element is not supported yet"))
    upper_se, _ = strel_split(CartesianIndex, se)
    fill!(out, zero(T))
    sets = DisjointMinSets{T}()
    sizehint!(sets.parents, floor(Int, sqrt(length(A))))
    @inbounds for i in CartesianIndices(A)
        val = A[i]
        val == bkg && continue
        label = typemax(T)    # sentinel value
        for Δi in upper_se
            ii = i + Δi
            checkbounds(Bool, A, ii) || continue
            if A[ii] == val
                newlabel = out[ii]
                label = if ((label == typemax(T)) | (label == newlabel))
                    newlabel
                else
                    union!(sets, label, newlabel)
                end
            end
        end
        if label == typemax(T)
            label = push!(sets)
        end
        out[i] = label
    end
    # Now parse sets to find the labels
    newlabel = minlabel(sets)
    @inbounds for i in eachindex(A, out)
        if A[i] != bkg
            out[i] = newlabel[find_root!(sets, out[i])]
        end
    end
    return out
end

function _maybe_build_symmetric_strel(iter)
    # NOTE(johnnychen94): label_components! (<v0.4) used to accept only the splited upper
    # half part as `se`. We should still support it with a deprecation warning.
    new_se = strel(CartesianIndex, strel(Bool, iter)) # remove the center position
    sym_se = union(.-new_se, new_se)
    if length(sym_se) == 2length(new_se)
        msg =
            "`label_components!` now requires a complete symmetric structuring element," *
            " but the provided structuring element is only the upper half." *
            " This is a deprecated behavior and will be removed in a future version."
        Base.depwarn(msg, :label_components!)
        return sym_se
    else
        return iter
    end
end

function throw_dmm(ax1, ax2)
    throw(DimensionMismatch("axes of input and output must match, got $ax1 and $ax2"))
end

# Copied directly from DataStructures.jl, but specialized
# to always make the parent be the smallest label
struct DisjointMinSets{T<:Integer}
    parents::Vector{T}

    DisjointMinSets{T}(n::Integer) where {T} = new([T(1):T(n);])
end
DisjointMinSets{T}() where {T<:Integer} = DisjointMinSets{T}(T(0))
DisjointMinSets() = DisjointMinSets{INt}()

Base.@propagate_inbounds function find_root!(sets::DisjointMinSets, m::Integer)
    p = sets.parents[m]   # don't use @inbounds here, it might not be safe
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    return p
end

# an unsafe variant of the above
function find_root_unsafe!(sets::DisjointMinSets, m::Int)
    @inbounds p = sets.parents[m]
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    return p
end

Base.@propagate_inbounds function union!(sets::DisjointMinSets, m::Integer, n::Integer)
    mp = find_root!(sets, m)
    np = find_root!(sets, n)
    if mp < np
        sets.parents[np] = mp
        return mp
    elseif np < mp
        sets.parents[mp] = np
        return np
    end
    return mp
end

function Base.push!(sets::DisjointMinSets)
    m = length(sets.parents) + 1
    m >= typemax(eltype(sets.parents)) && error("labels exhausted, use a larger integer type")
    push!(sets.parents, m)
    return m
end

function minlabel(sets::DisjointMinSets)
    out = Vector{Int}(undef, length(sets.parents))
    k = 0
    @inbounds for i in 1:length(sets.parents)
        if sets.parents[i] == i
            k += 1
        end
        out[i] = k
    end
    return out
end

"""
    boxes = component_boxes(labeled_array)

Calculates the minimal bounding boxes for each label including the background label. The
labels can be computed by [`label_components`](@ref).

Each bounding box is represented as a `CartesianIndices`. `boxes` is shifted to 0-based
indexing vector so that background region is `boxes[0]`.

```jldoctest; setup=:(using ImageMorphology)
julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{$Int}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{$Int}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> boxes = component_boxes(label) # get bounding boxes of all regions
5-element OffsetArray(::Vector{CartesianIndices{2, Tuple{UnitRange{$Int}, UnitRange{$Int}}}}, 0:4) with eltype CartesianIndices{2, Tuple{UnitRange{$Int}, UnitRange{$Int}}} with indices 0:4:
 CartesianIndices((2:5, 2:4))
 CartesianIndices((1:1, 1:5))
 CartesianIndices((2:5, 1:3))
 CartesianIndices((3:5, 3:5))
 CartesianIndices((2:3, 4:5))

julia> A[boxes[1]] # crop the image region with label 1
1×5 Matrix{$Int}:
 2  2  2  2  2

julia> A[boxes[4]] # crop the image region with label 4
2×2 Matrix{$Int}:
 0  1
 1  1
```
"""
function component_boxes(A::AbstractArray{T,N}) where {T<:Integer,N}
    mn, mx = extrema(A)
    if !(mn == 0 || mn == 1)
        throw(ArgumentError("The input labeled array should contain background label `0` as the minimum value"))
    end
    boxes = if mn == 1
        OffsetArray(Matrix{CartesianIndex{N}}(undef, mx, 2), 0, 0)
    elseif mn == 0
        OffsetArray(Matrix{CartesianIndex{N}}(undef, 1 + mx, 2), -1, 0)
    end
    R = CartesianIndices(A)
    boxes[:, 1] .= Ref(last(R))
    boxes[:, 2] .= Ref(first(R))

    @inbounds for i in R
        label_idx = A[i]
        # this actually does a few unnecessary comparison, but it's relatively harmless as
        # the main computation is limited by memory
        boxes[label_idx, 1] = min(i, boxes[label_idx, 1])
        boxes[label_idx, 2] = max(i, boxes[label_idx, 2])
    end

    # convert to CartesianIndices
    return boxes = [boxes[i, 1]:boxes[i, 2] for i in axes(boxes, 1)]
end

"""
    counts = component_lengths(labeled_array)

Count the number of each labels in the input labeled array. `counts` is shifted to 0-based
indexing vector so that the number of background pixels is `counts[0]`.

```jldoctest; setup=:(using ImageMorphology)
julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{$Int}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{$Int}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> component_lengths(label)
5-element OffsetArray(::Vector{$Int}, 0:4) with eltype $Int with indices 0:4:
 3
 5
 7
 7
 3
```

For gray images, labels can be computed by [`label_components`](@ref).
"""
function component_lengths(A::AbstractArray{T}) where {T<:Integer}
    mn, mx = extrema(A)
    if !(mn == 0 || mn == 1)
        throw(ArgumentError("The input labeled array should contain background label `0` as the minimum value"))
    end
    counts = zeros(Int, 0:mx)
    for i in eachindex(A)
        counts[A[i]] += 1
    end
    return counts
end

"""
    indices = component_indices([T], labeled_array)

Get the indices of each label in the input labeled array. `indices` is shifted to 0-based
indexing vector so that the indices of background pixels is `indices[0]`.

The optional type `T` can be either `Int`/`IndexLinear()` or
`CartesianIndex`/`IndexCartesian()` that is used to specify the type of the indices. The
default choice is `IndexStyle(labeled_array)`.

```jldoctest; setup=:(using ImageMorphology)
julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{$Int}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{$Int}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> indices = component_indices(label)
5-element OffsetArray(::Vector{Vector{$Int}}, 0:4) with eltype Vector{$Int} with indices 0:4:
 [8, 10, 17]
 [1, 6, 11, 16, 21]
 [2, 3, 4, 5, 7, 9, 12]
 [13, 14, 15, 19, 20, 24, 25]
 [18, 22, 23]

julia> indices = component_indices(CartesianIndex, label)
5-element OffsetArray(::Vector{Vector{CartesianIndex{2}}}, 0:4) with eltype Vector{CartesianIndex{2}} with indices 0:4:
 [CartesianIndex(3, 2), CartesianIndex(5, 2), CartesianIndex(2, 4)]
 [CartesianIndex(1, 1), CartesianIndex(1, 2), CartesianIndex(1, 3), CartesianIndex(1, 4), CartesianIndex(1, 5)]
 [CartesianIndex(2, 1), CartesianIndex(3, 1), CartesianIndex(4, 1), CartesianIndex(5, 1), CartesianIndex(2, 2), CartesianIndex(4, 2), CartesianIndex(2, 3)]
 [CartesianIndex(3, 3), CartesianIndex(4, 3), CartesianIndex(5, 3), CartesianIndex(4, 4), CartesianIndex(5, 4), CartesianIndex(4, 5), CartesianIndex(5, 5)]
 [CartesianIndex(3, 4), CartesianIndex(2, 5), CartesianIndex(3, 5)]
```

For gray images, labels can be computed by [`label_components`](@ref).
"""
component_indices(A) = component_indices(IndexStyle(A), A)
component_indices(::Type{T}, A) where {T<:Integer} = component_indices(IndexLinear(), A)
component_indices(::Type{T}, A) where {T<:CartesianIndex} = component_indices(IndexCartesian(), A)
component_indices(::IndexCartesian, A) = _component_indices(A, CartesianIndices(A))
component_indices(::IndexLinear, A) = _component_indices(A, LinearIndices(A))

function _component_indices(A::AbstractArray{<:Integer,N}, R::AbstractArray{T,N}) where {T,N}
    mn, mx = extrema(A)
    if !(mn == 0 || mn == 1)
        throw(ArgumentError("The input labeled array should contain background label `0` as the minimum value"))
    end
    indices = OffsetVector([T[] for _ in 0:mx], -1)
    @inbounds for i in R
        push!(indices[A[i]], i)
    end
    return indices
end

"""
    centroids = component_centroids(labeled_array)

Compute the centroid of each label in the input labeled array. `centroids` is shifted to
0-based indexing vector so that the centroid of background pixels is `centroids[0]`.

The centroid of a finite set `X`, also known as geometric center, is calculated using
`sum(X)/length(X)`. For label `i`, all (Cartesian) indices of pixels with label `i` are used
to build the set `X`

```jldoctest; setup=:(using ImageMorphology)
julia> A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
5×5 Matrix{$Int}:
 2  2  2  2  2
 1  1  1  0  1
 1  0  2  1  1
 1  1  2  2  2
 1  0  2  2  2

julia> label = label_components(A) # four disjoint components
5×5 Matrix{$Int}:
 1  1  1  1  1
 2  2  2  0  4
 2  0  3  4  4
 2  2  3  3  3
 2  0  3  3  3

julia> component_centroids(label)
5-element OffsetArray(::Vector{Tuple{Float64, Float64}}, 0:4) with eltype Tuple{Float64, Float64} with indices 0:4:
 (3.3333333333333335, 2.6666666666666665)
 (1.0, 3.0)
 (3.142857142857143, 1.5714285714285714)
 (4.285714285714286, 3.857142857142857)
 (2.6666666666666665, 4.666666666666667)
```

For gray images, labels can be computed by [`label_components`](@ref).
"""
function component_centroids(A::AbstractArray{<:Integer,N}) where {N}
    mn, mx = extrema(A)
    if !(mn == 0 || mn == 1)
        throw(ArgumentError("The input labeled array should contain background label `0` as the minimum value"))
    end
    n = fill(zero(CartesianIndex{N}), mn:mx)
    counts = fill(0, mn:mx)
    @inbounds for I in CartesianIndices(size(A))
        v = A[I]
        n[v] += I
        counts[v] += 1
    end
    return map((v, l) -> v.I ./ l, n, counts)
end

# deprecated
"`component_subscripts(labeled_array)` -> an array of pixels for each label, including the background label 0"
function component_subscripts(img::AbstractArray{Int})
    Base.depwarn("`component_subscripts` is deprecated. Use `component_indices(CartesianIndex, A)` instead.", :component_subscripts)
    n = [Tuple[] for i in 0:maximum(img)]
    s = CartesianIndices(size(img))
    for i in 1:length(img)
        push!(n[img[i] + 1], s[i])
    end
    return n
end


text:
include("Flux.jl")
include("CUDA.jl")
include("Zygote.jl")
include("ReinforcementLearningBase.jl")
include("ElasticArrays.jl")
include("Distributions.jl")


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
grp = MPI.Comm_group(comm)
@test MPI.Comm_size(comm) == MPI.Group_size(grp)
@test MPI.Comm_rank(comm) == MPI.Group_rank(grp)
@test MPI.Group_compare(grp, grp) == MPI.IDENT
grp2 = MPI.Group_union(grp, MPI.GROUP_EMPTY)
@test MPI.Group_compare(grp, grp2) == MPI.IDENT
grp3 = MPI.Group_difference(grp, grp2)
@test MPI.Group_compare(grp3, MPI.GROUP_EMPTY) == MPI.IDENT
grp4 = MPI.Group_intersection(grp, grp2)
@test MPI.Group_compare(grp4, MPI.GROUP_EMPTY) == MPI.UNEQUAL
grp5 = MPI.Group_excl(grp, Int32[0])
@test MPI.Group_size(grp5) == MPI.Group_size(grp) - 1
grp6 = MPI.Group_incl(grp, Int32[0])
@test MPI.Group_size(grp6) == 1

# Don't free the other groups
MPI.free(grp)

MPI.Finalize()
@test MPI.Finalized()


text:
using Documenter
using GeoRegions

makedocs(
    modules  = [GeoRegions],
    doctest  = false,
    format   = Documenter.HTML(
        collapselevel = 1,
        prettyurls    = false
    ),
    authors  = "Nathanael Wong",
    sitename = "GeoRegions.jl",
    pages    = [
        "Home"       => "index.md",
        "GeoRegion Basics" => [
            "What is a GeoRegion?"  => "georegions/intro.md",
            "Creating GeoRegions"   => "georegions/create.md",
            "Predefined GeoRegions" => "georegions/predefined.md",
            "Retrieving GeoRegions" => "georegions/read.md",
        ],
        "Using GeoRegions" => [
            "Is it in a GeoRegion?" => "using/isin.md",
            "RegionInfo Basics"     => "using/regiongrid.md",
            "Custom GeoRegions"     => "using/custom.md",
        ],
        "An Example with ERA5" => "example/era5.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaClimate/GeoRegions.jl.git",
    devbranch = "main"
)


text:
__precompile__()

module FNVHash

export fnv32, fnv16

if VERSION < v"0.6.0" const xor = $ end # FIXME

# References:
#   https://en.wikipedia.org/wiki/Fowler–Noll–Vo_hash_function
#   http://tools.ietf.org/html/draft-eastlake-fnv-03

# FNV1
#   http://www.isthe.com/chongo/tech/comp/fnv/#FNV-1
function fnv32(data::Vector{UInt8})
    hash = 0x811C9DC5
    for c in data
        hash = xor((hash * 0x01000193), c)
    end
    return hash
end

fnv32(s) = fnv32(Vector{UInt8}(s))

fnv16(hash::UInt32) = UInt16(hash >>> 16) $ UInt16(hash & 0xFFFF)
fnv16(data) = fnv16(fnv32(data))

# FNV1a
#   http://www.isthe.com/chongo/tech/comp/fnv/#FNV-1a
function fnv1a32(data::Vector{UInt8})
    hash = 0x811C9DC5
    for c in data
        hash = xor(hash, c) * 0x01000193
    end
    return hash
end

fnv1a32(s) = fnv1a32(Vector{UInt8}(s))

fnv1a16(hash::UInt32) = UInt16(hash >>> 16) $ UInt16(hash) & 0xFFFF
fnv1a16(data) = fnv1a16(fnv1a32(data))

end # module FNVHash


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module ForwardNets

using LightGraphs

export
    ForwardNet,
    Node,
    Layer,
    Variable,
    Activation,
    Affine,
    LSTM,
    Conv2d,
    BatchNorm1,
    BatchNorm3,
    Concatenator,
    Reshaper,
    ReLU,
    TanH,
    SoftPlus,
    ForwardPass,

    read_binary_vec,
    infer_shape,
    convert_to_column_major_array,

    forward!,
    add_node!,
    restore!,
    indexof,
    lastindex,
    get_name,
    get_output,

    sigmoid,
    sigmoid!,
    tanh!,

    zero!,

    calc_forwardpass,
    print_forward_pass

include("read_binary_vec.jl")

abstract Node

const DEFAULT_NAME = :unnamed
function forward!(n::Node)
    # do nothing by default
end
get_name(n::Node) = DEFAULT_NAME

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

type ForwardNet
    dag::DiGraph # graph over nodes
    nodes::Vector{Node}
    name_to_index::Dict{Symbol,Int} # Symbol → index in dag and nodes (note: not all Nodes have names)

    ForwardNet() = new(DiGraph(0), Node[], Dict{Symbol,Int}())
end

Base.getindex(net::ForwardNet, index::Int) = net.nodes[index]
Base.getindex(net::ForwardNet, name::Symbol) = net.nodes[net.name_to_index[name]]
indexof(net::ForwardNet, name::Symbol) = net.name_to_index[name]
lastindex(net::ForwardNet) = nv(net.dag)
function add_node!(net::ForwardNet, node::Node, parents::Vector{Int}=Int[])
    add_vertex!(net.dag)
    push!(net.nodes, node)
    i = length(net.nodes)

    name = get_name(node)
    if name != DEFAULT_NAME
        @assert(!haskey(net.name_to_index, name))
        net.name_to_index[name] = i
    end

    for parent in parents
        add_edge!(net.dag, parent, i)
    end

    net
end
add_node!(net::ForwardNet, node::Node, parent::Int) = add_node!(net, node, Int[parent])

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

abstract Layer <: Node
abstract Activation <: Node

type Variable{n} <: Node
    name::Symbol
    tensor::Array{Float32, n}
end
get_name(a::Variable) = a.name
get_output(a::Variable) = a.tensor

type Affine <: Layer
    name::Symbol
    W::Matrix{Float32} # o×i
    b::Vector{Float32} # o

    parent::Vector{Float32} # i
    child::Vector{Float32}  # o
end
get_name(a::Affine) = a.name
get_output(a::Affine) = a.child
function forward!(a::Affine)
    copy!(a.child, a.b) # y ← b
    Base.LinAlg.BLAS.gemv!('N', 1.0f0, a.W, a.parent, 1.0f0, a.child) # y ← W*x + y
    a
end
function add_node!(net::ForwardNet, ::Type{Affine},
    name::Symbol,
    parent_index::Int,
    output_dim::Int,
    )

    parent_node = net[parent_index]

    parent = get_output(parent_node)::Vector{Float32}

    W = Array(Float32, output_dim, length(parent))
    b = Array(Float32, output_dim)
    child = Array(Float32, output_dim)

    node = Affine(name, W, b, parent, child)
    add_node!(net, node, parent_index)
end
function restore!(a::Affine, filename_W::AbstractString, filename_b::AbstractString)
    input_dim = length(a.parent)

    vec = open(read_binary_vec, filename_W)
    shape = infer_shape(vec, (-1, input_dim))
    a.W[:] = convert_to_column_major_array(vec, (shape[2], shape[1]))'

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end

sigmoid(x::Real) = 1 / (1 + exp(-x))
function sigmoid{F<:Real}(x::Vector{F})
    retval = deepcopy(x)
    for i in 1 : length(x)
        retval[i] = sigmoid(x[i])
    end
    retval
end
function sigmoid!{F<:Real}(x::Vector{F})
    for i in 1 : length(x)
        x[i] = sigmoid(x[i])
    end
    x
end
function tanh!{F<:Real}(x::Vector{F})
    for i in 1 : length(x)
        x[i] = tanh(x[i])
    end
    x
end

type LSTM <: Layer
    name::Symbol
    W::Matrix{Float32} # 4H×(D+H)
        # consists of Wx and Wh
        #  Wx: Input-to-hidden weights, of shape (4H, D)
        #  Wh: Hidden-to-hidden weights, of shape (4H, H)
    b::Vector{Float32} # 4H
    forget_bias::Float32

    parent::Vector{Float32} # input, [D]
    state::Vector{Float32} # state (2H), contains [cell state, hidden state]

    # preallocated memory
    c::Vector{Float32} # [H] cell state
    m::Vector{Float32} # [H] member state
    cell_inputs::Vector{Float32}  # [D+H]
    lstm_mult_res::Vector{Float32} # [4H]
    i::Vector{Float32} # [H] input gate
    j::Vector{Float32} # [H] new input
    f::Vector{Float32} # [H] forget gate
    o::Vector{Float32} # [H] output gate
end
get_name(a::LSTM) = a.name
get_output(a::LSTM) = a.m
zero!(a::LSTM) = fill!(a.state, 0.0f0)
function add_node!(net::ForwardNet, ::Type{LSTM},
    name::Symbol,
    parent_index::Int,
    H::Int; # hidden layer size
    forget_bias::Float32 = 1.0f0
    )

    parent = get_output(net[parent_index])::Vector{Float32}
    D = length(parent)

    W = Array(Float32, 4H, D+H)
    b = Array(Float32, 4H)
    state = Array(Float32, 2H)

    c = Array(Float32, H)
    m = Array(Float32, H)
    cell_inputs = Array(Float32, D+H)
    lstm_mult_res = Array(Float32, 4H)
    i = Array(Float32, H)
    j = Array(Float32, H)
    f = Array(Float32, H)
    o = Array(Float32, H)

    node = LSTM(name, W, b, forget_bias, parent, state, c, m, cell_inputs, lstm_mult_res, i, j, f, o)
    add_node!(net, node, parent_index)
end
function forward!(a::LSTM)

    H = length(a.c)
    D = length(a.parent)

    copy!(a.c, 1, a.state, 1, H)
    copy!(a.m, 1, a.state, H+1, H)

    copy!(a.cell_inputs, 1, a.parent, 1, D)
    copy!(a.cell_inputs, D+1, a.m, 1, H)

    # lstm_matrix = a.W * a.cell_inputs + a.b
    copy!(a.lstm_mult_res, a.b) # y ← b
    Base.LinAlg.BLAS.gemv!('N', 1.0f0, a.W, a.cell_inputs, 1.0f0, a.lstm_mult_res) # y ← W*x + y

    copy!(a.i, 1, a.lstm_mult_res,    1, H)
    copy!(a.j, 1, a.lstm_mult_res,  H+1, H)
    copy!(a.f, 1, a.lstm_mult_res, 2H+1, H)
    copy!(a.o, 1, a.lstm_mult_res, 3H+1, H)

    # compute c_next and m_next
    for k in 1 : H
        a.f[k] += a.forget_bias
        c = sigmoid(a.f[k])*a.c[k] + sigmoid(a.i[k])*tanh(a.j[k]) # c_next
        a.state[k] = a.c[k] = c
        a.state[k+H] = a.m[k] = sigmoid(a.o[k])*tanh(c) # m_next
    end

    a
end
function restore!(a::LSTM, filename_W::AbstractString, filename_b::AbstractString)
    H = length(a.c)
    D = length(a.parent)

    vec = open(read_binary_vec, filename_W)
    a.W[:] = convert_to_column_major_array(vec, (D+H,4H))'

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end

type Conv2d <: Layer
    name::Symbol
    W::Array{Float32, 4} # filter_h × filter_w × in_c × out_c
    b::Vector{Float32} # out_c
    strides::Tuple{Int, Int} # width, height

    parent::Array{Float32, 3} # h × w × in_c
    child::Array{Float32, 3}  # out_h × out_w × out_c
end
get_name(a::Conv2d) = a.name
get_output(a::Conv2d) = a.child
function add_node!(net::ForwardNet, ::Type{Conv2d},
    name::Symbol,
    parent_index::Int,
    strides::Tuple{Int, Int},
    filter_h::Int,
    filter_w::Int,
    out_c::Int)

    parent = get_output(net[parent_index])::Array{Float32, 3}
    h, w, in_c = size(parent)

    W = Array(Float32, filter_h, filter_w, in_c, out_c)
    b = Array(Float32, out_c)

    # FOR VALID: (padding is always zero)
    # out_w = ceil(float(w - filter_w + 1) / float(strides[2]))
    # out_h = ceil(float(h - filter_h + 1) / float(strides[1]))

    # FOR SAME:
    out_w  =ceil(Int, w / strides[2])
    out_h = ceil(Int, h / strides[1])
    child = Array(Float32, out_h, out_w, out_c)

    node = Conv2d(name, W, b, strides, parent, child)
    add_node!(net, node, parent_index)
end
function forward!(a::Conv2d)

    h_in = size(a.parent, 1)
    h_out = size(a.child, 1)
    w_in = size(a.parent, 2)
    w_out = size(a.child, 2)
    hh, ww, f_in, f_out = size(a.W)

    s₁, s₂ = a.strides
    pad_h = (s₁ * (h_out - 1) + hh - h_in) / 2
    pad_w = (s₂ * (w_out - 1) + ww - w_in) / 2

    pad_h_lo = floor(Int, pad_h)
    pad_h_hi = ceil(Int, pad_h)
    pad_w_lo = floor(Int, pad_w)
    pad_w_hi = ceil(Int, pad_w)

    # println("parent_size: ", size(a.parent))
    # println("child_size:  ", size(a.child))

    # println("padding: ", pad_h_lo, "  ", pad_h_hi)
    # println("         ", pad_w_lo, "  ", pad_w_hi)

    i_out = 1
    for i_out in 1 : size(a.child, 1)
        i_in_lo = (i_out-1)*s₁ - pad_h_lo + 1
        i_in_hi = i_in_lo + hh - 1

        h_in_lo = 1
        if i_in_lo < 1
            h_in_lo += 1-i_in_lo
            i_in_lo = 1
        end
        h_in_hi = hh
        if i_in_hi > h_in
            h_in_hi -= i_in_hi - h_in
            i_in_hi = h_in
        end

        for j_out in 1 : size(a.child, 2)

            j_in_lo = (j_out-1)*s₂ - pad_w_lo + 1
            j_in_hi = j_in_lo + ww - 1

            # println("j_in_lo: ", j_in_lo)
            # println("j_in_hi: ", j_in_hi)

            w_in_lo = 1
            if j_in_lo < 1
                w_in_lo += 1-j_in_lo
                j_in_lo = 1
            end
            w_in_hi = ww
            if j_in_hi > w_in
                w_in_hi -= j_in_hi - w_in
                j_in_hi = w_in
            end

            # println("j_in_lo: ", j_in_lo)
            # println("j_in_hi: ", j_in_hi)
            # println("w_in_lo: ", w_in_lo)
            # println("w_in_hi: ", w_in_hi)

            for k_out in 1 : f_out

                # println("parent: ", i_in_lo:i_in_hi, ", ", j_in_lo:j_in_hi, ", ", 1:f_in)
                # println("W:      ", h_in_lo:h_in_hi, ", ", w_in_lo,w_in_hi, ", ", 1:f_in, ", ", k_out)

                a.child[i_out,j_out,k_out] = a.b[k_out]

                for (i, i_in) in enumerate(i_in_lo:i_in_hi)
                    h₁ = h_in_lo + i - 1
                    for (j, j_in) in enumerate(j_in_lo:j_in_hi)
                        w₁ = w_in_lo + j - 1
                        for k_in in 1 : f_in
                            W_contrib = a.W[h₁,w₁,k_in,k_out]
                            p_contrib = a.parent[i_in,j_in,k_in]
                            a.child[i_out,j_out,k_out] += W_contrib*p_contrib
                        end
                    end
                end
            end
        end
    end

    a
end
function restore!(a::Conv2d, filename_W::AbstractString, filename_b::AbstractString)

    filter_h, filter_w, in_c, out_c = size(a.W)

    vec = open(read_binary_vec, filename_W)
    a.W[:] = convert_to_column_major_array(vec, (filter_h, filter_w, in_c, out_c))

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end

type BatchNorm1 <: Layer
    # https://arxiv.org/pdf/1502.03167.pdf
    name::Symbol
    γ::Vector{Float32}
    β::Vector{Float32}
    ϵ::Float32
    μ::Vector{Float32} # train set mean
    ν::Vector{Float32} # train set variance

    parent::Vector{Float32}
    child::Vector{Float32}
end
get_name(a::BatchNorm1) = a.name
get_output(a::BatchNorm1) = a.child
function forward!(a::BatchNorm1)
    x = (a.parent - a.μ) ./ sqrt(a.ν + a.ϵ)
    a.child[:] = (a.γ .* x) .+ a.β
    a
end
function add_node!(net::ForwardNet, ::Type{BatchNorm1},
    name::Symbol,
    parent_index::Int,
    ϵ::Float32=Float32(1e-5)
    )

    parent_node = net[parent_index]

    parent = get_output(parent_node)::Vector{Float32}
    child = Array(Float32, length(parent))
    γ = Array(Float32, length(parent))
    β = Array(Float32, length(parent))
    μ = Array(Float32, length(parent))
    ν = Array(Float32, length(parent))

    node = BatchNorm1(name, γ, β, ϵ, μ, ν, parent, child)
    add_node!(net, node, parent_index)
end
function restore!(a::BatchNorm1, filename_gamma::AbstractString, filename_beta::AbstractString,
    filename_mean::AbstractString,
    filename_variance::AbstractString)

    copy!(a.γ, open(read_binary_vec, filename_gamma))
    copy!(a.β, open(read_binary_vec, filename_beta))
    copy!(a.μ, open(read_binary_vec, filename_mean))
    copy!(a.ν, open(read_binary_vec, filename_variance))

    a
end

type BatchNorm3 <: Layer
    # https://arxiv.org/pdf/1502.03167.pdf
    name::Symbol
    γ::Vector{Float32}
    β::Vector{Float32}
    ϵ::Float32
    μ::Vector{Float32} # train set mean
    ν::Vector{Float32} # train set variance

    parent::Array{Float32, 3}
    child::Array{Float32, 3}
end
get_name(a::BatchNorm3) = a.name
get_output(a::BatchNorm3) = a.child
function forward!(a::BatchNorm3)
    for i in 1 : size(a.child)[3]
        x = (a.parent[:, :, i] - a.μ[i]) / sqrt(a.ν[i] + a.ϵ)
        a.child[:, :, i] = a.γ[i] * x + a.β[i]
    end
    a
end
function add_node!(net::ForwardNet, ::Type{BatchNorm3},
    name::Symbol,
    parent_index::Int,
    ϵ::Float32=Float32(1e-5)
    )

    parent_node = net[parent_index]

    parent = get_output(parent_node)::Array{Float32, 3}
    child = Array(Float32, size(parent))
    _, _, in_c = size(parent)
    γ = Array(Float32, in_c)
    β = Array(Float32, in_c)
    μ = Array(Float32, in_c)
    ν = Array(Float32, in_c)

    node = BatchNorm3(name, γ, β, ϵ, μ, ν, parent, child)
    add_node!(net, node, parent_index)
end
function restore!(a::BatchNorm3, filename_gamma::AbstractString, filename_beta::AbstractString,
    filename_mean::AbstractString,
    filename_variance::AbstractString)

    copy!(a.γ, open(read_binary_vec, filename_gamma))
    copy!(a.β, open(read_binary_vec, filename_beta))
    copy!(a.μ, open(read_binary_vec, filename_mean))
    copy!(a.ν, open(read_binary_vec, filename_variance))

    a
end

type Concatenator <: Layer
    name::Symbol
    parents::Vector{Vector{Float32}}
    child::Vector{Float32}
end
get_name(a::Concatenator) = a.name
get_output(a::Concatenator) = a.child
function forward!(a::Concatenator)
    i = 0
    for parent in a.parents
        for j in 1 : length(parent)
            i += 1
            a.child[i] = parent[j]
        end
    end
    a
end
function add_node!(net::ForwardNet, ::Type{Concatenator},
    name::Symbol,
    parent_indeces::Vector{Int}
    )

    tot_len = 0
    parents = Array(Vector{Float32}, length(parent_indeces))
    for (i,parent_index) in enumerate(parent_indeces)
        parents[i] = get_output(net[parent_index])::Vector{Float32}
        tot_len += length(parents[i])
    end

    child = Array(Float32, tot_len)

    node = Concatenator(name, parents, child)
    add_node!(net, node, parent_indeces)
end

type Reshaper <: Layer
    parent::Array{Float32}
    child::Array{Float32}
end
get_output(a::Reshaper) = a.child
function forward!(a::Reshaper)
    # copy!(a.child, a.parent)
    dest_ind = 0
    for i in 1 : size(a.parent, 1)
        for j in 1 : size(a.parent, 2)
            for k in 1 : size(a.parent, 3)
                a.child[dest_ind+=1] = a.parent[i,j,k]
            end
        end
    end
    a
end
function add_node!(net::ForwardNet, ::Type{Reshaper},
    parent_index::Int,
    new_shape::Tuple{Vararg{Int}},
    )

    parent = get_output(net[parent_index])::Array{Float32}
    child = Array(Float32, new_shape...)

    node = Reshaper(parent, child)
    add_node!(net, node, parent_index)
end
function add_node!(net::ForwardNet, ::Type{Reshaper},
    parent_index::Int,
    )

    # this just makes it flat

    parent = get_output(net[parent_index])::Array{Float32}
    child = Array(Float32, length(parent))

    node = Reshaper(parent, child)
    add_node!(net, node, parent_index)
end

type TanH <: Activation
    member::Array{Float32}
end
get_output(a::TanH) = a.member
function forward!(a::TanH)
    for i in 1 : length(a.member)
        a.member[i] = tanh(a.member[i])
    end
    a
end

type ReLU <: Activation
    member::Array{Float32}
end
get_output(a::ReLU) = a.member
function forward!(a::ReLU)
    for i in 1 : length(a.member)
        if a.member[i] < 0.0
            a.member[i] = 0.0
        end
    end
    a
end

type SoftPlus <: Activation
    member::Array{Float32}
end
get_output(a::SoftPlus) = a.member
function forward!(a::SoftPlus)
    for i in 1 : length(a.member)
        a.member[i] = ln(1+exp(a.member[i]))
    end
end

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

immutable ForwardPass
    net::ForwardNet
    input::Array{Symbol}
    output::Array{Symbol}
    activation_order::Vector{Int} # in order
end
function Base.show(io::IO, forwardpass::ForwardPass)
    net = forwardpass.net

    if isempty(forwardpass.activation_order)
        print(io, "ForwardPass: (empty) ", forwardpass.input, " → ", forwardpass.output)
    else
        print(io, "ForwardPass: ", get_name(net[forwardpass.activation_order[1]]))
        for i in 2 : length(forwardpass.activation_order)
            print(io, " → ", get_name(net[forwardpass.activation_order[i]]))
        end
    end
end

function calc_forwardpass(net::ForwardNet, input::Array{Symbol}, output::Array{Symbol})
    #=
    1 - get a topological sort of net
    2 - run through nodes in topologogical order and activate those that are:
            # have at least one ancestor in input and one descendent in output
    3 - will need to run forward! on all active nodes, in topological order
    =#

    input_indeces = Set{Int}()
    input_dijkstra = Dict{Int, LightGraphs.DijkstraState{Int}}()
    for name in input
        index = indexof(net, name)
        push!(input_indeces, index)
        input_dijkstra[index] = dijkstra_shortest_paths(net.dag, index)
    end

    output_indeces = Set{Int}()
    for name in output
        push!(output_indeces, indexof(net, name))
    end

    activation_order = Int[]
    for i in topological_sort_by_dfs(net.dag)

        if !isa(net[i], Variable)

            add_it = false

            for parent in input_indeces

                if parent == i
                    add_it = true
                elseif !in(i, output_indeces) && input_dijkstra[parent].dists[i] != typemax(Int) # is descendent of parent
                    dijkstra = dijkstra_shortest_paths(net.dag, i)

                    for child in output_indeces
                        if dijkstra.dists[child] != typemax(Int)
                            # child is a descendent of our node
                            add_it = true
                        end
                    end

                end
            end

            if add_it
                push!(activation_order, i)
            end
        end
    end

    ForwardPass(net, input, output, activation_order)
end
function forward!(forwardpass::ForwardPass)
    for index in forwardpass.activation_order
        forward!(forwardpass.net.nodes[index])
    end
    forwardpass
end
function print_forward_pass(io::IO, forwardpass::ForwardPass)
    for (count,index) in enumerate(forwardpass.activation_order)
        node = net[index]
        @printf(io, "node%d = net[%d]::%s # %s \n", count, index, typeof(node), get_name(node))
        @printf(io, "forward!(node%d)\n", count)
    end
end
print_forward_pass(forwardpass::ForwardPass) = print_forward_pass(STDOUT, forwardpass)

end # module