text:
"""
`struct ProductEvaluator` : specifies a ProductEvaluator, which is basically defined
through a PIBasis and its coefficients. The n-correlations are evaluated directly
via a naive product of the atomic base.
"""
mutable struct ProductEvaluator{T, TPI <: PIBasis, REAL}
   pibasis::TPI        # AA basis from ACE papers
   coeffs::Vector{T}   # c̃ coefficients from ACE papers
   real::REAL          # the real operation stored in the SymmetricBasis
end

==(V1::ProductEvaluator, V2::ProductEvaluator) =
      (V1.pibasis == V2.pibasis) && (V1.coeffs == V2.coeffs)

# ----------- FIO

write_dict(ev::ProductEvaluator) = Dict( "__id__" => "ACE_ProductEvaluator" )

read_dict(::Val{:ACE_ProductEvaluator}, D::Dict, basis, c) = ProductEvaluator(basis, c)

# ------------------------------------------------------------
#   Initialisation and Parameter manipulation code
# ------------------------------------------------------------

ProductEvaluator(basis::SymmetricBasis, c::AbstractVector) =
      ProductEvaluator(basis.pibasis, _get_eff_coeffs(basis, c), basis.real)

# basic setter interface without any checks
function set_params!(ev::ProductEvaluator, c̃::AbstractVector)
   ev.coeffs[:] .= c̃
   return ev
end

# trivial setter when the parameters already come in the c̃ format (AA-basis)
function set_params!(ev::ProductEvaluator, basis::PIBasis, c̃::AbstractVector)
   @assert ev.pibasis === basis
   set_params!(ev, c̃)
end

# if parameters come in c (B-basis) format then they first need to be converted
# to c̃ format (AA basis)
function set_params!(ev::ProductEvaluator, basis::SymmetricBasis, c::AbstractVector)
   len_AA = length(ev.pibasis)
   @assert len_AA == size(basis.A2Bmap, 2)
   c̃ = _acquire_ctilde(basis, len_AA, c)
   _get_eff_coeffs!(c̃, basis, c)
   set_params!(ev, basis.pibasis, c̃)
   release!(basis.B_pool, c̃)
   return ev
end

_get_eff_coeffs!(c̃, basis::SymmetricBasis, c::AbstractVector) =
   genmul!(c̃, transpose(basis.A2Bmap), c, *)

function _get_eff_coeffs(basis::SymmetricBasis, c::AbstractVector)
   # c̃ = acquire_B!(basis, size(basis.A2Bmap, 2))
   c̃ = _alloc_ctilde(basis,c)
   return _get_eff_coeffs!(c̃, basis, c)
end

# TODO: we may need a second pool to allocate ctilde vectors...

struct _One <: Number
end

import Base: *
*(x, ::_One) = x
*(::_One, x) = x
*(x::AbstractProperty, ::_One) = x
*(::_One, x::AbstractProperty) = x
*(x::StaticArray, ::_One) = x
*(::_One, x::StaticArray) = x
*(::ACE._One, x::ACE.XState) = x
*(x::ACE.XState, ::ACE._One) = x

_acquire_ctilde(basis::SymmetricBasis, len_AA, c::AbstractVector{<: Number}) =
      zeros(promote_type(eltype(basis.A2Bmap), eltype(c)), len_AA)

_acquire_ctilde(basis::SymmetricBasis, len_AA, c::AbstractVector{<: SVector}) =
      zeros(SVector{length(c[1]),
                    promote_type(eltype(basis.A2Bmap), eltype(c[1]))
                   } , len_AA )

_alloc_ctilde(basis::SymmetricBasis, c::AbstractVector{<: SVector}) =
      zeros(SVector{length(c[1]),eltype(basis.A2Bmap)}, size(basis.A2Bmap, 2))

_alloc_ctilde(basis::SymmetricBasis, c::AbstractVector{<: Number}) =
      zeros(eltype(basis.A2Bmap), size(basis.A2Bmap, 2))

_alloc_dAco(dAAdA::AbstractVector, A::AbstractVector, c̃::AbstractArray, args...) =
            _alloc_dAco(dAAdA, A, c̃[1], args...)

function _alloc_dAco(dAAdA::AbstractVector, A::AbstractVector,
                     c̃::Union{TP, SVector{N, TP}}, dp = _One()
                     ) where {N, TP <: AbstractProperty}
   c̃_dp = contract(c̃, dp)
   zeros( promote_type(eltype(dAAdA), typeof(c̃_dp)), length(A) )
end

# ------------------------------------------------------------
#   Standard Evaluation code
# ------------------------------------------------------------

evaluate(::LinearACEModel, V::ProductEvaluator, cfg::AbstractConfiguration) =
      evaluate(V::ProductEvaluator, cfg)

# compute one "site energy"
function evaluate(V::ProductEvaluator, cfg::AbstractConfiguration)
   A = acquire_B!(V.pibasis.basis1p, cfg)
   evaluate!(A, V.pibasis.basis1p, cfg)
   spec = V.pibasis.spec
   pireal = V.pibasis.real
   symreal = V.real
   # initialize output with a sensible type
   val = symreal(zero(eltype(V.coeffs)) * pireal(zero(eltype(A))))

   # constant (0-order)
   if spec.orders[1] == 0
      val += V.coeffs[1]
      iAAinit = 2
   else
      iAAinit = 1
   end

   @inbounds for iAA = iAAinit:length(spec)
      aa = A[spec.iAA2iA[iAA, 1]]
      for t = 2:spec.orders[iAA]
         aa *= A[spec.iAA2iA[iAA, t]]
      end
      val += symreal(pireal(aa) * V.coeffs[iAA])
   end

   release_B!(V.pibasis.basis1p, A)
   return val
end

# compute one site energy gradient
grad_config!(g, m::LinearACEModel, V::ProductEvaluator,
                     cfg::AbstractConfiguration) =
      _rrule_evaluate!(g, _One(), m, V, cfg)

function _rrule_evaluate(dp, model::LinearACEModel, cfg::AbstractConfiguration)
   g = acquire_grad_config!(model, cfg, dp)
   return _rrule_evaluate!(g, dp, model, model.evaluator, cfg)
end

# NB - testing shows that pre-allocating everything gains about 10% for small
#      configs and ca 20% for larger, more realistic configs.
#      worth doing at some point, but not really an immediate priority!
function _rrule_evaluate!(g, dp, m::LinearACEModel, V::ProductEvaluator,
                          cfg::AbstractConfiguration)
   basis1p = V.pibasis.basis1p
   pireal = V.pibasis.real
   symreal = V.real
   A = acquire_B!(V.pibasis.basis1p, cfg)
   dA = acquire_dB!(V.pibasis.basis1p, cfg)    # MAJOR ALLOCATION!!
   dAAdA = _acquire_dAAdA!(V.pibasis)

   # stage 1: precompute all the A values
   evaluate_ed!(A, dA, basis1p, cfg)

   # stage 2: compute the coefficients for the ∇A_{nlm} = ∇ϕ_{nlm}
   # dAco[nlm] = coefficient of ∇A_{nlm} (via adjoints)
   c̃ = V.coeffs
   dAco =  _alloc_dAco(dAAdA, A, c̃, dp)        # TODO: ALLOCATION
   spec = V.pibasis.spec

   if spec.orders[1] == 0; iAAinit = 2; else iAAinit = 1; end

   fill!(dAco, zero(eltype(dAco)))
   @inbounds for iAA = iAAinit:length(spec)
      _AA_local_adjoints!(dAAdA, A, spec.iAA2iA, iAA, spec.orders[iAA], pireal)
      @fastmath for t = 1:spec.orders[iAA]
         dAco[spec.iAA2iA[iAA, t]] += dAAdA[t] * contract(dp, c̃[iAA])
      end
   end

   # stage 3: get the gradients
   fill!(g, zero(eltype(g)))
   for iX = 1:length(cfg), iA = 1:length(basis1p)
      g[iX] += symreal( coco_o_daa(dAco[iA], dA[iA, iX]) )
   end

   release_B!(V.pibasis.basis1p, A)
   release_dB!(V.pibasis.basis1p, dA)

   return g
end

function adjoint_EVAL_D(m::LinearACEModel, V::ProductEvaluator, cfg, w)
   basis1p = V.pibasis.basis1p
   TDX = gradtype(m.basis, cfg)
   _real = V.real
   A = acquire_B!(V.pibasis.basis1p, cfg)
   dA = acquire_dB!(V.pibasis.basis1p, cfg)
   dAAdA = _acquire_dAAdA!(V.pibasis)

   # dAw = acquire_B!(V.pibasis.basis1p, cfg)
   # dAAw = acquire_B!(V.pibasis, cfg)
   _dAw = contract(w[1], dA[1])
   dAw = zeros(typeof(_dAw), length(V.pibasis.basis1p))

   _dAAw = _real(_dAw * dAAdA[1])
   dAAw = zeros(typeof(_dAAw), length(V.pibasis))

   # [1] dA_t = ∑_j ∂ϕ_t / ∂X_j
   evaluate_ed!(A, dA, basis1p, cfg)
   # fill!(dAw, 0)
   for k = 1:length(basis1p), j = 1:length(w)
      dAw[k] += contract(w[j], dA[k, j])
   end

   # [2] dAA_k
   spec = V.pibasis.spec
   # fill!(dAAw, 0)
   if spec.orders[1] == 0; iAAinit=2; else; iAAinit=1; end
   @inbounds for iAA = iAAinit:length(spec)
      ord = spec.orders[iAA]
      _AA_local_adjoints!(dAAdA, A, spec.iAA2iA, iAA, ord, _real)
      @fastmath for t = 1:ord
         vt = spec.iAA2iA[iAA, t]
         dAAw[iAA] += _real(dAw[vt] * dAAdA[t])
      end
   end

   # dB = similar(m.c)
   # @show m.basis.A2Bmap[1]
   # @show dAAw[1]
   # @show m.basis.A2Bmap[1].val * dAAw[1]
   # genmul!(dB, m.basis.A2Bmap, dAAw, (a, x) -> a * x)
   dB = m.basis.A2Bmap * dAAw

   release_B!(V.pibasis.basis1p, A)
   release_dB!(V.pibasis.basis1p, dA)
   # release_B!(V.pibasis.basis1p, dAw)
   # release_B!(V.pibasis, dAAw)

   # [3] dB_k
   return dB
end

# #for multiple properties. dispatch on the pullback input being a matrix.
# #Basically the same code, except for some parts where we loop over all properties.
# #We generate a list of size "nprop" and keep the same objects as for a single property
# #inside the list.
# function adjoint_EVAL_D(m::LinearACEModel, V::ProductEvaluator, cfg, wt::Matrix)
#    _contract = ACE.contract

#    basis1p = V.pibasis.basis1p
#    dAAdA = zero(MVector{10, ComplexF64})   # TODO: VERY RISKY -> FIX THIS
#    A = zeros(ComplexF64, length(basis1p))
#    TDX = gradtype(m.basis, cfg)
#    dA = zeros(complex(TDX) , length(A), length(cfg))
#    _real = V.real
#    dAAw = [acquire_B!(V.pibasis, cfg) for _ in 1:length(m.c[1])]
#    dAw = [similar(A) for _ in 1:length(m.c[1])]
#    dB = similar(m.c)

#    # [1] dA_t = ∑_j ∂ϕ_t / ∂X_j
#    evaluate_ed!(A, dA, basis1p, cfg)
#    for i in 1:length(m.c[1])
#       fill!(dAw[i], 0)
#    end
#    for prop in 1:length(m.c[1])
#       w = wt[:,prop]
#       for k = 1:length(basis1p), j = 1:length(w)
#          dAw[prop][k] += _contract(w[j], dA[k, j])
#       end
#    end

#    # [2] dAA_k
#    spec = V.pibasis.spec
#    for i in 1:length(m.c[1])
#       fill!(dAAw[i], 0)
#    end
#    for prop in 1:length(m.c[1])
#       if spec.orders[1] == 0; iAAinit=2; else; iAAinit=1; end
#       @inbounds for iAA = iAAinit:length(spec)
#          _AA_local_adjoints!(dAAdA, A, spec.iAA2iA, iAA, spec.orders[iAA], _real)
#          @fastmath for t = 1:spec.orders[iAA]
#             vt = spec.iAA2iA[iAA, t]
#             dAAw[prop][iAA] += _real(dAw[prop][vt] * dAAdA[t])
#          end
#       end
#    end

#    adjointgenmul!(dB, m.basis.A2Bmap, dAAw, (a, x) -> a.val * x)

#    for i in 1:length(m.c[1])
#       release_B!(V.pibasis, dAAw[i]) #TODO check that this indeed releases
#    end

#    # [3] dB_k
#    return dB
# end


