text:
## CATEGORICAL ARRAY DECODER UTILITY

"""
    CategoricalDecoder(C::CategoricalArray; eltype=nothing)

Construct a decoder for transforming a `CategoricalArray{T}` object
into an ordinary array, and for re-encoding similar arrays back into a
`CategoricalArray{T}` object having the same `pool` (and, in
particular, the same levels) as `C`. If `eltype` is not specified then
the element type of the transformed array is `T`. Otherwise, the
element type is `eltype` and the elements are promotions of the
internal (integer) `ref`s of the `CategoricalArray`. One
must have `R <: eltype <: Real` where `R` is the reference type of the
`CategoricalArray` (usually `UInt32`).

    transform(decoder::CategoricalDecoder, C::CategoricalArray)

Transform `C` into an ordinary `Array`.

    inverse_transform(decoder::CategoricalDecoder, A::Array)

Transform an array `A` suitably compatible with `decoder` into a
`CategoricalArray` having the same `pool` as `C`.

### Example

````
julia> using CategoricalArrays
julia> C = categorical(["a" "b"; "a" "c"])
2×2 CategoricalArray{String,2,UInt32}:
 "a"  "b"
 "a"  "c"

julia> decoder = MLJBase.CategoricalDecoder(C, eltype=Float64);
julia> A = transform(decoder, C)
2×2 Array{Float64,2}:
 1.0  2.0
 1.0  3.0

julia> inverse_transform(decoder, A[1:1,:])
1×2 CategoricalArray{String,2,UInt32}:
 "a"  "b"

julia> levels(ans)
3-element Array{String,1}:
 "a"
 "b"
 "c"

"""
struct CategoricalDecoder{I<:Real,T,N,R<:Integer}  # I the output eltype
    pool::CategoricalPool{T,R} # abstract type, not optimal
    use_original_type::Bool
    CategoricalDecoder{I,T,N,R}(X::CategoricalArray{T,N,R}, use_original_type) where {I,T,N,R}  =
        new(X.pool, use_original_type)
end

function CategoricalDecoder(X::CategoricalArray{T,N,R}; eltype=nothing) where {T,N,R}
    if eltype ==  nothing
        eltype = R # any integer type will do here; not used
        use_original_type = true
    else
        use_original_type = false
    end
    return CategoricalDecoder{eltype,T,N,R}(X, use_original_type)
end

function transform(decoder::CategoricalDecoder{I,T,N,R}, C::CategoricalArray) where {I,T,N,R}
    if decoder.use_original_type
        return collect(C)
    else
        return broadcast(C.refs) do element
            ref = convert(I, element)
        end
    end
end

function inverse_transform(decoder::CategoricalDecoder{I,T,N,R}, A::Array{J}) where {I,T,N,R,J<:Union{I,T}}
    if decoder.use_original_type
        refs = broadcast(A) do element
            decoder.pool.invindex[element]
        end
    else
        refs = broadcast(A) do element
            round(R, element)
        end
    end
    return CategoricalArray{T,N}(refs, decoder.pool)

end

## UTILITY FOR CONVERTING TABULAR DATA INTO MATRIX

""""
    MLJBase.matrix(X)

Convert an iteratable table source `X` into an `Matrix`; or, if `X` is
a `Matrix`, return `X`.

"""
function matrix(X)
    TableTraits.isiterabletable(X) || error("Argument is not an iterable table.")

    df = @from row in X begin
        @select row
        @collect DataFrames.DataFrame
    end
    return convert(Matrix, df)

end

matrix(X::Matrix) = X

## TOOLS FOR INDEXING QUERYVERSE ITERABLE TABLES

struct Rows end
struct Cols end
struct Schema
    nrows::Int
    ncols::Int
    names::Vector{Symbol}
    eltypes::Vector{DataType}
end

# for accessing tabular data:
retrieve(X, args...) = retrieve(Val(TableTraits.isiterabletable(X)), X, args...)
retrieve(::Val{false}, X, args...) = error("Argument is not an iterable table.")

# Note: to `retrieve` from matrices and other non-iterabletable data,
# we overload the second method above.

# project named tuple onto a tuple with only specified `labels` or indices:
project(t::NamedTuple, labels::AbstractArray{Symbol}) =
    NamedTuple{tuple(labels...)}(tuple([getproperty(t, n) for n in labels]...))
project(t::NamedTuple, label::Symbol) = project(t, [label,])
project(t::NamedTuple, indices::AbstractArray{<:Integer}) =
    NamedTuple{tuple(keys(t)[indices]...)}(tuple([t[i] for i in indices]...))
project(t::NamedTuple, i::Integer) = project(t, [i,])

# to select columns `c` of any tabular data `X` with `retrieve(X, Cols, c)`:
# CURRENTLY RETURNS A DATAFRAME AND NOT THE ORIGINAL TYPE
function retrieve(::Val{true}, X::T, ::Type{Cols}, c::AbstractArray{I}) where {T,I<:Union{Symbol,Integer}}

    row_iterator = @from row in X begin
        @select project(row, c)
        @collect DataFrames.DataFrame # `T` does not work here unless typeof(X) has no type-parameters
    end

end

# to select a single column `c` of any tabular data `X` with
# `retrieve(X, Cols, c)`; note this will try to return a
# CategoricalArray, and does so for X a DataFrame and TypedTable.Table:
function retrieve(::Val{true}, X::T, ::Type{Cols}, c::I) where {T, I<:Union{Symbol,Integer}}

    row_iterator = @from row in X begin
        @select project(row, c)
        @collect DataFrames.DataFrame
    end

    return row_iterator[1]

end

# to select rows `r` of any tabular data `X` with `retrieve(X, Rows, c)`:
# CURRENTLY RETURNS A DATAFRAME AND NOT THE ORIGINAL TYPE
function retrieve(::Val{true}, X::T, ::Type{Rows}, r) where T

    row_iterator = @from row in X begin
        @select row
        @collect
    end

    return @from row in row_iterator[r] begin
        @select row
        @collect DataFrames.DataFrame # `T` does not work here unless typeof(X) has no type-parameters
    end

end

# to get the number of nrows, ncols, feature names and eltypes of an
# tabular data:
function retrieve(::Val{true}, X, ::Type{Schema})

    TableTraits.isiterabletable(X) || error("Argument is not an iterable table.")

    row_iterator = @from row in X begin
        @select row
        @collect
    end

    nrows = length(row_iterator)
    if nrows == 0
        ncols = 0
        _names = Symbol[]
        _eltypes = DataType[]
    else
        row = first(row_iterator) # a named tuple
        ncols = length(row)
        _names = keys(row) |> collect
        _eltypes = DataType[eltype(x) for x in row]
    end

    return Schema(nrows, ncols, _names, _eltypes)

end

# retrieve(df::JuliaDB.IndexedTable, ::Type{Rows}, r) = df[r]
# retrieve(df::JuliaDB.NextTable, ::Type{Cols}, c) = select(df, c)
# retrieve(df::JuliaDB.NextTable, ::Type{Names}) = getfields(typeof(df.columns.columns))
# retrieve(df::JuliaDB.NextTable, ::Type{NRows}) = length(df)

retrieve(::Val{false}, A::Matrix, ::Type{Rows}, r) = A[r,:]
retrieve(::Val{false}, A::Matrix, ::Type{Cols}, c) = A[:,c]
function retrieve(::Val{false}, A::Matrix{T}, ::Type{Schema}) where T
    nrows = size(A, 1)
    ncols = size(A, 2)
    _names = [Symbol(:x, j) for j in 1:ncols]
    _eltypes = fill(T, ncols)
    return Schema(nrows, ncols, _names, _eltypes)
end

retrieve(::Val{false}, v::Vector, ::Type{Rows}, r) = v[r]
retrieve(::Val{false}, v::Vector, ::Type{Cols}, c) = error("Abstract vectors are not column-indexable.")
retrieve(::Val{false}, v::Vector{T}, ::Type{Schema}) where T = Schema(length(v), 1, [:x], [T,])
retrieve(::Val{false}, v::CategoricalArray{T,1,S} where {T,S}, ::Type{Rows}, r) = @inbounds v[r]
retrieve(::Val{false}, v::CategoricalArray{T,1,S} where {T,S}, ::Type{Cols}, r) =
    error("Categorical vectors are not column-indexable.")
retrieve(::Val{false}, v::CategoricalArray{T,1,S} where {T,S}, ::Type{Schema}) =
    retrieve(Val(false), collect(v), Schema)


text:
"""
    expected_CpG(seq::String)

Return the fraction of CpG dinucleotides expected in a random DNA sequence.
"""
function expected_CpG(seq::String)

    target = uppercase(seq)
    G = sum( [char for char in target] .== 'G' )
    C = sum( [char for char in target] .== 'C' )
    N = sum( [char for char in target] .== 'N' )
    Gexp = G / (length(seq) - N) #freq G
    Cexp = C / (length(seq) - N) #freq C
    exp_count = 2 * ((Cexp * Gexp) * (length(seq) - N))
    return exp_count

end

ext: jl


text:
using Flux
println("NNPDE_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
using GalacticOptim
using Optim

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

## Example 1, 1D ode
@parameters t θ
@variables u(..)
@derivatives Dt'~t

# 1D ODE
eq = Dt(u(t,θ)) ~ t^3 + 2*t + (t^2)*((1+3*(t^2))/(1+t+(t^3))) - u(t,θ)*(t + ((1+3*(t^2))/(1+t+t^3)))

# Initial and boundary conditions
bcs = [u(0.,θ) ~ 1.0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,1.0)]
# Discretization
dt = 0.1
# Neural network
chain = FastChain(FastDense(1,12,Flux.σ),FastDense(12,1))

discretization = NeuralPDE.PhysicsInformedNN(dt,
                                             chain,
                                             init_params = nothing;
                                             phi = nothing,
                                             autodiff=false,
                                             derivative = nothing,
                                             strategy = NeuralPDE.GridTraining())

pde_system = PDESystem(eq,bcs,domains,[t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=1500)
phi = discretization.phi

analytic_sol_func(t) = exp(-(t^2)/2)/(1+t+t^3) + t^2
ts = [domain.domain.lower:dt/10:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(phi(t,res.minimizer)) for t in ts]

@test u_predict ≈ u_real atol = 0.2

# using Plots
# t_plot = collect(ts)
# plot(t_plot ,u_real)
# plot!(t_plot ,u_predict)

## Example 2, 2D Poisson equation
@parameters x y θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# 2D PDE
eq  = Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

# Initial and boundary conditions
bcs = [u(0,y,θ) ~ 0.f0, u(1,y,θ) ~ -sin(pi*1)*sin(pi*y),
       u(x,0,θ) ~ 0.f0, u(x,1,θ) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.05

# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN(dx,
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac = 0.4))

pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=700)
phi = discretization.phi

xs,ys = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains]
analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)

@test u_predict ≈ u_real atol = 3.0

# p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
# p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)

## Example 3, 3D PDE
@parameters x y t θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y
@derivatives Dt'~t

# 3D PDE
eq  = Dt(u(x,y,t,θ)) ~ Dxx(u(x,y,t,θ)) + Dyy(u(x,y,t,θ))
# Initial and boundary conditions
bcs = [u(x,y,0,θ) ~ exp(x+y)*cos(x+y) ,
      # u(x,y,2,θ) ~ exp(x+y)*cos(x+y+4*2) ,
      u(0,y,t,θ) ~ exp(y)*cos(y+4t),
      u(2,y,t,θ) ~ exp(2+y)*cos(2+y+4t) ,
      u(x,0,t,θ) ~ exp(x)*cos(x+4t),
      u(x,2,t,θ) ~ exp(x+2)*cos(x+2+4t)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
          y ∈ IntervalDomain(0.0,2.0),
          t ∈ IntervalDomain(0.0,2.0)]

