text:
#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")
#Pkg.add("PyCall")

using PyCall
@pyimport nltk.stem as nltk_stem
@pyimport nltk.stem.isri as nltk_stem_isri
@pyimport stopwords as stopwordsmod

PUNCTUACTION = ";:,.@#&\\-\"'/:*"
SYMBOLS = "()[]¿?¡!{}~<>|"
SKIP_SYMBOLS = string(PUNCTUACTION, SYMBOLS)
# SKIP_WORDS = set(["…", "..", "...", "...."])

type TextConfig
    lc::Bool
    del_diac::Bool
    del_dup::Bool

    del_punc::Bool
    del_users::Bool
    del_urls::Bool

    del_num::Bool
    del_stopwords::Bool
    del_nonascii::Bool
    stem::Bool

    min_length::Int
    max_length::Int
    ngrams::Int
    lang::String

    stopwords::Set{String}
    stemmer::Any
    TextConfig() = new(true, true, true, true, true, true, true, true, false, true, 2, 16, 1, "english", Set{String}(), -1)
end

function load(config::TextConfig, stopwords)
    length(config.lang) == 0 && return

    if config.del_stopwords
        for w in stopwordsmod.get_stopwords(config.lang)
            push!(config.stopwords, normchars(w, config) |> join)
        end
        if length(stopwords) > 0
            for w in readlines(stopwords)
                w = normchars(strip(w), config) |> join
                push!(config.stopwords, w)
            end
        end
    end

    if config.lang in ("en", "english")
        config.stemmer = nltk_stem.PorterStemmer()
    elseif config.lang in ("ar", "arabic")
        config.stemmer = nltk_stem_isri.ISRIStemmer()
    else
        config.stemmer = nltk_stem.SnowballStemmer(config.lang)
    end
end

function normchars(text, config::TextConfig)
    if config.del_users
        text = replace(text, r"@\w+\s?", "")
    end

    if config.del_urls
        text = replace(text, r"https?:\S+\s?", "")
        text = replace(text, r"ftp:\S+\s?", "")
    end

    if config.lc
        text = lowercase(text)
    end

    L = Vector{Char}()
    blank = ' '
    prev = blank

    for u in normalize_string(text, :NFD)
        if config.del_diac
            o = Int(u)
            0x300 <= o && o <= 0x036F && continue
        end

        if u in ('\n', '\r', ' ', '\t')
            u = blank
        elseif config.del_dup && prev == u
            continue
        elseif config.del_punc && u in SKIP_SYMBOLS
            prev = u
            continue
        elseif config.del_num && isdigit(u)
            continue
        end
        prev = u
        push!(L, u)
    end

    return L
end

function tokenize(text, config::TextConfig)
    if contains(text, "(@")
        text = split(text, "(@")[1]
    end
    text = join(normchars(text, config))
    tokens = String[]

    for word in split(text)
        len = length(word)
        if len <= config.min_length || len >= config.max_length
            continue
        end

        config.del_stopwords && word in config.stopwords && continue

        if config.stem
            word = try
                config.stemmer[:stem](word)
            catch
                word
            end

        end
        config.del_nonascii && !isascii(word) && continue

        push!(tokens, word)
    end
    if config.ngrams == 1
        return tokens
    end

    return [join(tokens[i:i+config.ngrams-1], ' ') for i in 1:(length(tokens)-config.ngrams+1)]
end

function normtext{T <: AbstractString}(text::T, config::TextConfig)
    join(tokenize(text, config), ' ')
end

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
module Fatou
using SyntaxTree,Reduce,LaTeXStrings,Requires,Base.Threads

#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

export fatou, juliafill, mandelbrot, newton, basin, plot

abstract type ComplexBundle end

struct Rectangle <: ComplexBundle
    ∂::Vector{Float64} # bounds
    n::UInt16 # number of grid points
    function Rectangle(∂::T=π/2,n::Integer=176) where T
        !(T <: Array) && (∂ = [-float(∂),∂,-∂,∂])
        length(∂) == 2 && (∂ = [∂[1],∂[2],∂[1],∂[2]])
        new(convert(Vector{Float64},∂),UInt16(n))
    end
end

struct ComplexRectangle <: ComplexBundle
    ∂::Rectangle
    Ω::Matrix{Complex{Float64}}
end

bounds(R::ComplexRectangle) = bounds(Rectangle(R))
bounds(Ω::Rectangle) = Ω.∂
Base.size(Ω::Rectangle) = (round(UInt16,(Ω.∂[4]-Ω.∂[3])/(Ω.∂[2]-Ω.∂[1])*Ω.n),Ω.n)
Base.size(R::ComplexRectangle) = size(R.Ω)

"""
    Fatou.Define(E::Any;                  # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      newt::Bool  = false,                # toggle Newton mode
      m::Number   = 0,                    # Newton multiplicity factor
      O::String   = F,                    # original Newton map
      mandel::Bool= false,                # toggle Mandelbrot mode
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` the metadata for a `Fatou.FilledSet`.
"""
struct Define{FT<:Function,QT<:Function,CT<:Function,M,N,P,D} <: ComplexBundle
    E::Any # input expression
    F::FT # primary map
    Q::QT # escape criterion
    C::CT # complex fixed point coloring
    Ω::Rectangle # bounded grid points
    N::UInt16 # number of iterations
    ϵ::Float64 # epsilon Limit criterion
    iter::Bool # toggle iteration mode
    p::Float64 # iteration color exponent
    newt::Bool # toggle Newton mode
    m::Number # newton multiplicity factor
    mandel::Bool # toggle Mandelbrot mode
    seed::Number # Mandelbrot seed value
    x0 # orbit starting point
    orbit::Int # orbit cobweb depth
    depth::Int # depth of function composition
    cmap::String # imshow color map
    plane::Bool # convert input disk to half-plane
    disk::Bool # convert output half-plane to disk
    function Define(E;
            Q=:(abs2(z)),
            C=:((angle(z)/(2π))*n^p),
            ∂=π/2,
            n::Integer=176,
            N::Integer=35,
            ϵ::Number=4,
            iter::Bool=false,
            p::Number=0,
            newt::Bool=false,
            m::Number=0,
            mandel::Bool=false,
            seed::Number=0.0+0.0im,
            x0=nothing,
            orbit::Int=0,
            depth::Int=1,
            cmap::String="",
            plane::Bool=false,
            disk::Bool=false)
        !newt ? (f = genfun(E,[:z,:c]); q = genfun(Q,[:z,:c])) :
        (f = genfun(newton_raphson(E,m),[:z,:c]); q = genfun(Expr(:call,:abs,E),[:z,:c]))
        c = genfun(C,[:z,:n,:p])
        e = typeof(E) == String ? parse(E) : E
        return new{typeof(f),typeof(q),typeof(c),mandel,newt,plane,disk}(e,f,q,c,Rectangle(∂,n),UInt16(N),float(ϵ),iter,float(p),newt,m,mandel,seed,x0,orbit,depth,cmap,plane,disk)
    end
end

"""
    Fatou.FilledSet(::Fatou.Define)

Compute the `Fatou.FilledSet` set using `Fatou.Define`.
"""
struct FilledSet{FT,QT,CT,M,N,P,D} <: ComplexBundle
    meta::Define{FT,QT,CT,M,N,P,D}
    set::ComplexRectangle
    iter::Matrix{UInt16}
    mix::Matrix{Float64}
    function FilledSet{FT,QT,CT,M,N,P,D}(K::Define{FT,QT,CT,M,N,P,D},Z::ComplexRectangle) where {FT,QT,CT,M,N,P,D}
        (i,s) = Compute(K,Z)
        return new{FT,QT,CT,M,N,P,D}(K,s,i,broadcast(K.C,s.Ω,broadcast(float,i./K.N),K.p))
    end
end

Rectangle(K::Define) = K.Ω
Rectangle(K::FilledSet) = Rectangle(ComplexRectangle(K))
Rectangle(R::ComplexRectangle) = R.∂

ComplexRectangle(K::FilledSet) = K.set
ComplexRectangle(Ω::Matrix{Complex{Float64}}) = ComplexRectangle(Rectangle([0,size(Ω)[2],0,size(Ω)[1]],size(Ω)[2]),Ω)

bounds(K::Define) = bounds(Rectangle(K))
bounds(K::FilledSet) = bounds(ComplexRectangle(K))

ranges(K::FilledSet) = ranges(K.meta)
ranges(K::Define) = ranges(Rectangle(K))
function ranges(Ω::Rectangle)
    yn,xn = size(Ω)
    x = range(Ω.∂[1]+0.0001,stop=Ω.∂[2],length=xn)
    y = range(Ω.∂[4],stop=Ω.∂[3],length=yn)
    return x,y
end

(K::Define)(Z) = fatou(K,Z)
(K::FilledSet)(Z) = fatou(K,Z)