text:
"""
	grdgradient(cmd0::String="", arg1=nothing, kwargs...)

Compute the directional derivative in a given direction, or to find the direction [and the magnitude]
of the vector gradient of the data.

Full option list at [`grdgradient`]($(GMTdoc)grdgradient.html)

Parameters
----------

- **A** | **azim** :: [Type => Str | Number]    ``Arg = azim[/azim2]``

    Azimuthal direction for a directional derivative.
    ($(GMTdoc)grdgradient.html#a)
- **D** | **find_dir** :: [Type => Str]      ``Arg = [a][c][o][n]``

    Find the direction of the positive (up-slope) gradient of the data.
    ($(GMTdoc)grdgradient.html#d)
- **G** | **save** | **write** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdgradient(....) form.
    ($(GMTdoc)grdgradient.html#g)
- **E** | **lambert** :: [Type => Str]    ``Arg = [m|s|p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine] ``

    Compute Lambertian radiance appropriate to use with grdimage and grdview.
    ($(GMTdoc)grdgradient.html#e)
- **N** | **norm** | **normalize** :: [Type => Str]     ``Arg = [e|t][amp][+ssigma][+ooffset]``

    Normalization. [Default is no normalization.] The actual gradients g are offset and scaled
    to produce normalized gradients.
    ($(GMTdoc)grdgradient.html#n)
- **Q** | **save_stats** :: [Type => Str]		``Arg = c|r|R``

    Controls how normalization via N is carried out.
    ($(GMTdoc)grdgradient.html#q)
- $(GMT.opt_R)
- **S** | **slopegrid** :: [Type => Str]

    Name of output grid file with scalar magnitudes of gradient vectors. Requires D but makes G optional.
    ($(GMTdoc)grdgradient.html#s)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdgradient(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdgradient", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:G :R :V_params :f :n])
	cmd  = parse_these_opts(cmd, d, [[:A :azim], [:D :find_dir], [:S :slopegrid]])
	cmd  = add_opt(d, cmd, "E", [:E :lambert],
	       (manip=("m", nothing, 1), simple=("s", nothing, 1), peucker=("p", nothing, 1), view=("", arg2str), ambient="+a", difuse="+d", specular="+p", shine="+s") )
	cmd  = add_opt(d, cmd, "N", [:N :norm :normalize],
           (laplace=("e", nothing, 1), cauchy=("t", nothing, 1), amp="", sigma="+s", offset="+o"))
	cmd = parse_Q_grdgrad(d, [:Q :save_stats], cmd)

	common_grd(d, cmd0, cmd, "grdgradient ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
function parse_Q_grdgrad(d::Dict, symbs::Array{<:Symbol}, cmd::String)
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "String")
    if ((val = find_in_dict(d, symbs)[1]) !== nothing)
		val = string(val)[1]
		(val == 's') && (val = 'c')
		(val == 'c' || val == 'r' || val == 'R') && (cmd *= " -Q" * val)
    end
    return cmd
end

# ---------------------------------------------------------------------------------------------------
grdgradient(arg1, cmd0::String=""; kw...) = grdgradient(cmd0, arg1; kw...)

ext: jl


text:
import Splittings:Geometry

@testset " Geometry " begin

geom1 = Geometry( 11, 21, 0.0, 0.0, 0.1, 0.1 )
geom2 = Geometry( 11, 21, (0.0, 1.0), (0.0, 2.0), :none )

@test geom1.n1     == geom2.n1     == 11
@test geom1.n2     == geom2.n2     == 21
@test geom1.delta1 == geom2.delta1 == 0.1
@test geom1.delta2 == geom2.delta2 == 0.1
@test geom1.x1min  == geom2.x1min  == 0.0
@test geom1.x2min  == geom2.x2min  == 0.0
@test geom1.x1grid[end] == geom2.x1grid[end]  == 1.0
@test geom1.x2grid[end] == geom2.x2grid[end]  == 2.0

end

@testset " Geometry perx" begin

geom1 = Geometry( 10, 21, 0.0, 0.0, 0.1, 0.1 )
geom2 = Geometry( 10, 21, (0.0, 1.0), (0.0, 2.0), :perx )

@test geom1.n1     == geom2.n1     == 10
@test geom1.n2     == geom2.n2     == 21
@test geom1.delta1 == geom2.delta1 == 0.1
@test geom1.delta2 == geom2.delta2 == 0.1
@test geom1.x1min  == geom2.x1min  == 0.0
@test geom1.x2min  == geom2.x2min  == 0.0

end

@testset " Geometry pery" begin

geom1 = Geometry( 11, 20, 0.0, 0.0, 0.1, 0.1 )
geom2 = Geometry( 11, 20, (0.0, 1.0), (0.0, 2.0), :pery )

@test geom1.n1     == geom2.n1     == 11
@test geom1.n2     == geom2.n2     == 20
@test geom1.delta1 == geom2.delta1 == 0.1
@test geom1.delta2 == geom2.delta2 == 0.1
@test geom1.x1min  == geom2.x1min  == 0.0
@test geom1.x2min  == geom2.x2min  == 0.0

end

@testset " Geometry perxy" begin

geom1 = Geometry( 10, 20, 0.0, 0.0, 0.1, 0.1 )
geom2 = Geometry( 10, 20, (0.0, 1.0), (0.0, 2.0), :perxy )

@test geom1.n1     == geom2.n1     == 10
@test geom1.n2     == geom2.n2     == 20
@test geom1.delta1 == geom2.delta1 == 0.1
@test geom1.delta2 == geom2.delta2 == 0.1
@test geom1.x1min  == geom2.x1min  == 0.0
@test geom1.x2min  == geom2.x2min  == 0.0

end


text:
function factorial(n::Integer)
    if n < 0
        return zero(n)
    end
    f = one(n)
    for i = 2:n
        f *= i
    end
    return f
end

# computes n!/k!
function factorial{T<:Integer}(n::T, k::T)
    if k < 0 || n < 0 || k > n
        return zero(T)
    end
    f = one(T)
    while n > k
        f *= n
        n -= 1
    end
    return f
end

function binomial{T<:Integer}(n::T, k::T)
    k < 0 && return zero(T)
    sgn = one(T)
    if n < 0
        n = -n + k -1
        if isodd(k)
            sgn = -sgn
        end
    end
    k > n && return zero(T)
    (k == 0 || k == n) && return sgn
    k == 1 && return sgn*n
    if k > (n>>1)
        k = (n - k)
    end
    x::T = nn = n - k + 1
    nn += 1
    rr = 2
    while rr <= k
        xt = div(widemul(x, nn), rr)
        x = xt
        x == xt || throw(OverflowError())
        rr += 1
        nn += 1
    end
    convert(T, copysign(x, sgn))
end

## other ordering related functions ##

function shuffle!(a::AbstractVector)
    for i = length(a):-1:2
        j = rand(1:i)
        a[i], a[j] = a[j], a[i]
    end
    return a
end

shuffle(a::AbstractVector) = shuffle!(copy(a))

function randperm(n::Integer)
    a = Array(typeof(n), n)
    a[1] = 1
    for i = 2:n
        j = rand(1:i)
        a[i] = a[j]
        a[j] = i
    end
    return a
end

function randcycle(n::Integer)
    a = Array(typeof(n), n)
    a[1] = 1
    for i = 2:n
        j = rand(1:i-1)
        a[i] = a[j]
        a[j] = i
    end
    return a
end

function nthperm!(a::AbstractVector, k::Integer)
    k -= 1 # make k 1-indexed
    n = length(a)
    f = factorial(oftype(k, n-1))
    for i=1:n-1
        j = div(k, f) + 1
        k = k % f
        f = div(f, n-i)

        j = j+i-1
        elt = a[j]
        for d = j:-1:i+1
            a[d] = a[d-1]
        end
        a[i] = elt
    end
    a
end
nthperm(a::AbstractVector, k::Integer) = nthperm!(copy(a),k)

function nthperm{T<:Integer}(p::AbstractVector{T})
    isperm(p) || error("argument is not a permutation")
    k, n = 1, length(p)
    for i = 1:n-1
        f = factorial(n-i)
        for j = i+1:n
            k += ifelse(p[j] < p[i], f, 0)
        end
    end
    return k
end

function invperm(a::AbstractVector)
    b = zero(a) # similar vector of zeros
    n = length(a)
    for i = 1:n
        j = a[i]
        ((1 <= j <= n) && b[j] == 0) ||
            error("argument is not a permutation")
        b[j] = i
    end
    b
end

function isperm(A)
    n = length(A)
    used = falses(n)
    for a in A
        (0 < a <= n) && (used[a] $= true) || return false
    end
    true
end

function permute!!{T<:Integer}(a, p::AbstractVector{T})
    count = 0
    start = 0
    while count < length(a)
        ptr = start = findnext(p, start+1)
        temp = a[start]
        next = p[start]
        count += 1
        while next != start
            a[ptr] = a[next]
            p[ptr] = 0
            ptr = next
            next = p[next]
            count += 1
        end
        a[ptr] = temp
        p[ptr] = 0
    end
    a
end

permute!(a, p::AbstractVector) = permute!!(a, copy(p))

function ipermute!!{T<:Integer}(a, p::AbstractVector{T})
    count = 0
    start = 0
    while count < length(a)
        start = findnext(p, start+1)
        temp = a[start]
        next = p[start]
        count += 1
        while next != start
            temp_next = a[next]
            a[next] = temp
            temp = temp_next
            ptr = p[next]
            p[next] = 0
            next = ptr
            count += 1
        end
        a[next] = temp
        p[next] = 0
    end
    a
end

ipermute!(a, p::AbstractVector) = ipermute!!(a, copy(p))

immutable Combinations{T}
    a::T
    t::Int
end

eltype(c::Combinations) = typeof(c.a)
eltype{T}(c::Combinations{UnitRange{T}}) = Array{T,1}
eltype{T}(c::Combinations{Range{T}}) = Array{T,1}

length(c::Combinations) = binomial(length(c.a),c.t)

function combinations(a, t::Integer)
    if t < 0
        # generate 0 combinations for negative argument
        t = length(a)+1
    end
    Combinations(a, t)
end

start(c::Combinations) = [1:c.t]
function next(c::Combinations, s)
    comb = c.a[s]
    if c.t == 0
        # special case to generate 1 result for t==0
        return (comb,[length(c.a)+2])
    end
    s = copy(s)
    for i = length(s):-1:1
        s[i] += 1
        if s[i] > (length(c.a) - (length(s)-i))
            continue
        end
        for j = i+1:endof(s)
            s[j] = s[j-1]+1
        end
        break
    end
    (comb,s)
end
done(c::Combinations, s) = !isempty(s) && s[1] > length(c.a)-c.t+1

immutable Permutations{T}
    a::T
end

eltype(c::Permutations) = typeof(c.a)
eltype{T}(c::Permutations{UnitRange{T}}) = Array{T,1}
eltype{T}(c::Permutations{Range{T}}) = Array{T,1}

length(c::Permutations) = factorial(length(c.a))

permutations(a) = Permutations(a)

start(p::Permutations) = [1:length(p.a)]
function next(p::Permutations, s)
    if length(p.a) == 0
        # special case to generate 1 result for len==0
        return (p.a,[1])
    end
    s = copy(s)
    perm = p.a[s]
    k = length(s)-1
    while k > 0 && s[k] > s[k+1];  k -= 1;  end
    if k == 0
        s[1] = length(s)+1   # done
    else
        l = length(s)
        while s[k] >= s[l];  l -= 1;  end
        s[k],s[l] = s[l],s[k]
        reverse!(s,k+1)
    end
    (perm,s)
end
done(p::Permutations, s) = !isempty(s) && s[1] > length(p.a)

# Integer Partitions

immutable IntegerPartitions
    n::Int
end

length(p::IntegerPartitions) = npartitions(p.n)

partitions(n::Integer) = IntegerPartitions(n)

start(p::IntegerPartitions) = Int[]
done(p::IntegerPartitions, xs) = length(xs) == p.n
next(p::IntegerPartitions, xs) = (xs = nextpartition(p.n,xs); (xs,xs))

function nextpartition(n, as)
    if isempty(as);  return Int[n];  end

    xs = similar(as,0)
    sizehint(xs,length(as)+1)

    for i = 1:length(as)-1
        if as[i+1] == 1
            x = as[i]-1
            push!(xs, x)
            n -= x
            while n > x
                push!(xs, x)
                n -= x
            end
            push!(xs, n)

            return xs
        end
        push!(xs, as[i])
        n -= as[i]
    end
    push!(xs, as[end]-1)
    push!(xs, 1)

    xs
end

let _npartitions = (Int=>Int)[]
    global npartitions
    function npartitions(n::Int)
        if n < 0
            0
        elseif n < 2
            1
        elseif (np = get(_npartitions, n, 0)) > 0
            np
        else
            np = 0
            sgn = 1
            for k = 1:n
                np += sgn * (npartitions(n-k*(3k-1)>>1) + npartitions(n-k*(3k+1)>>1))
                sgn = -sgn
            end
            _npartitions[n] = np
        end
    end
end

# Algorithm H from TAoCP 7.2.1.4
# Partition n into m parts
# in colex order (lexicographic by reflected sequence)

immutable FixedPartitions
    n::Int
    m::Int
end

length(f::FixedPartitions) = npartitions(f.n,f.m)

partitions(n::Integer, m::Integer) = (@assert 2 <= m <= n; FixedPartitions(n,m))

start(f::FixedPartitions) = Int[]
done(f::FixedPartitions, s::Vector{Int}) = !isempty(s) && s[1]-1 <= s[end]
next(f::FixedPartitions, s::Vector{Int}) = (xs = nextfixedpartition(f.n,f.m,s); (xs,xs))

function nextfixedpartition(n, m, bs)
    as = copy(bs)
    if isempty(as)
        # First iteration
        as = [n-m+1, ones(Int, m-1)]
    elseif as[2] < as[1]-1
        # Most common iteration
        as[1] -= 1
        as[2] += 1
    else
        # Iterate
        local j
        s = as[1]+as[2]-1
        for j = 3:m
            if as[j] < as[1]-1; break; end
            s += as[j]
        end
        x = as[j] += 1
        for k = j-1:-1:2
            as[k] = x
            s -= x
        end
        as[1] = s
    end

    return as
end

let _nipartitions = ((Int,Int)=>Int)[]
    global npartitions
    function npartitions(n::Int,m::Int)
        if n < m || m == 0
            0
        elseif n == m
            1
        elseif (np = get(_nipartitions, (n,m), 0)) > 0
            np
        else
            _nipartitions[(n,m)] = npartitions(n-1,m-1) + npartitions(n-m,m)
        end
    end
end

# Algorithm H from TAoCP 7.2.1.5
# Set partitions

immutable SetPartitions{T<:AbstractVector}
    s::T
end

length(p::SetPartitions) = nsetpartitions(length(p.s))

partitions(s::AbstractVector) = SetPartitions(s)

start(p::SetPartitions) = (n = length(p.s); (zeros(Int32, n), ones(Int32, n-1), n, 1))
done(p::SetPartitions, s) = !isempty(s) && s[1][1] > 0
next(p::SetPartitions, s) = nextsetpartition(p.s, s...)

function nextsetpartition(s::AbstractVector, a, b, n, m)
    function makeparts(s, a, m)
        temp = [ similar(s,0) for k = 0:m ]
        for i = 1:n
            push!(temp[a[i]+1], s[i])
        end
        filter!(x->!isempty(x), temp)
    end

    if isempty(s);  return ({s}, ([1], Int[], n, 1));  end

    part = makeparts(s,a,m)

    if a[end] != m
        a[end] += 1
    else
        local j
        for j = n-1:-1:1
            if a[j] != b[j]
                break
            end
        end
        a[j] += 1
        m = b[j] + (a[j] == b[j])
        for k = j+1:n-1
            a[k] = 0
            b[k] = m
        end
        a[end] = 0
    end

    return (part, (a,b,n,m))

end

let _nsetpartitions = (Int=>Int)[]
    global nsetpartitions
    function nsetpartitions(n::Int)
        if n < 0
            0
        elseif n < 2
            1
        elseif (wn = get(_nsetpartitions, n, 0)) > 0
            wn
        else
            wn = 0
            for k = 0:n-1
                wn += binomial(n-1,k)*nsetpartitions(n-1-k)
            end
            _nsetpartitions[n] = wn
        end
    end
end

immutable FixedSetPartitions{T<:AbstractVector}
    s::T
    m::Int
end

length(p::FixedSetPartitions) = nfixedsetpartitions(length(p.s),p.m)

partitions(s::AbstractVector,m::Int) = (@assert 2 <= m <= length(s); FixedSetPartitions(s,m))

start(p::FixedSetPartitions) = (n = length(p.s);m=p.m; (vcat(ones(Int, n-m),1:m), vcat(1,n-m+2:n), n))
# state consists of vector a of length n describing to which partition every element of s belongs and
# a vector b of length m describing the first index b[i] that belongs to partition i

done(p::FixedSetPartitions, s) = !isempty(s) && s[1][1] > 1
next(p::FixedSetPartitions, s) = nextfixedsetpartition(p.s,p.m, s...)

function nextfixedsetpartition(s::AbstractVector, m, a, b, n)
    function makeparts(s, a)
        part = [ similar(s,0) for k = 1:m ]
        for i = 1:n
            push!(part[a[i]], s[i])
        end
        return part
    end

    part = makeparts(s,a)

    if a[end] != m
        a[end] += 1
    else
        local j, k
        for j = n-1:-1:1
            if a[j]<m && b[a[j]+1]<j
                break
            end
        end
        if j>1
            a[j]+=1
            for p=j+1:n
                if b[a[p]]!=p
                    a[p]=1
                end
            end
        else
            for k=m:-1:2
                if b[k-1]<b[k]-1
                    break
                end
            end
            b[k]=b[k]-1
            b[k+1:m]=n-m+k+1:n
            a[1:n]=1
            a[b]=1:m
        end
    end

    return (part, (a,b,n))
end

function nfixedsetpartitions(n::Int,m::Int)
    numpart=0
    for k=0:m
        numpart+=(-1)^(m-k)*binomial(m,k)*(k^n)
    end
    numpart=div(numpart,factorial(m))
    return numpart
end

# For a list of integers i1, i2, i3, find the smallest
#     i1^n1 * i2^n2 * i3^n3 >= x
# for integer n1, n2, n3
function nextprod(a::Vector{Int}, x)
    if x > typemax(Int)
        error("unsafe for x bigger than typemax(Int)")
    end
    k = length(a)
    v = ones(Int, k)                  # current value of each counter
    mx = [nextpow(ai,x) for ai in a]  # maximum value of each counter
    v[1] = mx[1]                      # start at first case that is >= x
    p::widen(Int) = mx[1]             # initial value of product in this case
    best = p
    icarry = 1

    while v[end] < mx[end]
        if p >= x
            best = p < best ? p : best  # keep the best found yet
            carrytest = true
            while carrytest
                p = div(p, v[icarry])
                v[icarry] = 1
                icarry += 1
                p *= a[icarry]
                v[icarry] *= a[icarry]
                carrytest = v[icarry] > mx[icarry] && icarry < k
            end
            if p < x
                icarry = 1
            end
        else
            while p < x
                p *= a[1]
                v[1] *= a[1]
            end
        end
    end
    best = mx[end] < best ? mx[end] : best
    return int(best)  # could overflow, but best to have predictable return type
end

# For a list of integers i1, i2, i3, find the largest
#     i1^n1 * i2^n2 * i3^n3 <= x
# for integer n1, n2, n3
function prevprod(a::Vector{Int}, x)
    if x > typemax(Int)
        error("unsafe for x bigger than typemax(Int)")
    end
    k = length(a)
    v = ones(Int, k)                  # current value of each counter
    mx = [nextpow(ai,x) for ai in a]  # allow each counter to exceed p (sentinel)
    first = int(prevpow(a[1], x))     # start at best case in first factor
    v[1] = first
    p::widen(Int) = first
    best = p
    icarry = 1

    while v[end] < mx[end]
        while p <= x
            best = p > best ? p : best
            p *= a[1]
            v[1] *= a[1]
        end
        if p > x
            carrytest = true
            while carrytest
                p = div(p, v[icarry])
                v[icarry] = 1
                icarry += 1
                p *= a[icarry]
                v[icarry] *= a[icarry]
                carrytest = v[icarry] > mx[icarry] && icarry < k
            end
            if p <= x
                icarry = 1
            end
        end
    end
    best = x >= p > best ? p : best
    return int(best)
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
function _sqsum(a)
    m, n = size(a)
    ans = Array(Float64, n, 1)
    for i = 1:n
        ans[i] = sum(abs2(a[:, i]))
    end
    return(ans)
end

# The function below is adapted from Dahua Lin's blog post here: http://julialang.org/blog/2013/09/fast-numeric/
function vecdist(a1::Array{Float64, 2}, a2::Array{Float64, 2})
    m, n = size(a1)
    sa = _sqsum(a1)
    sb = _sqsum(a2)
    r = sa .+ reshape(sb, 1, n)

    #Update C as alpha*A*B + beta*C
    # or the other three variants according to tA (transpose A) and tB. Returns the updated C.
    gemm!('T', 'N', -2.0, a1, a2, 1.0, r)
    for i = 1:length(r)
        if r[i] < 0
            if r[i] > -1e-9
                r[i] = 0
            else
                error("Negative squared distance in vecdist")
            end
        else
            r[i] = sqrt(r[i])
        end
    end
    return(r)
end

### functions for attractor reconstruction and distance calculation
function calc_dists(shadowmat::Array{Float64,2})
    return vecdist(shadowmat, shadowmat)
end

function precalc_manif_dists(Evals::AbstractVector{Int}, tau_vals::AbstractVector{Int}, vector::AbstractVector)
    shadowmats = (Int=>Dict)[tt => (Int=>Array{Float64,2})[E => construct_shadow(vector, E, tt) for E in Evals] for tt in tau_vals]
    distmats   = (Int=>Dict)[tt => (Int=>Array{Float64,2})[E => calc_dists(shadowmats[tt][E]) for E in Evals] for tt in tau_vals]
    return shadowmats, distmats
end

function construct_shadow(vector::AbstractVector, E::Int, tau_s::Int=1)
    if tau_s < 1
        throw(ArgumentError("tau_s must be greater than 0!"))
    end
    n::Int   = length(vector)
    lag::Int = 0
    shadowmat  = fill(NaN, E, n)

    for ii in 1:E
        shadowmat[ii, (lag+1):n] = vector[1:(n-lag)]
        lag += tau_s
    end
    return shadowmat
end
#### End attractor reconstruction and distance calculation functions


text:
# This file is ugly. Don't look at it. Should be replaced by fixed-size arrays/tuples

using Polynomials

import Base: getindex, start, done, next, length, size, eltype,
             promote_rule, zero, one, zeros, ones, conj, copy

using ImmutableArrays
using Meshes
import Meshes: Vertex2

# Coefficient vector for a polynomial basis
immutable Coeffs{N,T}
    coeffs::Vector{T} #NTuple{N,Float64}
end
Coeffs{T}(N,x::Vector{T}) = (Coeffs{N,T}(x))

copy{N,T}(c::Coeffs{N,T}) = Coeffs{N,T}(copy(c.coeffs))

order(x::Coeffs) = length(x.coeffs) #N
order(x::Vector) = length(x)
eltype{N,T}(c::Coeffs{N,T}) = T
eltype{N,T}(::Type{Coeffs{N,T}}) = T

fapply(f::Poly,x) = polyval(f,x)
fapply(vec::Vector,x) = [fapply(f,x) for f in vec]
fapply(f,x) = f(x)
fapply(x::Number,y) = x

promote_rule{N,T}(::Type{Coeffs{N,T}},::Type{Coeffs{N,T}}) = Coeffs{N,T}
promote_rule{N,T,S<:Number}(::Type{Coeffs{N,T}},::Type{S}) = Coeffs{N,T}

+{N,T}(x::Coeffs{N,T}, y::Coeffs{N,T}) = Coeffs{N,T}(x.coeffs+y.coeffs)
+{N,T}(x::Coeffs{N,T},y::Vector{Float64}) = Coeffs{N,T}(x.coeffs+y)
-{N,T}(x::Coeffs{N,T},y::Vector{Float64}) = Coeffs{N,T}(x.coeffs-y)
*{N,T}(s::Number, x::Coeffs{N,T}) = Coeffs{N,T}(s*x.coeffs)
.*{N,T}(s, x::Coeffs{N,T}) = Coeffs{N,T}(s.*x.coeffs)
.^{N,T}(s::Coeffs{N,T}, y) = Coeffs{N,T}(s.coeffs.^y)
size(x::Coeffs) = size(x.coeffs)

for f in (:getindex, :start, :done, :next, :length)
    @eval ($f)(x::Coeffs,args...) = ($f)(x.coeffs,args...)
end

zero{N,T}(::Type{Coeffs{N,T}}) = Coeffs{N,T}(zeros(Float64,N))
one{N,T}(::Type{Coeffs{N,T}}) = Coeffs{N,T}(ones(Float64,N))
zeros{N,T}(::Type{Coeffs{N,T}},n) = [zero(Coeffs{N,T}) for _=1:n]
ones{N,T}(::Type{Coeffs{N,T}},n) = [one(Coeffs{N,T}) for _=1:n]
conj(x::Coeffs) = x

# Vertex2 (from Meshes)

zero(::Type{Vertex2}) = Vertex2(zero(Vector2{Float64}))

import Base: *, -, +, dot, ./, \

for f in (:(-),:(+))
    @eval ($f)(a::Vertex2,b::Vertex2) = Vertex2(($f)(a.coords,b.coords))
end
dot(a::Vertex2,b::Vertex2) = dot(a.coords,b.coords)
*(a::Float64,b::Vertex2) = Vertex2(a*b.coords.e1,a*b.coords.e2)
+(a::Vector2{Float64},b::Vertex2) = Vertex2(a+b.coords)
-(a::Vector2{Float64},b::Vertex2) = Vertex2(a-b.coords)
+(b::Vertex2,a::Vector2{Float64}) = a+b
./(a::Vertex2,b::Number) = Vertex2(a.coords./b)
\(a,b::Vertex2) = (r = a\[b.coords.e1,b.coords.e2]; Vertex2(r[1],r[2]))
\{T}(a::Array{T,2},b::Vector2{T}) = (r = a\[b.e1,b.e2]; Vector2{T}(r[1],r[2]))
*(a::Array{Float64,2}, b::Vertex2) = (r=a*[b.coords.e1,b.coords.e2]; Vertex2(r[1],r[2]))

# Manually inlined, because this function is hot in laxcf
dot(a::Vector2,b::Vertex2) = a.e1*b.coords.e1 + a.e2*b.coords.e2

function (*){T,S}(A::StridedMatrix{T}, x::Vector2{S})
    Vector2{S}((A*[x[1];x[2]])...)
end


text:
import FractionalCalculus.FracDiffAlg

"""
Caputo sense fractional derivative algorithms, please refer to [Caputo derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Caputo_fractional_derivative) for more details.
"""
abstract type Caputo <: FracDiffAlg end

"""
# Caputo sense fractional derivative.

    fracdiff(f::Function, α, start_point, end_point, step_size ::Caputo)