# Discretization
dx = 0.25; dy= 0.25; dt = 0.25
# Neural network
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dx,dy,dt],
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))
pde_system = PDESystem(eq,bcs,domains,[x,y,t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1), progress = false; cb = cb, maxiters=1000)
phi = discretization.phi

xs,ys,ts = [domain.domain.lower:dx:domain.domain.upper for (dx,domain) in zip([dx,dy,dt],domains)]
analytic_sol_func(x,y,t) = exp(x+y)*cos(x+y+4t)
u_real = [reshape([analytic_sol_func(x,y,t) for x in xs  for y in ys], (length(xs),length(ys)))  for t in ts ]
u_predict = [reshape([first(phi([x,y,t],res.minimizer)) for x in xs  for y in ys], (length(xs),length(ys)))  for t in ts ]
@test u_predict ≈ u_real atol = 200.0

# p1 =plot(xs, ys, u_predict, st=:surface);
# p2 = plot(xs, ys, u_real, st=:surface);
# plot(p1,p2)

## Example 4, high-order ode
@parameters x θ
@variables u(..)
@derivatives Dxxx'''~x
@derivatives Dx'~x

# ODE
eq = Dxxx(u(x,θ)) ~ cos(pi*x)

# Initial and boundary conditions
bcs = [u(0.,θ) ~ 0.0,
       u(1.,θ) ~ cos(pi),
       Dx(u(1.,θ)) ~ 1.0]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.05

# Neural network
chain = FastChain(FastDense(1,8,Flux.σ),FastDense(8,1))

discretization = NeuralPDE.PhysicsInformedNN(dx,chain,strategy= NeuralPDE.StochasticTraining(include_frac=0.5))
pde_system = PDESystem(eq,bcs,domains,[x],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob, ADAM(0.01), progress = false; cb = cb, maxiters=2000)
phi = discretization.phi

analytic_sol_func(x) = (π*x*(-x+(π^2)*(2*x-3)+1)-sin(π*x))/(π^3)

xs = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(x) for x in xs]
u_predict  = [first(phi(x,res.minimizer)) for x in xs]

@test u_predict ≈ u_real atol = 1.0

# x_plot = collect(xs)
# plot(x_plot ,u_real)
# plot!(x_plot ,u_predict)

## Example 5, system of pde
@parameters x, y, θ
@variables u1(..), u2(..)
@derivatives Dx'~x
@derivatives Dy'~y

# System of pde
eqs = [Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
      Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0,
      3*u1(x,0,θ) ~ 2*u2(x,0,θ)]

# Initial and boundary conditions
bcs = [[u1(x,0,θ) ~ 2x, u2(x,0,θ) ~ 3x]]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.1

# Neural network
chain = FastChain(FastDense(2,8,Flux.σ),FastDense(8,2))
discretization = NeuralPDE.PhysicsInformedNN(dx,chain)
pde_system = PDESystem(eqs,bcs,domains,[x,y],[u1,u2])
prob = NeuralPDE.discretize(pde_system,discretization)

res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=300)
phi = discretization.phi

analytic_sol_func(x,y) =[1/3*(6x - y), 1/2*(6x - y)]
xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
u_real  = [[analytic_sol_func(x,y)[i] for x in xs  for y in ys] for i in 1:2]
u_predict  = [[phi([x,y],res.minimizer)[i] for x in xs  for y in ys] for i in 1:2]

@test u_predict ≈ u_real atol = 10.0

# p1 =plot(xs, ys, u_predict, st=:surface);
# p2 = plot(xs, ys, u_real, st=:surface);
# plot(p1,p2)

## Example 6, 2d wave equation, neumann boundary condition
#here we use low level api for build solution
@parameters x, t, θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dtt''~t
@derivatives Dt'~t

#2D PDE
C=1
eq  = Dtt(u(x,t,θ)) ~ C^2*Dxx(u(x,t,θ))

# Initial and boundary conditions
bcs = [u(0,t,θ) ~ 0.,# for all t > 0
       u(1,t,θ) ~ 0.,# for all t > 0
       u(x,0,θ) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(x,0,θ)) ~ 0. ] #for all  0 < x < 1]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           t ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1

# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
strategy=NeuralPDE.GridTraining()
discretization = NeuralPDE.PhysicsInformedNN(dx,chain,strategy=strategy)

phi = discretization.phi
derivative = discretization.derivative
initθ = discretization.initθ

indvars = [x,t]
depvars = [u]
dim = length(domains)

expr_pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars)
expr_bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars) for bc in bcs]

train_sets = NeuralPDE.generate_training_sets(domains,dx,bcs,indvars,depvars)

train_domain_set,train_bound_set,train_set= train_sets

pde_loss_function = NeuralPDE.get_loss_function(eval(expr_pde_loss_function),
                                      train_domain_set,
                                      phi,
                                      derivative,
                                      strategy)
bc_loss_function = NeuralPDE.get_loss_function(eval.(expr_bc_loss_functions),
                                     train_bound_set,
                                     phi,
                                     derivative,
                                     strategy)

function loss_function(θ,p)
    return pde_loss_function(θ) + bc_loss_function(θ)
end

f = OptimizationFunction(loss_function, initθ, GalacticOptim.AutoZygote())
prob = GalacticOptim.OptimizationProblem(f, initθ)

res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=400)

xs,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
analytic_sol_func(x,t) =  sum([(8/(k^3*pi^3)) * sin(k*pi*x)*cos(C*k*pi*t) for k in 1:2:50000])

u_predict = reshape([first(phi([x,t],res.minimizer)) for x in xs for t in ts],(length(xs),length(ts)))
u_real = reshape([analytic_sol_func(x,t) for x in xs for t in ts], (length(xs),length(ts)))

@test u_predict ≈ u_real atol = 10.0

# diff_u = abs.(u_predict .- u_real)
# p1 = plot(xs, ts, u_real, linetype=:contourf,title = "analytic");
# p2 =plot(xs, ts, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(xs, ts, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)


text:
# -----------------------------------------------------------
#
# Constructing the matrix elements of the derivative
# operator in the 1-D DG Basis
#
# -----------------------------------------------------------

# Efficiency criticality: LOW
# Computations only performed once

# Accuracy criticality: HIGH
# Critical for accurate PDE evolution

# -----------------------------------------------------
# 1-D Derivatives
# -----------------------------------------------------

function dLegendreP(k,x)
	k<=K_max || throw(DomainError())
	return array2poly(symbolic_diff(leg_coeffs[k+1]),x)
end

function dh(k,mode,x)
	mode<=k || throw(DomainError())
	return array2poly(symbolic_diff((dg_coeffs[k])[mode]),x)
end

function dleg{T<:Real}(mode::Int, x::T)
	return sqrt(2.0)*dLegendreP(mode-1, 2*x-1) *2
end

function dbasis{T<:Real}(level::Int, cell::Int, mode::Int, x::T)
	return dleg(mode, (1<<level)*x - (cell-1)) * (2.0)^(level/2) * (1<<level)
end

function dbasis(level::Int, cell::Int, mode::Int)
	return x->dbasis(level, cell, mode, x)
end

# -----------------------------------------------------
# Shifted and scaled derivatives: v
# -----------------------------------------------------

function dv(k::Int, level::Int, cell::Int, mode::Int, x::Real)
	if level==0 # At the base level, we are not discontinuous, and we simply
				# use Legendre polynomials up to degree k-1 as a basis
		return 2*dLegendreP(mode-1,2*x-1)*sqrt(2.0)
	else
		return (1<<(level))*dh(k, mode, (1<<level)*x - (2*cell-1)) * (2.0)^(level/2)
		# Otherwise we use an appropriately shifted, scaled, and normalized
		# DG function
	end
end

function dv(k::Int, level::Int, cell::Int, mode::Int)
	return (xs::Real -> dv(k,level,cell,mode,xs))
end

# -----------------------------------------------------
# Numerical Inner Product
# -----------------------------------------------------

# TODO: There is a way to do this all symbolically, with no use for numerics
function inner_product1D(f::Function, g::Function, lvl::Int, cell::Int;
								rel_tol = REL_TOL, abs_tol = ABS_TOL, max_evals=MAX_EVALS)
	xmin = (cell-1)/(1<<max(0, lvl-1))
	xmax = (cell)/(1<<max(0, lvl-1))
	h = (x-> f(x)*g(x))
	(val, err) = hquadrature(h, xmin, xmax; reltol=rel_tol, abstol=abs_tol, maxevals=max_evals)
	return val
end

# -----------------------------------------------------
# Taking Derivative of Array Representing Polynomial
# -----------------------------------------------------

function symbolic_diff{T<:Real}(v::Array{T})
	n=length(v)
	k=div(n,2)
	ans = zeros(T,n)
	for i in 1:n
		if i<k
			ans[i] = i*v[i+1]
		elseif i > k && i<2k
			ans[i] = (i-k) * v[i+1]
		else
			ans[i]=0
		end
	end
	return ans
end

# -----------------------------------------------------
# < f_1 | D | f_2 > matrix elements
# -----------------------------------------------------

# Matrix element of legendre basis on [-1, 1]
function legendreDlegendre(mode1::Int, mode2::Int)
	return inner_product(leg_coeffs[mode1], symbolic_diff(leg_coeffs[mode2]))
end

# Matrix element of position basis elements
function legvDv(level, cell1, mode1, cell2, mode2;
				rel_tol = REL_TOL, abs_tol=ABS_TOL, max_evals=MAX_EVALS)
	if cell1 == cell2
		return (1<<(level+1))*legendreDlegendre(mode1, mode2)
	end
	return 0.0
end

# Matrix element of DG basis on [-1, 1]
function hDh(k::Int, mode1::Int, mode2::Int)
	return inner_product(dg_coeffs[k][mode1], symbolic_diff(dg_coeffs[k][mode2]))
end

# Matrix element of hierarchical basis elements
function vDv(k::Int, lvl1::Int, cell1::Int, mode1::Int, lvl2::Int, cell2::Int, mode2::Int;
					rel_tol = REL_TOL, abs_tol = ABS_TOL, max_evals=MAX_EVALS)
	if lvl1 == lvl2
		if lvl1 == 0
			return 2*legendreDlegendre(mode1, mode2)
		end
		if cell1 == cell2
			return hDh(k, mode1, mode2)*(1<<max(0, lvl1))
		end
		return 0.0
	end
	if lvl1 < lvl2
		if lvl1 == 0
			return inner_product1D(v(k,0,1,mode1), dv(k,lvl2,cell2,mode2), lvl2, cell2;
									rel_tol = rel_tol, abs_tol = abs_tol, max_evals=max_evals)
		end
		if (1<<(lvl2-lvl1))*cell1 >= cell2 && (1<<(lvl2-lvl1))*(cell1-1) < cell2
			return inner_product1D(v(k,lvl1,cell1,mode1), dv(k,lvl2,cell2,mode2), lvl2, cell2;
									rel_tol = rel_tol, abs_tol = abs_tol, max_evals=max_evals)
		end
		return 0.0
	end
	return 0.0

end

# -----------------------------------------------------
# < f_1 | D | f_2 > for a specific f_2 over all f_1
# -----------------------------------------------------

function diff_basis_DG(k::Int, level::Int, cell::Int, mode::Int)
	dcoeffs = Array{Real}((level+1, k))
	p = cell
	for l in level:-1:0
		for f_n in 1:k
			dcoeffs[l+1,f_n]=vDv(k, l, p, f_n, level, cell, mode)
		end
		p = Int(ceil(p/2))
	end
	return dcoeffs
end


text:
# test edge lists

using Graphs
using Test

## simple edge list

edge_pairs = [(1,2), (1,3), (2,3), (2,4), (3,5), (4,5), (2,5)]
eds = Edge{Int}[Edge(i,p[1],p[2]) for (i,p) in enumerate(edge_pairs)]

gd = simple_edgelist(5, eds)
gu = simple_edgelist(5, eds; is_directed=false)

@test implements_vertex_list(gd)    == true
@test implements_edge_list(gd)      == true
@test implements_vertex_map(gd)     == true
@test implements_edge_map(gd)       == true

@test implements_adjacency_list(gd) == false
@test implements_incidence_list(gd) == false
@test implements_bidirectional_adjacency_list(gd) == false
@test implements_bidirectional_incidence_list(gd) == false
@test implements_adjacency_matrix(gd) == false

# test for gd

@test vertex_type(gd) == Int
@test edge_type(gd) == Edge{Int}
@test is_directed(gd)

@test num_vertices(gd) == 5
@test vertices(gd) == 1:5

@test num_edges(gd) == 7
@test edges(gd) === eds

# test for gu

@test vertex_type(gu) == Int
@test edge_type(gu) == Edge{Int}
@test !is_directed(gu)

@test num_vertices(gu) == 5
@test vertices(gu) == 1:5

@test num_edges(gu) == 7
@test edges(gu) === eds

## edge list (based on vector of vertices)

for T in [ ExVertex, String ]
    g = edgelist(T[], ExEdge{T}[])

    vs = [ add_vertex!(g, "a"),
           add_vertex!(g, "b"),
           add_vertex!(g, "c") ]

    es = [ add_edge!(g, vs[1], vs[2]),
           add_edge!(g, vs[1], vs[3]) ]

    @test vertex_type(g) == T
    @test edge_type(g) == ExEdge{T}
    @test is_directed(g)

    @test num_vertices(g) == 3
    @test vertices(g) == vs

    @test num_edges(g) == 2
    @test edges(g) == es
end


text:
using QuadraticFormsMGHyp
import QuadraticFormsMGHyp: lklam
using FinancialToolbox
using LinearAlgebra

portfolio = 1 # between 1 and 16; these are the portfolios from Broda (2012)
xvec2 = 3.5:.01:17.5 # results in VaR levels up to 10% for portfolio 1 for NIG with chi=psi=1

lam = -.5
chi = 1
psi = 1
m = 10 # number of factors=number of options
S0 = ones(m) * 100.0 # initial prices
K = S0 # strikes
r = .05 # risk-free rate
sigma = ones(m) * .3 # annual return volatility
put=[false, false, false, false, false, true, true, true, true, true] # true for put, false for call
h = 1 / 252 # forecast horizon of x days = x/250 years

if portfolio < 9
    T = ones(m) * 126 / 252 # time to maturity in years
else
    T = ones(m) * 21/252 # time to maturity in years
end
if mod(portfolio, 2) == 1
    number = -ones(m) # number of options (negative for short positions)
else
    number = ones(m) # number of options (negative for short positions)
end
if portfolio == 3 || portfolio == 4 || portfolio == 7 || portfolio == 8 || portfolio == 11 || portfolio == 12 || portfolio == 15 || portfolio == 16
    deltahedged = true # include stocks to delta-hedge
else
    deltahedged = false
end
if (4 < portfolio && portfolio < 9) || portfolio > 12
    equicorrelated = true # equicorrelated or uncorrelated assets
else
    equicorrelated = false
end
if equicorrelated
    R = .5 * ones(m, m) + .5 * I # the correlation matrix
else
    R = I
end
# get initial option prices
price = blsprice.(S0, K, r, T, sigma, 0.0, put)
delta = blsdelta.(S0, K, r, T, sigma, 0.0, put)
gamma = blsgamma.(S0, K, r, T, sigma, 0.0, put)
theta = blstheta.(S0, K, r, T, sigma, 0.0, put)
if deltahedged
    stocks = -number .* delta
else
    stocks = zeros(m)
end
a = - delta .* number - stocks
A = diagm(0=>-.5 * gamma .* number)
portfolio_value = sum(price .* number) + sum(stocks .* S0)
sigs = diagm(0=>sigma .* S0 * sqrt(h)) # price volatilities 250*h days ahead
V = exp(lklam(lam+1, chi, psi) - lklam(lam, chi, psi)) # variance, assuming gamma=0
Sigma = sigs * R * sigs / V # the dispersion matrix
C = ((Sigma + Sigma') / 2) ^ .5 # and its square root
a0 = -h * sum(number.*theta)
ccdf, pm = qfmgh(xvec2, a0, a, A, C, zeros(m), zeros(m), lam, chi, psi; do_spa=false)
ccdfspa, pmspa = qfmgh(xvec2, a0, a, A, C, zeros(m), zeros(m), lam, chi, psi; do_spa=true, order=2)


text:
module Pluto
export Notebook, Cell, run

import Pkg

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
include_dependency(joinpath(PKG_ROOT_DIR, "Project.toml"))
const PLUTO_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(PKG_ROOT_DIR, "Project.toml"))["version"])
const PLUTO_VERSION_STR = 'v' * string(PLUTO_VERSION)
const JULIA_VERSION_STR = 'v' * string(VERSION)
const ENV_DEFAULTS = Dict(
    "PLUTO_WORKSPACE_USE_DISTRIBUTED" => "true",
    "PLUTO_RUN_NOTEBOOK_ON_LOAD" => "true",
    "PLUTO_ROOT_URL" => "/",
)

if get(ENV, "PLUTO_SHOW_BANNER", "true") == "true"
@info """\n
    Welcome to Pluto $(PLUTO_VERSION_STR) 🎈
    Start a notebook server using:

  julia> Pluto.run(1234)

    Have a look at the FAQ:
    https://github.com/fonsp/Pluto.jl/wiki
\n"""
end

include("./react/ExploreExpression.jl")
include("./runner/PlutoRunner.jl")
using .ExploreExpression

include("./notebookserver/Cell.jl")
include("./notebookserver/Notebook.jl")
include("./notebookserver/Client.jl")

include("./react/Parse.jl")
include("./react/WorkspaceManager.jl")
include("./react/Errors.jl")
include("./react/React.jl")
include("./react/Run.jl")

include("./webserver/WebServer.jl")
include("./webserver/Static.jl")
include("./webserver/Dynamic.jl")
include("./webserver/REPLTools.jl")

end


text:
using Cthulhu
using REPL
using InteractiveUtils
using Test
using Random
using StaticArrays
using Revise

@test isempty(detect_ambiguities(Cthulhu))

include("utils.jl")

function testf_simple()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

function empty_func(::Bool) end

isordered(::Type{T}) where {T<:AbstractDict} = false

@testset "Callsites" begin
    callsites = find_callsites_by_ftt(testf_simple)
    @test length(callsites) >= 4

    callsites = find_callsites_by_ftt(testf_simple; optimize=false)
    @test length(callsites) == 4

    callsites = find_callsites_by_ftt(empty_func, Tuple{Bool}; optimize=true)
    @test isempty(callsites)

    # handle pure
    callsites = find_callsites_by_ftt(iterate, Tuple{SVector{3,Int}, Tuple{SOneTo{3}}}; optimize=false)
    @test occursin("< pure >", string(callsites[1]))
    @test occursin(r"< constprop > getindex\(::.*Const.*,::.*Const\(1\)\)::.*Const\(1\)", string(callsites[2]))
    callsites = @eval find_callsites_by_ftt(; optimize=false) do
        length($(QuoteNode(Core.svec(0,1,2))))
    end
    @test occursin("< pure >", string(callsites[1]))

    # Some weird methods get inferred
    callsites = find_callsites_by_ftt(iterate, (Base.IdSet{Any}, Union{}); optimize=false)
    @test callsites[1].info isa Cthulhu.ConstPropCallInfo

    # Broken stuff in Julia
    @test_broken find_callsites_by_ftt(Core.Compiler._limit_type_size, Tuple{Any, Type{Any}, Core.SimpleVector, Int, Int})  # ssair/ir.jl bug
end

@testset "Expr heads" begin
    # Check that the Expr head (:invoke or :call) is preserved
    @noinline twice(x::Real) = 2x
    calltwice(c) = twice(c[1])

    callsites = find_callsites_by_ftt(calltwice, Tuple{Vector{Float64}})
    @test length(callsites) == 1 && callsites[1].head === :invoke
    io = IOBuffer()
    print(io, callsites[1])
    @test occursin("invoke twice(::Float64)::Float64", String(take!(io)))

    callsites = find_callsites_by_ftt(calltwice, Tuple{Vector{AbstractFloat}})
    @test length(callsites) == 1 && callsites[1].head === :call
    io = IOBuffer()
    print(io, callsites[1])
    @test occursin("call twice(::AbstractFloat)", String(take!(io)))

    # Note the failure of `callinfo` to properly handle specialization
    @test_broken Cthulhu.callinfo(Tuple{typeof(twice), AbstractFloat}, AbstractFloat) isa Cthulhu.MultiCallInfo
end

# https://github.com/JuliaDebug/Cthulhu.jl/issues/196
let
    callsites = @find_callsites_by_ftt tryparse(VersionNumber, "v2.1.3")
    @test !isempty(callsites)
    callsites = @find_callsites_by_ftt optimize=false tryparse(VersionNumber, "v2.1.3")
    @test !isempty(callsites)
end

# Check that we see callsites that are the rhs of assignments
@noinline bar_callsite_assign() = nothing
function foo_callsite_assign()
    x = bar_callsite_assign()
    x
end
let callsites = find_callsites_by_ftt(foo_callsite_assign, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

@eval function call_rt()
    S = $(Core.Compiler.return_type)(+, Tuple{Int, Int})
end
let callsites = find_callsites_by_ftt(call_rt, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

function call_by_apply(args...)
    identity(args...)
end
let callsites = find_callsites_by_ftt(call_by_apply, Tuple{Tuple{Int}}; optimize=false)
    @test length(callsites) == 1
end

# # TODO run this testset in a separate process
# # julia --check-bounds=auto --code-coverage=none
# @testset "DCE & boundscheck" begin
#     M = Module()
#     @eval M begin
#         Base.@propagate_inbounds function f(x)
#             @boundscheck error()
#         end
#         g(x) = @inbounds f(x)
#         h(x) = f(x)
#     end
#
#     let
#         (_,CI, _, _, _, _) = process(M.g, Tuple{Vector{Float64}})
#         @test all(CI.stmts.inst) do stmt
#             isa(stmt, Core.GotoNode) || (isa(stmt, Core.ReturnNode) && isdefined(stmt, :val))
#         end
#     end
#
#     let
#         (_,CI, _, _, _, _) = process(M.h, Tuple{Vector{Float64}})
#         @test count(!isnothing, CI.stmts.inst) == 2
#         stmt = CI.stmts.inst[end]
#         @test isa(stmt, Core.ReturnNode) && !isdefined(stmt, :val)
#     end
# end

# Something with many methods, but enough to be under the match limit
g_matches(a::Int, b::Int) = a+b
g_matches(a::Float64, b::Float64) = a+b
f_matches(a, b) = g_matches(a, b)
let callsites = find_callsites_by_ftt(f_matches, Tuple{Any, Any}; optimize=false)
    @test length(callsites) == 1
    callinfo = callsites[1].info
    @test callinfo isa Cthulhu.MultiCallInfo
    io = IOBuffer()
    Cthulhu.show_callinfo(io, callinfo)
    @test occursin(r"→ g_matches\(::Any, ?::Any\)::Union{Float64, ?Int\d+}", String(take!(io)))
end

@testset "wrapped callinfo" begin
    let
        m = Module()
        @eval m begin
            # mutually recursive functions
            f(a) = g(a)
            g(a) = somecode::Bool ? h(a) : a
            h(a) = f(Type{a})
        end

        # make sure we form `UncachedCallInfo` so that we won't try to retrieve non-existing cache
        callsites = @find_callsites_by_ftt m.f(Int)
        @test length(callsites) == 1
        ci = first(callsites).info
        @test isa(ci, Cthulhu.UncachedCallInfo)
        @test Cthulhu.is_callsite(ci, ci.wrapped.mi)
        io = IOBuffer()
        show(io, first(callsites))
        @test occursin("< uncached >", String(take!(io)))
        # TODO do some test with `LimitedCallInfo`, but they happen at deeper callsites
    end
end

@testset "union-split constant-prop'ed callsites" begin
    # constant prop' on all the splits
    let callsites = (@eval Module() begin
            struct F32
                val::Float32
                _v::Int
            end
            struct F64
                val::Float64
                _v::Int
            end

            $find_callsites_by_ftt((Union{F32,F64},); optimize = false) do f
                f.val
            end
        end)
        @test length(callsites) == 1
        callinfo = callsites[1].info
        @test isa(callinfo, Cthulhu.MultiCallInfo)
        callinfos = callinfo.callinfos
        @test length(callinfos) == 2
        @test all(ci->isa(ci, Cthulhu.ConstPropCallInfo), callinfos)
        io = IOBuffer()
        Cthulhu.show_callinfo(io, callinfos[1])
        @test startswith(String(take!(io)), "getproperty")
        io = IOBuffer()
        print(io, callsites[1])
        @test occursin("= call → getproperty", String(take!(io)))
    end

    # const and non-const splits
    let callsites = find_callsites_by_ftt((Bool,Vector{Any},); optimize=false) do cond, t
            if cond
                t = (1, nothing)
            end
            t[1]
        end
        @test length(callsites) == 1                                        # getindex(::Union{Vector{Any}, Const(tuple(1,nothing))}, ::Const(1))
        callinfo = callsites[1].info
        @test isa(callinfo, Cthulhu.MultiCallInfo)
        callinfos = callinfo.callinfos
        @test length(callinfos) == 2
        @test count(ci->isa(ci, Cthulhu.MICallInfo), callinfos) == 1        # getindex(::Vector{Any}, ::Const(1))
        @test count(ci->isa(ci, Cthulhu.ConstPropCallInfo), callinfos) == 1 # getindex(::Const(tuple(1,nothing)), ::Const(1))
    end

    let callsites = (@eval Module() begin
            struct F32
                val::Float32
                _v::Int
            end

            $find_callsites_by_ftt((F32,); optimize = false) do f
                f.val
            end
        end)
        io = IOBuffer()
        print(io, callsites[1])
        @test occursin("= < constprop > getproperty(", String(take!(io)))
    end
end

struct SingletonPureCallable{N} end

@testset "PureCallInfo" begin
    s = sprint(Cthulhu.show_callinfo, Cthulhu.PureCallInfo(Any[typeof(sin), Float64], Float64))
    @test s == "sin(::Float64)::Float64"

    s = sprint(Cthulhu.show_callinfo, Cthulhu.PureCallInfo(Any[SingletonPureCallable{1}, Float64], Float64))
    @test s == "SingletonPureCallable{1}()(::Float64)::Float64"
end

@testset "ReturnTypeCallInfo" begin
    only_ints(::Integer) = 1

    callsites = find_callsites_by_ftt(; optimize=false) do
            t1 = Base._return_type(only_ints, Tuple{Int})     # successful `return_type`
            t2 = Base._return_type(only_ints, Tuple{Float64}) # failed `return_type`
            t1, t2
        end
    @test length(callsites) == 2
    callinfo1 = callsites[1].info
    @test callinfo1 isa Cthulhu.ReturnTypeCallInfo
    @test callinfo1.vmi isa Cthulhu.MICallInfo
    io = IOBuffer()
    Cthulhu.show_callinfo(io, callinfo1)
    @test String(take!(io)) == "only_ints(::$Int)::$Int"
    io = IOBuffer()
    print(io, callsites[1])
    @test occursin("return_type < only_ints(::$Int)::$Int >", String(take!(io)))

    callinfo2 = callsites[2].info
    @test callinfo2 isa Cthulhu.ReturnTypeCallInfo
    @test callinfo2.vmi isa Cthulhu.FailedCallInfo
    io = IOBuffer()
    Cthulhu.show_callinfo(io, callinfo2)
    @test String(take!(io)) == "only_ints(::Float64)::Union{}"
    io = IOBuffer()
    print(io, callsites[2])
    @test occursin("return_type < only_ints(::Float64)::Union{} >", String(take!(io)))
end

# tasks
ftask() = @sync @async show(io, "Hello")
let callsites = find_callsites_by_ftt(ftask, Tuple{})
    task_callsites = filter(c->c.info isa Cthulhu.TaskCallInfo, callsites)
    @test !isempty(task_callsites)
    @test filter(c -> c.info.ci isa Cthulhu.FailedCallInfo, task_callsites) == []
    io = IOBuffer()
    show(io, first(task_callsites))
    @test occursin("= task < #", String(take!(io)))
end

@testset "invoke" begin
    m = Module()
    @eval m begin
        f(a::Integer) = :Integer
        f(a::Int) = :Int
    end

    io = IOBuffer()
    let
        callsites = @eval m begin
            $find_callsites_by_ftt(; optimize=false) do
                Base.@invoke f(0::Integer)
            end
        end
        @test any(callsites) do callsite
            info = callsite.info
            isa(info, Cthulhu.InvokeCallInfo) && print(io, callsite)
            isa(info, Cthulhu.InvokeCallInfo) && info.ci.rt === Core.Compiler.Const(:Integer)
        end
    end
    @test occursin("= invoke < f(::Int", String(take!(io)))

    let
        callsites = @eval m begin
            $find_callsites_by_ftt(; optimize=false) do
                Base.@invoke f(0::Int)
            end
        end
        @test any(callsites) do callsite
            info = callsite.info
            isa(info, Cthulhu.InvokeCallInfo) && info.ci.rt === Core.Compiler.Const(:Int)
        end
    end
end

##
# Union{f, g}
##
callf(f::F, x) where F = f(x)
let callsites = find_callsites_by_ftt(callf, Tuple{Union{typeof(sin), typeof(cos)}, Float64})
    @test !isempty(callsites)
    if length(callsites) == 1
        @test first(callsites).info isa Cthulhu.MultiCallInfo
        callinfos = first(callsites).info.callinfos
        @test !isempty(callinfos)
        mis = map(Cthulhu.get_mi, callinfos)
        @test any(mi->mi.def.name == :cos, mis)
        @test any(mi->mi.def.name == :sin, mis)
    else
        @test all(cs->cs.info isa Cthulhu.MICallInfo, callsites)
    end
end

function toggler(toggle)
    if toggle
        g = sin
    else
        g = cos
    end
    g(0.0)
end
let callsites = find_callsites_by_ftt(toggler, Tuple{Bool})
    @test !isempty(callsites)
    if length(callsites) == 1
        @test first(callsites).info isa Cthulhu.MultiCallInfo
        callinfos = first(callsites).info.callinfos
        @test !isempty(callinfos)
        mis = map(Cthulhu.get_mi, callinfos)
        @test any(mi->mi.def.name == :cos, mis)
        @test any(mi->mi.def.name == :sin, mis)
    else
        @test all(cs->cs.info isa Cthulhu.MICallInfo, callsites)
    end
end

@testset "Varargs and printing" begin
    @test Cthulhu.headstring(Float32) == "Float32"
    @test Cthulhu.headstring(Vararg{Float32}) == "Float32"
    @test Cthulhu.headstring(Union{Char,Float32}) == "Union{Char, Float32}"
    @test Cthulhu.headstring(Union{}) == "Union{}"
    @test Cthulhu.headstring(Vector{Int}) == "Array"
    @test Cthulhu.headstring(Vector{T} where T) == "Array"
    @test Cthulhu.headstring(TypeVar(:F)) == "F"

    function checklim(f, n, strchar::Union{AbstractString,AbstractChar})
        io = IOBuffer()
        iolim = Cthulhu.TextWidthLimiter(io, n)
        print(iolim, strchar)
        return f(String(take!(iolim)))
    end
    @test checklim(str -> str == "α", 3, 'α')
    @test checklim(str -> str == "", 1, 'a')
    @test checklim(str ->  occursin("…", str), 80, "abcd"^21)
    @test checklim(str -> !occursin("…", str), 80, "abcd"^19)

    function checklim(f, n, info::Cthulhu.CallInfo)
        io = IOBuffer()
        iolim = Cthulhu.TextWidthLimiter(io, n)
        Cthulhu.show_callinfo(iolim, info)
        return f(String(take!(iolim)))
    end
    function fsplat(::Type{Int}, a...)
        z = zero(Int)
        for v in a
            z += v
        end
        return z
    end
    gsplat1(T::Type, a...) = fsplat(T, a...)   # does not force specialization
    hsplat1(A) = gsplat1(eltype(A), A...)
    for (Atype, haslen) in ((Tuple{Tuple{Int,Int,Int}}, true),
                            (Tuple{Vector{Int}}, false),
                            (Tuple{SVector{3,Int}}, true))
        let callsites = find_callsites_by_ftt(hsplat1, Atype; optimize=false)
            @test !isempty(callsites)
            cs = callsites[end]
            @test cs isa Cthulhu.Callsite
            mi = cs.info.mi
            @test mi.specTypes.parameters[end] === (haslen ? Int : Vararg{Int})
            @test checklim(str -> occursin("...", str) != haslen, 80, cs.info)
        end
    end
    callsites = find_callsites_by_ftt(hsplat1, Tuple{NTuple{10,Int}}; optimize=false)
    cs = callsites[end]
    @test checklim(str -> occursin("gsplat1(…,…,…,…,…,…,…,…,…,…,…)::Int", str), 80, cs.info)
    callsites = find_callsites_by_ftt(hsplat1, Tuple{NTuple{50,Int}}; optimize=false)
    cs = callsites[end]
    @test checklim(str -> occursin("gsplat1(…)::Int", str), 80, cs.info)

    # foo(x::Vector{Vector{Vector{Vector{Char}}}}) = -1
    foo(x::Vector) = -1
    bar() = foo([[[['c']]]])
    callsites = find_callsites_by_ftt(bar, Tuple{}; optimize=false)
    cs = callsites[end]
    checklim(80, cs.info) do str
        @test !occursin("Array{…}", str)
        @test occursin("::Core.Const(-1)", str)
    end
    checklim(55, cs.info) do str
        @test !occursin("Array{…}", str)
        @test !occursin("::Core.Const(-1)", str)
    end
    checklim(40, cs.info) do str
        @test occursin("Array{…}", str)
        @test occursin("::Core.Const(-1)", str)
    end
    checklim(25, cs.info) do str
        @test occursin("Array{…}", str)
        @test !occursin("::Core.Const(-1)", str)
    end
    checklim(8, cs.info) do str
        @test str == "foo(…)"
    end
    checklim(4, cs.info) do str
        @test str == "…"
    end
end

@testset "MaybeUndef" begin
    function undef(b::Bool)
        b || @goto final_step
        str = randstring(8)
        @label final_step
        return str*"end"
    end
    @test isa(undef(true), String)
    @test_throws UndefVarError undef(false)
    cs = find_callsites_by_ftt(undef, Tuple{Bool})[end]
    @test cs.head === :invoke
    @test cs.info.mi.def == which(string, (String,String))
end

like_cat(dims, xs::AbstractArray{T}...) where T = like_cat_t(T, xs...; dims=dims)
like_cat_t(::Type{T}, xs...; dims) where T = T

let callsites = find_callsites_by_ftt(like_cat, Tuple{Val{3}, Vararg{Matrix{Float32}}})
    @test length(callsites) == 1
    cs = callsites[1]
    @test cs isa Cthulhu.Callsite
    mi = cs.info.mi
    @test mi.specTypes.parameters[4] === Type{Float32}
end

@testset "warntype variables" begin
    src, rettype = code_typed(identity, (Any,); optimize=false)[1]
    io = IOBuffer()
    ioctx = IOContext(io, :color=>true)
    Cthulhu.cthulhu_warntype(ioctx, :none, src, rettype, nothing)
    str = String(take!(io))
    @test occursin("x\e[91m\e[1m::Any\e[22m\e[39m", str)
end

@testset "Limit printing (issue #94)" begin
    m = Module()
    @eval m begin
        const x = collect(1:1000)
        f1() = x
        function f2()
            y = x
            return sum(y)
        end
    end
    function doprint(f)
        interp, mi = Cthulhu.mkinterp(f, ())
        src, rt = Cthulhu.lookup(interp, mi, true)
        io = IOBuffer()
        Cthulhu.cthulhu_typed(io, :none, src, rt, mi; iswarn=false)
        return String(take!(io))
    end
    @test occursin("invoke f1()::…\n", doprint(getfield(m, :f1)))
    str = doprint(getfield(m, :f2))
    @test occursin("x::Core.Const([1, 2, 3", str)
    @test !occursin("500,", str)
end

@testset "Issue #132" begin
    f132(w, dim) = [i == dim ? w[i]/2 : w[i]/1 for i in eachindex(w)]
    interp, mi = Cthulhu.mkinterp(f132, (Vector{Int}, Int))
    @test isa(mi, Core.MethodInstance)   # just check that the above succeeded
end

## Functions for "backedges & treelist"
# The printing changes when the functions are defined inside the testset
fbackedge1() = 1
fbackedge2(x) = x > 0 ? fbackedge1() : -fbackedge1()
fst1(x) = backtrace()
@inline fst2(x) = fst1(x)
@noinline fst3(x) = fst2(x)
@inline fst4(x) = fst3(x)
fst5(x) = fst4(x)

@testset "backedges and treelist" begin
    @test fbackedge2(0.2) == 1
    @test fbackedge2(-0.2) == -1
    mspec = @which(fbackedge1()).specializations
    mi = isa(mspec, Core.SimpleVector) ? mspec[1] : mspec.func
    root = Cthulhu.treelist(mi)
    @test Cthulhu.count_open_leaves(root) == 2
    @test root.data.callstr == "fbackedge1()"
    @test root.children[1].data.callstr == " fbackedge2(::Float64)"

    # issue #114
    unspecva(@nospecialize(i::Int...)) = 1
    @test unspecva(1, 2) == 1
    mi = something(first(methods(unspecva)).specializations...)
    root = Cthulhu.treelist(mi)
    @test occursin("Vararg", root.data.callstr)

    # Test highlighting and other printing
    mi = Cthulhu.get_specialization(:, Tuple{T, T} where T<:Integer)
    root = Cthulhu.treelist(mi)
    @test occursin("\e[31m::T\e[39m", root.data.callstr)
    mi = Cthulhu.get_specialization(Vector{Int}, Tuple{typeof(undef), Int})
    io = IOBuffer()
    @test Cthulhu.callstring(io, mi) == "Vector{$Int}(::UndefInitializer, ::$Int)"
    mi = Cthulhu.get_specialization(similar, Tuple{Type{Vector{T}}, Dims{1}} where T)
    @test occursin(r"31m::Type", Cthulhu.callstring(io, mi))

    # treelist for stacktraces
    tree = Cthulhu.treelist(fst5(1.0))
    @test match(r"fst1 at .*:\d+ => fst2 at .*:\d+ => fst3\(::Float64\) at .*:\d+", tree.data.callstr) !== nothing
    @test length(tree.children) == 1
    child = tree.children[1]
    @test match(r" fst4 at .*:\d+ => fst5\(::Float64\) at .*:\d+", child.data.callstr) !== nothing

    # issue #184
    tree = Cthulhu.treelist(similar(fst5(1.0), 0))
    @test isempty(tree.data.callstr)
    @test isempty(Cthulhu.callstring(io, similar(stacktrace(fst5(1.0)), 0)))
    @test Cthulhu.instance(similar(stacktrace(fst5(1.0)), 0)) === Core.Compiler.Timings.ROOTmi
end

@testset "ascend" begin
    # This tests only the non-interactive "look up the caller" portion
    callee(x) = 2x
    function caller(x)
        val = 0.0
        val += callee(3); line1 = @__LINE__
        val += callee(3.0); line2 = @__LINE__
        val += callee(x); line3 = @__LINE__
        val = sum([val])     # FIXME: without this line, `lookup` fails because codeinst.inferred === nothing
        return val, line1, line2, line3
    end
    _, line1, line2, line3 = caller(7)
    micaller = Cthulhu.get_specialization(caller, Tuple{Int})
    micallee_Int = Cthulhu.get_specialization(callee, Tuple{Int})
    micallee_Float4 = Cthulhu.get_specialization(callee, Tuple{Float64})
    info, lines = only(Cthulhu.find_caller_of(micallee_Int, micaller))
    @test info == (:caller, Symbol(@__FILE__), 0) && lines == [line1, line3]
    info, lines = only(Cthulhu.find_caller_of(micallee_Float4, micaller))
    @test info == (:caller, Symbol(@__FILE__), 0) && lines == [line2]

    # Detection in optimized (post-inlining) code
    @noinline nicallee(x) = 2x
    midcaller(x) = nicallee(x), @__LINE__
    function outercaller(x)
        val, line2 = midcaller(x); line1 = @__LINE__
        val = sum([val])     # FIXME: without this line, `lookup` fails because codeinst.inferred === nothing
        return val, line1, line2
    end
    _, line1, line2 = outercaller(7)
    micaller = Cthulhu.get_specialization(outercaller, Tuple{Int})
    micallee = Cthulhu.get_specialization(nicallee, Tuple{Int})
    callerinfo = Cthulhu.find_caller_of(micallee, micaller)
    @test length(callerinfo) == 2
    info, lines = callerinfo[1]
    @test info == (:outercaller, Symbol(@__FILE__), 0)
    @test lines == [line1]
    info, lines = callerinfo[2]
    @test info == (:midcaller, Symbol(@__FILE__), 1)
    @test lines == [line2]
end

@testset "ascend interface" begin
    m = Module()
    @eval m begin
        using Cthulhu
        struct FunnyMI end
        struct HasName
            name::Symbol
        end
        Cthulhu.method(::FunnyMI) = HasName(:funny)
        funny(c::Char) = "haha"
        Cthulhu.specTypes(::FunnyMI) = Tuple{typeof(funny),Char}
    end

    io = IOBuffer()
    @test Cthulhu.callstring(io, m.FunnyMI()) == "funny(::Char)"
end

##
# Cthulhu config test
##
let config = Cthulhu.CthulhuConfig(enable_highlighter=false)
    for lexer in ["llvm", "asm"]
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT", lexer, config)
        end == "INPUT\n"
    end
end

let config = Cthulhu.CthulhuConfig(
    highlighter = `I_am_hoping_this_command_does_not_exist`,
    enable_highlighter = true,
)
    for lexer in ["julia"]
        @test begin
            @test_logs (:warn, r"Highlighter command .* does not exist.") begin
                sprint() do io
                    Cthulhu.highlight(io, "INPUT", lexer, config)
                end == "INPUT\n"
            end
        end
    end
end

let config = Cthulhu.CthulhuConfig(
    # Implementing `cat` in Julia:
    highlighter = `$(Base.julia_cmd()) -e "write(stdout, read(stdin))" --`,
    enable_highlighter = true,
)
    for lexer in ["llvm", "asm"]
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT", lexer, config)
        end == "INPUT\n"
        @test sprint() do io
            Cthulhu.highlight(io, "INPUT\n", lexer, config)
        end == "INPUT\n"
    end
end

# https://github.com/JuliaDebug/Cthulhu.jl/issues/152
let callsites = @find_callsites_by_ftt optimize=false replace("s", "a"=>"b")
    @test !isempty(callsites)
end
function issue152_another(t)
    s = 0
    for a in t
        s += a
    end
    return s
end
let callsites = find_callsites_by_ftt(issue152_another, (Tuple{Float64,Vararg{Float64}},); optimize=false)
    @test !isempty(callsites)
end

include("codeview.jl")

@testset "Bookmarks" begin
    (interp, mi) = Cthulhu.mkinterp(sqrt, Tuple{Float64})
    b = Cthulhu.Bookmark(mi, interp)

    @testset "code_typed(bookmark)" begin
        ci, rt = code_typed(b)
        @test ci isa Core.Compiler.CodeInfo
        @test rt isa Type
    end

    @testset "code_typed(bookmark; optimize = false)" begin
        ci, rt = code_typed(b; optimize = false)
        @test ci isa Core.Compiler.CodeInfo
        @test rt isa Type
    end

    @testset "show(io, bookmark)" begin
        str = sprint(io -> show(io, "text/plain", b))
        @test occursin("Cthulhu.Bookmark (world: ", str)
    end

    @testset "show(io, [bookmark])" begin
        # Test that it does not print the full IR:
        str = sprint(io -> show(io, "text/plain", [b]))
        @test occursin("world:", str)
        @test !occursin("Cthulhu.Bookmark (world: ", str)
    end

    @testset "Smoke tests" begin
        @test code_warntype(devnull, b) isa Any
        @test code_llvm(devnull, b) isa Any
        @test code_native(devnull, b) isa Any
    end
end

include("terminal.jl")

using Cthulhu: MultiCallInfo, show_callinfo, CallInfo

@testset "printing of MultiCallInfo" begin
    # PR 226
    m = Module()
    @eval m begin
        struct Foo
            x
        end
        struct Bar end
    end
    ci = MultiCallInfo(Tuple{m.Foo, Float64}, Float64, CallInfo[])
    @test sprint(show_callinfo, ci) == "→ (::Foo)(::Float64)::Float64"
    ci = MultiCallInfo(Tuple{Union{m.Foo, m.Bar}, Float64}, Float64, CallInfo[])
    @test sprint(show_callinfo, ci) == "→ (::Union{Bar, Foo})(::Float64)::Float64"
end


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

"""
    LibGit2.tag_list(repo::GitRepo) -> Vector{String}

Get a list of all tags in the git repository `repo`.
"""
function tag_list(repo::GitRepo)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_tag_list, :libgit2), Cint,
                 (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, repo.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    LibGit2.tag_delete(repo::GitRepo, tag::AbstractString)

Remove the git tag `tag` from the repository `repo`.
"""
function tag_delete(repo::GitRepo, tag::AbstractString)
    @check ccall((:git_tag_delete, :libgit2), Cint,
                  (Ptr{Void}, Cstring), repo.ptr, tag)
end

"""
    LibGit2.tag_create(repo::GitRepo, tag::AbstractString, commit; kwargs...)

Create a new git tag `tag` (e.g. `"v0.5"`) in the repository `repo`, at
the commit `commit`.

The keyword arguments are:
  * `msg::AbstractString=""`: the message for the tag.
  * `force::Bool=false`: if `true`, existing references will be overwritten.
  * `sig::Signature=Signature(repo)`: the tagger's signature.
"""
function tag_create(repo::GitRepo, tag::AbstractString, commit::Union{AbstractString,AbstractGitHash};
                    msg::AbstractString = "",
                    force::Bool = false,
                    sig::Signature = Signature(repo))
    oid_ptr = Ref(GitHash())
    with(GitCommit(repo, commit)) do commit_obj
        commit_obj === nothing && return oid_ptr[] # return empty oid
        with(convert(GitSignature, sig)) do git_sig
            @check ccall((:git_tag_create, :libgit2), Cint,
                 (Ptr{GitHash}, Ptr{Void}, Cstring, Ptr{Void}, Ptr{SignatureStruct}, Cstring, Cint),
                  oid_ptr, repo.ptr, tag, commit_obj.ptr, git_sig.ptr, msg, Cint(force))
        end
    end
    return oid_ptr[]
end

"""
    LibGit2.name(tag::GitTag)

The name of `tag` (e.g. `"v0.5"`).
"""
function name(tag::GitTag)
    Base.@gc_preserve tag begin
        str_ptr = ccall((:git_tag_name, :libgit2), Cstring, (Ptr{Void},), tag.ptr)
        str_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
        str = unsafe_string(str_ptr)
    end
    return str
end

# should we return the actual object? i.e. git_tag_target?
"""
    LibGit2.target(tag::GitTag)

The `GitHash` of the target object of `tag`.
"""
function target(tag::GitTag)
    Base.@gc_preserve tag begin
        oid_ptr = ccall((:git_tag_target_id, :libgit2), Ptr{GitHash}, (Ptr{Void},), tag.ptr)
        oid_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
        str = unsafe_load(oid_ptr)
    end
    return str
end

Base.show(io::IO, tag::GitTag) = print(io, "GitTag:\nTag name: $(name(tag)) target: $(target(tag))")

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
module SerializingChannels

using TypeSerializers

export SerializingChannel, put!, take!, isready, fetch, close, isopen, bind, iterate

struct SerializingChannel{T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    channel::TChannel
end

function Base.put!(channel::SerializingChannel{T, TChannel, TSerializer}, value::T) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    stream = serialize(TSerializer, value)
    data = take!(stream)
    put!(channel.channel, data)
end

function Base.take!(channel::SerializingChannel{T, TChannel, TSerializer})::T where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    message = take!(channel.channel)
    stream = IOBuffer(message)
    return deserialize(TSerializer, stream)
end

function Base.isready(channel::SerializingChannel{T, TChannel, TSerializer}) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    return isready(channel.channel)
end

function Base.fetch(channel::SerializingChannel{T, TChannel, TSerializer}) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    message = fetch(channel.channel)
    stream = IOBuffer(message)
    return deserialize(TSerializer, stream)
end

function Base.close(channel::SerializingChannel{T, TChannel, TSerializer}) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    return close(channel.channel)
end

function Base.isopen(channel::SerializingChannel{T, TChannel, TSerializer}) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    return isopen(channel.channel)
end

function Base.bind(channel::SerializingChannel{T, TChannel, TSerializer}, task::Task) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    return bind(channel.channel, task)
end

function Base.iterate(channel::SerializingChannel{T, TChannel, TSerializer}, state::Any=nothing) where {T, TChannel <: AbstractChannel{Vector{UInt8}}, TSerializer <: AbstractTypeSerializer{T}}
    if (isopen(channel))
        return (take!(channel), nothing)
    end
    return nothing
end

end # module


text:
for filename in (
        "compressedio.jl",
        "imgcondition.jl",
        "spacetime.jl",
        "wavelet.jl")
    include(filename)
end


text:
"""
    expectation_step(X, z, dists)

Computes the expectation step with:
- `X`: (m×n) multi-variate observations (rows are observations, columns are variables)
- `dists`: distributions, one for each category (unique label)
- `prior`: prior probability of each label
Returns `Z`: (m×n) matrix of probability distributions of each label for each observation
"""
function expectation_step(X, dists, prior)
    n = size(X, 1)
    Z = Matrix{Float64}(undef, n, length(prior))
    for k in eachindex(prior)
        for i in 1:n
            Z[i,k] = prior[k] * pdf(dists[k], X[i,:]) /
                     sum(
                        prior[j] * pdf(dists[j], X[i,:])
                        for j in eachindex(prior)
                     )
        end
    end
    return Z
end

"""
    maximization_step(::Type{<:D}, X, Z)

Maximization step, parameterized by the type of distribution used.
Returns a tuple (dists, π, loglike) with:
- `dists`: the created vector of distributions
- `prior` vector of updated prior probabilities of each label
"""
function maximization_step(::Type{<:MvNormal}, X, Z)
    n = size(X, 1)
    Nk = size(Z, 2)
    µ = map(1:Nk) do k
        num = sum(Z[i,k] .* X[i,:] for i in 1:n)
        den = sum(Z[i,k] for i in 1:n)
        num / den
    end

    Σ = map(1:Nk) do k
        num = zeros(size(X, 2), size(X, 2))
        for i in 1:n
            r = X[i,:] .- µ[k]
            num .= num .+ Z[i,k] .* (r * r')
        end
        den = sum(Z[i,k] for i in 1:n)
        num ./ den
    end
    prior = [inv(n) * sum(Z[i,k] for i in 1:n)
            for k in 1:Nk
        ]
    dists = map(1:Nk) do k
        MvNormal(µ[k], Σ[k] + 10e-7I)
    end
    return (dists, prior)
end

"""
    expectation_maximization(D::Type{<:Distribution}, X, Nk, maxiter::Integer = 500, loglike_diff = 10e-5; all_dists = nothing)

Complete EM algorithm, takes a type of distribution `D` used to
determine the maximization step, observations `X` and number of classes `Nk`.
"""
function expectation_maximization(D::Type{<:Distribution}, X, Nk, maxiter::Integer = 500, loglike_diff = 10e-5; all_dists = nothing)
    # initialize classes
    n = size(X,1)
    Z = zeros(n, Nk)
    for i in 1:n
        j0 = mod(i,Nk)+1
        j1 = j0 > 1 ? j0-1 : 2
        Z[i,j0] = 0.75
        Z[i,j1] = 0.25
    end
    (dists, prior) = maximization_step(D, X, Z)
    all_dists isa Vector && push!(all_dists, (dists, prior))
    l = loglike_mixture(X, dists, prior)
    lprev = 0.0
    iter = 0
    # EM iterations
    while iter < maxiter && abs(lprev-l) > loglike_diff
        Z = expectation_step(X, dists, prior)
        (dists, prior) = maximization_step(D, X, Z)
        all_dists isa Vector && push!(all_dists, (dists, prior))
        lprev = l
        l = loglike_mixture(X, dists, prior)
        iter += 1
    end
    return (dists, prior, Z, l, iter)
end

function loglike_mixture(X, dists, prior)
    l = zero(eltype(X))
    n = size(X,1)
    for i in 1:n
        l += log(
            sum(prior[k] * pdf(dists[k], X[i,:]) for k in eachindex(prior))
        )
    end
    return l
end

let
    Random.seed!(42)
    X = [randn(1000, 2); 3.0 .+ randn(1000, 2)]
    n = size(X, 1)
    for maxiter in (1, 5, 10, 100, 150, 250)
        (dists, prior, Z, l, iter) = expectation_maximization(MvNormal, X, 2, maxiter)
        contour_mixture(x1, x2) = maximum([pdf(dists[k], [x1,x2]) for k in 1:2])
        colors = map(1:n) do i
            (_, j) = findmax(Z[i,:])
            j == 1 ? :red : :blue
        end
        plt = Plots.contour(-2.8:0.05:6.0, -2.8:0.05:6.0, contour_mixture, contours = false, legend = nothing)
        Plots.scatter!(plt, X[:,1],X[:,2], label = "", color = colors)
        Plots.savefig(joinpath(base_img, "mixture_iter$(maxiter)_$iter.pdf"))
    end
end


text:
using FinancialMonteCarlo, Test
@show "Test Parameters"
S0 = 100.0;
K = 100.0;
r = 0.02;
Tneg = -1.0;
T = -Tneg;
d = 0.01;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
McConfig = MonteCarloConfiguration(Nsim, Nstep);
toll = 0.8;
rfCurve2 = ZeroRate([0.00, 0.02], T);
rfCurve = ZeroRate(r);

@show "Test Geometric Brownian Motion Parameters"
drift = 0.0
@test_throws(ErrorException, GeometricBrownianMotion(-sigma, drift, 1.0))
@test_throws(ErrorException, GeometricBrownianMotion(-sigma, rfCurve2.r, 1.0))


text:
module Branch_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

S = BnBSolver()
B = BnBModel([Interval(1.0,2.0)])
tL = 3.0
tU = 5.0
X1 = [Interval(1.0,2.0)]
X2 = [Interval(3.0,5.0)]
X = [Interval(3.0,5.0)]
pos = 4
EAGOBranchBound.BM_breadth!(S,B,tL,tU,X1,X2,pos)
@test (B.pos[1:3] == [5; 5; 1])
@test (B.LBD[1:3] == [3.0; 3.0; -Inf])
@test (B.UBD[1:3] == [5.0; 5.0; Inf])
@test (B.id[1:3] == [2; 3; 1])

EAGOBranchBound.BM_Single!(S,B,tL,tU,X,pos)
@test (B.pos[4] == 4)
@test (B.LBD[4] == 3.0)
@test (B.UBD[4] == 5.0)
@test (B.id[4] == 3)

end


text:
###
### logistic and logit
###

logistic(x::Real) = inv(one(x) + exp(-x))

function logistic_logjac(x::Real)
    mx = -abs(x)
    mx - 2*log1p(exp(mx))
end

logit(x::Real) = log(x / (one(x) - x))

logit_logjac(y) = -log(y) - log1p(-y)

###
### calculations
###

"""
    $SIGNATURES

Number of elements (strictly) above the diagonal in an ``n×n`` matrix.
"""
unit_triangular_dimension(n::Int) = n * (n-1) ÷ 2

###
### type calculations
###

"""
    $(SIGNATURES)

Extend element type of argument so that it is closed under the algebra used by this package.

Pessimistic default for non-real types.
"""
function extended_eltype(::Type{S}) where S
    T = eltype(S)
    T <: Real ? typeof(√(one(T))) : Any
end

extended_eltype(x::T) where T = extended_eltype(T)

###
### macros
###

"""
$(SIGNATURES)

Workaround for https://github.com/JuliaLang/julia/issues/14919 to make
transformation types callable.

TODO: remove when this issue is closed, also possibly remove MacroTools as a
dependency if not used elsewhere.
"""
macro calltrans(ex)
    if @capture(ex, struct T1_ fields__ end)
        @capture T1 (T2_ <: S_|T2_)
        @capture T2 (T3_{params__}|T3_)
        quote
            Base.@__doc__ $(esc(ex))
            (t::$(esc(T3)))(x) = transform(t, x)
        end
    else
        throw(ArgumentError("can't find anything to make callable in $(ex)"))
    end
end

####
#### random values
####

"Shared part of docstrings for keyword arguments of or passed to [`random_reals`](@ref)."
const _RANDOM_REALS_KWARGS_DOC = """
A standard multivaritate normal or Cauchy is used, depending on `cauchy`, then scaled with
`scale`. `rng` is the random number generator used.
"""

_random_reals_scale(rng::AbstractRNG, scale::Real, cauchy::Bool) =
    cauchy ? scale / abs2(randn(rng)) : scale * 1.0

"""
$(SIGNATURES)

Random real number.

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_real(; scale::Real = 1, cauchy::Bool = false, rng::AbstractRNG = GLOBAL_RNG) =
    randn(rng) * _random_reals_scale(rng, scale, cauchy)

"""
$(SIGNATURES)

Random vector in ``ℝⁿ``.

$(_RANDOM_REALS_KWARGS_DOC)
"""
random_reals(n::Integer; scale::Real = 1, cauchy::Bool = false,
             rng::AbstractRNG = GLOBAL_RNG) =
                 randn(rng, n) .* _random_reals_scale(rng, scale, cauchy)


text:
module FastMarching

using LinearAlgebra

include("libmsfm.jl")
include("pointmin.jl")
include("rk4.jl")
include("shortestpath.jl")

end # module


text:
# 2D test distributions for generative models with the change of variables formula
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export sample_banana, sample_swirl

# Banana distribution
"""
    S = sample_banana(batchsize; c=[1f0, 4f0])

Generate samples from a 2D banana distribution. Samples S have
dimensions `1 x 1 x 2 x batchsize`.
"""
function sample_banana(batchsize; c=[1f0, 4f0])
    x = randn(Float32, 2, batchsize)
    y = zeros(Float32, 1, 1, 2, batchsize)
    y[1,1,1,:] = x[1,:] ./ c[1]
    y[1,1,2,:] = x[2,:].*c[1] + c[1].*c[2].*(x[1,:].^2 .+ c[1]^2)
    return y
end

# Swirl distribution
"""
    S = sample_swirl(batchsize; noise=1f0)

Generate samples from a 2D swirl distribution. Samples S have
dimensions `1 x 1 x 2 x batchsize`.
"""
function sample_swirl(batchsize; noise=1f0)
    n = sqrt.(rand(Float32, batchsize)) .* 800f0 .* (2f0*pi) ./ 360f0
    d1x = -cos.(n).*n .+ rand(Float32, batchsize) .* noise
    d1y = sin.(n).*n .+ rand(Float32, batchsize) .* noise
    return reshape(cat(d1x', d1y'; dims=1), 1, 1, 2, batchsize)
end


text:
abstract type QuadratureMethod end
struct DefaultQuadrature <: QuadratureMethod end
struct Analytic <: QuadratureMethod end

struct GaussHermite <: QuadratureMethod
    n_points::Int
end
GaussHermite() = GaussHermite(20)

struct MonteCarlo <: QuadratureMethod
    n_samples::Int
end
MonteCarlo() = MonteCarlo(20)

_default_quadrature(_) = GaussHermite()

"""
    expected_loglik(quadrature::QuadratureMethod, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik)

This function computes the expected log likelihood:

```math
    ∫ q(f) log p(y | f) df
```
where `p(y | f)` is the process likelihood. This is described by `lik`, which should be a callable that takes `f` as input and returns a Distribution over `y` that supports `loglikelihood(lik(f), y)`.

`q(f)` is an approximation to the latent function values `f` given by:
```math
    q(f) = ∫ p(f | u) q(u) du
```
where `q(u)` is the variational distribution over inducing points (see
[`elbo`](@ref)). The marginal distributions of `q(f)` are given by `q_f`.

`quadrature` determines which method is used to calculate the expected log
likelihood - see [`elbo`](@ref) for more details.

# Extended help

`q(f)` is assumed to be an `MvNormal` distribution and `p(y | f)` is assumed to
have independent marginals such that only the marginals of `q(f)` are required.
"""
expected_loglik(quadrature, y, q_f, lik)

"""
    expected_loglik(::DefaultQuadrature, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik)

The expected log likelihood.
Defaults to a closed form solution if it exists, otherwise defaults to
Gauss-Hermite quadrature.
"""
function expected_loglik(
    ::DefaultQuadrature, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik
)
    quadrature = _default_quadrature(lik)
    return expected_loglik(quadrature, y, q_f, lik)
end

function expected_loglik(
    mc::MonteCarlo, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik
)
    # take `n_samples` reparameterised samples
    f_μ = mean.(q_f)
    fs = f_μ .+ std.(q_f) .* randn(eltype(f_μ), length(q_f), mc.n_samples)
    lls = loglikelihood.(lik.(fs), y)
    return sum(lls) / mc.n_samples
end

function expected_loglik(
    gh::GaussHermite, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik
)
    # Compute the expectation via Gauss-Hermite quadrature
    # using a reparameterisation by change of variable
    # (see e.g. en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature)
    xs, ws = gausshermite(gh.n_points)
    # size(fs): (length(y), n_points)
    fs = √2 * std.(q_f) .* xs' .+ mean.(q_f)
    lls = loglikelihood.(lik.(fs), y)
    return sum((1 / √π) * lls * ws)
end

ChainRulesCore.@non_differentiable gausshermite(n)

function expected_loglik(::Analytic, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik)
    return error(
        "No analytic solution exists for ",
        typeof(lik),
        ". Use `DefaultQuadrature()`, `GaussHermite()` or `MonteCarlo()` instead.",
    )
end

# The closed form solution for independent Gaussian noise
function expected_loglik(
    ::Analytic,
    y::AbstractVector{<:Real},
    q_f::AbstractVector{<:Normal},
    lik::GaussianLikelihood,
)
    return sum(
        -0.5 * (log(2π) .+ log.(lik.σ²) .+ ((y .- mean.(q_f)) .^ 2 .+ var.(q_f)) / lik.σ²)
    )
end

_default_quadrature(::GaussianLikelihood) = Analytic()

# The closed form solution for a Poisson likelihood with an exponential inverse link function
function expected_loglik(
    ::Analytic,
    y::AbstractVector{<:Real},
    q_f::AbstractVector{<:Normal},
    ::PoissonLikelihood{ExpLink},
)
    f_μ = mean.(q_f)
    return sum((y .* f_μ) - exp.(f_μ .+ (var.(q_f) / 2)) - loggamma.(y .+ 1))
end

_default_quadrature(::PoissonLikelihood{ExpLink}) = Analytic()

# The closed form solution for an Exponential likelihood with an exponential inverse link function
function expected_loglik(
    ::Analytic,
    y::AbstractVector{<:Real},
    q_f::AbstractVector{<:Normal},
    ::ExponentialLikelihood{ExpLink},
)
    f_μ = mean.(q_f)
    return sum(-f_μ - y .* exp.((var.(q_f) / 2) .- f_μ))
end

_default_quadrature(::ExponentialLikelihood{ExpLink}) = Analytic()

# The closed form solution for a Gamma likelihood with an exponential inverse link function
function expected_loglik(
    ::Analytic,
    y::AbstractVector{<:Real},
    q_f::AbstractVector{<:Normal},
    lik::GammaLikelihood{<:Any,ExpLink},
)
    f_μ = mean.(q_f)
    return sum(
        (lik.α - 1) * log.(y) .- y .* exp.((var.(q_f) / 2) .- f_μ) .- lik.α * f_μ .-
        loggamma(lik.α),
    )
end

_default_quadrature(::GammaLikelihood{<:Any,ExpLink}) = Analytic()


text:
using Test;
using aoc_2019_01;

@testset "Should correctly calculate part one per-module cost" begin
    calculate = aoc_2019_01.calculate_fuel_cost

    example_one   = calculate(12)
    @test example_one   == 2

    example_two   = calculate(14)
    @test example_two   == 2

    example_three = calculate(1969)
    @test example_three == 654

    example_four  = calculate(100756)
    @test example_four  == 33583
end

@testset "Should correctly solve part one" begin
    solve  = aoc_2019_01.solve_part_one
    answer = solve("../src/input.txt")
    @test answer == 3_393_938
end

@testset "Should correctly calculate part two per-module cost" begin
    calculate = aoc_2019_01.calculate_fuel_cost_recursive

    example_one   = calculate(14)
    @test example_one   == 2

    example_two   = calculate(1969)
    @test example_two   == 966

    example_three = calculate(100_756)
    @test example_three == 50_346
end

@testset "Should correctly solve part two" begin
    solve  = aoc_2019_01.solve_part_two
    answer = solve("../src/input.txt")
    @test answer == 5_088_037
end

max_stars_repo_name: BsideAlgorithms/CodingDojo

max_issues_repo_name: BsideAlgorithms/CodingDojo

max_forks_repo_name: BsideAlgorithms/CodingDojo


text:
function readMat(data::Vector{UInt8})
    n = Int(length(data)/3073)
    x = Array(Float64, 3072, n)
    y = Array(Int, n)
    for i = 1:n
        k = (i-1) * 3073 + 1
        y[i] = Int(data[k])
        x[:,i] = data[k+1:k+3072] / 255
    end
    x = reshape(x, 32, 32, 3, n)
    x, y
end
function readLine(dir, files)
  lines = open(readlines, joinpath(dir,files[1]))
  x=map(l -> Vector{String}(split(chomp(l))), lines)
  y=copy(x) #Needs improvement
  x, y
end


text:
# test file to test testing

# Test @test
@test true
@test 1 == 1
@test 1 != 2
@test strip("\t  hi   \n") == "hi"
@test strip("\t  this should fail   \n") != "hi"

a = Array(Float64, 2, 2, 2, 2, 2)
a[1,1,1,1,1] = 10
@test a[1,1,1,1,1] == 10
@test a[1,1,1,1,1] != 2

@test rand() != rand()

# Test with_handler
successflag = false
failureflag = false
errorflag = false
test_handler(r::Test.Success) = !successflag
test_handler(r::Test.Failure) = !failureflag
test_handler(r::Test.Error) = !errorflag

Test.with_handler(test_handler) do
    @test true
    @test successflag
    @test !failureflag
    @test !errorflag
    successflag = false
    @test false
    @test !successflag
    @test failureflag
    @test !errorflag
    failureflag = false
    @test error("throw error")
    @test !successflag
    @test !failureflag
    @test errorflag
end

# Test @test_throws
domainerror_thrower() = throw(DomainError())
boundserror_thrower() = throw(BoundsError())
error_thrower() = error("An error happened")
@test_throws DomainError domainerror_thrower()
@test_throws BoundsError boundserror_thrower()

failureflag = false
successflag = false
Test.with_handler(test_handler) do
    @test_throws DomainError boundserror_thrower()
    @test failureflag
    @test_throws DomainError domainerror_thrower()
    @test successflag
end

# Test @test_approx_eq
# TODO
@test isapprox(.1+.1+.1, .3)
@test !isapprox(.1+.1+.1, .4)

# Test @test_approx_eq_eps
# TODO

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
@testset "Functions" begin
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    @testset "Copy" begin
        @testset "VectorOfVariables" begin
            f = MOI.VectorOfVariables([x, y])
            f_copy = copy(f)
            f_copy.variables[2] = z
            @test f.variables[2] == y
        end
    end
end


text:
const _DEPRECATED_ERROR_MESSAGE = """
The C API of CPLEX.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl/47814

Here is a brief summary of the changes.

* Function names have changed. For example, `CPLEX.close_CPLEX(env)` is now
    `CPXcloseCPLEX(env)`.
* For users of `CPLEX.Optimizer`, `model.inner` has been replaced by the fields
    `model.env` and `model.lp`, which correspond to the environment and problem
    pointers at the C API level. For example:
    ```julia
    stat = CPLEX.get_status(model.inner)
    ```
    is now:
    ```julia
    stat = CPXgetstat(model.env, model.lp)
    ```
* Querying functionality has changed. For example:
    ```julia
    is_point = CPLEX.cbcandidateispoint(cb_data)
    ```
    is now:
    ```julia
    is_point_P = Ref{Cint}()
    CPXcallbackcandidateispoint(cb_data, is_point_P)
    if ret != 0
        # Do something because the call failed
    end
    is_point = is_point_P[]
    ```

The new API is more verbose, but the names and function arguments are now
identical to the C API, documentation for which is available at:
https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/refcallablelibrary/groups/homepagecallable.html

To revert to the old API, use:

    import Pkg
    Pkg.add(Pkg.PackageSpec(name = "CPLEX", version = v"0.6"))

Then restart Julia for the change to take effect.
"""

add_constr!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_constrs!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_constrs_t!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_diag_qpterms!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_indicator_constraint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_qconstr!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_qpterms!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_rangeconstrs!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_rangeconstrs_t!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_sos!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_var!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

add_vars!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_addrows(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgbds(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgcoef(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgctype(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgname(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgobj(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgobjoffset(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgobjsen(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgrhs(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_chgsense(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_delcols(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_delqconstrs(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_delrows(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_delsos(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getax(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getbaritcnt(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getconflict(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getdj(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getitcnt(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getlb(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getnumcols(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getnumrows(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getobj(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getobjoffset(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getobjsen(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getobjval(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getpi(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getqconstr(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getqconstrslack(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getquad(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getrhs(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getrows(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getsos(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getstat(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getstatstring(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getub(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getx(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_getxqxax(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

c_api_solninfo(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

callback_wrapper(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

callbackgetcandidatepoint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

callbackgetrelaxationpoint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbabort(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbaddboundbranchdown!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbaddboundbranchdown!

cbaddboundbranchup!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbaddboundbranchup!

cbaddconstrbranch!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbaddusercuts(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbbranch(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbbranchconstr(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbcandidateispoint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbcut(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbcutlocal(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbgetcandidatepoint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbgetdetstarttime(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetdetstarttime

cbgetdettimestamp(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetdettimestamp

cbgetfeasibility(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetfeasibility

cbgetgap(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetgap

cbgetintfeas(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetintfeas

cbgetmipiterations(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetmipiterations

cbgetnodelb(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetnodelb

cbgetnodeobjval(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetnodeobjval

cbgetnodesleft(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetnodesleft

cbgetnodeub(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetnodeub

cbgetrelaxationpoint(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbgetstarttime(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgetstarttime

cbgettimestamp(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export cbgettimestamp

cblazy(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cblazylocal(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbpostheursoln(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbprocessincumbent!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbrejectcandidate(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cbsetfunc(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cchar(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

check_moi_callback_validity(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

close_CPLEX(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

compute_conflict(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cplex_model(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cvec(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

cvecx(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

default_moi_callback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

eval(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

free_problem(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

fvec(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

fvecx(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_basis(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_best_bound(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constrLB(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constrUB(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constr_duals(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constr_matrix(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constr_senses(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_constr_solution(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_error_msg(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_infeasibility_ray(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_nnz(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_node_count(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_num_cuts(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_num_sos(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_obj(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_objval(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_param(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_param_type(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_prob_type(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_reduced_costs(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_rel_gap(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_rhs(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_sense(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_solution(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_status(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_status_code(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_unbounded_ray(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_varLB(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_varUB(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

get_vartype(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

getdettime(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

include(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

intervalize(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

is_valid(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

ivec(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

masterbranchcallback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

mastercallback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

masterheuristiccallback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

masterincumbentcallback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

masterinfocallback(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

newlongannotation(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export newlongannotation

notify_freed_model(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

notify_new_model(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

num_constr(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

num_qconstr(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

num_var(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

optimize!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

read_model(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

return_status_or_throw(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_branching_priority(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_constrLB!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_constrUB!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_constr_senses!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_logfile(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_obj!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_param!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_prob_type!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_rhs!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_sense!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_terminate(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_varLB!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_varUB!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_varname!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_vartype!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

set_warm_start!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setbranchcallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setcallbackcut(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setcallbackcutlocal(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setincumbentcallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setlongannotations(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export setlongannotations

setmathprogbranchcallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export setmathprogbranchcallback!

setmathprogcutcallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setmathprogheuristiccallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setmathprogincumbentcallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setmathproginfocallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

setmathproglazycallback!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

terminate(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

toggleproblemtype!(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

tune_param(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

version(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

write_model(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)

CplexSolver(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)
export CplexSolver


text:
function write_input(	file_path::String,
					particle_type::String,
					R::Array{Float64, 2},
					Lx::Float64,
					Ly::Float64,
					Lz::Float64,
					phi::Float64,
					phi_target::Float64,
					position_constraint_axis::String,
					position_constraint_lower::Float64,
					position_constraint_upper::Float64,
					orientation_axis::Array{Float64, 1},
					orientation_constraint_axis::Array{Float64, 1},
					orientation_constraint_lower::Float64,
					orientation_constraint_upper::Float64,
					sigma_translation_max::Float64,
					sigma_rotation_max::Float64,
					sigma_ratio::Float64,
					number_of_equilibration_sweeps::Int64,
					delta_phi::Float64,
					output_file_path::String)

	file_stream::IOStream = open(file_path, "w")

	@printf(file_stream, "%s", "<input_generation>\n")

	write_key(file_stream, "particle_type", particle_type)
	write_key(file_stream, "R", R[:])
	write_key(file_stream, "domain_size_x", Lx)
	write_key(file_stream, "domain_size_y", Ly)
	write_key(file_stream, "domain_size_z", Lz)
	write_key(file_stream, "phi", phi)
	write_key(file_stream, "phi_target", phi_target)
	write_key(file_stream, "position_constraint_axis", position_constraint_axis)
	write_key(file_stream, "position_constraint_lower", position_constraint_lower)
	write_key(file_stream, "position_constraint_upper", position_constraint_upper)
	write_key(file_stream, "orientation_axis", orientation_axis)
	write_key(file_stream, "orientation_constraint_axis", orientation_constraint_axis)
	write_key(file_stream, "orientation_constraint_lower", orientation_constraint_lower)
	write_key(file_stream, "orientation_constraint_upper", orientation_constraint_upper)
	write_key(file_stream, "sigma_translation_max", sigma_translation_max)
	write_key(file_stream, "sigma_rotation_max", sigma_rotation_max)
	write_key(file_stream, "sigma_ratio", sigma_ratio)
	write_key(file_stream, "number_of_equilibration_sweeps", number_of_equilibration_sweeps)
	write_key(file_stream, "delta_phi", delta_phi)
	write_key(file_stream, "output_file_path", output_file_path)

	@printf(file_stream, "%s", "</input_generation>")

	close(file_stream)

	nothing
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
module ARules

# export # common types and functions
# FPTree,
# FPEmptyNode,
# FPTreeNode,
# haschild,
# growchild!,
# climb_down!,
# remove_node,
# getnode

include("common.jl")
include("apriori.jl")
include("fp_tree.jl")
# include("fp_growth.jl")

end # module end


text:
BLPName(name::Symbol) = BLPName(blpapi_Name_create(String(name)))
has_name(blp_name::BLPName, name::Symbol) = blp_name.symbol == name
has_name(blp_name::BLPName, name::AbstractString) = has_name(blp_name, Symbol(name))
Base.Symbol(name::BLPName) = name.symbol

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
nlevels(c::CategoricalValue) = length(levels(c.pool))
nlevels(c::CategoricalString) = length(levels(c.pool))

abstract type Found end
    struct Unknown <: Found end
    abstract type Known <: Found end
        abstract type Infinite <: Known end
           struct Continuous <: Infinite end
           struct Count <: Infinite end
        abstract type Finite{N} <: Known end
            struct Multiclass{N} <: Finite{N} end
            struct OrderedFactor{N} <: Finite{N} end

# aliases:
const Other = Unknown # TODO: depreciate:
const Binary = Multiclass{2}

"""
    scitype(x)

Return the scientific type that an object `x` can represent, when
appearing as an element of a table or vector used as input or target
in fitting MLJ models.

    julia> scitype(4.5)
    Continous

    julia> scitype("book")
    Unknown

    julia> using CategoricalArrays
    julia> v = categorical([:m, :f, :f])
    julia> scitype(v[1])
    Multiclass{2}

Note that `scitype` "commutes" with the formation of tuples or arrays,
as these examples illustrate:

```julia
scitype((42, float(π), "Julia"))
```

```julia
Tuple{Count,Continuous,Unknown}
```

```julia
scitype(rand(7,3))
```

```julia
AbstractArray{Continuous,2}
```

For getting the union of the scitypes of all elements of an iterable,
use `scitype_union`.

"""
scitype(::Any) = Unknown
scitype(::Missing) = Missing
scitype(::AbstractFloat) = Continuous
scitype(::Integer) = Count
scitype(c::CategoricalValue) =
    c.pool.ordered ? OrderedFactor{nlevels(c)} : Multiclass{nlevels(c)}
scitype(c::CategoricalString) =
    c.pool.ordered ? OrderedFactor{nlevels(c)} : Multiclass{nlevels(c)}

scitype(t::Tuple) = Tuple{scitype.(t)...}
MLJBase.scitype(A::B) where {T,N,B<:AbstractArray{T,N}} = AbstractArray{scitype(first(A)),N}

"""
    scitype_union(A)

Return the type union, over all elements `x` generated by the iterable
`A`, of `scitype(x)`.

"""
scitype_union(A) = reduce((a,b)->Union{a,b}, (scitype(el) for el in A))

"""
    scitypes(X)

Returns a named tuple keyed on the column names of the table `X` with
values the corresponding scitype unions over a column's entries.

"""
function scitypes(X)
    container_type(X) in [:table, :sparse] ||
        throw(ArgumentError("Container should be a table or sparse table. "))
    names =    schema(X).names
    return NamedTuple{names}(scitype_union(selectcols(X, c)) for c in names)
end


text:
using HorizonSideRobots

function find_marker(r::Robot)
    side = Ost
    spiral_side = 1
    while ismarker(r) == false
        for i in 1:2 # Каждые два поворота нужно увеличить длину стороны спирали на 1
            move_in_spiral(r, side, spiral_side)
            side = HorizonSide(mod(Int(side) - 1, 4)) # Следующее направление по часовой стрелке
        end
        spiral_side += 1
    end
end

function move_in_spiral(r::Robot, side::HorizonSide, n::Int)
    for i in 1:n
        if ismarker(r)
            return 0
        end
        move!(r, side)
    end
end


text:
"""
	update_core(params::ParametersOrthophoto, laspoint::LasIO.LasPoint, h::LasIO.LasHeader, s::Union{Nothing,IOStream}, n::Int64)
"""
function processPoint(params::ParametersOrthophoto, point::Point)
    position = point.position
    rgb = point.color
    pt = params.coordsystemmatrix * position
    xcoord = map(Int ∘ trunc, (pt[1] - params.refX) / params.GSD) + 1
    ycoord = map(Int ∘ trunc, (params.refY - pt[2]) / params.GSD) + 1

    # color for image
    if params.rasterquote[ycoord, xcoord] < pt[3]
        params.rasterquote[ycoord, xcoord] = pt[3]
        params.RGBtensor[1, ycoord, xcoord] = rgb[1]
        params.RGBtensor[2, ycoord, xcoord] = rgb[2]
        params.RGBtensor[3, ycoord, xcoord] = rgb[3]
    end

    # point for point cloud
    if params.pc
        Common.update_boundingbox!(
            params.tightBB,
            vcat(Common.apply_matrix(params.ucs, position)...),
        )
        plas = newPointRecord(
            point,
            LasIO.LasPoint2,
            params.mainHeader;
            affineMatrix = params.ucs,
        )
        # plas = FileManager.newPointRecord(laspoint,header,LasIO.LasPoint2,params.mainHeader; affineMatrix = params.ucs)
        write(params.stream_tmp, plas) # write this record on temporary file
    end

    params.numPointsProcessed += 1

end

"""
	updateWithControl!(params::Union{ParametersOrthophoto,ParametersExtraction}, file::String, s::Union{Nothing,IOStream}, n::Int64)

Process all points, in file, falling in region of interest.
"""
function updateWithControl!(params::ParametersOrthophoto, file::String)
    header, laspoints = FileManager.read_LAS_LAZ(file) # read file
    for laspoint in laspoints # read each point
        # point = FileManager.xyz(laspoint,header)
        point = Point(laspoint, header)
        if Common.inmodel(params.model)(point.position) # if point in model
            # processPoint(params,laspoint,header)
            processPoint(params, point)
        end
    end
end

"""
	updateWithoutControl!(params::Union{ParametersOrthophoto,ParametersExtraction}, file::String, s::Union{Nothing,IOStream}, n::Int64)

Process points in file without further checks.
"""
function updateWithoutControl!(params::ParametersOrthophoto, file::String)
    header, laspoints = FileManager.read_LAS_LAZ(file) # read file
    for laspoint in laspoints # read each point
        point = Point(laspoint, header)
        processPoint(params, point)
        # processPoint(params,laspoint,header)
    end
end