"""
      fatou(::Fatou.Define)

Compute the `Fatou.FilledSet` set using `Fatou.Define`.

# Examples
```Julia
julia> fatou(K)
```
"""
fatou(K::Define{FT,QT,CT,M,N,P,D},Z::ComplexRectangle) where {FT,QT,CT,M,N,P,D} = FilledSet{FT,QT,CT,M,N,P,D}(K,Z)
fatou(K::Define,Z::Rectangle=Rectangle(K)) = fatou(K,fatou(Z))
fatou(K::Define,Z::FilledSet) = fatou(K,ComplexRectangle(Z))
fatou(K::Define,Z::Define) = fatou(K,fatou(Z))
fatou(K::FilledSet,Z=ComplexRectangle(K)) = fatou(K.meta,Z)
function fatou(K::Rectangle) # generate coordinate grid
    x, y = ranges(K)
    ComplexRectangle(K, x' .+ im*y)
end # apply Newton-Orbit function element-wise to coordinate grid

"""
    juliafill(::Expr;                     # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` filled Julia basin in `Fatou`

# Exmaples
```Julia
julia> juliafill("z^2-0.06+0.67im",∂=[-1.5,1.5,-1,1],N=80,n=1501,cmap="RdGy")
```
"""
function juliafill(E;
        Q=:(abs2(z)),
        C=:((angle(z)/(2π))*n^p),
        ∂=π/2,
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=4,
        iter::Bool=false,
        p::Number=0,
        newt::Bool=false,
        m::Number=0,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    return Define(E,Q=Q,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=newt,m=m,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

"""
    mandelbrot(::Expr;                    # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :(exp(-abs(z))*n^p),  # coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      m::Number   = 0,                    # Newton multiplicity factor
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` Mandelbrot basin in `Fatou`

# Examples
```Julia
mandelbrot(:(z^2+c),n=800,N=20,∂=[-1.91,0.51,-1.21,1.21],cmap="nipy_spectral")
```
"""
function mandelbrot(E;
        Q=:(abs2(z)),
        ∂=π/2,
        C=:(exp(-abs(z))*n^p),
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=4,
        iter::Bool=false,
        p::Number=0,
        newt::Bool=false,
        m::Number=0,
        seed::Number=0.0+0.0im,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    m ≠ 0 && (newt = true)
    return Define(E,Q=Q,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=newt,m=m,mandel=true,seed=seed,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

"""
    newton(::Expr;                        # primary map, (z, c) -> F
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      m::Number   = 0,                    # Newton multiplicity factor
      mandel::Bool= false,                # toggle Mandelbrot mode
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` Newton basin in `Fatou`

# Examples
```Julia
julia> newton("z^3-1",n=800,cmap="brg")
```
"""
function newton(E;
        C=:((angle(z)/(2π))*n^p),
        ∂=π/2,
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=0.01,
        iter::Bool=false,
        p::Number=0,
        m::Number=1,
        mandel::Bool=false,
        seed::Number=0.0+0.0im,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    return Define(E,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=true,m=m,mandel=mandel,seed=seed,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

# load additional functionality
include("internals.jl"); include("orbitplot.jl")

"""
    basin(::Fatou.Define, ::Integer)

Output the `j`-th basin of `Fatou.Define` as LaTeX.
Each subsequent iteration of the Newton-Raphson method will yield a more complicated set.

# Examples
```Julia
julia> basin(newton("z^3-1"),2)
L"\$\\displaystyle D_2(\\epsilon) = \\left\\{z\\in\\mathbb{C}:\\left|\\,z - \\frac{\\left(z - \\frac{z^{3} - 1}{3 z^{2}}\\right)^{3} - 1}{3 \\left(z - \\frac{z^{3} - 1}{3 z^{2}}\\right)^{2}} - \\frac{z^{3} - 1}{3 z^{2}} - r_i\\,\\right|<\\epsilon,\\,\\forall r_i(\\,f(r_i)=0 )\\right\\}\$"
```
"""
basin(K::Define,j) = K.newt ? nrset(K.E,K.m,j) : jset(K.E,j)

plane(z::Complex) = (2z.re/(z.re^2+(1-z.im)^2))+im*(1-z.re^2-z.im^2)/(z.re^2+(1-z.im)^2)
disk(z::Complex) = (2z.re/(z.re^2+(1+z.im)^2))+im*(z.re^2+z.im^2-1)/(z.re^2+(1+z.im)^2)

# define function for computing orbit of a z0 input
function orbit(K::Define{FT,QT,CT,M,N,P,D},z0::Complex{Float64}) where {FT,QT,CT,M,N,P,D}
    M ? (z = K.seed) : (z = P ? plane(z0) : z0)
    zn = 0x0000
    while (N ? (K.Q(z,z0)::Float64>K.ϵ)::Bool : (K.Q(z,z0)::Float64<K.ϵ))::Bool && K.N>zn
        z = K.F(z,z0)::Complex{Float64}
        zn+=0x0001
    end; #end
    # return the normalized argument of z or iteration count
    return (zn::UInt16,(D ? disk(z) : z)::Complex{Float64})
end

"""
    Compute(::Fatou.Define)::Union{Matrix{UInt8},Matrix{Float64}}

`Compute` the `Array` for `Fatou.FilledSet` as specefied by `Fatou.Define`.
"""
function Compute(K::Define{FT,QT,CT,M,N,D},Z::ComplexRectangle)::Tuple{Matrix{UInt16},ComplexRectangle} where {FT,QT,CT,M,N,D}
    yn,xn = size(Z)
    (matU,matF) = (Matrix{UInt16}(undef,yn,xn),Matrix{Complex{Float64}}(undef,yn,xn))
    @time @threads for j = 1:yn; for k = 1:xn;
        (matU[j,k],matF[j,k]) = orbit(K,Z.Ω[j,k])::Tuple{UInt16,Complex{Float64}}
    end; end
    return (matU,ComplexRectangle(Rectangle(Z),matF))
end

# determine if plot is Iteration, Roots, or Limit
typeplot(K::FilledSet) = typeof(K.meta.iter ? K.iter : K.mix) == Matrix{UInt16} ? "iter." : K.meta.m==1 ? "roots" : "limit"

function Base.String(K::FilledSet) # annotate title using LaTeX
    text,t = "f : z ↦ $(K.meta.E),",typeplot(K)
    K.meta.newt ? "$text m = $(K.meta.m), $t" : "$text $t"
end

function __init__()
    println("Fatou detected $(Threads.nthreads()) julia threads.")
    @require ColorSchemes="35d6a980-a343-548e-a6ea-1d62b119f2f4" begin
        nonan(x) = isnan(x) ? 0.0 : x
        function (C::ColorSchemes.ColorScheme)(K::FilledSet)
            S = size(K.iter)
            H = zeros(ColorSchemes.RGB{Float64},S...)
            if K.meta.iter
                M = length(C)/(maximum(K.iter)+1)
                for x ∈ 1:S[1], y ∈ 1:S[2]
                    H[x,y] = C[round(Int,M*(K.iter[x,y]+1),RoundUp)]
                end
            else
                for x ∈ 1:S[1], y ∈ 1:S[2]
                    H[x,y] = get(C,nonan(K.mix[x,y]))
                end
            end
            return H
        end
    end
    @require ImageInTerminal="d8c32880-2388-543b-8c61-d9f865259254" begin
        import ColorSchemes
        function Base.show(io::IO,K::FilledSet;c::String="",bare::Bool=false)
            isempty(c) && (c = K.meta.cmap)
            display(getproperty(ColorSchemes, isempty(c) ? :balance : Symbol(c))(K))
            !bare && print(io,String(K))
        end
    end
    @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" include("makie.jl")
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("pyplot.jl")
    @require UnicodePlots="b8865327-cd53-5732-bb35-84acbb429228" include("uniplots.jl")
end

end # module

ext: jl


text:
using LowRankModels
import StatsBase: sample, Weights

# test ordistic loss

## generate data
srand(1);
m,n,k = 100,100,3;
kfit = k+1
nlevels = 5; # number of levels
d = nlevels-1 # embedding dimension
D = n*d;
# coordinates of covariates
X_real = randn(m,k)
# directions of observations
Y_real = randn(k,n)
# measurement thresholds
T_real = k*randn(d,n) # notice x^T y has variance k; so making this bigger makes the problem easier
for j=1:n # this scheme doesn't work to ensure uniform sampling
	T_real[:,j] = sort(T_real[:,j])
end

signedsums = @compat Array{Float64}(d, nlevels)
for i=1:d
    for j=1:nlevels
        signedsums[i,j] = i<j ? 1 : -1
    end
end

XY = X_real*Y_real;
XYplusT = zeros(Float64, (m,D))
A = zeros(Int, (m, n))
for i=1:m
	for j=1:n
		u = XY[i,j] + T_real[:,j]
		XYplusT[i,(j-1)*d+(1:d)] = u
		diffs = u'*signedsums
		wv = Weights(Float64[exp(-diffs[l]) for l in 1:nlevels])
		l = sample(wv)
		A[i,j] = l
	end
end
# loss is insensitive to shifts; regularizer should pick this shift
XYplusT = XYplusT .- mean(XYplusT, 2)

# and the model
losses = BvSLoss(nlevels)
rx, ry = lastentry1(QuadReg(.01)), OrdinalReg(QuadReg(.01)) #lastentry_unpenalized(QuadReg(10));
glrm = GLRM(A,losses,rx,ry,kfit, scale=false, offset=false, X=randn(kfit,m), Y=randn(kfit,D));

# fit w/o initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y#[:,1:d:D];
@show ch.objective
println("After fitting, parameters differ from true parameters by $(vecnorm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs,(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)\% of entries would be wrong.)\n")

# initialize
init_svd!(glrm)
XYh = glrm.X' * glrm.Y
println("After initialization with the svd, parameters differ from true parameters by $(vecnorm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After initialization with the svd, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After initialization with the svd, imputed entries are off by $(sum(abs,(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)\% of entries would be wrong.)\n")

# fit w/ initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y#[:,1:d:D];
@show ch.objective
println("After fitting, parameters differ from true parameters by $(vecnorm(XYplusT - XYh)/sqrt(prod(size(XYplusT)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs,(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, $((nlevels-1)/nlevels*100)\% of entries would be wrong.)\n")

# test scaling
scale!(glrm.ry, 3)


text:
# update_weights.jl

export lω_to_v, v_to_lω;

function lω_to_v(lω::Array{T, 1}) where T <: Real
    H = length(lω)
    H > 2 || throw(ArgumentError("H must be at least 3."))

    v = Array{T, 1}(undef, H-1)

    v[1] = exp(lω[1])
    cumsuml1minusv = log(1.0 - v[1])

    for h = 2:(H-1)
        v[h] = exp(lω[h] - cumsuml1minusv)
        cumsuml1minusv += log(1.0 - v[h])
    end
    cumsuml1minusv ≈ lω[H] || throw(error("lω_to_v unstable."))

    return v
end

function v_to_lω(v::Array{T, 1}) where T <: Real
    H = length(v) + 1
    H > 2 || throw(ArgumentError("H must be at least 3."))

    lω = Array{T, 1}(undef, H)

    lω[1] = log(v[1])
    cumsuml1minusv = log(1.0 - v[1])

    for h = 2:(H-1)
        lω[h] = log(v[h]) + cumsuml1minusv
        cumsuml1minusv += log(1.0 - v[h])
    end

    lω[H] = cumsuml1minusv

    return lω
end

function lfc_v(v::Array{T, 1}, a_v::Array{T, 1}, b_v::Array{T, 1},
    lNX_mat::Array{T, 2}) where T <: Real
    ### log full (block) conditional for v

    lω = v_to_lω(v)

    lωNX_vec = lωNXvec(lω, lNX_mat)

    ## replaced by lωXvec
    # lωNX_mat = broadcast(+, permutedims(lω), lNX_mat) # n by H
    # lωNX_vec = BayesInference.logsumexp(lωNX_mat, 2) # n vector

    ldb = [ logpdf( Beta(a_v[h], b_v[h]), v[h] ) for h = 1:length(v) ]

    return sum(ldb) - sum(lωNX_vec), lωNX_vec
end

function mvSlice_v(v::Array{T, 1}, a_v::Array{T, 1}, b_v::Array{T, 1}, lNX_mat::Array{T, 2},
    w_slice::T=1.0, maxtries::Int=1000) where T <: Real

    ### Multivariate slice sampler Figure 8 from Neal (2003)

    HH = length(v)

    ### Step A: Define the slice

    lf0, trash = lfc_v(v, a_v, b_v, lNX_mat)
    ee = rand(Exponential(1.0))
    zz = lf0 - ee

    ### Step B: Randomly position hyperrectangle

    U = rand(HH)
    L = max.(v .- w_slice .* U, 0.0 + eps(T))
    R = min.(L .+ w_slice, 1.0 - eps(T))

    ### Step C: Sample from hyperrectangle, shrinking when rejected

    keeptrying = true
    tries = 1
    cand = Array{T, 1}(undef, HH) # define outside scope of while loop
    lωNX_vec_out = Array{T, 1}(undef, length(trash))

    while keeptrying

        cand = L .+ rand(HH) .* (R .- L)
        lfcand, lωNX_vec_out = lfc_v(cand, a_v, b_v, lNX_mat)

        keeptrying = zz > lfcand

        if keeptrying

            tries += 1
            tries <= maxtries || throw(error("Exceeded maximum slice attempts."))

            for h = 1:HH
                if cand[h] < v[h]
                    L[h] = deepcopy(cand[h])
                else
                    R[h] = deepcopy(cand[h])
                end
                ## equivalent to
                # cand[h] < v[h] ? L[h] = copy(cand[h]) : R[h] = copy(cand[h])
            end
        end

    end

    return cand, lωNX_vec_out
end

function update_vlω_mvSlice!(model::Model_BNP_WMReg_Joint) where T <: Real

    M = StatsBase.counts(model.state.S, 1:model.H)
    a_v = 1.0 .+ M[1:(model.H-1)]
    b_v = reverse( model.state.α .+ cumsum( reverse( M[2:model.H] ) ) )

    if ( model.state.γδc == Inf || isnothing(model.state.γδc) ) && sum(model.state.γ) == 0
        v_new = [ rand( Beta(a_v[h], b_v[h]) ) for h = 1:(model.H-1) ]
    else
        v_new, lωNX_vec_new = mvSlice_v(model.state.v, a_v, b_v, model.state.lNX)
        model.state.lωNX_vec = lωNX_vec_new
    end

    lω_new = v_to_lω(v_new)

    model.state.v = v_new
    model.state.lω = lω_new

    return nothing
end


text:
using FastSpike
using LinearAlgebra: I
using Plots
using ProgressMeter

const time = 100
const N = 20
const input_dim = 5
# Define Network
net = Network(LIF(1), 2, STDP(1e-2, 1e-4, 20, 10; softbound=true, min_weight=0, max_weight=1))
input = add_group!(net, input_dim)
neurons = add_group!(net, N)
# input to neurons
adjacency = rand(0:1, input_dim, N)
weights = rand(input_dim, N) .* adjacency
connect!(net, input, neurons, weights, adjacency)
# connections in neural group
adjacency = ones(N, N) - I
weights = rand(N, N) .* adjacency
connect!(net, neurons, neurons, weights, adjacency)
net.weight
# Define Monitor to record network activities
monitor = Monitor()
#plot weights histogram
plot_weights(activity, active_neurons) = display(
    plot(
        histogram(collect(Iterators.flatten(net.weight[input.idx, neurons.idx])), bins=25),
        histogram(collect(Iterators.flatten(net.weight[neurons.idx, neurons.idx])), bins=25),
        plot(activity),
        plot(scatter(active_neurons, markersize=2, c=:black)),
        layout=(4, 1)
    )
)
# Define training function
function train(time)
    # Generate input spikes
    input_spikes = convert(Matrix{Bool}, [rand(0:1, time, input_dim) zeros(time, N)])
    active_neurons = []
    @showprogress 1 "training " for t = 1:time
        run!(net, input_spikes=reshape(input_spikes[t, :], 1, :))
        record!(monitor, net)
        append!(active_neurons, [(t, i[2]) for i in findall(net.spikes)])
        # plot_weights(sum(monitor.spikes[:, 1, neurons.idx], dims = 2), Tuple.(active_neurons))
    end
end
# Run
plot_weights([], [])
train(time)


text:
@testset "2D p-FEM" begin
    W = JacobiWeight(1,1) .* Jacobi(1,1)
    x = axes(W,1)
    D = Derivative(x)

    D2 = -((D*W)'*(D*W))
    M = W'W
    A = KronTrav(D2,M)
    N = 30;
    V = view(A,Block(N,N));
    @time MemoryLayout(arguments(V)[2]) isa LazyBandedMatrices.ApplyBandedLayout{typeof(*)}

    Δ = KronTrav(D2,M) + KronTrav(M,D2-M)
    N = 100; @time L = Δ[Block.(1:N+2),Block.(1:N)];
    r = KronTrav(M,M)[Block.(1:N+2),1]
    @time F = qr(L);
    @time u = F \ r;

    # u = Δ \ [1; zeros(∞)];
end


text:
module NoPool

# dummy allocator that passes through any requests, calling into the GC if that fails.

using ..CUDA
import ..CUDA: actual_alloc, actual_free

using ..PoolUtils

init() = return

function alloc(dev, sz)
    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 alloc" begin
            block = actual_alloc(dev, sz, phase==3)
        end
        block === nothing || break
    end

    return block
end

function free(dev, block)
    actual_free(dev, block)
    return
end

reclaim(dev, target_bytes::Int=typemax(Int)) = return 0

cached_memory(dev=device()) = 0

end


text:
using Base.Test
using Attributes

g = 5;
threshold = 0.55
v = 3
b = OrderedAttributes(g, threshold, v)
nodes = 3;
attr = get_attributes(b, nodes)

@test min(attr...) == 1 #test rozmiaru
@test max(attr...) == v #test wartosci powinny byc +-1

attr = ones(Int, nodes, g)
attr[1, :] = 3
attr[2, 1:2] = 3
attr[2, 3:5] = 2
#attr is:
# [3  3  3  3  3
#  3  3  2  2  2
#  1  1  1  1  1]
#This should give following weights:
# [0.0  0.5  -0.9
#  0.0  0.0  -0.3
#  0.0  0.0   0.0]
exp_weights = zeros(3, 3)
exp_weights[1, 2] = 0.5;
exp_weights[1, 3] = -0.9;
exp_weights[2, 3] = -0.3;

weights = get_attribute_layer_weights(b, attr)

@test size(weights) == (nodes, nodes)
@test weights ≈ exp_weights

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Examples

using ..OpticSim
using ..OpticSim.Geometry
using ..OpticSim.Repeat
using ..OpticSim.Emitters
using ..OpticSim.ParaxialAnalysis
using ..OpticSim.ParaxialAnalysis.HeadEye

using StaticArrays
using Colors
import Makie

emitterarray(emitters, transform = identitytransform()) = Emitters.Sources.CompositeSource(transform, emitters)

function example1()

    t = identitytransform()
    t = Transform(OpticSim.Vec4(t[:,1]), OpticSim.Vec4(t[:,2]), OpticSim.Vec4(t[:,3])*-1.0, OpticSim.Vec4(t[:,4]))
    h = HeadEye.Head(t)

    # HeadEye.lookAt!(HeadEye.eye(h, :left), SVector(-0.5, 0.0, 1.0))

    # -----------------------------------------------------------
    # mla stuff
    # -----------------------------------------------------------

    # compute center of mla - in front of left eye - oposite the tirection of the head
    left_eye_tr = Transform(origin(HeadEye.tr(h, :left_eye)), forward(HeadEye.tr(h, :head)))
    pos = local2world(left_eye_tr) * SVector(0.0, 0.0, 30.0)
    mla_rt = Transform(pos, forward(HeadEye.tr(h, :head)) * -1.0)
    # eli = Ellipse(7.0, 7.0, unitZ3(), pos)
    # Vis.draw!(eli, transparency=true, color=RGBA(1.0, 0.2, 0.2, 0.4))

    csg = HeadEye.csg_sphere(radius=20.0)
    # csg = HeadEye.csg_cylinder(radius = 20.0, added_rotation = rotationX(π/2.0))
    # csg = HeadEye.csg_plane()

    # shapes_2d = HeadEye.get_shapes(HexBasis1(), resolution=(8,5), size=1.0)

     shapes_2d = HeadEye.get_shapes(RectangularBasis(), resolution=(5,5), size=1.0)

    shapes_3d = HeadEye.project(shapes_2d, csg)

    # create the planar polygons and the lenses
    paraxial_lenses = Vector{ParaxialLens}(undef, 0)
    emitters = Vector{Emitters.Sources.Source}(undef, 0)
    for shape in shapes_3d

        # sanity test to change the lens optical center according to the shape's coordinates (relative to the center)
        # this one is for rectangles - need better algorithm for hexagons
        coord = HeadEye.coordinates(shape)
        center_point = SVector(-3.0*coord[1], 0.0)
        center_point = SVector(-3.0*coord[1], -3.0*coord[2])
        # center_point = SVector(0.0, -3.0*coord[2])
        # center_point = SVector(0.0, 0.0)

        lens = HeadEye.build_paraxial_lens(
            shape,
            local_center_point = center_point,
            parent_transform = mla_rt,
            scale=0.99,
            focal_length=30.0
        )
        push!(paraxial_lenses, lens)

        emitter = HeadEye.build_emitter(lens; distance=1.0, size=(1.0, 1.0))
        push!(emitters, emitter)
    end

    mla = LensAssembly(paraxial_lenses...)
    display = emitterarray(emitters)

    pupil = HeadEye.pupil(HeadEye.eye(h, :left))
    pupil_tr = HeadEye.tr(h, :left_pupil)
    detector = Circle(HeadEye.pupilsize(pupil) / 2.0, forward(pupil_tr), origin(pupil_tr), interface = opaqueinterface())
    sys = CSGOpticalSystem(mla, detector)

    # -----------------------------------------------------------
    # rendering
    # -----------------------------------------------------------
    Vis.set_current_mode(:normal)

    if (true)
        Vis.drawtracerays(sys; raygenerator=display, trackallrays = true, colorbynhits = true, test = true, numdivisions = 100, drawgen = false)
        # Vis.draw!(sys)

        Vis.draw!(h, draw_head=false)

        Vis.draw!(display, debug=false)
    end
end
export example1

end # module Test


text:
using Cardio
using Test
using Random

Random.seed!(1776)

#testfiles = ["test_adaptiveHRVFilter.jl", "test_detectRPeaks.jl", "test_util.jl"]
(root, dirs, files) = first(walkdir(dirname(@__FILE__)))

for testfile in files[occursin.("test_", files)]
    eval(:(@testset $testfile begin include($testfile) end))
end


text:
using Dates
using TimeZones

const JULIAN_EPOCH_USDA = 2415078.5 # 1990-03-01
const JULIAN_EPOCH_UNIX = 2440587.5 # 1970-01-01

datetime_from_julian_day_WEA(year, jday, time::Time, tz::TimeZone, occurrence) =
    zoned_datetime(Date(year) + (Day(jday) - Day(1)) + time, tz, occurrence)
datetime_from_julian_day_WEA(year, jday, tz::TimeZone) = datetime_from_julian_day_WEA(year, jday, "00:00", tz)
#HACK: handle different API for Fixed/VariableTimeZone
zoned_datetime(dt::DateTime, tz::TimeZone, occurrence=1) = ZonedDateTime(dt, tz)
zoned_datetime(dt::DateTime, tz::VariableTimeZone, occurrence=1) = ZonedDateTime(dt, tz, occurrence)

#julian_day_from_datetime(clock::Dates.AbstractDateTime) = dayofyear(clock)

#round_datetime(clock::Dates.AbstractDateTime) = round(clock, Minute)

datetime_from_julian_day_2DSOIL(jday, jhour=0) = begin
    d = (jday + jhour) + (JULIAN_EPOCH_USDA - JULIAN_EPOCH_UNIX)
    #HACK prevent degenerate timestamps due to precision loss
    t = ZonedDateTime(1970, 1, 1, tz"UTC") + Day(d)
    round(t, Minute)
end

julian_day_from_datetime_2DSOIL(clock::ZonedDateTime; hourly=false) = begin
    s = clock - ZonedDateTime(1970, 1, 1, tz"UTC") |> Second
    j = s - (JULIAN_EPOCH_USDA - JULIAN_EPOCH_UNIX)
    #FIXME: type instability
    hourly ? j : Int(round(j))
end

julian_hour_from_datetime_2DSOIL(clock::ZonedDateTime) =
    julian_day_from_datetime_2DSOIL(clock; hourly=true) - julian_day_from_datetime_2DSOIL(clock; hourly=false)

date(year, month::Int, day::Int; tz) = ZonedDateTime(year, month, day, tz)
date(year, doy::Int; tz) = ZonedDateTime(year, 1, 1, tz) + Dates.Day(doy - 1)
date(year, d::ZonedDateTime; tz) = d
date(year, d::DateTime; tz) = ZonedDateTime(d, tz)
date(year, ::Nothing; tz) = nothing

storagedays(t::ZonedDateTime) = (t - ZonedDateTime(year(t), 6, 30, timezone(t))) |> Day |> Dates.value

using CSV
using DataFrames: DataFrames, DataFrame

loadwea(filename, timezone; indexkey=:index) = begin
    df = CSV.File(filename) |> DataFrame
    df[!, indexkey] = map(r -> begin
        occurrence = 1
        i = DataFrames.row(r)
        if i > 1
            r0 = parent(r)[i-1, :]
            r0.time == r.time && (occurrence = 2)
        end
        datetime_from_julian_day_WEA(r.year, r.jday, r.time, timezone, occurrence)
    end, eachrow(df))
    df
end

datapath() = joinpath(@__DIR__,  "../../data/")
datapath(filename) = joinpath(datapath(), filename)


text:
#using Images
# Converts img to greyscale matrix representation of the image
#
# returns: mxn Float64 array
#function img_to_array(img_name)
#	img = load(img_name)
#	img = Gray.(img)
#	return convert(Array{Float64,2}, img)
#end

function load_array(file_name)
	value_arrays = []
	open(file_name) do file
		for ln in eachline(file)
			indices = [ c != "," for c in ln]
			values = split(ln, ",")
			push!(value_arrays, values)
		end
	end

	size = (parse(Int,value_arrays[1][1]),parse(Int,value_arrays[1][2]))
	vals = parse.(Int,value_arrays[2][1:end-1])

	v = reshape(vals,size[1],size[2])

end

# Used to keep pixel intensity sums lower
switcheroo(x) = x = 1-x

man_norm(x) = sum(x.^2)^(1/2)

# Partitions image into num_cols * num_rows rectangles and returns the
# the partitions as t_images and the sum of the region's pixel intensity
# as sums
function get_regions(img,s,num_cols=2,num_rows=3)
	# initialize returned values
	sums = Float64[]
	t_images = []

	for col in 1:num_cols
		for row in 1:num_rows
			# Reference the (col * row)th partition of img
			region = @view img[max(1,Int(round((row-1)*s[1]/num_rows))):Int(round(row*s[1]/num_rows)),
												max(1,Int(round((col-1)*s[2]/num_cols))):Int(round(col*s[2]/num_cols))]

			push!(t_images, region)

			push!(sums, man_norm(region))
		end
	end
	return sums, t_images
end

# Take the sums matrix and output a matrix that represents the
# braille representation
function get_output_vector(sums,tol_div=4, zero_tol=5)
	largest = maximum(sums)
	tolerance = largest/tol_div
	vect = [ abs(largest - sum) < tolerance for sum in sums]
	if largest < zero_tol
		vect = zeros(6)
	end
	return Int.(reshape(vect,3,2))
end

# Given an image of a single char we return the regions of the
# image, the sums of those regions and the output vector
#
function process_char(img)
	s = size(img)
	sums, t_images = get_regions(img,s)
	output_vector = get_output_vector(sums)
	return output_vector
end

# Preliminary version of processing the image to feed the
# single characters to the process_char function
#
function process_chars(img,s,num_chars)
	chars = [ @view img[:,max(1,Int(round((col-1)*s[2]/num_chars))):Int(round(col*s[2]/num_chars))]
						for col in 1:num_chars ]

	outputs = [ process_char(char) for char in chars ]
	return outputs
end

function row_to_num_chars(l_row)
	last_col = 0
	distance = 0
	distances = Int64[]
	dot_size = 0
	dot_sizes = Int64[]

	for col in l_row
		distance += 1
		dot_size += 1
		if col > 0 && last_col == 0
			push!(distances, distance)
			dot_size = 0

		elseif col == 0 && last_col > 0
			push!(dot_sizes, dot_size)
			distance = 0
		end

		last_col = col
	end

	return distances, maximum(dot_sizes)
end

function optimize_x(t_sum, dot_size, inner_dist, outer_dist)
	x = 2
	y = 1
	z = 0

	while dot_size*x + y*inner_dist + z*outer_dist < t_sum
		x += 2
		y = x/2
		z = x/2 - 1
	end

	return (x, y, z)
end

function get_distances(img,s)
	l_row = img[1,:]
	l_row_sum = 0
	for row_num in 1:s[1]
		row = @view img[row_num, :]
		row_sum = sum(row)
		if row_sum >= l_row_sum
			l_row = row
			l_row_sum = row_sum
		end
	end
	distances, dot_size = row_to_num_chars(l_row)
	return distances, dot_size
end

function determine_chars(distances, dot_size)

	t_sum = sum([distance + dot_size for distance in distances[2:end]])

	inner_dist = minimum(distances)
	outer_dist = maximum(distances)
	opt_x = optimize_x(t_sum, dot_size, inner_dist, outer_dist)

	return t_sum, opt_x
end

function get_min_space(distances, tolerance=3)
	s_dists = sort(distances, rev=true)
	min_space = s_dists[1]
	for col in 2:size(s_dists)[1]
		gap = abs(min_space - s_dists[col])
		if gap > min_space/tolerance
			return min_space
		end
		min_space = s_dists[col]
	end
	return Inf
end

function get_words(img, distances, dot_size, indices)
	cur_region_strt = 1
	cur_region_end = cur_region_strt

	dot_ticker = 1

	words = []
	for col in 2:size(distances)[1]
		if dot_ticker > 0
			cur_region_end += dot_size
		end
		dot_ticker = mod(dot_ticker + 1,3)
		cur_region_end += distances[col]
		if indices[col]
			word = @view img[:,cur_region_strt:cur_region_end]
			push!(words,word)
			cur_region_strt = cur_region_end
		end
	end
	word = @view img[:,cur_region_strt:end]
	push!(words,word)
	return words
end

function vert_trim(word)
	while true
		if sum(word[1,:]) > 0
			break
		end
		word = @view word[2:end,:]
	end
	while true
		if sum(word[end,:]) > 0
			break
		end
		word = @view word[1:end-1,:]
	end
	return word
end

function horiz_trim(word)
	while true
		if sum(word[:,1]) > 0
			break
		end
		word = @view word[:,2:end]
	end
	while true
		if sum(word[:,end]) > 0
			break
		end
		word = @view word[:,1:end-1]
	end

	return word
end

function add_padding(word)
	s = size(word)
	word = [ zeros(5,s[2]); word; zeros(5,s[2])]

	word = [ zeros(s[1]+10,10) word zeros(s[1]+10,10) ]
	return word
end

function trim_word(word)
	word = vert_trim(word)
	word = horiz_trim(word)
	word = add_padding(word)
	return word
end

function image_to_braille(file_name)
	img = load_array(file_name)
	img = map(switcheroo,img)

	s = size(img)
	distances, dot_size = get_distances(img,s)
	braille_chars = []

	min_space = get_min_space(distances)
	indices = [ dist >= min_space for dist in distances]

	words = get_words(img, distances, dot_size, indices)
	for word in words
		word = trim_word(word)
		distances, dot_size = get_distances(word,size(word))
		t_sum, opt_x = determine_chars(distances,dot_size)
		num_chars = opt_x[1]/2
		outputs = process_chars(word, size(word), num_chars)
		if sum(sum(outputs)) > 0
			push!(braille_chars, outputs, [[0 0; 0 0; 0 0]])
		end
	end
	return braille_chars
end


text:
"""
	$(SIGNATURES)

Utility from consumption.
"""
function utility(uS :: CRRA, cM :: AbstractArray{F1}) where F1
    if uS.dbg
        @assert all(x -> x .> zero(F1), cM)  "Negative consumption"
    end

    if uS.sigma == 1.0
        utilM = log.(cM);
    else
        utilM = utility_crra(cM, 1.0 - uS.sigma);
    end
    return utilM
end

# utility_log(cM) = log.(cM);
utility_crra(cM, oneMinusSigma) = (cM .^ oneMinusSigma) ./ oneMinusSigma .- 1.0;

function utility(uS :: CRRA, c :: T1) where T1 <: AbstractFloat
    if uS.sigma == 1.0
        utilM = log(c);
    else
        # oneMinusSigma = 1.0 - uS.sigma;
        utilM = utility_crra(c, 1.0 - uS.sigma);
    end
    return utilM
end

"""
	$(SIGNATURES)

Marginal utility.
"""
function marginal_utility(uS :: CRRA, cM :: AbstractArray{T1}) where
    T1 <: AbstractFloat

    if uS.dbg
        @assert all(x -> x > 0.0, cM)  "Negative consumption"
    end

    if uS.sigma == 1.0
        muM = one(T1) ./ cM;
    else
        muM = cM .^ (-uS.sigma);
    end
    return muM
end

function marginal_utility(uS :: CRRA, c :: T1) where T1 <: AbstractFloat
    if uS.sigma == 1.0
        muM = one(T1) / c;
    else
        muM = c ^ (-uS.sigma);
    end
    return muM
end

"""
    $(SIGNATURES)

Euler equation deviation.
"""
function euler_deviation(uS :: UtilityOneArg, beta :: T1, R :: T1,
    cV :: Vector{T1}) where T1 <: AbstractFloat

    muV = marginal_utility(uS, cV);
    T = length(muV);
    return muV[2 : T] .* beta .* R ./ muV[1 : (T-1)] .- 1.0;
end

"""
	$(SIGNATURES)

Consumption growth factor. Constant for CRRA. Otherwise this would have to depend on consumption levels.
"""
function cons_growth(uS :: CRRA,  betaR :: T1) where
    T1 <: AbstractFloat

    return betaR ^ (1.0 / uS.sigma)
end

# """
# 	$(SIGNATURES)

# Growth rate of utility.
# Input is consumption growth factor.
# """
# util_growth(uS :: CRRA, cGrowthFactor :: T1) where T1 <: AbstractFloat =
#     cGrowthFactor ^ (one(T1) - uS.sigma);

"""
    $(SIGNATURES)

Present value of consumption as multiple of c1.
Periods 1...T. First not discounted.
"""
function pv_consumption(uS :: CRRA,
    beta :: T1, R :: T1, T :: T2) where
    {T1 <: AbstractFloat, T2 <: Integer}

    g = cons_growth(uS, beta * R);
    pv = const_growth_sum(g/R, T);
    return pv
    # return ((g/R) ^ T - one(T1)) / (g/R - one(T1))
end

"""
	$(SIGNATURES)

Consumption at age 1 for given lifetime income.
"""
function cons_age1(uS :: CRRA, beta :: T1, R :: T1, T :: T2,
    ltIncome) where {T1 <: AbstractFloat, T2 <: Integer}

    c1 = ltIncome ./ pv_consumption(uS, beta, R, T);
    return c1
end

"""
    $(SIGNATURES)

Consumption path.
"""
function cons_path(uS :: CRRA, beta :: T1, R :: T1, T :: T2,
    ltIncome :: T1) where {T1 <: AbstractFloat, T2 <: Integer}

    c1 = cons_age1(uS, beta, R, T, ltIncome);
    g = cons_growth(uS, beta .* R);
    cV = const_growth_vector(c1, g, T);
    return cV
    # return c1 .* (g .^ (0 : (T-1)))
end

"""
    $(SIGNATURES)

Lifetime utility, given present value of income.

Note that the growth rate of utility is not constant (because of the -1 in the utility function). So lifetime utility has to be computed from the consumption or utility path.
"""
function lifetime_utility(uS :: CRRA, beta :: T1, R :: T1, T :: T2,
    ltIncome :: T1) where {T1 <: AbstractFloat, T2 <: Integer}

    if uS.dbg
        @assert ltIncome .> 0.0  "Negative income"
    end

    # cGrowthFactor = cons_growth(uS, beta * R);
    # c1 = cons_age1(uS, beta, R, T, ltIncome);

    # if uS.sigma == 1
    #   uFct = c -> log.(c);
    # else
    #   uFct = c -> utility_crra(c, 1.0 - uS.sigma);
    # end
    # util1 = uFct(c1);
    # # util1 = utility(uS, c1);

    # ltu = util1;
    # ct = c1;
    # betaFactor = one(T1);
    # for t = 2 : T
    #     ct *= cGrowthFactor;
    #     betaFactor *= beta;
    #     ltu += betaFactor * uFct(ct);
    #     # utility(uS, ct);
    # end
    # return ltu
    utilV = utility(uS, cons_path(uS, beta, R, T, ltIncome));
    return lifetime_utility(utilV, beta, T)
end

"""
	$(SIGNATURES)

Marginal utility of initial assets given present value of income.
"""
function mu_wealth(uS :: CRRA, beta :: T1, R :: T1, T :: T2,
    ltIncome) where {T1 <: AbstractFloat, T2 <: Integer}

    c1 = cons_age1(uS, beta, R, T, ltIncome);
    return marginal_utility(uS, c1)
end

"""
	$(SIGNATURES)

Lifetime utility for a given utility sequence.
"""
function lifetime_utility(utilV :: Vector{T1},  beta :: T1,  T :: Integer)  where T1 <: AbstractFloat
    return discounted_sum(beta, utilV);
end

# Same for a constant utility flow
function lifetime_utility(util :: T1, beta :: T1, T :: Integer) where T1 <: AbstractFloat

    return util * (beta ^ T - 1.0) / (beta - 1.0)
end

# -----------


text:
# some convenient type-piracy for scalars
AD.gradient(ad::AD.AbstractBackend, f, x::Real) = AD.derivative(ad, f, x)
AD.hessian(ad::AD.AbstractBackend, f, x::Real) = first.(AD.hessian(ad, f∘first, [x]))

# ForwardDiffBackend is built-in ot AutomaticDifferentiation
@init @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
    @eval const ForwardDiffBackend = AD.ForwardDiffBackend
end

# Zygote not yet in AutomaticDifferentiation, but remove this once it is
@init @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    struct ZygoteBackend <: AD.AbstractBackend end
    AD.gradient(ad::ZygoteBackend, f, x::Real) = first(Zygote.gradient(f, x))
    AD.gradient(ad::ZygoteBackend, f, x) = Zygote.gradient(f, x)
    function AD.value_and_gradient(::ZygoteBackend, f, x)
        y, back = Zygote.pullback(f, x)
        y, back(Zygote.sensitivity(y))
    end
end

function optim_only_fg!(func, autodiff)
    Optim.only_fg!() do F, G, z
        if G != nothing
            f, g = AD.value_and_gradient(autodiff, func, z)
            G .= first(g)
            return f
        end
        if F != nothing
            return neglogp(z)
        end
    end
end


text:
# Glue functions for natively using EmpiricalBayes with NetworkInference structures.
# Will not load if EmpiricalBayes does not exist.

# Check for EmpiricalBayes package
EB_EXISTS = in("EmpiricalBayes",keys(Pkg.installed())) ? true : false

# Only load if package exists
if EB_EXISTS

import EmpiricalBayes.empirical_bayes # for overloading

"""
    to_index(label1::String, label2::String)

Convert any ordering of two strings into a unique tuple for that pair.
"""
function to_index(label1::AbstractString, label2::AbstractString)
    if label1 > label2
        return (label1, label2)
    else
        return (label2, label1)
    end
end

"""
    to_index(node1::Node, node2::Node)

Convert any ordering of two Node objects into a unique tuple for that pair.
"""
function to_index(node1::Node, node2::Node)
    return to_index(node1.label, node2.label)
end

"""
    to_index(nodes::Array{Node, 1})

Convert any ordering of two Node objects into a unique tuple for that pair.
"""
function to_index(nodes::Array{Node, 1})
    n1, n2 = nodes
    return to_index(n1.label, n2.label)
end

"""
    make_priors(filepath::String, weight = 1.0, skiplines = 1)

Convert a file containing prior data into a dictionary of (id, prior) pairs.

# Arguments
* `filepath` : path to prior file. File should contain data in three columns as
  follows: node1, node2, prior.
* `weight` : weight to apply to the prior.
* `skiplines = 1` : number of initial lines to skip in the file. Defaults to 1, to
  skip a header line.
"""
function make_priors(filepath::String, weight = 1.0, skiplines = 1)
    prior_mat = readdlm(filepath, skipstart = skiplines)

    num_priors = size(prior_mat, 1)

    prior_dict = Dict()
    for i in 1:num_priors
        n1, n2, p = prior_mat[i, :]
        n1 = string(n1)
        n2 = string(n2)
        index = to_index(n1, n2)
        if haskey(prior_dict, index)
            prior_dict[index] = max(p * weight, prior_dict[index])
        else
            prior_dict[index] = p * weight
        end
    end

    return prior_dict
end

"""
    make_priors(filepaths::Array{String}, weights = ones(length(filepaths)), skiplines = 1)

Convert a file containing prior data into a dictionary of (id, prior) pairs.

# Arguments
* `filepaths` : paths to prior files. Each file should contain data in three columns as
  follows: node1, node2, prior.
* `weights` : weights to apply to each prior, in the same order. Defaults to array of ones.
* `skiplines = 1` : number of initial lines to skip in the file. Defaults to 1, to
  skip a header line.
"""
function make_priors(filepaths::Array{String}, weights = ones(length(filepaths)), skiplines = 1)
    num_filepaths = length(filepaths)

    # Check there is a weight for each prior file; if not, fall back to default
    weights = length(weights) == num_filepaths ? weights : ones(num_filepaths)

    prior_dicts = make_priors.(filepaths, weights, skiplines)
    prior_dict = merge(+, prior_dicts...)

    return prior_dict
end

"""
    empirical_bayes(network::InferredNetwork, priors::Dict, num_bins, distr::Symbol, proportion_to_keep = 1.0, key_func = to_index)

Calculate the empirical Bayes posteriors of the input statistics using the priors.

# Arguments
* `network::InferredNetwork` : network to apply priors to.
* `priors::Dict` : dictionary of priors such that looking up a pair of nodes
  returns the prior value for the edge between them.
* `num_bins` : number of uniform width bins to discretize into.
* `distr` : form of the null distribution to be fitted.
* `proportion_to_keep = 1.0` : Proportion of lowest test statistics to
   keep when calculating null distribution.
* `key_func = to_index` : a function mapping an Edge object to a key useable in
   the `priors` dictionary.
* `tail` : Whether the test is two-tailed (:two) or one-tailed (:lower or :upper).
- `w0` : Default constant for the prior calculation
"""
function empirical_bayes(network::InferredNetwork, priors::Dict, num_bins, distr::Symbol; proportion_to_keep = 1.0, key_func = to_index,
    tail = :upper, w0 = 2.2)

    edge_list = network.edges
    test_statistics = [e.weight for e in edge_list]
    prior_list = [ get(priors, key_func(e.nodes), 0) for e in edge_list ]

    eb_edges = Array{Edge}(undef, length(edge_list))

    posteriors = empirical_bayes(test_statistics, prior_list, num_bins, distr, proportion_to_keep = proportion_to_keep,
        tail = tail, w0 = w0)

    for i in 1:length(edge_list)
        nodes = edge_list[i].nodes
        eb_edges[i] = Edge(nodes, posteriors[i])
    end

    # Remove infinite values
    eb_edges = filter(x->isfinite(x.weight), eb_edges)

    sort!(eb_edges, rev = true, by = x->x.weight)

    return InferredNetwork(network.nodes, eb_edges)
end

"""
  empirical_bayes(network::InferredNetwork num_bins, distr::Symbol, proportion_to_keep = 1.0, key_func = to_index)

Calculate the empirical Bayes posteriors of the input statistics with no priors.

# Arguments
* `network::InferredNetwork` : network to apply priors to.
* `num_bins` : number of uniform width bins to discretize into.
* `distr` : form of the null distribution to be fitted.
* `proportion_to_keep = 1.0` : proportion of lowest test statistics to
   keep when calculating null distribution.
* `key_func = to_index` : a function mapping an Edge object to a key useable in
   the `priors` dictionary.
* `tail` : whether the test is two-tailed (:two) or one-tailed (:lower or :upper).
- `w0` : Default constant for the prior calculation
"""
function empirical_bayes(network::InferredNetwork, num_bins, distr::Symbol; proportion_to_keep = 1.0, key_func = to_index, tail = :upper, w0 = 2.2)
    return empirical_bayes(network, Dict(), num_bins, distr, proportion_to_keep = proportion_to_keep, key_func = key_func, tail = tail, w0 = w0)
end

end


text:
#####
##### Tendency specification
#####

#####
##### Sources
#####

eq_tends(pv::AbstractPrognosticVariable, m::AtmosModel, tt::Source) =
    (m.source[pv]..., eq_tends(pv, turbconv_model(m), tt)...)

#####
##### First order fluxes
#####

eq_tends(::Mass, ::Anelastic1D, ::Flux{FirstOrder}) = ()

eq_tends(::Mass, ::Compressible, ::Flux{FirstOrder}) = (Advect(),)

# Mass
eq_tends(pv::Mass, atmos::AtmosModel, tt::Flux{FirstOrder}) =
    (eq_tends(pv, compressibility_model(atmos), tt))

# Momentum
eq_tends(pv::Momentum, ::Compressible, ::Flux{FirstOrder}) =
    (PressureGradient(),)

eq_tends(pv::Momentum, ::Anelastic1D, ::Flux{FirstOrder}) = ()

eq_tends(pv::Momentum, m::AtmosModel, tt::Flux{FirstOrder}) =
    (Advect(), eq_tends(pv, compressibility_model(m), tt)...)

# Energy
eq_tends(::Energy, m::TotalEnergyModel, tt::Flux{FirstOrder}) =
    (Advect(), Pressure())

eq_tends(::ρθ_liq_ice, m::θModel, tt::Flux{FirstOrder}) = (Advect(),)

# TODO: make radiation aware of which energy formulation is used:
# eq_tends(pv::PV, m::AtmosModel, tt::Flux{FirstOrder}) where {PV <: AbstractEnergyVariable} =
#     (eq_tends(pv, energy_model(m), tt)..., eq_tends(pv, energy_model(m), radiation_model(m), tt)...)
eq_tends(pv::AbstractEnergyVariable, m::AtmosModel, tt::Flux{FirstOrder}) = (
    eq_tends(pv, energy_model(m), tt)...,
    eq_tends(pv, radiation_model(m), tt)...,
)

# AbstractMoistureVariable
eq_tends(::AbstractMoistureVariable, ::AtmosModel, ::Flux{FirstOrder}) =
    (Advect(),)

# AbstractPrecipitationVariable
eq_tends(
    pv::AbstractPrecipitationVariable,
    m::AtmosModel,
    tt::Flux{FirstOrder},
) = (eq_tends(pv, precipitation_model(m), tt)...,)

# Tracers
eq_tends(pv::Tracers{N}, ::AtmosModel, ::Flux{FirstOrder}) where {N} =
    (Advect(),)

#####
##### Second order fluxes
#####

eq_tends(
    ::Union{Mass, Momentum, AbstractMoistureVariable},
    ::DryModel,
    ::Flux{SecondOrder},
) = ()
eq_tends(
    ::Union{Mass, Momentum, AbstractMoistureVariable},
    ::AbstractMoistureModel,
    ::Flux{SecondOrder},
) = (MoistureDiffusion(),)

# Mass
eq_tends(pv::Mass, m::AtmosModel, tt::Flux{SecondOrder}) =
    (eq_tends(pv, moisture_model(m), tt)...,)

# Momentum
eq_tends(pv::Momentum, m::AtmosModel, tt::Flux{SecondOrder}) = (
    ViscousStress(),
    eq_tends(pv, moisture_model(m), tt)...,
    eq_tends(pv, turbconv_model(m), tt)...,
    eq_tends(pv, hyperdiffusion_model(m), tt)...,
)

# Energy
eq_tends(::Energy, m::TotalEnergyModel, tt::Flux{SecondOrder}) =
    (ViscousFlux(), DiffEnthalpyFlux())

eq_tends(::ρθ_liq_ice, m::θModel, tt::Flux{SecondOrder}) = (ViscousFlux(),)

eq_tends(pv::AbstractEnergyVariable, m::AtmosModel, tt::Flux{SecondOrder}) = (
    eq_tends(pv, energy_model(m), tt)...,
    eq_tends(pv, turbconv_model(m), tt)...,
    eq_tends(pv, hyperdiffusion_model(m), tt)...,
)

# AbstractMoistureVariable
eq_tends(pv::AbstractMoistureVariable, m::AtmosModel, tt::Flux{SecondOrder}) = (
    eq_tends(pv, moisture_model(m), tt)...,
    eq_tends(pv, turbconv_model(m), tt)...,
    eq_tends(pv, hyperdiffusion_model(m), tt)...,
)

# AbstractPrecipitationVariable
eq_tends(
    pv::AbstractPrecipitationVariable,
    m::AtmosModel,
    tt::Flux{SecondOrder},
) = (eq_tends(pv, precipitation_model(m), tt)...,)

# Tracers
eq_tends(pv::Tracers{N}, m::AtmosModel, tt::Flux{SecondOrder}) where {N} =
    (eq_tends(pv, tracer_model(m), tt)...,)


text:
using FromFile
@from "test_params.jl" import maximum_residual

using Distributed
procs = addprocs(4)
using Test, Pkg
project_path = splitdir(Pkg.project().path)[1]
@everywhere procs begin
    Base.MainInclude.eval(quote
        using Pkg
        Pkg.activate($$project_path)
    end)
end
@everywhere using SymbolicRegression
@everywhere _inv(x::Float32)::Float32 = 1f0/x
X = rand(Float32, 5, 100) .+ 1
y = 1.2f0 .+ 2 ./ X[3, :]

options = SymbolicRegression.Options(
    binary_operators=(+, *),
    unary_operators=(_inv,),
    npopulations=8
)
hallOfFame = EquationSearch(X, y, niterations=8, options=options, procs=procs)
rmprocs(procs)

dominating = calculateParetoFrontier(X, y, hallOfFame, options)
best = dominating[end]
# Test the score
@test best.score < maximum_residual / 10


text:
export BasicBlock, unsafe_delete!,
       terminator, name,
       move_before, move_after

@checked struct BasicBlock <: Value
    ref::reftype(Value)
end
identify(::Type{Value}, ::Val{API.LLVMBasicBlockValueKind}) = BasicBlock

# forward declarations
@checked struct Function <: GlobalObject
    ref::reftype(GlobalObject)
end

BasicBlock(ref::API.LLVMBasicBlockRef) = BasicBlock(API.LLVMBasicBlockAsValue(ref))
blockref(bb::BasicBlock) = API.LLVMValueAsBasicBlock(ref(bb))

BasicBlock(f::Function, name::String) =
    BasicBlock(API.LLVMAppendBasicBlock(ref(f), name))
BasicBlock(f::Function, name::String, ctx::Context) =
    BasicBlock(API.LLVMAppendBasicBlockInContext(ref(ctx), ref(f), name))
BasicBlock(bb::BasicBlock, name::String) =
    BasicBlock(API.LLVMInsertBasicBlock(blockref(bb), name))
BasicBlock(bb::BasicBlock, name::String, ctx::Context) =
    BasicBlock(API.LLVMInsertBasicBlockInContext(ref(ctx), blockref(bb), name))

unsafe_delete!(::Function, bb::BasicBlock) = API.LLVMDeleteBasicBlock(blockref(bb))
Base.delete!(::Function, bb::BasicBlock) =
    API.LLVMRemoveBasicBlockFromParent(blockref(bb))

parent(bb::BasicBlock) = Function(API.LLVMGetBasicBlockParent(blockref(bb)))

terminator(bb::BasicBlock) = Instruction(API.LLVMGetBasicBlockTerminator(blockref(bb)))

name(bb::BasicBlock) =
    unsafe_string(API.LLVMGetBasicBlockName(blockref(bb)))

move_before(bb::BasicBlock, pos::BasicBlock) =
    API.LLVMMoveBasicBlockBefore(blockref(bb), blockref(pos))
move_after(bb::BasicBlock, pos::BasicBlock) =
    API.LLVMMoveBasicBlockAfter(blockref(bb), blockref(pos))

## instruction iteration

export instructions

struct BasicBlockInstructionSet
    bb::BasicBlock
end

instructions(bb::BasicBlock) = BasicBlockInstructionSet(bb)

Base.eltype(::BasicBlockInstructionSet) = Instruction

function Base.iterate(iter::BasicBlockInstructionSet,
                      state=API.LLVMGetFirstInstruction(blockref(iter.bb)))
    state == C_NULL ? nothing : (Instruction(state), API.LLVMGetNextInstruction(state))
end

Base.last(iter::BasicBlockInstructionSet) =
    Instruction(API.LLVMGetLastInstruction(blockref(iter.bb)))

Base.isempty(iter::BasicBlockInstructionSet) =
    API.LLVMGetLastInstruction(blockref(iter.bb)) == C_NULL

Base.IteratorSize(::BasicBlockInstructionSet) = Base.SizeUnknown()


text:
module AWSBatch

using AWS
using AutoHashEquals
using OrderedCollections: OrderedDict
using Dates
using Memento
using Mocking

@service Batch
@service Cloudwatch_Logs

export BatchJob, ComputeEnvironment, BatchEnvironmentError, BatchJobError
export JobQueue, JobDefinition, JobState, LogEvent
export run_batch, describe, status, status_reason, wait, log_events, isregistered, register, deregister
export list_job_queues, list_job_definitions, create_compute_environment, create_job_queue

const logger = getlogger(@__MODULE__)
# Register the module level logger at runtime so that folks can access the logger via `getlogger(MyModule)`
# NOTE: If this line is not included then the precompiled `MyModule.logger` won't be registered at runtime.
__init__() = Memento.register(logger)

include("exceptions.jl")
include("log_event.jl")
include("compute_environment.jl")
include("job_queue.jl")
include("job_state.jl")
include("job_definition.jl")
include("batch_job.jl")

"""
    run_batch(;
        name::AbstractString="",
        queue::AbstractString="",
        region::AbstractString="",
        definition::Union{AbstractString, JobDefinition, Nothing}=nothing,
        image::AbstractString="",
        vcpus::Integer=1,
        memory::Integer=-1,
        role::AbstractString="",
        cmd::Cmd=``,
        num_jobs::Integer=1,
        parameters::Dict{String, String}=Dict{String, String}(),
    ) -> BatchJob

Handles submitting a BatchJob based on various potential defaults.
For example, default job fields can be inferred from an existing job definition or an
existing job (if currently running in a batch job).

Order of priority from highest to lowest:

1. Explict arguments passed in via `kwargs`.
2. Inferred environment (e.g., `AWS_BATCH_JOB_ID` environment variable set)
3. Job definition parameters

If no valid job definition exists (see [`AWSBatch.job_definition_arn`](@ref) then a new job
definition will be created and registered based on the job parameters.
"""
function run_batch(;
    name::AbstractString="",
    queue::AbstractString="",
    region::AbstractString="",
    definition::Union{AbstractString, JobDefinition, Nothing}=nothing,
    image::AbstractString="",
    vcpus::Integer=1,
    memory::Integer=-1,
    role::AbstractString="",
    cmd::Cmd=``,
    num_jobs::Integer=1,
    parameters::Dict{String, String}=Dict{String, String}(),
    allow_job_registration::Bool=true,
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    if isa(definition, AbstractString)
        definition = isempty(definition) ? nothing : definition
    end

    # Determine if the job definition already exists and update the default job parameters
    if definition !== nothing
        response = describe_job_definition(definition; aws_config=aws_config)
        if !isempty(response["jobDefinitions"])
            details = first(response["jobDefinitions"])

            container = details["containerProperties"]
            isempty(image) && (image = container["image"])
            isempty(role) && (role = container["jobRoleArn"])

            # Update container override parameters
            vcpus == 1 && (vcpus = container["vcpus"])
            memory < 0 && (memory = container["memory"])
            isempty(cmd) && (cmd = Cmd(Vector{String}(container["command"])))
        end
    end

    # Get inferred environment parameters
    if haskey(ENV, "AWS_BATCH_JOB_ID")
        # Environmental variables set by the AWS Batch service. They were discovered by
        # inspecting the running AWS Batch job in the ECS task interface.
        job_id = ENV["AWS_BATCH_JOB_ID"]
        job_queue = ENV["AWS_BATCH_JQ_NAME"]

        # if not specified, get region from the aws_config
        isempty(region) && (region = aws_config.region)

        # Requires permissions to access to "batch:DescribeJobs"
        response = @mock Batch.describe_jobs([job_id]; aws_config=aws_config)

        # Use the job's description to only update fields that are using the default
        # values since explict arguments passed in via `kwargs` have higher priority
        if length(response["jobs"]) > 0
            details = first(response["jobs"])

            # Update the job's required parameters
            isempty(name) && (name = details["jobName"])
            definition === nothing && (definition = details["jobDefinition"])
            isempty(queue) && (queue = job_queue)

            # Update the container parameters
            container = details["container"]
            isempty(image) && (image = container["image"])
            isempty(role) && (role = container["jobRoleArn"])

            # Update container overrides
            vcpus == 1 && (vcpus = container["vcpus"])
            memory < 0 && (memory = container["memory"])
            isempty(cmd) && (cmd = Cmd(Vector{String}(container["command"])))
        else
            warn(logger, "No jobs found with id: $job_id.")
        end
    end

    # Error if required parameters were not explicitly set and cannot be inferred
    if isempty(name) || isempty(queue) || memory < 0
        throw(BatchEnvironmentError(
            "Unable to perform AWS Batch introspection when not running within " *
            "an AWS Batch job. Current job parameters are: " *
            "\nname=$name" *
            "\nqueue=$queue" *
            "\nmemory=$memory"
        ))
    end

    # Reuse a previously registered job definition if available.
    if isa(definition, AbstractString)
        reusable_job_definition_arn = job_definition_arn(definition; image=image, role=role)

        if reusable_job_definition_arn !== nothing
            definition = JobDefinition(reusable_job_definition_arn)
        end
    elseif definition === nothing
         # Use the job name as the definiton name since the definition name was not specified
        definition = name
    end

    # If no job definition exists that can be reused, a new job definition is created
    # under the current job specifications.
    if isa(definition, AbstractString)
        if allow_job_registration
            definition = register(
                definition;
                image=image,
                role=role,
                vcpus=vcpus,
                memory=memory,
                cmd=cmd,
                parameters=parameters,
                aws_config=aws_config,
            )
        else
            throw(BatchEnvironmentError(string(
                "Attempting to register job definition \"$definition\" but registering ",
                "job definitions is disallowed. Current job definition parameters are: ",
                "\nimage=$image",
                "\nrole=$role",
                "\nvcpus=$vcpus",
                "\nmemory=$memory",
                "\ncmd=$cmd",
                "\nparameters=$parameters",
            )))
        end
    end

    # Parameters that can be overridden are `memory`, `vcpus`, `command`, and `environment`
    # See https://docs.aws.amazon.com/batch/latest/APIReference/API_ContainerOverrides.html
    container_overrides = Dict(
        "vcpus" => vcpus,
        "memory" => memory,
        "command" => cmd.exec,
    )

    return submit(
        name,
        definition,
        JobQueue(queue);
        container=container_overrides,
        parameters=parameters,
        num_jobs=num_jobs,
    )
end

end  # AWSBatch


text:
abstract type Selector end

matches(sel::Selector, x) = false

struct SelectCond <: Selector
    f
end
matches(sel::SelectCond, x) = sel.f(x)

struct SelectAll <: Selector end
matches(::SelectAll, x) = true

struct SelectNode <: Selector end
matches(::SelectNode, ::XNode) = true
matches(::SelectNode, _) = false

struct SelectLeaf <: Selector end
matches(::SelectLeaf, ::XLeaf) = true
matches(::SelectLeaf, _) = false

struct SelectTag <: Selector
    tag::Symbol
end

matches(sel::SelectTag, x::XNode) = sel.tag == tag(x)

struct SelectOr{T<:Tuple} <: Selector
    selectors::T
end

Base.:(|)(sel1::Selector, sel2::Selector) = SelectOr((sel1, sel2))
Base.:(|)(selor::SelectOr, sel::Selector) = SelectOr((selor.seletors..., sel))

matches(sel::SelectOr, x) = any(matches(s, x) for s in sel.selectors)

struct SelectAnd{T<:Tuple} <: Selector
    selectors::T
end

Base.:(&)(sel1::Selector, sel2::Selector) = SelectAnd((sel1, sel2))
Base.:(&)(seland::SelectAnd, sel::Selector) = SelectAnd((seland.selectors..., sel))

matches(sel::SelectAnd, x) = all(matches(s, x) for s in sel.selectors)

struct SelectNot <: Selector
    sel::Selector
end
matches(sel::SelectNot, x) = !matches(sel.sel, x)

Base.:(!)(sel::Selector) = SelectNot(sel)

"""
    SelectAttrEq(attr, val)

Select an x-expression with attribute `attr == val`
"""
struct SelectAttrEq{T} <: Selector
    attr::Symbol
    val::T
end

matches(sel::SelectAttrEq, x::XNode) = get(attributes(x), sel.attr, nothing) == sel.val

struct SelectHasAttr <: Selector
    attr::Symbol
end

matches(sel::SelectHasAttr, x::XNode) = haskey(attributes(x), sel.attr)

## API

select(xtree::XTree, sel::Selector) = (x for x in PostOrderDFS(xtree) if matches(sel, x))

function selectfirst(xtree::XTree, sel::Selector)
    for x in PostOrderDFS(xtree)
        if matches(sel, x)
            return x
        end
    end
    return nothing
end


text:
@testset "gaussian state bhd" begin
    dim = 100
    ρ = SqueezedState(0.8, π, Matrix, dim=dim)
    d = GaussianStateBHD(ρ)

    θ, x = rand(d)
    @test 0 ≤ θ < 2π

    points = rand(d, 4096)
    @test all(0 ≤ θ < 2π for θ in points[1, :])
end

@testset "operator" begin
    T = Float64
    dim = 100

    ρ = Matrix{T}(I, dim, dim)
    c = QSB.Creation(T, dim=dim)
    a = QSB.Annihilation(T, dim=dim)

    @test QSD.mean_of_create(ρ) ≈ tr(c * ρ)
    @test QSD.mean_of_create²(ρ) ≈ tr(c^2 * ρ)
    @test QSD.mean_of_annihilate(ρ) ≈ tr(a * ρ)
    @test QSD.mean_of_annihilate²(ρ) ≈ tr(a^2 * ρ)
    @test QSD.mean_of_create_annihilate(ρ) ≈ tr(c*a * ρ)

    θ = π/8
    @test QSD.mean_of_π̂ₓ²(ρ, θ) ≈ (
        exp(-2im*θ) * tr(a^2 * ρ) +
        exp(2im*θ) * tr(c^2 * ρ) +
        1 + 2tr(c*a * ρ)
    ) / 4
    @test QSD.mean_of_π̂ₓ(ρ, θ) ≈ (
        exp(-im*θ) * tr(a * ρ) +
        exp(im*θ) * tr(c * ρ)
    ) / 2

    θs = LinRange(0, 2π, 4096)
    @test QSD.mean_of_π̂ₓ²(ρ, θs) ≈ (
        exp.(-2im*θs) .* tr(a^2 * ρ) .+
        exp.(2im*θs) .* tr(c^2 * ρ) .+
        1 .+ 2tr(c*a * ρ)
    ) ./ 4
    @test QSD.mean_of_π̂ₓ(ρ, θs) ≈ (
        exp.(-im*θs) .* tr(a * ρ) .+
        exp.(im*θs) .* tr(c * ρ)
    ) ./ 2
end


text:
# Module with functions for time-domain modeling and inversion using OPESCI/devito
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January, 2017
#

module TimeModeling
using PyCall, JOLI, SeisIO
@pyimport numpy as np
@pyimport Jcontainers as ct
@pyimport JAcoustic_codegen as ac
@pyimport JTTI_codegen as tc
@pyimport devito as dvt
import Base.*, Base./, Base.+, Base.-, Base.ctranspose, Base.conj, Base.vcat, Base.vec, Base.dot, Base.norm, Base.abs, Base.getindex

#############################################################################
# Containers
include("ModelStructure.jl")	# model container
include("InfoStructure.jl")	# basic information required by all operators
include("GeometryStructure.jl")	# source or receiver setup, recording time and sampling
include("OptionsStructure.jl")
include("auxiliaryFunctions.jl")

#############################################################################
# Abstract vectors
include("joRHS.jl")	# RHS to be multiplied with linear operator
include("joData.jl")	# Julia data container

#############################################################################
# PDE solvers
include("time_modeling_serial.jl")	# forward/adjoint linear/nonlinear modeling
include("time_modeling_serial_rho.jl")	# forward/adjoint linear/nonlinear modeling w/ density
include("time_modeling_parallel.jl")	# parallelization for modeling
include("fwi_objective_serial.jl")	# FWI objective function value and gradient
include("fwi_objective_serial_rho.jl")  # with density
include("fwi_objective_parallel.jl")	# parallelization for FWI gradient

#############################################################################
# Linear operators
include("joModeling.jl")	# nonlinear modeling operator F (no projection operators)
include("joProjection.jl")	# source/receiver projection operator
include("joPDEfull.jl")	# modeling operator with source and receiver projection: P*F*P'
include("joPDE.jl")	# modeling operator with lhs projection only: P*F
include("joJacobian.jl")	# linearized modeling operator J

end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
const DUMMY_PRECISION = 1e-12

Base.isfinite(a::AbstractVec) = all(isfinite.(a))
Base.isinf(a::AbstractVec) = any(isinf.(a))

function StaticArrays.similar_type(::Type{V}, ::Type{F}, size::Size) where {V<:AbstractVec, F<:AbstractFloat}
    # TODO: implement this for other types
    if size == Size(V) && eltype(V) == F
        return V
    else
        error("StaticArrays.similar_type is not yet implemented for Vec except for Float64 and matching size.")
    end
end

Base.abs(a::AbstractVec) = error("abs(v::AbstractVec) has been removed. Use norm(v) to get the norm; abs.(v) to get the element-wise absolute value.")