### Example:

```julia-repl
julia> fracdiff(x->x^5, 0.5, 0, 2.5, 0.0001, Caputo_Direct())
```

Returns a tuple (**result**, **error**), which means the value of this derivative is 141.59714979764541, and the error estimate is 1.1532243848672914e-6.

Refer to [Caputo derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Caputo_fractional_derivative)

!!! note
    *Caputo Direct* algorithms belong to direct computing, precise are ensured, but maybe cause more memory allocation and take more compilation time.

Using the direct mathematic expression:

```math
^CD_t^α=\\frac{1}{\\Gamma(n-α)}\\int_0^t\\frac{f^{(n)}(τ)}{(t-τ)^{α+1-n}}
```

As for the derivative inside the integral, we use the **Complex Step Differentiation** to obtain a more accurate value.
"""
struct Caputo_Direct <: Caputo end

"""
# Caputo sense fractional derivative with first derivative known.

    fracdiff(fd, α, start_point, end_point, Caputo_Direct_First_Diff_Known())

If the first order derivative of a function is already known, we can use this method to compute the fractional order derivative more precisely.

The inout function should be the first order derivative of a function

### Example

```julia-repl
julia> fracdiff(x->5*x^4, 0.5, 0, 2.5, Caputo_Direct_First_Diff_Known())
```
Return the semi-derivative of ``f(x)=x^5`` at ``x=2.5``.

Compared with **Caputo_Direct** method, this method don't need to specify step size, more precision are guaranteed.
"""
struct Caputo_Direct_First_Diff_Known <: Caputo end

"""
# Caputo sense fractional derivative with first and second derivative known.

    fracdiff(fd1, fd2, α, start_point, end_point, Caputo_Direct_First_Second_Diff_Known)

If the first and second order derivative of a function is already known, we can use this method to compute the fractional order derivative more precisely.

### Example

```julia-repl
julia> fracdiff(x->5*x^4, x->20*x^3, 0.5, 0, 2.5, Caputo_Direct_First_Second_Diff_known())
```
Return the semi-derivative of ``f(x)=x^5`` at ``x=2.5``.

Compared with **Caputo_Direct** method, this method don't need to specify step size, more precision are guaranteed.
"""
struct Caputo_Direct_First_Second_Diff_Known <: Caputo end

"""

"""

"""
# Caputo sense Piecewise algorithm

    fracdiff(f, α, end_point, h, Caputo_Piecewise())

Using piecewise linear interpolation function to approximate input function and combining Caputo derivative then implement summation.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.001, Caputo_Piecewise())
```

Return the fractional derivative of ``f(x)=x^5`` at point ``x=2.5``.

```tex
@article{LI20113352,
title = {Numerical approaches to fractional calculus and fractional ordinary differential equation},
author = {Changpin Li and An Chen and Junjie Ye},
}
```

"""
struct Caputo_Piecewise <: Caputo end

"""
# Caputo sense Diethelm computation

    fracdiff(f, α, end_point, h, Caputo_Diethelm())

Using quadrature weights(derived from product trapezoidal rule) to approximate the derivative.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.007, Caputo_Diethelm())
1.128378318687192
```

!!! info "Scope"
    0 < α < 1

Diethelm's method for computingn Caputo sense fractional derivative using [product trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) and [Richardsin extrapolation](https://en.wikipedia.org/wiki/Richardson_extrapolation)

```tex
@article{10.1093/imanum/17.3.479,
author = {DIETTELM, KAI},
title = "{Generalized compound quadrature formulae for finite-part integrals}",
journal = {IMA Journal of Numerical Analysis},
doi = {10.1093/imanum/17.3.479},
}
```

"""
struct Caputo_Diethelm <: Caputo end

"""
# Caputo sense fractioal derivative with p-th order precision.

    fracdiff(f, α, t, Caputo_High_Precision())

Use the high precision algorithm to compute the Caputo sense fractional derivative.
"""
struct Caputo_High_Precision <: Caputo end

################################################################
###                    Type defination done                  ###
################################################################

function fracdiff(f::Union{Function, Number}, α::Float64, start_point, end_point, h::Float64, ::Caputo_Direct)
    #checks(f, α, start_point, end_point)

    #Using complex step differentiation to calculate the first order differentiation
    g(τ) = imag(f(τ+1*im*h) ./ h) ./ ((end_point-τ) .^α)
    result = quadgk(g, start_point, end_point) ./gamma(1-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, start_point, end_point::AbstractArray, h::Float64, ::Caputo_Direct)::Vector
    #=
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, h, Caputo_Direct()))
    end
    return ResultArray
    =#
    result = map(x->fracdiff(f, α, start_point, x, h, Caputo_Direct()), end_point)
    return result
end

function fracdiff(fd::Function, α, start_point, end_point, ::Caputo_Direct_First_Diff_Known)

    g(τ) = (end_point-τ) .^ (-α) * fd(τ)
    result = quadgk(g, start_point, end_point) ./gamma(1-α)
    return result
end

function fracdiff(f::Function, α, start_point, end_point::AbstractArray, step_size, ::Caputo_Direct_First_Diff_Known)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, step_size, Caputo_Direct_First_Diff_Known()))
    end
    return ResultArray
end

function fracdiff(fd1::Function, fd2, α, start_point, end_point, ::Caputo_Direct_First_Second_Diff_Known)

    temp1 = fd1(start_point) .* (end_point-start_point) .^ (1-α) ./ gamma(2-α)
    g(τ) = fd2(τ) .* ((end_point-τ) .^ (1-α))
    temp2 = quadgk(g, start_point, end_point) ./ gamma(2-α)
    result = temp1 .+ temp2
    return result
end

function fracdiff(f::Function, α::Float64, start_point, end_point::AbstractArray, h, ::Caputo_Direct_First_Second_Diff_Known)::Vector
    ResultArray = Float64[]
    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracdiff(f, α, start_point, value, h, Caputo_Direct_First_Second_Diff_Known()))
    end
    return ResultArray
end

function fracdiff(f::Union{Function, Number}, α::Float64, end_point, h, ::Caputo_Piecewise)

    m = floor(α)+1

    summation = 0
    n = Int64(floor(end_point/h))

    @fastmath @inbounds @simd for i ∈ 0:n
        summation += W₅(i, n, m, α)*first_order(f, i*h, h)
    end

    result=summation*h^(m-α)/gamma(m-α+2)
    return result
end
function W₅(i, n, m, α)
    if i == 0
        return n^(m-α)*(m-α+1-n) + (n-1)^(m-α+1)
    elseif i == n
        return 1
    else
        (n-i-1)^(m-α+1) + (n-i+1)^(m-α+1) - 2*(n-i)^(m-α+1)
    end
end
function first_order(f, point, h)
    return (f(point+h)-f(point-h))/(2*h)
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::Caputo_Piecewise)::Vector
    result = map(x->fracdiff(f, α, x, h, Caputo_Piecewise()), end_point)
    return result
end

#=
Caputo Diethelm algorithm
=#
function fracdiff(f::Union{Function, Number}, α::Float64, end_point::Number, h, ::Caputo_Diethelm)
    #checks(f, α, 0, end_point)

    N = Int64(floor(end_point/h))

    result = 0

    @fastmath @inbounds @simd for i ∈ 0:N
        result += quadweights(i, N, α)*(f(end_point-i*h) - f(0))
    end

    return result*h^(-α)/gamma(2-α)
end
function quadweights(n, N, α)
    if n == 0
        return 1
    elseif  0 < n < N
        return (n+1)^(1-α) - 2*n^(1-α) + (n-1)^(1-α)
    elseif n == N
        return (1-α)*N^(-α) - N^(1-α) + (N-1)^(1-α)
    end
end

#=
Caputo sense high precision algorithm
=#
function fracdiff(y, α, t, p, ::Caputo_High_Precision)
    h = t[2]-t[1]
    t = t[:]
    n = length(t)
    y = y.(t)
    q = Int64(ceil(α))
    r = Int64(max(p, q))
    R = reverse(Vandermonde(collect(0:(r-1)).*h))
    c = inv(R)*y[1:r]
    u = 0
    du = 0
    for i = 1:r
        u = u.+c[i]*t.^(i-1)
    end
    if q < r
        for i = (q+1):p
            du = du.+c[i]*t.^(i-1-α)*gamma(1)/gamma(i-α)
        end
    end

    v = y[:] - u[:]
    dv = fracdiff(v, α, t, p, GL_High_Precision())
    dy = dv[:] + du[:]
    return dy
end


text:
using Genie, Genie.Router, Genie.Renderer.Html, Stipple

# define (reactive) model
Base.@kwdef mutable struct Name <: ReactiveModel
    name::R{String} = "World!"
end

# initialize the model
model = Stipple.init(Name)

# define UI
function ui()
    page(
    root(model), class="container", [
      h1([
        "Hello "
        span("", @text(:name); style="color:red")
      ])

      p([
        "What is your name? "
        input("", placeholder="Type your name", @bind(:name))
      ])
    ], title="Basic Stipple"
  ) |> html
end

# bind root "/" to ui
route("/", ui)

# start web server
# up()

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
import ProximalOperators: gradient! # this can be removed when moved to Prox

export PrecomposeNonlinear

struct PrecomposeNonlinear{P <: ProximableFunction,
                           T <: AbstractOperator,
                           D, C
                          } <: ProximableFunction
	g::P
    G::T
    bufD::D
    bufC::C
    bufC2::C
end

function PrecomposeNonlinear(g::P, G::T) where {P, T}
    bufD  = blockzeros(domainType(G),  size(G,2))
    bufC  = blockzeros(codomainType(G),size(G,1))
    bufC2 = blockzeros(codomainType(G),size(G,1))
    PrecomposeNonlinear{P, T, typeof(bufD), typeof(bufC)}(g, G, bufD, bufC, bufC2)
end

is_smooth(f::PrecomposeNonlinear) = is_smooth(f.g)

function (f::PrecomposeNonlinear)(x)
    return f.g(f.G*x)
end

function gradient!(y::D, f::PrecomposeNonlinear{P,T,D,C}, x::D) where {P,T,D,C}
    mul!(f.bufC, f.G, x)
    v = gradient!(f.bufC2, f.g, f.bufC)
    J = Jacobian(f.G, x)
    y = mul!(y, J', f.bufC2)
    return v
end


text:
using Distributions: Normal, Exponential, pdf
using LinearAlgebra

struct BeamPDF{T<:AbstractFloat} <: AbstractBeamPDF
    a_hit::T
    a_short_lin::T
    a_short_exp::T
    a_max::T
    a_rand::T
    sigma_hit::T
    lambda_short::T
    zmax::T

    function BeamPDF(a_hit, a_short_lin, a_short_exp, a_max, a_rand, sigma_hit, lambda_short, zmax)
        eta = a_short_lin + a_short_exp + a_max + a_rand + a_hit
        a_hit /= eta
        a_short_lin /= eta
        a_short_exp /= eta
        a_max /= eta
        a_rand /= eta
        new{Float64}(a_hit, a_short_lin, a_short_exp, a_max, a_rand, sigma_hit, lambda_short, zmax)
    end
end

function (m::BeamPDF)(z_W, zexp_W)
    # Not a valid measurement, but we don't want to collapse the distribution
    if z_W > m.zmax || z_W < 0 return m.a_rand/m.zmax end
    # p_hit
    p = pdf(Normal(zexp_W, m.sigma_hit), z_W) * m.a_hit
    # p_short_exp
    p += 0<=z_W<=zexp_W ? pdf(Exponential(1/m.lambda_short), z_W) * m.a_short_exp : 0
    # p_short_lin
    p += 0<=z_W<zexp_W ? 2 * m.a_short_lin * (zexp_W-z_W) / zexp_W : 0
    # p_max
    p += z_W == m.zmax ? m.a_max : 0
    # p_rand
    p += m.a_rand/m.zmax
    @assert !iszero(p) && !isnan(p)
    p
end

struct DiscBeamPDF{P, S, M} <: AbstractBeamPDF
    p_z_zexp_W::P # p_z_zexp_W[z, z_exp] = p(z | zexp_W)
    scale::S
    zmax::M
    arrsize::Int

    function DiscBeamPDF(scale_WM, zmax_W; a_hit=0.75, a_short_lin=0.1, a_short_exp=0.1, a_max=0.05, a_rand=0.1, sigma_hit=1.0,
        lambda_short=1/5)
        # +1 b/c z in inclusive range [0,zmax]
        dists_W = collect(LinRange(0, zmax_W, ceil(Int, zmax_W/scale_WM) + 1))
        p_z_zexp_W = zeros(Float64, length(dists_W), length(dists_W))
        model = BeamPDF(a_hit, a_short_lin, a_short_exp, a_max, a_rand, sigma_hit, lambda_short, zmax_W)
        for (i_zexp_W, zexp_W) in enumerate(dists_W)
            norm = 0.0
            for (i_z_W, z_W) in enumerate(dists_W)
                p = model(z_W, zexp_W)
                p_z_zexp_W[i_z_W, i_zexp_W] = p
                norm += p
            end
            p_z_zexp_W[:, i_zexp_W] ./= norm
        end
        @assert !any(isnan, p_z_zexp_W)
        @assert !any(iszero, p_z_zexp_W)
        new{typeof(p_z_zexp_W), typeof(scale_WM), typeof(zmax_W)}(p_z_zexp_W, 1/scale_WM, zmax_W, size(p_z_zexp_W,1))
    end
end

# In world units
function (bm::DiscBeamPDF)(z_W::Real, zexp_W::Real)
    z_W = round(Int, z_W * bm.scale + 1)
    zexp_W = round(Int, zexp_W * bm.scale + 1)
    z_W = clamp(z_W, 1, bm.arrsize)
    zexp_W = clamp(zexp_W, 1, bm.arrsize)
    @assert 1 <= z_W <= bm.arrsize
    @assert 1 <= zexp_W <= bm.arrsize
    @inbounds bm.p_z_zexp_W[z_W, zexp_W]
end


text:
export ignorefirstn, this, keepall, tonothing

# FIXME: Move these to MH common

"`burnin(n)` keep function that ignores first `n` samples"
ignorefirstn(n) = i -> i > n

"`thin(n)` keep function that accept only every `n` samples"
thin(n) = i -> i % n == 0

@inline keepall(i) = true

"`a &ₚ b` Pointwise logical `x->a(x) & b(x)`"
a &ₚ b = (x...)->a(x...) &ₚ b(x...)

tonothing(args...) = nothing

struct Capture{T}
  x::T
end

(c::Capture)(data) = push!(c.x, data)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using DataStructures: Deque

"""
    TimeWindowAssociativeOp{Value,Op,Op!,Time}(window::TimeDiff)
    TimeWindowAssociativeOp{Value,Op,Time}(window::TimeDiff)

State necessary for accumulation over a rolling window of fixed size, in terms of time.

When presented with a new time t', we guarantee that all times t remaining in the window
satisfy:

    t > t' - w

That is, at time t' this window represents the open-closed time interval (t' - w, t']

We require that `window` be of a type that, when added to a `Time`, gives a `Time`.

# Fields
- `window_state::WindowedAssociativeOp{Value}`: The underlying general-window state.
- `window::TimeDiff`: The window, as a difference between two times.
- `times::Deque{Time}`: The same length as the values stored in `window_state`, and
    representing the times of those observations.
- `window_full::Bool`: For internal use - will be set to true once a point has dropped out
    of the window.
"""
mutable struct TimeWindowAssociativeOp{Value,Op,Op!,Time,TimeDiff}
    window_state::WindowedAssociativeOp{Value,Op,Op!}
    window::TimeDiff
    times::Deque{Time}
    window_full::Bool

    function TimeWindowAssociativeOp{Value,Op,Op!,Time}(
        window::TimeDiff
    ) where {Value,Op,Op!,Time,TimeDiff}
        # Verify that TimeDiff and Time are compatible.
        ret_types = Base.return_types(+, (Time, TimeDiff))
        isempty(ret_types) && throw(ArgumentError("Incompatible: $Time and $TimeDiff"))
        only(ret_types) == Time || throw(ArgumentError("Incompatible: $Time and $TimeDiff"))

        if window <= zero(TimeDiff)
            throw(ArgumentError("Got window $window, but it must be positive."))
        end
        return new{Value,Op,Op!,Time,TimeDiff}(
            WindowedAssociativeOp{Value,Op,Op!,}(), window, Deque{Time}(), false
        )
    end
end

function TimeWindowAssociativeOp{T,Op,Time}(window::TimeDiff) where {T,Op,Time,TimeDiff}
    return TimeWindowAssociativeOp{T,Op,Op,Time}(window)
end

"""
    update_state!(state::TimeWindowAssociativeOp, time, value) -> state

Add the specified `value` to the state with associated `time`, and drop any values that
are no longer in the time window.

# Arguments
- `state::TimeWindowAssociativeOp`:
- `time`: The time to which `value` corresponds.
- `value`: The value to add to the window.

# Returns
- `::TimeWindowAssociativeOp`: `state`, which has been mutated.
"""
function update_state!(state::TimeWindowAssociativeOp, time, value)
    if !isempty(state.times) && time <= last(state.times)
        throw(
            ArgumentError(
                "Got out-of-order time $time. Previous time was $(last(state.times))"
            ),
        )
    end

    push!(state.times, time)

    # Drop off times from the front of the deque, keeping track of how many values we need
    # to remove from the window state.
    # This is a linear search
    num_dropped_from_window = 0
    while (first(state.times) + state.window) <= time
        popfirst!(state.times)
        num_dropped_from_window += 1
    end

    update_state!(state.window_state, value, num_dropped_from_window)

    if !state.window_full && num_dropped_from_window > 0
        # The window has now filled; record this for use in `window_full` below.
        state.window_full = true
    end

    return state
end

window_value(state::TimeWindowAssociativeOp) = window_value(state.window_state)
window_size(state::TimeWindowAssociativeOp)::Int = window_size(state.window_state)

"""
    window_full(state::TimeWindowAssociativeOp)::Bool

Returns true iff the given `state` has had at least one value drop out of the window,
indicating that the window is now full.
"""
window_full(state::TimeWindowAssociativeOp)::Bool = state.window_full


text:
export
	ConstraintVals

@with_kw mutable struct ConstraintParams{T}
	ϕ::T = 10.0  	  # penalty scaling parameter
	μ0::T = 1.0 	  # initial penalty parameter
	μ_max::T = 1e8    # max penalty parameter
	λ_max::T = 1e8    # max Lagrange multiplier
end

""" $(TYPEDEF)
Struct that stores all of the values associated with a particular constraint.
Importantly, `ConstraintVals` stores the list of knotpoints to which the constraint
is applied. This type should be fairly transparent to the user, and only needs to be
directly dealt with when writing solvers or setting fine-tuned updates per constraint
(via the `.params` field).
"""
struct ConstraintVals{T,W,C,P,N}
	con::C
	inds::UnitRange{Int}
	vals::Vector{SVector{P,T}}
	vals_prev::Vector{SVector{P,T}}
	∇c::Vector{SizedMatrix{P,N,T,2}}
	λ::Vector{SVector{P,T}}
	μ::Vector{SVector{P,T}}
	active::Vector{SVector{P,Bool}}
	c_max::Vector{T}
	params::ConstraintParams{T}

	function ConstraintVals(con::AbstractConstraint{S, W},
			inds::UnitRange{Int}, vals::V, vals_prev,
			∇c::Vector{<:SizedMatrix{P,N}}, λ::V, μ::V,
			active::Vector{SVector{P,Bool}}, c_max::Vector{T},
			params::ConstraintParams) where {S,W,T,P,N,V}
		new{T,W,typeof(con),P,N}(con,inds,vals,vals_prev,∇c,λ,μ,
			active,c_max, params)
	end
end

function ConstraintVals(con::C, inds::UnitRange, kwargs...) where C
	p = length(con)
	w = width(con)
	P = length(inds)
	λ    = [@SVector zeros(p) for k = 1:P]
	μ    = [@SVector ones(p)  for k = 1:P]
	atv  = [@SVector ones(Bool,p) for k = 1:P]
	vals = [@SVector zeros(p) for k = 1:P]
	∇c = [SizedMatrix{p,w}(zeros(Float64,p,w)) for k = 1:P]
	params = ConstraintParams(;kwargs...)
	ConstraintVals(con, inds, vals, deepcopy(vals), ∇c, λ, μ, atv, zeros(P),
		params)
end

get_params(con::ConstraintVals)::ConstraintParams = con.params

function _index(con::ConstraintVals, k::Int)
	if k ∈ con.inds
		return k - con.inds[1] + 1
	else
		return 0
	end
end

Base.length(::ConstraintVals{T,W,C,P}) where {T,W,C,P} = P
Base.length(con::ConstraintVals, k::Int) = k ∈ con.inds ? length(con) : 0
state_dim(con::ConstraintVals) = state_dim(con.con)
control_dim(con::ConstraintVals) = control_dim(con.con)
check_dims(con::ConstraintVals,n,m) = check_dims(con.con,n,m)
constraint_type(con::ConstraintVals{T,W,C}) where {T,W,C} = C
is_bound(con::ConstraintVals) = is_bound(con.con)
duals(con::ConstraintVals) = con.λ
duals(con::ConstraintVals, k::Int) = con.λ[_index(con,k)]
penalty(con::ConstraintVals) = con.μ
penalty(con::ConstraintVals, k::Int) = con.μ[_index(con,k)]
penalty_matrix(con::ConstraintVals, i::Int) = Diagonal(con.active[i] .* con.μ[i])
lower_bound(con::ConstraintVals) = lower_bound(con.con)
upper_bound(con::ConstraintVals) = upper_bound(con.con)
contype(con::ConstraintVals) = contype(con.con)
sense(con::ConstraintVals) = sense(con.con)

evaluate!(con::ConstraintVals, Z::Traj) = evaluate!(con.vals, con.con, Z, con.inds)
jacobian!(con::ConstraintVals, Z::Traj) = jacobian!(con.∇c, con.con, Z, con.inds)

function update_active_set!(con::ConstraintVals{T,W,C}, ::Val{tol}) where
		{T,W,C<:AbstractConstraint{Inequality},tol}
	for i in eachindex(con.vals)
		con.active[i] = @. (con.vals[i] >= -tol) | (con.λ[i] > 0)
	end
	return nothing
end

update_active_set!(con::ConstraintVals{T,W,C}, ::Val{tol}) where
	{T,W,C<:AbstractConstraint{Equality},tol} = nothing

function viol_ineq(v::T, a)::T where T
	for i in eachindex(a)
		v = max(v, max(a[i], 0.0))
	end
	return v
end

function viol_eq(v::T, a)::T where T
	for i in eachindex(a)
		v = max(v, abs(a[i]))
	end
	return v
end

function violation(::Type{Inequality}, vals::SVector)
	return max.(vals, 0)
end

function violation(::Type{Equality}, vals::SVector)
	return abs.(vals)
end

function max_violation!(con::ConstraintVals{T,W,C}) where
		{T,W,C<:AbstractConstraint{Inequality}}
	for i in eachindex(con.c_max)
		con.c_max[i] = viol_ineq(0.0, con.vals[i])
	end
	return nothing
end

function max_violation!(con::ConstraintVals{T,W,C}) where
		{T,W,C<:AbstractConstraint{Equality}}
	for i in eachindex(con.c_max)
		# con.c_max[i] = norm(con.vals[i],Inf)
		con.c_max[i] = viol_eq(0.0, con.vals[i])
	end
	return nothing
end

function max_penalty!(con::ConstraintVals)
	for i in eachindex(con.c_max)
		con.c_max[i] = maximum(con.μ[i])
	end
	return nothing
end

get_c_max(con::ConstraintVals) = maximum(con.c_max)

function reset!(con::ConstraintVals{T,W,C,P}) where {T,W,C,P}
	reset_duals!(con)
	reset_penalties!(con)
end

function reset_duals!(con::ConstraintVals)
	λ = con.λ
	for i in eachindex(con.inds)
		λ[i] *= 0.0
	end
end

function reset_penalties!(con::ConstraintVals{T,W,C,P}) where {T,W,C,P}
	c = con.vals
	μ = con.μ
	for i in eachindex(con.inds)
		μ[i] = con.params.μ0 * @SVector ones(T,P)
		c[i] *= 0.0
	end
end

function shift_fill!(con::ConstraintVals)
	shift_fill!(con.μ)
	shift_fill!(con.λ)
end

function cost!(J, con::ConstraintVals, Z)
	for (i,k) in enumerate(con.inds)
		c = con.vals[i]
		λ = con.λ[i]
		Iμ = penalty_matrix(con, i)
		J[k] += λ'c + 0.5*c'Iμ*c
	end
end

@generated function cost_expansion!(E::Vector{<:AbstractExpansion}, con::ConstraintVals{T,W},
		Z::Vector{<:KnotPoint{T,N,M}}) where {T,W<:Stage,N,M}
	if W <: State
		expansion = quote
			cx = SMatrix(con.∇c[i])
			E[k].xx .+= cx'Iμ*cx
			E[k].x .+= cx'g
		end
	elseif W <: Control
		expansion = quote
			cu = SMatrix(con.∇c[i])
			E[k].uu .+= cu'Iμ*cu
			E[k].u .+= cu'g
		end
	else
		expansion = quote
			cx = SMatrix(con.∇c[i])[:,ix]
			cu = SMatrix(con.∇c[i])[:,iu]

			E[k].xx .+= cx'Iμ*cx
			E[k].uu .+= cu'Iμ*cu
			E[k].ux .+= cu'Iμ*cx

			E[k].x .+= cx'g
			E[k].u .+= cu'g
		end
	end
	quote
		ix,iu = Z[1]._x, Z[1]._u
		@inbounds for i in eachindex(con.inds)
			k = con.inds[i]
			c = con.vals[i]
			λ = con.λ[i]
			μ = con.μ[i]
			a = con.active[i]
			Iμ = Diagonal( a .* μ )
			g = Iμ*c + λ

			$expansion
		end
	end
end


text:
module FastDEQ

using CUDA, DiffEqBase, DiffEqCallbacks, DiffEqSensitivity, Flux, FluxExperimental, LinearAlgebra, LinearSolve,
      OrdinaryDiffEq, SciMLBase, Statistics, SteadyStateDiffEq, UnPack, Zygote

abstract type AbstractDeepEquilibriumNetwork end

function Base.show(io::IO, l::AbstractDeepEquilibriumNetwork)
    return print(io, string(typeof(l).name.name), "() ", string(length(l.p)), " Trainable Parameters")
end

Flux.trainable(d::AbstractDeepEquilibriumNetwork) = (d.p,)

Base.deepcopy(op::DiffEqSensitivity.ZygotePullbackMultiplyOperator) = op

include("solve.jl")
include("utils.jl")

include("solvers/broyden.jl")
include("solvers/limited_memory_broyden.jl")

include("models/basics.jl")

include("layers/jacobian_stabilization.jl")
include("layers/utils.jl")
include("layers/deq.jl")
include("layers/sdeq.jl")
include("layers/mdeq.jl")
include("layers/smdeq.jl")

include("models/chain.jl")

include("losses.jl")

# DEQ Solvers
export ContinuousDEQSolver, DiscreteDEQSolver, BroydenSolver, LimitedMemoryBroydenSolver

# Utils
export NormalInitializer, SteadyStateAdjoint, get_and_clear_nfe!, compute_deq_jacobian_loss, DeepEquilibriumSolution, SupervisedLossContainer

# Layers
export MultiParallelNet

# DEQ Layers
export DeepEquilibriumNetwork, SkipDeepEquilibriumNetwork, MultiScaleDeepEquilibriumNetwork, MultiScaleSkipDeepEquilibriumNetwork
export DEQChain

end


text:
g = NEAT.Global(config)
g.cf.feedforward = false
p = NEAT.Population(g)
@test all(ch-> ch.node_gene_type == NEAT.Recurrent(), p.population)

g.cf.feedforward = true
p = NEAT.Population(g)
@test all(ch-> ch.node_gene_type == NEAT.FeedForward(), p.population)

g.cf.fully_connected = true
g.cf.feedforward = true
g.cf.hidden_nodes = 0
p = NEAT.Population(g)
@test all(ch-> length(ch.connection_genes) == 3, p.population)

g.cf.fully_connected = true
g.cf.feedforward = true
g.cf.hidden_nodes = 2
p = NEAT.Population(g)
@test all(ch-> length(ch.connection_genes) == 11, p.population)

g.cf.fully_connected = true
g.cf.feedforward = false
g.cf.hidden_nodes = 2
p = NEAT.Population(g)
@test all(ch-> length(ch.connection_genes) == 17, p.population)

g.cf.fully_connected = false
g.cf.feedforward = false
g.cf.hidden_nodes = 2
p = NEAT.Population(g)
@test all(ch-> length(ch.connection_genes) == 14, p.population)

g.cf.fully_connected = false
g.cf.feedforward = true
g.cf.hidden_nodes = 2
p = NEAT.Population(g)
@test all(ch-> length(ch.connection_genes) == 8, p.population)

# remove tests
p = NEAT.Population(g)
map(x->NEAT.remove(p, p.population[1]), 1:5)
@test length(p.population) == g.cf.pop_size - 5

g = NEAT.Global(config)
p = NEAT.Population(g)
map(ch->ch.fitness = rand(), p.population)
NEAT.speciate(g, p, true)
NEAT.compute_spawn_levels(g, p)
totalSpawns = sum([s.spawn_amount for s in p.species])
@test length(p.population)-1 <= totalSpawns <= length(p.population)+1

g = NEAT.Global(config)
g.cf.prob_addnode = 0.5
p = NEAT.Population(g)
p.evaluate = (chs)->map(ch->ch.fitness=max(0,randn()*0.3+ch.fitness*0.7), chs)
@time NEAT.epoch(g, p, 1, false)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

@testset "string normalization" begin
    # normalize_string (Unicode normalization etc.):
    @test normalize_string("\u006e\u0303", :NFC) == "\u00f1"
    @test "\u006e\u0303" == normalize_string("\u00f1", :NFD)
    @test normalize_string("\ufb00", :NFC) != "ff"
    @test normalize_string("\ufb00", :NFKC) == "ff"
    @test normalize_string("\u006e\u0303\ufb00", :NFKC) == "\u00f1"*"ff"
    @test normalize_string("\u00f1\ufb00", :NFKD) == "\u006e\u0303"*"ff"
    @test normalize_string("\u006e\u0303", compose=true) == "\u00f1"
    @test "\u006e\u0303" == normalize_string("\u00f1", decompose=true)
    @test normalize_string("\u006e\u0303\u00b5",compat=true) == "\u00f1\u03bc"
    @test normalize_string("Σσς",casefold=true) == "σσσ"
    @test normalize_string("∕⁄", lump=true) == "//"
    @test normalize_string("\ua\n\r\r\ua", newline2lf=true) == "\ua\ua\ua\ua"
    @test normalize_string("\ua\n\r\r\ua", newline2ls=true) == "\u2028\u2028\u2028\u2028"
    @test normalize_string("\ua\n\r\r\ua", newline2ps=true) == "\u2029\u2029\u2029\u2029"
    @test normalize_string("\u00f1", stripmark=true) == "n"
    @test isempty(normalize_string("\u00ad", stripignore=true))
    @test normalize_string("\t\r", stripcc=true) == "  "
    @test normalize_string("\t\r", stripcc=true, newline2ls=true) == " \u2028"
end

@testset "unicode sa#15" begin
    #Tests from Unicode SA#15, "Unicode normalization forms"
    #http://www.unicode.org/reports/tr15/

    @testset "canonical equivalence" begin
        let ==(a::Array{Char},b::Array{Char}) = normalize_string(string(a...), :NFC)==normalize_string(string(b...), :NFC)
            ==(a,b) = Base.:(==)(a,b)
            @test ['C', '̧'] == ['Ç']
            @test ['q', '̇', '̣'] == ['q', '̣', '̇']
            @test ['가'] == ['ᄀ', 'ᅡ']
            @test ['Ω'] == ['Ω']
        end
    end

    @testset "compatibility equivalence" begin
        let ==(a::Array{Char},b::Array{Char}) = normalize_string(string(a...), :NFKC)==normalize_string(string(b...), :NFKC)
            ==(a,b) = Base.:(==)(a,b)
            @test ['ℌ'] == ['ℍ'] == ['H']
            @test ['ﻨ'] == ['ﻧ'] == ['ﻦ'] == ['ﻥ']
            @test ['①'] == ['1']
            @test ['ｶ'] == ['カ']
            @test ['︷'] == ['{']
            @test ['⁹'] == ['₉']
            @test ['㌀'] == ['ア', 'パ', 'ー', 'ト']
            @test ['¼'] == ['1', '⁄', '4']
            @test ['ǆ'] == ['d', 'ž']
        end
    end

    @testset "singletons" begin
        @test normalize_string("\U212b", :NFD) == "A\U030a"
        @test normalize_string("\U212b", :NFC) == "\U00c5"
        @test normalize_string("\U2126", :NFC) == normalize_string("\U2126", :NFD) == "\U03a9"
    end

    @testset "canonical composites" begin
        @test normalize_string("\U00c5", :NFC) == "\U00c5"
        @test normalize_string("\U00c5", :NFD) == "A\U030a"
        @test normalize_string("\U00f4", :NFC) == "\U00f4"
        @test normalize_string("\U00f4", :NFD) == "o\U0302"
    end

    @testset "multiple combining marks" begin
        @test normalize_string("\U1e69", :NFD) == "s\U0323\U0307"
        @test normalize_string("\U1e69", :NFC) == "\U1e69"
        @test normalize_string("\U1e0b\U0323", :NFD) == "d\U0323\U0307"
        @test normalize_string("\U1e0b\U0323", :NFC) == "\U1e0d\U0307"
        @test normalize_string("q\U0307\U0323", :NFC) == "q\U0323\U0307"
        @test normalize_string("q\U0307\U0323", :NFD) == "q\U0323\U0307"
    end

    @testset "compatibility composites" begin
        @test normalize_string("\Ufb01", :NFD) == normalize_string("\Ufb01", :NFC) == "\Ufb01"
        @test normalize_string("\Ufb01", :NFKD) == normalize_string("\Ufb01", :NFKC) == "fi"
        @test normalize_string("2\U2075", :NFD) == normalize_string("2\U2075", :NFC) == "2\U2075"
        @test normalize_string("2\U2075", :NFKD) == normalize_string("2\U2075", :NFKC) == "25"
        @test normalize_string("\U1e9b\U0323", :NFD) == "\U017f\U0323\U0307"
        @test normalize_string("\U1e9b\U0323", :NFC) == "\U1e9b\U0323"
        @test normalize_string("\U1e9b\U0323", :NFKD) == "s\U0323\U0307"
        @test normalize_string("\U1e9b\U0323", :NFKC) == "\U1e69"
    end
end

@testset "#5939 uft8proc character predicates" begin
    alower=['a', 'd', 'j', 'y', 'z']
    ulower=['α', 'β', 'γ', 'δ', 'ф', 'я']
    for c in vcat(alower,ulower)
        @test islower(c) == true
        @test isupper(c) == false
        @test isdigit(c) == false
        @test isnumber(c) == false
    end

    aupper=['A', 'D', 'J', 'Y', 'Z']
    uupper= ['Δ', 'Γ', 'Π', 'Ψ', 'ǅ', 'Ж', 'Д']

    for c in vcat(aupper,uupper)
        @test islower(c) == false
        @test isupper(c) == true
        @test isdigit(c) == false
        @test isnumber(c) == false
    end

    nocase=['א','ﺵ']
    alphas=vcat(alower,ulower,aupper,uupper,nocase)

    for c in alphas
        @test isalpha(c) == true
        @test isnumber(c) == false
    end

    anumber=['0', '1', '5', '9']
    unumber=['٣', '٥', '٨', '¹', 'ⅳ' ]

    for c in anumber
        @test isdigit(c) == true
        @test isnumber(c) == true
    end
    for c in unumber
        @test isdigit(c) == false
        @test isnumber(c) == true
    end

    alnums=vcat(alphas,anumber,unumber)
    for c in alnums
        @test isalnum(c) == true
        @test ispunct(c) == false
    end

    asymbol = ['(',')', '~', '$' ]
    usymbol = ['∪', '∩', '⊂', '⊃', '√', '€', '¥', '↰', '△', '§']

    apunct =['.',',',';',':','&']
    upunct =['‡', '؟', '჻' ]

    for c in vcat(apunct,upunct)
        @test ispunct(c) == true
        @test isalnum(c) == false
    end

    for c in vcat(alnums,asymbol,usymbol,apunct,upunct)
        @test isprint(c) == true
        @test isgraph(c) == true
        @test isspace(c) == false
        @test iscntrl(c) == false
    end

    NBSP = Char(0x0000A0)
    ENSPACE = Char(0x002002)
    EMSPACE = Char(0x002003)
    THINSPACE = Char(0x002009)
    ZWSPACE = Char(0x002060)

    uspace = [ENSPACE, EMSPACE, THINSPACE]
    aspace = [' ']
    acntrl_space = ['\t', '\n', '\v', '\f', '\r']
    for c in vcat(aspace,uspace)
        @test isspace(c) == true
        @test isprint(c) == true
        @test isgraph(c) == false
    end

    for c in vcat(acntrl_space)
        @test isspace(c) == true
        @test isprint(c) == false
        @test isgraph(c) == false
    end

    @test isspace(ZWSPACE) == false # zero-width space

    acontrol = [ Char(0x001c), Char(0x001d), Char(0x001e), Char(0x001f)]
    latincontrol = [ Char(0x0080), Char(0x0085) ]
    ucontrol = [ Char(0x200E), Char(0x202E) ]

    for c in vcat(acontrol, acntrl_space, latincontrol)
        @test iscntrl(c) == true
        @test isalnum(c) == false
        @test isprint(c) == false
        @test isgraph(c) == false
    end

    for c in ucontrol  #non-latin1 controls
        if c!=Char(0x0085)
            @test iscntrl(c) == false
            @test isspace(c) == false
            @test isalnum(c) == false
            @test isprint(c) == false
            @test isgraph(c) == false
        end
    end

    @test  all(isspace,"  \t   \n   \r  ")
    @test !all(isgraph,"  \t   \n   \r  ")
    @test !all(isprint,"  \t   \n   \r  ")
    @test !all(isalpha,"  \t   \n   \r  ")
    @test !all(isnumber,"  \t   \n   \r  ")
    @test !all(ispunct,"  \t   \n   \r  ")

    @test !all(isspace,"ΣβΣβ")
    @test  all(isalpha,"ΣβΣβ")
    @test  all(isgraph,"ΣβΣβ")
    @test  all(isprint,"ΣβΣβ")
    @test !all(isupper,"ΣβΣβ")
    @test !all(islower,"ΣβΣβ")
    @test !all(isnumber,"ΣβΣβ")
    @test !all(iscntrl,"ΣβΣβ")
    @test !all(ispunct,"ΣβΣβ")

    @test  all(isnumber,"23435")
    @test  all(isdigit,"23435")
    @test  all(isalnum,"23435")
    @test !all(isalpha,"23435")
    @test  all(iscntrl,string(Char(0x0080)))
    @test  all(ispunct, "‡؟჻")

    @test  isxdigit('0')
    @test  isxdigit('a')
    @test !isxdigit('x')
    @test !isxdigit('g')
end

@testset "utf8proc" begin
    # check utf8proc handling of CN category constants
    let c_ll = 'β', c_cn = '\u038B'
        @test Base.UTF8proc.category_code(c_ll) == Base.UTF8proc.UTF8PROC_CATEGORY_LL
        # check codepoint with category code CN
        @test Base.UTF8proc.category_code(c_cn) == Base.UTF8proc.UTF8PROC_CATEGORY_CN
    end
end

@testset "graphemes" begin
    let grphtest = (("b\u0300lahβlahb\u0302láh", ["b\u0300","l","a","h",
                                                  "β","l","a","h",
                                                  "b\u0302","l","á","h"]),
                    ("", String[]),
                    ("x\u0302", ["x\u0302"]),
                    ("\U1d4c1\u0302", ["\U1d4c1\u0302"]),
                    ("\U1d4c1\u0302\U1d4c1\u0300", ["\U1d4c1\u0302",
                                                    "\U1d4c1\u0300"]),
                    ("x",["x"]),
                    ("abc",["a","b","c"]))
        for T in (String,GenericString)
            for nf in (:NFC, :NFD)
                for (s, g) in grphtest
                    s_ = T(normalize_string(s, nf))
                    g_ = map(s -> normalize_string(s, nf), g)
                    # #9261
                    if length(s_) > 0
                        @test typeof(first(graphemes(s_))) == SubString{typeof(s_)}
                    end
                    grph = collect(graphemes(s_))
                    @test eltype(grph) == SubString{typeof(s_)}
                    @test grph == g_
                    @test length(graphemes(s_)) == length(grph)
                end
                S = [T(normalize_string(s)) for (s,g) in grphtest]
                G = map(graphemes, S)
                @test map(graphemes, sort!(S)) == sort!(G)
            end
        end
    end
end

@testset "#3721, #6939 up-to-date character widths" begin
    @test charwidth('\U1f355') == 2
    @test strwidth("\U1f355") == 2
    @test strwidth(GenericString("\U1f355")) == 2
    @test strwidth("\U1f355\u0302") == 2
    @test strwidth(GenericString("\U1f355\u0302")) == 2
end

@testset "#10958 handling of embedded NUL chars" begin
    @test length("\0w") == length("\0α") == 2
    @test strwidth("\0w") == strwidth("\0α") == 1
    @test normalize_string("\0W", casefold=true) == "\0w"
end

@testset "ut8proc_map with GenericString" begin
    @test normalize_string(GenericString("\u006e\u0303"), :NFC) == "\u00f1"
end

@testset "normalize_string keywords" begin
    @test_throws ArgumentError normalize_string("\u006e\u0303", compose=false, compat=true)
    @test_throws ArgumentError normalize_string("\u006e\u0303", compose=false, stripmark=true)
end

@testset "fastplus" begin
    @test lowercase('A') == 'a'
    @test uppercase('a') == 'A'

    @test is_assigned_char('A')
end

@testset "isspace" begin
    @test isspace(' ')
    @test isspace('\t')
    @test isspace('\r')
    @test isspace('\u85')
    @test isspace('\ua0')
    @test !isspace('\ufffd')
    @test !isspace('\U10ffff')
end

@testset "grapheme iterators" begin
    let str = ascii("This is a test")
        g = graphemes(str)
        h = hash(str)
        @test hash(g) == h
        @test convert(GenericString, g) == str
        io = IOBuffer()
        show(io, g)
        check = "length-14 GraphemeIterator{String} for \"$str\""
        @test String(take!(io)) == check
    end
end

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
using Knet, AutoGrad
using Knet: cudnnSoftmaxForward, cudnnSoftmaxBackward

#=

"""
    logp(x, dims=1)

Treat entries in `x` as as unnormalized log probabilities and return
normalized log probabilities.
`dims` is an optional argument, if not specified the normalization is
over the first dimension of `x`, otherwise the normalization is performed over the
given dimensions.  In particular, if `x` is a matrix, `dims=1`
normalizes columns of `x`, `dims=2` normalizes rows of `x` and
`dims=(1,2)` normalizes the whole matrix.
Calls to `logsoftmax` are equivalent to `logp`, and `softmax(x,dims)` is
equivalent to `exp.(logp(x,dims)`.
"""
logp(x, dims=1) = _logp(x, dims) # generic fallback

function logp(x::A, dims=1) where A <: Union{<:KnetArray, Rec{<:KnetArray}}
    d = sort(union(dims))
    if ndims(x) == length(d)
        n = length(x)
        if n > 20000
            _logp(x, dims)
        else
            sz = size(x)
            x = cudnnSoftmaxForward(reshape(x,(1,1,n,1)),algo=2)
            reshape(x, sz)
        end
    elseif d == [1]
        sz = size(x)
        x = cudnnSoftmaxForward(reshape(x, (1,1,sz[1],:)), algo=2)
        reshape(x, sz)
    elseif ndims(x) == 2 && d == [2]
        logp(x.', 1).'
    else
        _logp(x, dims)
    end
end

# Math for the cross-entropy loss: x is unnormalized input, p is
# target probabilities, q is estimated probabilities. Read left column
# down, right column (loss gradients) back up.

# x			dx = -p + qz/z = -p + exp(logq)
# xmax  = max(x,1)	-sum(db)=0
# logqz = x .- xmax	-p + qz/z
# qz    = exp(logqz)	rep(1/z)
# z     = sum(qz,1)	1/z
# logz  = log(z)	sum(p)=1
# logq  = logqz.-logz	-p
# plogq = p .* logq	-1
# loss  = -sum(plogq)	1

# We keep the old implementation _logp for CPU arrays, slow cases and
# cases of d not handled by cudnn.

function _logp(x, dims=1)
    xval = getval(x)
    if isa(xval,Number)
        return zero(xval)
    elseif isempty(xval)
        return xval
    else
        x = x .- maximum(x, dims)
        return (x .- log.(sum(exp.(x), dims)))
        # Expanding for profiling:
        # x1 = maximum(x,dims)
        # x2 = x .- x1
        # x3 = exp.(x2)
        # x4 = sum(x3,d...)
        # x5 = log.(x4)
        # x6 = x2 .- x5
        # return x6
    end
end

function _logpback(x,y,dy,dims)
    xval = getval(x)
    if isa(xval,Number)
        return zero(xval)
    elseif isempty(xval)
        return xval
    else
        return dy .- exp.(y).*sum(dy, dims)
        # Expanding for profiling:
        # dx1 = sum(dy,dims)
        # dx2 = exp.(y)
        # dx3 = dx2 .* dx1
        # dx4 = dy - dx3
        # return dx4
    end
end

# dy should be -p and y=logq so this should give us -p+q
@primitive  _logp(x,dims),dy,y  _logpback(x,y,dy,dims)

"""
    logsoftmax(x, dims=1)

Equivalent to `logp(x, dims)`. See also `sotfmax`.
"""
const logsoftmax = logp

=#

"""
    softmax(x, dims=1; algo=1)

The softmax function typically used in classification.
Gives the same results as to `exp.(logp(x, dims))`.

If `algo=1` computation is more accurate, if `algo=0` it is
faster.

See also `logsoftmax`.
"""
softmax(x, dims=1; algo=1) = _softmax(x, dims; algo=algo) # generic fallback

function softmax(x::A, dims=1; algo=1) where A <: Union{<:KnetArray, Rec{<:KnetArray}}
    @assert algo ∈ [0, 1]
    d = sort(union(dims))
    if ndims(x) == length(d)
        sz = size(x)
        x = cudnnSoftmaxForward(reshape(x,(1,1,:,1)),algo=algo)
        reshape(x, sz)
    elseif d == [1]
        sz = size(x)
        x = cudnnSoftmaxForward(reshape(x, (1,1,sz[1],:)), algo=algo)
        reshape(x, sz)
    elseif ndims(x) == 2 && d == [2]
        softmax(x.', 1, algo=algo).'
    else
        _softmax(x, dims)
    end
end

function _softmax(x, dims; algo=1)
    @assert algo ∈ [0, 1]
    if algo == 1
        x = x .- maximum(x, dims)
    end
    x = exp.(x)
    return x ./ sum(x, dims)
end

function _softback(x,y,dy,dims)
    return y .* dy .- y .* sum(y .* dy, dims)
end

#=
typedef enum
{
    CUDNN_SOFTMAX_FAST     = 0,         /* straightforward implementation */
    CUDNN_SOFTMAX_ACCURATE = 1,         /* subtract max from every point to avoid overflow */
    CUDNN_SOFTMAX_LOG      = 2
} cudnnSoftmaxAlgorithm_t;
typedef enum
{
    CUDNN_SOFTMAX_MODE_INSTANCE = 0,   /* compute the softmax over all C, H, W for each N */
    CUDNN_SOFTMAX_MODE_CHANNEL = 1     /* compute the softmax over all C for each H, W, N */
} cudnnSoftmaxMode_t;
=#

#=
function cudnnSoftmaxForward{T}(x::KnetArray{T}; algo=0, mode=0, alpha=1, handle=cudnnhandle())
    beta = 0 # nonzero beta does not make sense when we create y
    y = similar(x)
    @cuda(cudnn, cudnnSoftmaxForward,
          (Cptr, Cint, Cint, Ptr{T}, Cptr, Ptr{T}, Ptr{T}, Cptr, Ptr{T}),
          handle, algo, mode, Ref(T(alpha)), TD(x), x, Ref(T(beta)), TD(y), y)
    return y
end

function cudnnSoftmaxBackward{T}(y::KnetArray{T}, dy::KnetArray{T}; algo=0, mode=0, alpha=1, handle=cudnnhandle())
    beta = 0
    dx = similar(dy)
    @cuda(cudnn, cudnnSoftmaxBackward,
          (Cptr, Cint, Cint, Ptr{T}, Cptr, Ptr{T}, Cptr, Ptr{T}, Ptr{T}, Cptr, Ptr{T}),
          handle, algo, mode, Ref(T(alpha)), TD(y), y, TD(dy), dy, Ref(T(beta)), TD(dx), dx)
    return dx
end

@primitive cudnnSoftmaxForward(x;o...),dy,y cudnnSoftmaxBackward(y,dy;o...)
@zerograd cudnnSoftmaxBackward(y,dy;o...)

"""
    logsumexp(x,[dims])

Compute `log(sum(exp(x),dims))` in a numerically stable manner.
`dims` is an optional argument, if not specified the summation is over
the whole `x`, otherwise the summation is performed over the given
dimensions.  In particular if `x` is a matrix, `dims=1` sums columns
of `x` and `dims=2` sums rows of `x`.
"""
function logsumexp(x,d...)
    xmax = maximum(x,d...)
    xmax + log.(sum(exp.(x .- xmax),d...))
end

@primitive logsumexp(x,d...),dy,y  (dy .* exp.(x .- y))

"""
    nll(scores, answers, d=1; average=true)

Given an unnormalized `scores` matrix and an `Integer` array of
correct `answers`, return the per-instance negative log
likelihood. `d=1` means instances are in columns, `d=2` means
instances are in rows.  Use `average=false` to return the sum instead
of per-instance average.
"""
function nll{T<:Integer}(y,a::Array{T},d=1; average=true)
    indices = findindices(y,a,d)
    lp = logp(y,d)[indices]
    average ? -mean(lp) : -sum(lp)
end

"""
    accuracy(scores, answers, d=1; average=true)

Given an unnormalized `scores` matrix and an `Integer` array of
correct `answers`, return the ratio of instances where the correct
answer has the maximum score. `d=1` means instances are in columns,
`d=2` means instances are in rows. Use `average=false` to return
the number of correct answers instead of the ratio.
"""
function accuracy{T<:Integer}(y,a::Array{T},d=1; average=true)
    indices = findindices(y,a,d)
    (maxval,maxind) = findmax(Array(y),d)
    correct = (vec(maxind) .== indices)
    average ? mean(correct) : sum(correct)
end

function findindices{T<:Integer}(y,a::Array{T},d=1)
    n = length(a)
    indices = Vector{Int}(n)
    if d == 1                   # instances in first dimension
        y1 = size(y,1)
        y2 = div(length(y),y1)
        if n != y2; throw(DimensionMismatch()); end
        @inbounds for j=1:n
            indices[j] = (j-1)*y1 + a[j]
        end
    elseif d == 2               # instances in last dimension
        y2 = size(y,ndims(y))
        y1 = div(length(y),y2)
        if n != y1; throw(DimensionMismatch()); end
        @inbounds for j=1:n
            indices[j] = (a[j]-1)*y1 + j
        end
    else
        error("findindices only supports d = 1 or 2")
    end
    return indices
end

=#

# # The xentloss interface is no good because of double normalization.

# """
# xentloss(x, p [,dims])

# Compute cross-entropy loss for unnormalized log probability estimates
# x and normalized probabilities p normalizing over the given
# dimensions.  By default normalization is over the whole array, dims=1
# normalizes over the columns, dims=2 normalizes over the rows of a 2-D
# array.

# """
# function xentloss(x,p,d...)
#     x = x .- maximum(x,d...)
#     z = log(sum(exp(x),d...))
#     return sum(p .* (x .- z)) / length(z)
# end

# function xentback(x,p,d...)
#     x = x .- maximum(x,d...)
#     x = exp(x)
#     z = sum(x,d...)
#     return x./z - p
# end

# @primitive xentloss(x,p,d...),dy,y  (dy.*xentback(x,p,d...))

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
import ParserCombinator: PInt64, Eos, set_fix, Drop, Star, Space, @with_pre, Parse, @p_str, Pattern, Alt, Delayed, parse_one, @E_str
import Dates

space = Drop(Star(Space()))

PNonZeroInt64() = Parse(p"-?[1-9][0-9]*", Int64)
PPosInt64() = Parse(p"[1-9][0-9]*", Int64)
PPosZeroInt64() = Parse(p"[0-9][0-9]*", Int64)
PCalendarNames() = p"[a-zA-Z0-9-\/\s|]+"

@with_pre space begin
    sum = Delayed()

    weekday_short = Alt(map(x -> Pattern(uppercase(x)), Dates.ENGLISH.days_of_week_abbr)...)
    month_short = Alt(map(x -> Pattern(uppercase(x)), Dates.ENGLISH.months_abbr)...)

    brackets = E"(" + space + sum + space + E")"
    repeat = E"Repeat(" + space + sum + space + E")" > rd -> Repeat(rd)
    rdate_term = Alt()

    next = E"Next(" + (space + sum) + (E"," + space + sum)[0:end] + space + E")" |> parts -> Next(parts)
    next_inc = E"Next!(" + (space + sum) + (E"," + space + sum)[0:end] + space + E")" |> parts -> Next(parts, true)
    previous = E"Previous(" + (space + sum) + (E"," + space + sum)[0:end] + space + E")" |> parts -> Previous(parts)
    previous_inc = E"Previous!(" + (space + sum) + (E"," + space + sum)[0:end] + space + E")" |> parts -> Previous(parts, true)

    rdate_expr = rdate_term | brackets | repeat | next | next_inc | previous | previous_inc

    neg = Delayed()
    neg.matcher = rdate_expr | (E"-" + neg > -)
    cal_adj = neg + (E"@" + PCalendarNames() + E"[" + Alt(map(Pattern, collect(keys(HOLIDAY_ROUNDING_MAPPINGS)))...) + E"]")[0:1] |> xs -> length(xs) == 1 ? xs[1] : CalendarAdj(map(String, split(xs[2], "|")), xs[1], HOLIDAY_ROUNDING_MAPPINGS[xs[3]])

    mult = Delayed()
    mult = cal_adj | ((PPosInt64() + space + E"*" + space + cal_adj) > (c,rd) -> multiply(rd, c)) | ((cal_adj + space + E"*" + space + PPosInt64()) > (rd,c) -> multiply(rd, c))

    add = E"+" + mult
    sub = E"-" + mult > -
    sum.matcher = (sub | mult) + (add | mult)[0:end] |> Base.sum

    entry = sum + Eos()
end

function register_grammar!(term)
    # Handle the spacing correctly
    push!(rdate_term.matchers[2].matchers, term)
end

macro rd_str(arg::String)
    val = parse_one(arg, entry)[1]
    isa(val, RDate) || error("Unable to parse $(arg) as RDate")
    return val
end

function rdate(arg::String)
    val = parse_one(arg, entry)[1]
    isa(val, RDate) || error("Unable to parse $(arg) as RDate")
    return val
end


text:
"""
    function vectorize_(abm::Agent,code::Expr)

Function to subtitute all the declared symbols of the agents in the expression into vector form.
optional arguments base and update add some intermediate names to the vectorized variables and to updated ones.
"""
function vectorize_(abm::Agent,code::Expr,p::Program_)

    #Vectorisation changes
    for (i,v) in enumerate(abm.declaredSymbols["Local"])

        bs = :localV

        code = postwalk(x->@capture(x, vAux_) && vAux == v ? :($bs[ic1_,$i]) : x, code)
        v2 = Meta.parse(string(v,"_i"))
        code = postwalk(x->@capture(x, vAux_) && vAux == v2 ? :($bs[ic1_,$i]) : x, code)
        v2 = Meta.parse(string(v,"_j"))
        code = postwalk(x->@capture(x, vAux_) && vAux == v2 ? :($bs[nnic2_,$i]) : x, code)

        if "Local" in keys(p.update)
            if v in keys(p.update["Local"])
                cc = findfirst(abm.declaredSymbols["Local"].==v)
                pos = p.update["Local"][v]
                #Code that may be optimized hopefully
                code = postwalk(x->@capture(x, localV[g_,c_] = v1__) && cc==c ? :(localVCopy[$g,$pos] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] += v1__) && cc==c ? :(localVCopy[$g,$pos] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] -= v1__) && cc==c ? :(localVCopy[$g,$pos] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] *= v1__) && cc==c ? :(localVCopy[$g,$pos] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] /= v1__) && cc==c ? :(localVCopy[$g,$pos] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] \= v1__) && cc==c ? :(localVCopy[$g,$pos] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] ÷= v1__) && cc==c ? :(localVCopy[$g,$pos] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] %= v1__) && cc==c ? :(localVCopy[$g,$pos] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] ^= v1__) && cc==c ? :(localVCopy[$g,$pos] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] &= v1__) && cc==c ? :(localVCopy[$g,$pos] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] |= v1__) && cc==c ? :(localVCopy[$g,$pos] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] ⊻= v1__) && cc==c ? :(localVCopy[$g,$pos] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] >>>= v1__) && cc==c ? :(localVCopy[$g,$pos] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] >>= v1__) && cc==c ? :(localVCopy[$g,$pos] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, localV[g_,c_] <<= v1__) && cc==c ? :(localVCopy[$g,$pos] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["Identity"])

        bs = :identityV

        code = postwalk(x->@capture(x, vAux_) && vAux == v ? :($bs[ic1_,$i]) : x, code)
        v2 = Meta.parse(string(v,"_i"))
        code = postwalk(x->@capture(x, vAux_) && vAux == v2 ? :($bs[ic1_,$i]) : x, code)
        v2 = Meta.parse(string(v,"_j"))
        code = postwalk(x->@capture(x, vAux_) && vAux == v2 ? :($bs[nnic2_,$i]) : x, code)

        if "Identity" in keys(p.update)
            if v in keys(p.update["Identity"])
                cc = findfirst(abm.declaredSymbols["Identity"].==v)
                pos = p.update["Identity"][v]
                #Code that may be optimized hopefully
                code = postwalk(x->@capture(x, identityV[g_,c_] = v1__) && cc==c ? :(identityVCopy[$g,$pos] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] += v1__) && cc==c ? :(identityVCopy[$g,$pos] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] -= v1__) && cc==c ? :(identityVCopy[$g,$pos] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] *= v1__) && cc==c ? :(identityVCopy[$g,$pos] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] /= v1__) && cc==c ? :(identityVCopy[$g,$pos] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] \= v1__) && cc==c ? :(identityVCopy[$g,$pos] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] ÷= v1__) && cc==c ? :(identityVCopy[$g,$pos] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] %= v1__) && cc==c ? :(identityVCopy[$g,$pos] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] ^= v1__) && cc==c ? :(identityVCopy[$g,$pos] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] &= v1__) && cc==c ? :(identityVCopy[$g,$pos] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] |= v1__) && cc==c ? :(identityVCopy[$g,$pos] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] ⊻= v1__) && cc==c ? :(identityVCopy[$g,$pos] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] >>>= v1__) && cc==c ? :(identityVCopy[$g,$pos] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] >>= v1__) && cc==c ? :(identityVCopy[$g,$pos] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, identityV[g_,c_] <<= v1__) && cc==c ? :(identityVCopy[$g,$pos] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["Global"])

        bs = :globalV

        code = postwalk(x->@capture(x, $v) ? :($bs[$i]) : x, code)

        if "Global" in keys(p.update)
            if v in keys(p.update["Global"])
                cc = findfirst(abm.declaredSymbols["Global"].==v)
                pos = p.update["Global"][v]
                #Code that may be optimized hopefully
                code = postwalk(x->@capture(x, globalV[c_] = v1__) && cc==c ? :(globalVCopy[$pos] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] += v1__) && cc==c ? :(globalVCopy[$pos] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] -= v1__) && cc==c ? :(globalVCopy[$pos] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] *= v1__) && cc==c ? :(globalVCopy[$pos] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] /= v1__) && cc==c ? :(globalVCopy[$pos] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] \= v1__) && cc==c ? :(globalVCopy[$pos] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ÷= v1__) && cc==c ? :(globalVCopy[$pos] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] %= v1__) && cc==c ? :(globalVCopy[$pos] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ^= v1__) && cc==c ? :(globalVCopy[$pos] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] &= v1__) && cc==c ? :(globalVCopy[$pos] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] |= v1__) && cc==c ? :(globalVCopy[$pos] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ⊻= v1__) && cc==c ? :(globalVCopy[$pos] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] >>>= v1__) && cc==c ? :(globalVCopy[$pos] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] >>= v1__) && cc==c ? :(globalVCopy[$pos] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] <<= v1__) && cc==c ? :(globalVCopy[$pos] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["GlobalArray"])

        if "GlobalArray" in keys(p.update)
            if v in keys(p.update["GlobalArray"])
                name = Meta.parse(string(v,GLOBALARRAYCOPY))
                code = postwalk(x->@capture(x, cc_[c__] = v1__) && cc==v ? :($name[$(c...)] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] += v1__) && cc==v ? :($name[$(c...)] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] -= v1__) && cc==v ? :($name[$(c...)] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] *= v1__) && cc==v ? :($name[$(c...)] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] /= v1__) && cc==v ? :($name[$(c...)] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] \= v1__) && cc==v ? :($name[$(c...)] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ÷= v1__) && cc==v ? :($name[$(c...)] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] %= v1__) && cc==v ? :($name[$(c...)] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ^= v1__) && cc==v ? :($name[$(c...)] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] &= v1__) && cc==v ? :($name[$(c...)] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] |= v1__) && cc==v ? :($name[$(c...)] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ⊻= v1__) && cc==v ? :($name[$(c...)] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] >>>= v1__) && cc==v ? :($name[$(c...)] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] >>= v1__) && cc==v ? :($name[$(c...)] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] <<= v1__) && cc==v ? :($name[$(c...)] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["Medium"])

        bs = :mediumV

        if abm.dims == 1
            code = postwalk(x->@capture(x, $v) ? :($bs[idMediumV[ic1_,1],$i]) : x, code)
        elseif abm.dims == 2
            code = postwalk(x->@capture(x, $v) ? :($bs[idMediumV[ic1_,1],idMediumV[ic1_,2],$i]) : x, code)
        elseif abm.dims == 3
            code = postwalk(x->@capture(x, $v) ? :($bs[idMediumV[ic1_,1],idMediumV[ic1_,2],idMediumV[ic1_,3],$i]) : x, code)
        end
    end

    return code
end

"""
    function vectorizeMedium_(abm::Agent,code::Expr)

Function to subtitute all the declared symbols of the medium in the expression into vector form.
"""
function vectorizeMedium_(abm::Agent,code::Union{Expr,Symbol,<:Number},p::Program_)

    for (i,v) in enumerate(abm.declaredSymbols["Global"])

        bs = :globalV

        code = postwalk(x->@capture(x, $v) ? :($bs[$i]) : x, code)

        if "Global" in keys(p.update)
            if v in keys(p.update["Global"])
                cc = findfirst(abm.declaredSymbols["Global"].==v)
                pos = p.update["Global"][v]
                #Code that may be optimized hopefully
                code = postwalk(x->@capture(x, globalV[c_] = v1__) && cc==c ? :(globalVCopy[$pos] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] += v1__) && cc==c ? :(globalVCopy[$pos] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] -= v1__) && cc==c ? :(globalVCopy[$pos] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] *= v1__) && cc==c ? :(globalVCopy[$pos] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] /= v1__) && cc==c ? :(globalVCopy[$pos] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] \= v1__) && cc==c ? :(globalVCopy[$pos] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ÷= v1__) && cc==c ? :(globalVCopy[$pos] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] %= v1__) && cc==c ? :(globalVCopy[$pos] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ^= v1__) && cc==c ? :(globalVCopy[$pos] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] &= v1__) && cc==c ? :(globalVCopy[$pos] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] |= v1__) && cc==c ? :(globalVCopy[$pos] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] ⊻= v1__) && cc==c ? :(globalVCopy[$pos] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] >>>= v1__) && cc==c ? :(globalVCopy[$pos] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] >>= v1__) && cc==c ? :(globalVCopy[$pos] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, globalV[c_] <<= v1__) && cc==c ? :(globalVCopy[$pos] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["GlobalArray"])

        if "GlobalArray" in keys(p.update)
            if v in keys(p.update["GlobalArray"])
                name = Meta.parse(string(v,GLOBALARRAYCOPY))
                code = postwalk(x->@capture(x, cc_[c__] = v1__) && cc==v ? :($name[$(c...)] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] += v1__) && cc==v ? :($name[$(c...)] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] -= v1__) && cc==v ? :($name[$(c...)] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] *= v1__) && cc==v ? :($name[$(c...)] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] /= v1__) && cc==v ? :($name[$(c...)] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] \= v1__) && cc==v ? :($name[$(c...)] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ÷= v1__) && cc==v ? :($name[$(c...)] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] %= v1__) && cc==v ? :($name[$(c...)] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ^= v1__) && cc==v ? :($name[$(c...)] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] &= v1__) && cc==v ? :($name[$(c...)] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] |= v1__) && cc==v ? :($name[$(c...)] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] ⊻= v1__) && cc==v ? :($name[$(c...)] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] >>>= v1__) && cc==v ? :($name[$(c...)] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] >>= v1__) && cc==v ? :($name[$(c...)] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, cc_[c__] <<= v1__) && cc==v ? :($name[$(c...)] <<= $(v1...)) : x, code)
            end
        end
    end

    for (i,v) in enumerate(abm.declaredSymbols["Medium"])

        bs = :mediumV

        if abm.dims == 1
            code = postwalk(x->@capture(x, $v) ? :($bs[ic1_,$i]) : x, code)
        elseif abm.dims == 2
            code = postwalk(x->@capture(x, $v) ? :($bs[ic1_,ic2_,$i]) : x, code)
        elseif abm.dims == 3
            code = postwalk(x->@capture(x, $v) ? :($bs[ic1_,ic2_,ic3_,$i]) : x, code)
        end

        if "Medium" in keys(p.update)
            if v in keys(p.update["Medium"])
                cc = findfirst(abm.declaredSymbols["Medium"].==v)
                pos = p.update["Medium"][v]
                #Code that may be optimized hopefully
                code = postwalk(x->@capture(x, mediumV[g__,c_] = v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] = $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] += v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] += $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] -= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] -= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] *= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] *= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] /= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] /= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] \= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] \= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] ÷= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] ÷= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] %= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] %= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] ^= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] ^= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] &= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] &= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] |= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] |= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] ⊻= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] ⊻= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] >>>= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] >>>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] >>= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] >>= $(v1...)) : x, code)
                code = postwalk(x->@capture(x, mediumV[g__,c_] <<= v1__) && cc==c ? :(mediumVCopy[$(g...),$pos] <<= $(v1...)) : x, code)
            end
        end
    end

    return code

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia