text:
using BenchmarkTools
using Flux
using CUDA
using Zygote: pullback, ignore

fw(m, x) = m(x)
bw(back) = back(1f0)
fwbw(m, ps, x) = gradient(() -> sum(fw(m, x)), ps)
pb(m, ps, x) = pullback(() -> sum(fw(m, x)), ps)

function run_benchmark(model, x; cuda=true)

    if cuda
        model = model |> gpu
        x = x |> gpu
    end

    ps = Flux.params(model)
    y, back =  pb(model, ps, x)

    if cuda
        CUDA.allowscalar(false)
        # CUDA.device!(3)
        println("  forward")
        fw(model, x); GC.gc(); CUDA.reclaim(); #warmup
        @btime CUDA.@sync(fw($model, $x)) teardown=(GC.gc(); CUDA.reclaim())

        println("  backward")
        bw(back); GC.gc(); CUDA.reclaim(); #warmup
        @btime CUDA.@sync(bw($back)) teardown=(GC.gc(); CUDA.reclaim())

        println("  forw and back")
        fwbw(model, ps, x); GC.gc(); CUDA.reclaim(); #warmup
        @btime CUDA.@sync(fwbw($model, $ps, $x)) teardown=(GC.gc(); CUDA.reclaim())
    else
        println("  forward")
        fw(model, x)  #warmup
        @btime fw($model, $x)

        println("  backward")
        bw(back)  #warmup
        @btime bw($back)

        println("  forw and back")
        fwbw(model, ps, x) # warmup
        @btime fwbw($model, $ps, $x)
    end
end


text:
function draw_hg(s::Sampling)
    xs = Float64[]
    ys = Float64[]
    yerrs = Float64[]
    for (n, hg) in s.hg
        N = n^2
        x = log2(N)
        hg = (hg * x) / sqrt(n)
        y = Statistics.mean(hg)
        yerr = Statistics.std(hg)
        push!(xs, x)
        push!(ys, y)
        push!(yerrs, yerr)
    end
    Plots.scatter(xs, ys, yerror=yerrs, grid=false, legend=false, mc=:black, shape=:hline,
                  bg=:transparent, fg=:black, yguidefontsize=14, xguidefontsize=14)
    Plots.xlabel!(L"\log_2{|V|}")
    Plots.ylabel!(L"H_G")
    Plots.hline!([0, 0], color=:gray, linestyle=:dash)
end

function draw_gg(s::Sampling)
    xs = Float64[]
    ys = Float64[]
    yerrs = Float64[]
    for (n, ngs) in s.ng
        N = n^2
        x = log2(N)
        model = geodesic_model(ngs)
        y = GLM.coef(model)[2]
        yerr = GLM.stderror(model)[2]
        push!(xs, x)
        push!(ys, y)
        push!(yerrs, yerr)
    end
    Plots.scatter(xs, ys, yerror=yerrs, grid=false, legend=false, mc=:black, shape=:hline,
                  bg=:transparent, fg=:black, yguidefontsize=14, xguidefontsize=14)
    Plots.xlabel!(L"\log_2{|V|}")
    Plots.ylabel!(L"\gamma")
    Plots.hline!([1, 1], color=:gray, linestyle=:dash)
end

function draw_tg(s::Sampling)
    xs = Float64[]
    ys = Float64[]
    yerrs = Float64[]
    for (n, tgs) in s.tg
        N = n^2
        x = log2(N)
        ts = collect(Iterators.flatten(values(tgs)))
        y = Statistics.mean(ts)
        yerr = Statistics.std(ts)
        push!(xs, x)
        push!(ys, y)
        push!(yerrs, yerr)
    end
    Plots.scatter(xs, ys, yerror=yerrs, grid=false, legend=false, mc=:black, shape=:hline,
                  bg=:transparent, fg=:black, yguidefontsize=14, xguidefontsize=14)
    Plots.xlabel!(L"\log_2{|V|}")
    Plots.ylabel!(L"T_G")
    Plots.hline!([0, 0], color=:gray, linestyle=:dash)
end

function draw_ag(s::Sampling)
    xs = Float64[]
    ys = Float64[]
    yerrs = Float64[]
    for (n, ags) in s.ag
        N = n^2
        x = log2(N)
        ag = collect(Iterators.flatten(values(ags)))
        y = Statistics.mean(ag)
        yerr = Statistics.std(ag)
        push!(xs, x)
        push!(ys, y)
        push!(yerrs, yerr)
    end
    Plots.scatter(xs, ys, yerror=yerrs, grid=false, legend=false, mc=:black, shape=:hline,
                  bg=:transparent, fg=:black, yguidefontsize=14, xguidefontsize=14)
    Plots.xlabel!(L"\log_2{|V|}")
    Plots.ylabel!(L"A_G")
    Plots.hline!([0, 0], color=:gray, linestyle=:dash)
end

function draw_rand(s::Sampling)
    draw_hg(s)
    Plots.savefig("figures/rand_hg.svg")
    draw_gg(s)
    Plots.savefig("figures/rand_gg.svg")
    draw_tg(s)
    Plots.savefig("figures/rand_tg.svg")
    draw_ag(s)
    Plots.savefig("figures/rand_ag.svg")
end

function draw_gamma(s::Sampling)
    draw_hg(s)
    Plots.savefig("figures/gamma_hg.svg")
    draw_gg(s)
    Plots.savefig("figures/gamma_gg.svg")
    draw_tg(s)
    Plots.savefig("figures/gamma_tg.svg")
    draw_ag(s)
    Plots.savefig("figures/gamma_ag.svg")
end

function draw_gamma_disordered(s::Sampling)
    draw_gg(s)
    Plots.savefig("figures/gamma_disordered_gg.svg")
end

ext: jl

max_stars_repo_name: moshelooks/eugrid

max_issues_repo_name: moshelooks/eugrid

max_forks_repo_name: moshelooks/eugrid


text:
# ### Diagonal ops

# we use Base.Diagonal(f) for diagonal operators so very little specific code is
# actually needed here.

simulate(rng::AbstractRNG, D::DiagOp; Nbatch=nothing) =
    sqrt(D) * white_noise(similar(diag(D), (isnothing(Nbatch) || Nbatch==1 ? () : (Nbatch,))...), rng)
global_rng_for(D::DiagOp) = global_rng_for(diag(D))

# automatic basis conversion (and NaN-zeroing)
(*)(D::DiagOp{<:Field{B}}, f::Field) where {B} = diag(D) .* B(f)
(\)(D::DiagOp{<:Field{B}}, f::Field) where {B} = nan2zero.(diag(D) .\ B(f))

# the generic versions of these kind of suck so we need these specializized
# versions:
(*)(x::Adjoint{<:Any,<:Field}, D::DiagOp) = (D*parent(x))'
(*)(x::Adjoint{<:Any,<:Field}, D::DiagOp, y::Field) = x*(D*y)
diag(D::DiagOp) = D.diag
(^)(D::DiagOp, p::Integer) = Diagonal(diag(D).^p)
pinv(D::DiagOp) = Diagonal(pinv.(diag(D)))
(≈)(D₁::DiagOp, D₂::DiagOp) = diag(D₁) ≈ diag(D₂)

# use getproperty here to ensure no basis conversion is done
getindex(D::DiagOp, s::Symbol) = Diagonal(getproperty(diag(D),s))

# the generic version of this is prohibitively slow so we need this
hash(D::DiagOp, h::UInt64) = foldr(hash, (typeof(D), D.diag), init=h)

# adapting
get_storage(D::DiagOp) = get_storage(diag(D))

# allows broadcasting over DiagOps
struct DiagOpStyle{S} <: AbstractArrayStyle{2} end
BroadcastStyle(::Type{Diagonal{T,F}}) where {T,F<:Field} = DiagOpStyle{typeof(BroadcastStyle(F))}()
BroadcastStyle(::DiagOpStyle{S₁}, ::DiagOpStyle{S₂}) where {S₁,S₂} = DiagOpStyle{typeof(result_style(S₁(), S₂()))}()
materialize(bc::Broadcasted{DiagOpStyle{S}}) where {S} = Diagonal(materialize(convert(Broadcasted{S}, preprocess(DiagOpStyle{S}(), bc))))
materialize!(dest::DiagOp, bc::Broadcasted{DiagOpStyle{S}}) where {S} =
    (materialize!(diag(dest), convert(Broadcasted{S}, preprocess(DiagOpStyle{S}(), bc))); dest)
preprocess(::DiagOpStyle, D::Diagonal) = diag(D)

### Derivative ops

# These ops just store the coordinate with respect to which a derivative is
# being taken, and each Field type F implements how to actually take the
# derivative.
#
# ∇diag(coord, covariance, prefactor) is a Field which represents the diagonal
# entries of the Fourier representation of the gradient operation, with `coord`
# specifiying the coordinate for the gradient direction (1 or 2), `covariance`
# being COVARIANT or CONTRAVARIANT, and multiplied by a `prefactor` (generally
# 1 or -1 to keep track of if an adjoint has been taken). ∇[i] returns one of
# these wrapped in a Diagonal so that multipliying by these does the necessary
# basis conversion. ∇ (along with ∇ⁱ, and ∇ᵢ) are StaticVectors which can also
# be used with FieldVector algebra, e.g. ∇*f.
#
# Also note: We define the components of vectors, including ∇, to be with respect to
# the _unnormalized_ covariant or contravariant basis vectors, hence ∇ⁱ = d/dxᵢ
# and ∇ᵢ = d/dxⁱ. This is different than with respect to the _normalized)
# covariant basis vectors, which, e.g., in spherical coordinates, gives the more
# familiar ∇ = (d/dθ, 1/sinθ d/dϕ), (but whose components are neither covariant
# nor contravariant).

Base.@enum Covariance COVARIANT CONTRAVARIANT

struct ∇diag <: ImplicitField{DerivBasis,Complex{Bottom}}
    coord :: Int
    covariance :: Covariance
    prefactor :: Int
end

struct ∇²diag <: Field{DerivBasis,Complex{Bottom}} end

get_g_metadata(::∇diag) = nothing

# adjoint(D::Diagonal{<:∇diag}) calls conj(D.diag) and here we keep
# track of that a conjugate was taken
conj(d::∇diag) = ∇diag(d.coord, d.covariance, -d.prefactor)
-(d::DiagOp{∇diag}) = Diagonal(conj(d))

# Gradient vector which can be used with FieldVector algebra.
# struct ∇Op <: StaticVector{2,Diagonal{Float32,∇diag}} end
# getindex(::∇Op, i::Int) = Diagonal(∇diag{i,covariance,prefactor}())
# -(::∇Op{covariance,prefactor}) where {covariance,prefactor} = ∇Op{covariance,-prefactor}()
const ∇ⁱ = @SVector[Diagonal(∇diag(coord, CONTRAVARIANT, 1)) for coord=1:2]
const ∇ᵢ = @SVector[Diagonal(∇diag(coord, COVARIANT,     1)) for coord=1:2]
const ∇ = ∇ⁱ # ∇ is contravariant by default if not specified
const ∇² = Diagonal(∇²diag())

@doc doc"""
    gradhess(f)

Compute the gradient $g^i = \nabla^i f$, and the hessian, $H_j^{\,i} = \nabla_j \nabla^i f$.
"""
function gradhess(f)
    g = ∇ⁱ*f
    H = @SMatrix[∇ᵢ[1]*g[1] ∇ᵢ[2]*g[1]; ∇ᵢ[1]*g[2] ∇ᵢ[2]*g[2]]
    (;g, H)
end

# ### FuncOp

# An Op which applies some arbitrary function to its argument.
# Transpose and/or inverse operations which are not specified will return an error.
@kwdef struct FuncOp <: ImplicitOp{Bottom}
    op   = nothing
    opᴴ  = nothing
    op⁻¹ = nothing
    op⁻ᴴ = nothing
end
SymmetricFuncOp(;op=nothing, op⁻¹=nothing) = FuncOp(op,op,op⁻¹,op⁻¹)
FuncOp(op::Function) = FuncOp(op=op)
SymmetricFuncOp(op::Function) = SymmetricFuncOp(op=op)
*(L::FuncOp, f::Field) = L.op   != nothing ? L.op(f)   : error("op*f not implemented")
\(L::FuncOp, f::Field) = L.op⁻¹ != nothing ? L.op⁻¹(f) : error("op\\f not implemented")
*(f::Adjoint{<:Any,Field}, L::FuncOp) = L.opᴴ  != nothing ? L.opᴴ(f)  : error("opᴴ*f not implemented")
adjoint(L::FuncOp) = FuncOp(L.opᴴ,L.op,L.op⁻ᴴ,L.op⁻¹)
inv(L::FuncOp) = FuncOp(L.op⁻¹,L.op⁻ᴴ,L.op,L.opᴴ)
adapt_structure(to, L::FuncOp) = FuncOp(adapt(to, fieldvalues(L))...)
hash(L::FuncOp, h::UInt64) = foldr(hash, (typeof(L), fieldvalues(L)...), init=h)

### BandPassOp

# An op which applies some bandpass, like a high or low-pass filter. This object
# stores the bandpass weights, Wℓ, and each Field type F should implement
# preprocess((b,metadata), ::BandPassOp) to describe how this is actually
# applied.

abstract type HarmonicBasis <: Basislike end

struct BandPass{W<:InterpolatedCℓs} <: ImplicitField{HarmonicBasis,Bottom}
    Wℓ::W
end
BandPassOp(ℓ,Wℓ) = Diagonal(BandPass(InterpolatedCℓs(promote(collect(ℓ),collect(Wℓ))...)))
BandPassOp(Wℓ::InterpolatedCℓs) = Diagonal(BandPass(Wℓ))
cos_ramp_up(length) = @. (cos($range(π,0,length=length))+1)/2
cos_ramp_down(length) = 1 .- cos_ramp_up(length)
HighPass(ℓ; Δℓ=50) = BandPassOp(ℓ:20000, [cos_ramp_up(Δℓ); ones(20000-ℓ-Δℓ+1)])
LowPass(ℓ; Δℓ=50) = BandPassOp(0:ℓ, [ones(ℓ-Δℓ+1); cos_ramp_down(Δℓ)])
MidPass(ℓmin, ℓmax; Δℓ=50) = BandPassOp(ℓmin:ℓmax, [cos_ramp_up(Δℓ); ones(ℓmax-ℓmin-2Δℓ+1); cos_ramp_down(Δℓ)])
MidPasses(ℓedges; Δℓ=10) = [MidPass(ℓmin-Δℓ÷2,ℓmax+Δℓ÷2; Δℓ=Δℓ) for (ℓmin,ℓmax) in zip(ℓedges[1:end-1],ℓedges[2:end])]
hash(b::BandPass, h::UInt64) = foldr(hash, (typeof(b), b.Wℓ), init=h)

### ParamDependentOp

# A FieldOp which depends on some parameters, θ.
# L(;θ...) recomputes the operator at a given set of parameters, but the
# operator can also be used as-is in which case it is evaluated at a fiducial θ
# (which is stored inside the operator when it is first constructed).

@doc doc"""
    ParamDependentOp(recompute_function::Function)

Creates an operator which depends on some parameters $\theta$ and can be
evaluated at various values of these parameters.

`recompute_function` should be a function which accepts keyword arguments for
$\theta$ and returns the operator. Each keyword must have a default value; the
operator will act as if evaluated at these defaults unless it is explicitly
evaluated at other parameters.

Example:

```julia
Cϕ₀ = Diagonal(...) # some fixed Diagonal operator
Cϕ = ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀) # create ParamDependentOp

Cϕ(Aϕ=1.1) * ϕ   # Cϕ(Aϕ=1.1) is equal to 1.1*Cϕ₀
Cϕ * ϕ           # Cϕ alone will act like Cϕ(Aϕ=1) because that was the default above
```

Note: if you are doing parallel work, global variables referred to in the
`recompute_function` need to be distributed to all workers. A more robust
solution is to avoid globals entirely and instead ensure all variables are
"closed" over (and hence will automatically get distributed). This will happen
by default if defining the `ParamDependentOp` inside any function, or can be
forced at the global scope by wrapping everything in a `let`-block, e.g.:

```julia
Cϕ = let Cϕ₀=Cϕ₀
    ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀)
end
```

After executing the code above, `Cϕ` is now ready to be (auto-)shipped to any workers
and will work regardless of what global variables are defined on these workers.
"""
struct ParamDependentOp{T, L<:FieldOp{T}, F<:Function} <: ImplicitOp{T}
    op :: L
    recompute_function :: F
    parameters :: Vector{Symbol}
end
function ParamDependentOp(recompute_function::Function)
    kwarg_names = filter(!=(Symbol("_...")), get_kwarg_names(recompute_function))
    if endswith(string(kwarg_names[end]), "...")
        # an "kwargs..."-like keyword argument (except for "_..."
        # which is filtered out above) indicates this operator depends
        # on everything, which is indicated by an empty kwarg_names
        empty!(kwarg_names)
    end
    # invokelatest here allows creating a ParamDependent op which calls a
    # BinRescaledOp (eg this is the case for the mixing matrix G which depends
    # on Cϕ) from inside function. this would otherwise fail due to
    # BinRescaledOp eval'ed function being too new
    ParamDependentOp(Base.invokelatest(recompute_function), recompute_function, kwarg_names)
end
function (L::ParamDependentOp)(θ::NamedTuple)
    if depends_on(L,θ)
        # if L got adapt'ed to CuArray since this op was created,
        # L.op will be GPU-backed, but depending on how
        # recompute_function is written, recompute_function may
        # still return something CPU-backed. in that case, copy it
        # to GPU here
        storage = get_storage(L.op)
        Lθ = L.recompute_function(;θ...)
        basetype(storage) == basetype(get_storage(Lθ)) ? Lθ : adapt(storage, Lθ)
    else
        L.op
    end
end
(L::ParamDependentOp)(;θ...) = L((;θ...))
@init @require ComponentArrays="b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" begin
    using ComponentArrays
    (L::ParamDependentOp)(θ::ComponentArray) = L(convert(NamedTuple, θ))
    (L::Union{FieldOp,UniformScaling})(::ComponentArray) = L
end

@auto_adjoint *(L::ParamDependentOp, f::Field) = L.op * f
@auto_adjoint \(L::ParamDependentOp, f::Field) = L.op \ f
for F in (:inv, :pinv, :sqrt, :adjoint, :Diagonal, :diag, :simulate, :zero, :one, :logdet, :global_rng_for)
    @eval $F(L::ParamDependentOp) = $F(L.op)
end
getindex(L::ParamDependentOp, x) = getindex(L.op, x)
simulate(rng::AbstractRNG, L::ParamDependentOp) = simulate(rng, L.op)
depends_on(L::ParamDependentOp, θ) = depends_on(L, keys(θ))
depends_on(L::ParamDependentOp, θ::Tuple) = isempty(L.parameters) || any(L.parameters .∈ Ref(θ))
depends_on(L,                   θ) = false

typealias_def(::Type{<:ParamDependentOp{T,L}}) where {T,L} = "ParamDependentOp{$(typealias(L))}"
function Base.summary(io::IO, L::ParamDependentOp)
    dependent_params = isempty(L.parameters) ? "..." : join(map(string, L.parameters), ",")
    print(io, join(size(L.op), "×"), " (", dependent_params, ")-dependent ")
    Base.showarg(io, L, true)
end

# we have to include recompute_function in the hash, but its hash
# might change if shipped to a distributed worker, despite being the
# same function. not sure if there's any way around this...
hash(L::ParamDependentOp, h::UInt64) = foldr(hash, (typeof(L), L.op, L.recompute_function), init=h)

adapt_structure(to, L::ParamDependentOp) =
    ParamDependentOp(adapt(to, L.op), adapt(to, L.recompute_function), L.parameters)

### LazyBinaryOp

# we use LazyBinaryOps to create new operators composed from other
# operators which don't actually evaluate anything until they've been
# multiplied by a field.
# L::LazyBinaryOp{λ} lazily represents λ(L.X,L.Y)
struct LazyBinaryOp{λ} <: ImplicitOp{Bottom}
    X :: FieldOpScal
    Y :: FieldOpScal
    LazyBinaryOp(λ, X::Union{FieldOp,Scalar}, Y::Union{FieldOp,Scalar}) = new{λ}(X, Y)
end

# creating LazyBinaryOps
for λ in (:+, :-, :*)
    @eval begin
        function ($λ)(
            X :: Union{ImplicitOp, Adjoint{<:Any,<:ImplicitOp}, DiagOp{<:Field{B₁}}},
            Y :: Union{ImplicitOp, Adjoint{<:Any,<:ImplicitOp}, DiagOp{<:Field{B₂}}}
        ) where {B₁,B₂}
            LazyBinaryOp($λ, X, Y)
        end
        function ($λ)(
            X :: DiagOp{<:Field{B}},
            Y :: DiagOp{<:Field{B}}
        ) where {B}
            Diagonal(broadcast($λ, diag(X), diag(Y)))
        end
    end
end
(*)(X::ImplicitOp,              Y::Scalar)     = LazyBinaryOp(*, X, Y)
(*)(X::Scalar,                  Y::ImplicitOp) = LazyBinaryOp(*, X, Y)
(/)(X::ImplicitOp,              Y::Real)       = LazyBinaryOp(/, X, Y)
(^)(X::ImplicitOp,              Y::Integer)    = LazyBinaryOp(^, X, Y)
(^)(X::DiagOp{<:ImplicitField}, Y::Integer)    = LazyBinaryOp(^, X, Y)
(-)(L::ImplicitOp)                             = LazyBinaryOp(*, -1, L)
pinv(L::LazyBinaryOp{*})                       = LazyBinaryOp(*, pinv(L.Y), pinv(L.X))
adjoint(L::LazyBinaryOp{*})                    = LazyBinaryOp(*, adjoint(L.Y), adjoint(L.X))

# evaluating LazyBinaryOps
for λ in (:+, :-)
    @eval (*)(L::LazyBinaryOp{$λ}, f::Field) = ($λ)(L.X * f, L.Y * f)
    @eval diag(L::LazyBinaryOp{$λ}) = ($λ)(diag(L.X), diag(L.Y))
    @eval adjoint(L::LazyBinaryOp{$λ}) = LazyBinaryOp(($λ), adjoint(L.X), adjoint(L.Y))
end
@auto_adjoint (*)(L::LazyBinaryOp{/}, f::Field) = (L.X * f) / L.Y
@auto_adjoint (*)(L::LazyBinaryOp{*}, f::Field) = L.X * (L.Y * f)
@auto_adjoint (\)(L::LazyBinaryOp{*}, f::Field) = L.Y \ (L.X \ f)
@auto_adjoint (*)(L::LazyBinaryOp{^}, f::Field) = foldr((L.Y>0 ? (*) : (\)), fill(L.X, abs(L.Y::Integer)), init=f)
adapt_structure(to, L::LazyBinaryOp{λ}) where {λ} = LazyBinaryOp(λ, adapt(to,L.X), adapt(to,L.Y))
hash(L::LazyBinaryOp, h::UInt64) = foldr(hash, (typeof(L), L.X, L.Y), init=h)


text:
module YASSAD

include("utils.jl")
include("compiler.jl")
include("lib.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "librealsense"
version = v"2.16.0"

# Collection of sources required to build librealsense
sources = [
    "https://github.com/IntelRealSense/librealsense.git" =>
    "b55069d177f69daea3b9437532f4b98fbe97a718",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd librealsense/
mkdir build
cd build
CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_EXAMPLES=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_GRAPHICAL_EXAMPLES=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_UNIT_TESTS=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_WITH_OPENMP=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DENFORCE_METADATA=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_PYTHON_BINDINGS=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_NODEJS_BINDINGS=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DFORCE_LIBUVC=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DTRACE_API=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DHWM_OVER_XU=true"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_SHARED_LIBS=true"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=Release"
CMAKE_FLAGS="${CMAKE_FLAGS} -DENABLE_ZERO_COPY=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_EASYLOGGINGPP=true"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_CV_EXAMPLES=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_PCL_EXAMPLES=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_WITH_TM2=false"
CMAKE_FLAGS="${CMAKE_FLAGS} -DBUILD_WITH_STATIC_CRT=true"
cmake ${CMAKE_FLAGS} -DCMAKE_SYSTEM_VERSION=0 ..
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "librealsense2", :librealsense2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPerception/LibusbBuilder/releases/download/v1.0.22-1/build_libusb.v1.0.22.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using PanelCountData
using Test

@testset "PanelCountData.jl" begin
    # Write your tests here.
    @test welcome() == "welcome to Panel Count."
end


text:
using CxxrtlCosim

cd(@__DIR__)

if !isfile("blinky.so")
    println("compiling blinky.v")
    CxxrtlCosim.synth("blinky.cpp", "blinky.v")
    CxxrtlCosim.compile("blinky.so", "blinky.cpp")
end

top = CxxrtlCosim.cxxrtl_design_create("./blinky.so")

blink = CxxrtlCosim.cxxrtl_create(top)

clk = CxxrtlCosim.cxxrtl_get(blink, "clk")
led = CxxrtlCosim.cxxrtl_get(blink, "led")

CxxrtlCosim.cxxrtl_reset(blink)

CxxrtlCosim.cxxrtl_step(blink)
prev_led = 0
for i = 1:10000
    unsafe_store!(unsafe_load(clk).next, UInt32(1))
    CxxrtlCosim.cxxrtl_step(blink)
    unsafe_store!(unsafe_load(clk).next, UInt32(0))
    CxxrtlCosim.cxxrtl_step(blink)
    global ledval = unsafe_load(unsafe_load(led).curr)
    if prev_led != ledval
        println("$i: $ledval")
        global prev_led = ledval
    end
end

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


text:
using BenchmarkTools
using Random
using CovarianceEstimation

const SUITE = BenchmarkGroup()

Random.seed!(1234)

const X10_100  = randn(15, 100)
const X100_100  = randn(100, 100)
const X100_1000  = randn(100, 1000)
const X100_10000  = randn(100, 10000)
const X1000_100  = randn(1000, 100)
const X1000_10000  = randn(1000, 10000)

matrices = [X10_100, X100_100, X100_1000, X100_10000, X1000_10000]
estimators = Dict(
    "Simple" => (ce = SimpleCovariance(), maxfeatures = 10000),
    "Linear constant correlation" => (ce = LinearShrinkage(ConstantCorrelation()), maxfeatures = 1000),
    "Linear diagonal unit variance" => (ce = LinearShrinkage(DiagonalUnitVariance()), maxfeatures = 1000),
    "Linear diagonal common variance" => (ce = LinearShrinkage(DiagonalCommonVariance()), maxfeatures = 1000),
    "Linear common covariance" => (ce = LinearShrinkage(CommonCovariance()), maxfeatures = 1000),
    "Linear diagonal unequal variance" => (ce = LinearShrinkage(DiagonalUnequalVariance()), maxfeatures = 1000),
    "Linear perfect positive correlation" => (ce = LinearShrinkage(PerfectPositiveCorrelation()), maxfeatures = 1000),
    "Linear RBLW" => (ce = LinearShrinkage(DiagonalCommonVariance(), :rblw), maxfeatures = 10000),
    "Linear OAS" => (ce = LinearShrinkage(DiagonalCommonVariance(), :oas), maxfeatures = 10000),
    "Analytical nonlinear shrinkage" => (ce = AnalyticalNonlinearShrinkage(), maxfeatures = 1000)
)

for (k, v) in estimators
    SUITE[k] = BenchmarkGroup()
    for m in matrices
        if size(m, 2) ≤ v[:maxfeatures]
            SUITE[k]["Size $(size(m))"] = @benchmarkable cov($v[:ce], $m)
        end
    end
end


text:
# This file was generated, do not modify it. # hide
cm = X |> Matrix |> cor
round.(cm, sigdigits=1)


text:
module ModuleCarriers
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using StaticArrays
export ∂ₓu, uₓ
include("carriers.jl")
end


text:
###############################################################################
#
#   acb_poly.jl : Polynomials over arb
#
###############################################################################

export AcbPolyRing, acb_poly, isreal, derivative, integral, evaluate,
       evaluate2, compose, from_roots, evaluate_iter, evaluate_fast, evaluate,
       interpolate_newton, interpolate_barycentric, interpolate_fast,
       interpolate, roots

###############################################################################
#
#   Basic manipulation
#
###############################################################################

parent_type(::Type{acb_poly}) = AcbPolyRing

elem_type(::Type{AcbPolyRing}) = acb_poly

length(x::acb_poly) = ccall((:acb_poly_length, :libarb), Int,
                                   (Ref{acb_poly},), x)

set_length!(x::acb_poly, n::Int) = ccall((:_acb_poly_set_length, :libarb), Void,
                                   (Ref{acb_poly}, Int), x, n)

degree(x::acb_poly) = length(x) - 1

function coeff(a::acb_poly, n::Int)
  n < 0 && throw(DomainError())
  t = parent(a).base_ring()
  ccall((:acb_poly_get_coeff_acb, :libarb), Void,
              (Ref{acb}, Ref{acb_poly}, Int), t, a, n)
  return t
end

zero(a::AcbPolyRing) = a(0)

one(a::AcbPolyRing) = a(1)

function gen(a::AcbPolyRing)
   z = acb_poly()
   ccall((:acb_poly_set_coeff_si, :libarb), Void,
        (Ref{acb_poly}, Int, Int), z, 1, 1)
   z.parent = a
   return z
end

# todo: write a C function for this
function isgen(a::acb_poly)
   return isequal(a, gen(parent(a)))
end

#function iszero(a::acb_poly)
#   return length(a) == 0
#end

#function isone(a::acb_poly)
#   return isequal(a, one(parent(a)))
#end

function deepcopy_internal(a::acb_poly, dict::ObjectIdDict)
   z = acb_poly(a)
   z.parent = parent(a)
   return z
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::AcbPolyRing)
  print(io, "Univariate Polynomial Ring in ")
  print(io, string(var(x)))
  print(io, " over ")
  show(io, x.base_ring)
end

function show(io::IO, f::acb_poly)
  if length(f) == 0
    print(io, "0")
  else
    print(io, "[ ")
    for i in 0:degree(f)-1
      print(io, coeff(f,i))
      print(io, ", ")
    end
    print(io, coeff(f,degree(f)))
    print(io, " ]")
  end
end

###############################################################################
#
#   Comparisons
#
###############################################################################

function isequal(x::acb_poly, y::acb_poly)
   return ccall((:acb_poly_equal, :libarb), Bool,
                                      (Ref{acb_poly}, Ref{acb_poly}), x, y)
end

doc"""
    overlaps(x::acb_poly, y::acb_poly)
> Return `true` if the coefficient boxes of $x$ overlap the coefficient boxes
> of $y$, otherwise return `false`.
"""
function overlaps(x::acb_poly, y::acb_poly)
   return ccall((:acb_poly_overlaps, :libarb), Bool,
                                      (Ref{acb_poly}, Ref{acb_poly}), x, y)
end

doc"""
    contains(x::acb_poly, y::acb_poly)
> Return `true` if the coefficient boxes of $x$ contain the corresponding
> coefficient boxes of $y$, otherwise return `false`.
"""
function contains(x::acb_poly, y::acb_poly)
   return ccall((:acb_poly_contains, :libarb), Bool,
                                      (Ref{acb_poly}, Ref{acb_poly}), x, y)
end

doc"""
    contains(x::acb_poly, y::fmpz_poly)
> Return `true` if the coefficient boxes of $x$ contain the corresponding
> exact coefficients of $y$, otherwise return `false`.
"""
function contains(x::acb_poly, y::fmpz_poly)
   return ccall((:acb_poly_contains_fmpz_poly, :libarb), Bool,
                                      (Ref{acb_poly}, Ref{fmpz_poly}), x, y)
end

doc"""
    contains(x::acb_poly, y::fmpq_poly)
> Return `true` if the coefficient boxes of $x$ contain the corresponding
> exact coefficients of $y$, otherwise return `false`.
"""
function contains(x::acb_poly, y::fmpq_poly)
   return ccall((:acb_poly_contains_fmpq_poly, :libarb), Bool,
                                      (Ref{acb_poly}, Ref{fmpq_poly}), x, y)
end

function ==(x::acb_poly, y::acb_poly)
    if length(x) != length(y)
        return false
    end
    for i = 0:degree(x)
        if !(coeff(x, i) == coeff(y, i))
            return false
        end
    end
    return true
end

function !=(x::acb_poly, y::acb_poly)
    for i = 0:max(degree(x), degree(y))
        if coeff(x, i) != coeff(y, i)
            return true
        end
    end
    return false
end

doc"""
    unique_integer(x::acb_poly)
> Return a tuple `(t, z)` where $t$ is `true` if there is a unique integer
> contained in the (constant) polynomial $x$, along with that integer $z$
> in case it is, otherwise sets $t$ to `false`.
"""
function unique_integer(x::acb_poly)
  z = FmpzPolyRing(var(parent(x)))()
  unique = ccall((:acb_poly_get_unique_fmpz_poly, :libarb), Int,
    (Ref{fmpz_poly}, Ref{acb_poly}), z, x)
  return (unique != 0, z)
end

doc"""
    isreal(x::acb_poly)
> Return `true` if all the coefficients of $x$ are real, i.e. have exact zero
> imaginary parts.
"""
function isreal(x::acb_poly)
  return ccall((:acb_poly_is_real, :libarb), Cint, (Ref{acb_poly}, ), x) != 0
end

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::acb_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:acb_poly_shift_left, :libarb), Void,
      (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, len)
   return z
end

function shift_right(x::acb_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:acb_poly_shift_right, :libarb), Void,
       (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, len)
   return z
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::acb_poly)
  z = parent(x)()
  ccall((:acb_poly_neg, :libarb), Void, (Ref{acb_poly}, Ref{acb_poly}), z, x)
  return z
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(x::acb_poly, y::acb_poly)
  z = parent(x)()
  ccall((:acb_poly_add, :libarb), Void,
              (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
              z, x, y, prec(parent(x)))
  return z
end

function *(x::acb_poly, y::acb_poly)
  z = parent(x)()
  ccall((:acb_poly_mul, :libarb), Void,
              (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
              z, x, y, prec(parent(x)))
  return z
end

function -(x::acb_poly, y::acb_poly)
  z = parent(x)()
  ccall((:acb_poly_sub, :libarb), Void,
              (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
              z, x, y, prec(parent(x)))
  return z
end

function ^(x::acb_poly, y::Int)
  y < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:acb_poly_pow_ui, :libarb), Void,
              (Ref{acb_poly}, Ref{acb_poly}, UInt, Int),
              z, x, y, prec(parent(x)))
  return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

for T in [Integer, fmpz, fmpq, Float64, BigFloat, arb, acb, fmpz_poly, fmpq_poly]
   @eval begin
      +(x::acb_poly, y::$T) = x + parent(x)(y)

      +(x::$T, y::acb_poly) = y + x

      -(x::acb_poly, y::$T) = x - parent(x)(y)

      -(x::$T, y::acb_poly) = parent(y)(x) - y

      *(x::acb_poly, y::$T) = x * parent(x)(y)

      *(x::$T, y::acb_poly) = y * x
   end
end

+(x::acb_poly, y::Rational{T}) where T <: Union{Int, BigInt} = x + parent(x)(y)

+(x::Rational{T}, y::acb_poly) where T <: Union{Int, BigInt} = y + x

-(x::acb_poly, y::Rational{T}) where T <: Union{Int, BigInt} = x - parent(x)(y)

-(x::Rational{T}, y::acb_poly) where T <: Union{Int, BigInt} = parent(y)(x) - y

*(x::acb_poly, y::Rational{T}) where T <: Union{Int, BigInt} = x * parent(x)(y)

*(x::Rational{T}, y::acb_poly) where T <: Union{Int, BigInt} = y * x

###############################################################################
#
#   Scalar division
#
###############################################################################

for T in [Integer, fmpz, fmpq, Float64, BigFloat, arb, acb]
   @eval begin
      divexact(x::acb_poly, y::$T) = x * inv(base_ring(parent(x))(y))

      //(x::acb_poly, y::$T) = divexact(x, y)
   end
end

divexact(x::acb_poly, y::Rational{T}) where {T <: Integer} = x * inv(base_ring(parent(x))(y))

//(x::acb_poly, y::Rational{T}) where {T <: Integer} = divexact(x, y)

###############################################################################
#
#   Euclidean division
#
###############################################################################

function divrem(x::acb_poly, y::acb_poly)
   iszero(y) && throw(DivideError())
   q = parent(x)()
   r = parent(x)()
   if (ccall((:acb_poly_divrem, :libarb), Int,
         (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
               q, r, x, y, prec(parent(x))) == 1)
      return (q, r)
   else
      throw(DivideError())
   end
end

function mod(x::acb_poly, y::acb_poly)
   return divrem(x, y)[2]
end

function divexact(x::acb_poly, y::acb_poly)
   return divrem(x, y)[1]
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(a::acb_poly, n::Int)
   n < 0 && throw(DomainError())
   if length(a) <= n
      return a
   end
   # todo: implement set_trunc in arb
   z = deepcopy(a)
   ccall((:acb_poly_truncate, :libarb), Void,
                (Ref{acb_poly}, Int), z, n)
   return z
end

function mullow(x::acb_poly, y::acb_poly, n::Int)
   n < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:acb_poly_mullow, :libarb), Void,
         (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int, Int),
            z, x, y, n, prec(parent(x)))
   return z
end

###############################################################################
#
#   Reversal
#
###############################################################################

#function reverse(x::acb_poly, len::Int)
#   len < 0 && throw(DomainError())
#   z = parent(x)()
#   ccall((:acb_poly_reverse, :libarb), Void,
#                (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, len)
#   return z
#end

###############################################################################
#
#   Evaluation
#
###############################################################################

function evaluate(x::acb_poly, y::acb)
   z = parent(y)()
   ccall((:acb_poly_evaluate, :libarb), Void,
                (Ref{acb}, Ref{acb_poly}, Ref{acb}, Int),
                z, x, y, prec(parent(y)))
   return z
end

doc"""
    evaluate2(x::acb_poly, y::acb)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::acb)
   z = parent(y)()
   w = parent(y)()
   ccall((:acb_poly_evaluate2, :libarb), Void,
                (Ref{acb}, Ref{acb}, Ref{acb_poly}, Ref{acb}, Int),
                z, w, x, y, prec(parent(y)))
   return z, w
end

function evaluate(x::acb_poly, y::Union{Int,Float64,fmpq,arb})
    return evaluate(x, base_ring(parent(x))(y))
end

function evaluate(x::acb_poly, y::fmpz)
    return evaluate(x, base_ring(parent(x))(y))
end

doc"""
    evaluate2(x::acb_poly, y::Integer)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::Integer)
    return evaluate2(x, base_ring(parent(x))(y))
end

doc"""
    evaluate2(x::acb_poly, y::Float64)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::Float64)
    return evaluate2(x, base_ring(parent(x))(y))
end

doc"""
    evaluate2(x::acb_poly, y::fmpq)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::fmpz)
    return evaluate2(x, base_ring(parent(x))(y))
end

doc"""
    evaluate2(x::acb_poly, y::fmpq)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::fmpq)
    return evaluate2(x, base_ring(parent(x))(y))
end

doc"""
    evaluate2(x::acb_poly, y::arb)
> Return a tuple $p, q$ consisting of the polynomial $x$ evaluated at $y$ and
> its derivative evaluated at $y$.
"""
function evaluate2(x::acb_poly, y::arb)
    return evaluate2(x, base_ring(parent(x))(y))
end

###############################################################################
#
#   Composition
#
###############################################################################

function compose(x::acb_poly, y::acb_poly)
   z = parent(x)()
   ccall((:acb_poly_compose, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
                z, x, y, prec(parent(x)))
   return z
end

###############################################################################
#
#   Derivative and integral
#
###############################################################################

function derivative(x::acb_poly)
   z = parent(x)()
   ccall((:acb_poly_derivative, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, prec(parent(x)))
   return z
end

function integral(x::acb_poly)
   z = parent(x)()
   ccall((:acb_poly_integral, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Int), z, x, prec(parent(x)))
   return z
end

###############################################################################
#
#   Multipoint evaluation and interpolation
#
###############################################################################

function acb_vec(n::Int)
   return ccall((:_acb_vec_init, :libarb), Ref{acb_struct}, (Int,), n)
end

function acb_vec(b::Array{acb, 1})
   v = ccall((:_acb_vec_init, :libarb), Ref{acb_struct}, (Int,), length(b))
   for i=1:length(b)
       ccall((:acb_set, :libarb), Void, (Ref{acb_struct}, Ref{acb}),
           v + (i-1)*sizeof(acb_struct), b[i])
   end
   return v
end

function array(R::AcbField, v::Ref{acb_struct}, n::Int)
   r = Array{acb}(n)
   for i=1:n
       r[i] = R()
       ccall((:acb_set, :libarb), Void, (Ref{acb}, Ref{acb_struct}),
           r[i], v + (i-1)*sizeof(acb_struct))
   end
   return r
end

function acb_vec_clear(v::Ref{acb_struct}, n::Int)
   ccall((:_acb_vec_clear, :libarb), Void, (Ref{acb_struct}, Int), v, n)
end

doc"""
    from_roots(R::AcbPolyRing, b::Array{acb, 1})
> Construct a polynomial in the given polynomial ring from a list of its roots.
"""
function from_roots(R::AcbPolyRing, b::Array{acb, 1})
   z = R()
   tmp = acb_vec(b)
   ccall((:acb_poly_product_roots, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_struct}, Int, Int), z, tmp, length(b), prec(R))
   acb_vec_clear(tmp, length(b))
   return z
end

function evaluate_iter(x::acb_poly, b::Array{acb, 1})
   return acb[evaluate(x, b[i]) for i=1:length(b)]
end

function evaluate_fast(x::acb_poly, b::Array{acb, 1})
   tmp = acb_vec(b)
   ccall((:acb_poly_evaluate_vec_fast, :libarb), Void,
                (Ref{acb_struct}, Ref{acb_poly}, Ref{acb_struct}, Int, Int),
            tmp, x, tmp, length(b), prec(parent(x)))
   res = array(base_ring(parent(x)), tmp, length(b))
   acb_vec_clear(tmp, length(b))
   return res
end

function interpolate_newton(R::AcbPolyRing, xs::Array{acb, 1}, ys::Array{acb, 1})
   length(xs) != length(ys) && error()
   z = R()
   xsv = acb_vec(xs)
   ysv = acb_vec(ys)
   ccall((:acb_poly_interpolate_newton, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_struct}, Ref{acb_struct}, Int, Int),
            z, xsv, ysv, length(xs), prec(R))
   acb_vec_clear(xsv, length(xs))
   acb_vec_clear(ysv, length(ys))
   return z
end

function interpolate_barycentric(R::AcbPolyRing, xs::Array{acb, 1}, ys::Array{acb, 1})
   length(xs) != length(ys) && error()
   z = R()
   xsv = acb_vec(xs)
   ysv = acb_vec(ys)
   ccall((:acb_poly_interpolate_barycentric, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_struct}, Ref{acb_struct}, Int, Int),
            z, xsv, ysv, length(xs), prec(R))
   acb_vec_clear(xsv, length(xs))
   acb_vec_clear(ysv, length(ys))
   return z
end

function interpolate_fast(R::AcbPolyRing, xs::Array{acb, 1}, ys::Array{acb, 1})
   length(xs) != length(ys) && error()
   z = R()
   xsv = acb_vec(xs)
   ysv = acb_vec(ys)
   ccall((:acb_poly_interpolate_fast, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_struct}, Ref{acb_struct}, Int, Int),
            z, xsv, ysv, length(xs), prec(R))
   acb_vec_clear(xsv, length(xs))
   acb_vec_clear(ysv, length(ys))
   return z
end

# todo: cutoffs for fast algorithm
function interpolate(R::AcbPolyRing, xs::Array{acb, 1}, ys::Array{acb, 1})
   return interpolate_newton(R, xs, ys)
end

# todo: cutoffs for fast algorithm
function evaluate(x::acb_poly, b::Array{acb, 1})
   return evaluate_iter(x, b)
end

###############################################################################
#
#   Root finding
#
###############################################################################

doc"""
    roots(x::acb_poly; target=0, isolate_real=false, initial_prec=0, max_prec=0, max_iter=0)
> Attempts to isolate the complex roots of the complex polynomial $x$ by
> iteratively refining balls in which they lie.
>
> This is done by increasing the working precision, starting at `initial_prec`.
> The maximal number of iterations can be set using `max_iter` and the maximal
> precision can be set using `max_prec`.
>
> If `isolate_real` is set and $x$ is strictly real, then the real roots will
> be isolated from the non-real roots. Every root will have either zero,
> positive or negative real part.
>
> It is assumed that $x$ is squarefree.
"""
function roots(x::acb_poly; target=0, isolate_real=false, initial_prec=0, max_prec=0, max_iter=0)
    deg = degree(x)
    if deg <= 0
        return Array{acb}(0)
    end

    initial_prec = (initial_prec >= 2) ? initial_prec : 32
    max_prec = (max_prec >= 2) ? max_prec : 3 * prec(parent(x))

    isolated = 0
    wp = initial_prec
    roots = acb_vec(deg)

    while true
        in_roots = (wp == initial_prec) ? C_NULL : roots
        step_max_iter = (max_iter >= 1) ? max_iter : min(max(deg, 32), wp)
        isolated = ccall((:acb_poly_find_roots, :libarb), Int,
            (Ref{acb_struct}, Ref{acb_poly}, Ref{acb_struct}, Int, Int),
                roots, x, in_roots, step_max_iter, wp)

        wp = wp * 2

        if isolated == deg
            ok = true
            if target > 0
                for i = 0 : deg-1
                    re = ccall((:acb_real_ptr, :libarb), Ref{arb_struct},
                        (Ref{acb}, ), roots + i * sizeof(acb_struct))
                    im = ccall((:acb_imag_ptr, :libarb), Ref{arb_struct},
                        (Ref{acb}, ), roots + i * sizeof(acb_struct))
                    t = ccall((:arb_rad_ptr, :libarb), Ref{mag_struct}, (Ref{arb}, ), re)
                    u = ccall((:arb_rad_ptr, :libarb), Ref{mag_struct}, (Ref{arb}, ), im)
                    ok = ok && (ccall((:mag_cmp_2exp_si, :libarb), Cint,
                        (Ref{mag_struct}, Int), t, -target) <= 0)
                    ok = ok && (ccall((:mag_cmp_2exp_si, :libarb), Cint,
                        (Ref{mag_struct}, Int), u, -target) <= 0)
                end
            end

            if isreal(x)
                real_ok = ccall((:acb_poly_validate_real_roots, :libarb),
                    Bool, (Ref{acb_struct}, Ref{acb_poly}, Int), roots, x, wp)

                if isolate_real && !real_ok
                    ok = false
                end

                if real_ok
                    for i = 0 : deg - 1
                        im = ccall((:acb_imag_ptr, :libarb), Ref{arb_struct},
                            (Ref{acb}, ), roots + i * sizeof(acb_struct))
                        if ccall((:arb_contains_zero, :libarb), Bool, (Ref{arb_struct}, ), im)
                            ccall((:arb_zero, :libarb), Void, (Ref{arb_struct}, ), im)
                        end
                    end
                end
            end

            if ok
                break
            end
        end

        if wp > max_prec
            break
        end
    end

    if isolated == deg
        ccall((:_acb_vec_sort_pretty, :libarb), Void,
            (Ref{acb_struct}, Int), roots, deg)
        res = array(base_ring(parent(x)), roots, deg)
    end

    acb_vec_clear(roots, deg)

    if isolated == deg
        return res
    else
        error("unable to isolate all roots (insufficient precision, or there is a multiple root)")
    end
end

###############################################################################
#
#   Root bounds
#
###############################################################################

doc"""
    roots_upper_bound(f::acb_poly) -> arb

> Returns an upper bound for the absolute value of all complex roots of $f$.
"""
function roots_upper_bound(x::acb_poly)
   z = ArbField(prec(base_ring(x)))()
   p = prec(base_ring(x))
   t = ccall((:arb_rad_ptr, :libarb), Ref{mag_struct}, (Ref{arb}, ), z)
   ccall((:acb_poly_root_bound_fujiwara, :libarb), Void,
         (Ref{mag_struct}, Ref{acb_poly}), t, x)
   s = ccall((:arb_mid_ptr, :libarb), Ref{arf_struct}, (Ref{arb}, ), z)
   ccall((:arf_set_mag, :libarb), Void, (Ref{arf_struct}, Ref{mag_struct}), s, t)
   ccall((:arf_set_round, :libarb), Void,
         (Ref{arf_struct}, Ref{arf_struct}, Int, Cint), s, s, p, ARB_RND_CEIL)
   ccall((:mag_zero, :libarb), Void, (Ref{mag_struct},), t)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::acb_poly)
   ccall((:acb_poly_zero, :libarb), Void, (Ref{acb_poly},), z)
   return z
end

function fit!(z::acb_poly, n::Int)
   ccall((:acb_poly_fit_length, :libarb), Void,
                    (Ref{acb_poly}, Int), z, n)
   return nothing
end

function setcoeff!(z::acb_poly, n::Int, x::fmpz)
   ccall((:acb_poly_set_coeff_fmpz, :libarb), Void,
                    (Ref{acb_poly}, Int, Ref{fmpz}), z, n, x)
   return z
end

function setcoeff!(z::acb_poly, n::Int, x::acb)
   ccall((:acb_poly_set_coeff_acb, :libarb), Void,
                    (Ref{acb_poly}, Int, Ref{acb}), z, n, x)
   return z
end

function mul!(z::acb_poly, x::acb_poly, y::acb_poly)
   ccall((:acb_poly_mul, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
                    z, x, y, prec(parent(z)))
   return z
end

function addeq!(z::acb_poly, x::acb_poly)
   ccall((:acb_poly_add, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
                    z, z, x, prec(parent(z)))
   return z
end

function add!(z::acb_poly, x::acb_poly, y::acb_poly)
   ccall((:acb_poly_add, :libarb), Void,
                (Ref{acb_poly}, Ref{acb_poly}, Ref{acb_poly}, Int),
                    z, x, y, prec(parent(z)))
   return z
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{acb_poly}, ::Type{Float64}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{Complex{Float64}}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{T}) where {T <: Integer} = acb_poly

promote_rule(::Type{acb_poly}, ::Type{Rational{T}}) where T <: Union{Int, BigInt} = acb_poly

promote_rule(::Type{acb_poly}, ::Type{Complex{Int}}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{fmpz}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{fmpq}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{arb}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{acb}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{fmpz_poly}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{fmpq_poly}) = acb_poly

promote_rule(::Type{acb_poly}, ::Type{arb_poly}) = acb_poly

################################################################################
#
#  Parent object call overloads
#
################################################################################

function (a::AcbPolyRing)()
   z = acb_poly()
   z.parent = a
   return z
end

for T in [Integer, fmpz, fmpq, Float64, Complex{Float64},
          Complex{Int}, arb, acb]
  @eval begin
    function (a::AcbPolyRing)(b::$T)
      z = acb_poly(base_ring(a)(b), a.base_ring.prec)
      z.parent = a
      return z
    end
  end
end

(a::AcbPolyRing)(b::Rational{T}) where {T <: Integer} = a(fmpq(b))

function (a::AcbPolyRing)(b::Array{acb, 1})
   z = acb_poly(b, a.base_ring.prec)
   z.parent = a
   return z
end

for T in [fmpz, fmpq, Float64, Complex{Float64}, Complex{Int}, arb]
  @eval begin
    (a::AcbPolyRing)(b::Array{$T, 1}) = a(map(base_ring(a), b))
  end
end

(a::AcbPolyRing)(b::Array{T, 1}) where {T <: Integer} = a(map(base_ring(a), b))

(a::AcbPolyRing)(b::Array{Rational{T}, 1}) where {T <: Integer} = a(map(base_ring(a), b))

function (a::AcbPolyRing)(b::fmpz_poly)
   z = acb_poly(b, a.base_ring.prec)
   z.parent = a
   return z
end

function (a::AcbPolyRing)(b::fmpq_poly)
   z = acb_poly(b, a.base_ring.prec)
   z.parent = a
   return z
end

function (a::AcbPolyRing)(b::arb_poly)
   z = acb_poly(b, a.base_ring.prec)
   z.parent = a
   return z
end

function (a::AcbPolyRing)(b::acb_poly)
   z = acb_poly(b, a.base_ring.prec)
   z.parent = a
   return z
end

################################################################################
#
#  PolynomialRing constructor
#
################################################################################

function PolynomialRing(R::AcbField, s::AbstractString; cached = true)
  S = Symbol(s)
  parent_obj = AcbPolyRing(R, S, cached)
  return parent_obj, parent_obj(fmpz_poly([fmpz(0), fmpz(1)]))
end


text:
#
# A vector of composite data stored in a column-oriented format.
#

abstract AbstractCompositeVector{R} <: AbstractVector{R}

@inline column(cv::AbstractCompositeVector, k::Int) =
    columns(cv)[k]

@inline offsets(cv::AbstractCompositeVector, k::Int) =
    offsets(column(cv, k))

@inline values(cv::AbstractCompositeVector, k::Int) =
    values(column(cv, k))

guesslength(cols::Columns) =
    !isempty(cols) ? length(cols[1]) : 0

guesstype(cols::Columns) =
    Tuple{map(datatype, cols)...}

guesstype(::Type{Any}, cols::Columns) =
    guesstype(cols)

guesstype{T}(::Type{T}, cols::Columns) = T

linearindexing{T<:AbstractCompositeVector}(::Type{T}) = Base.LinearFast()

Base.array_eltype_show_how(::AbstractCompositeVector) = (true, "")
Base.array_eltype_show_how{T,N,P<:AbstractCompositeVector,I,L}(::SubArray{T,N,P,I,L}) = (true, "")

Base.summary(cv::AbstractCompositeVector) =
    "$(length(cv))-element composite vector of $(typesummary(cv))"

function typesummary(cv::AbstractCompositeVector)
    names = fieldnames(eltype(cv))
    cols = columns(cv)
    fields =
        if eltype(names) == Symbol && length(names) == length(cols)
            map((n, c) -> "$n::$(typesummary(c))", names, cols)
        else
            map(typesummary, cols)
        end
    return string("{", join(fields, ", "), "}")
end

# Composite vector with an open signature.

immutable TupleVector{R,C<:ColumnTuple} <: AbstractCompositeVector{R}
    len::Int
    cols::C

    function TupleVector(len::Int, cols::C)
        @boundscheck for col in cols
            length(col) == len || error("unexpected column length")
        end
        new(len, cols)
    end
end

typealias RecordVector{R<:AbstractRecord} TupleVector{R}

TupleVector(len::Int, cols::ColumnTuple) =
    TupleVector{guesstype(cols),typeof(cols)}(len, cols)

TupleVector(len::Int, cols::Tuple{AbstractVector, Vararg{AbstractVector}}) =
    TupleVector(len, map(Column, cols))

TupleVector(cols::ColumnTuple) =
    TupleVector(guesslength(cols), cols)

TupleVector(cols::Tuple{AbstractVector, Vararg{AbstractVector}}) =
    TupleVector(map(Column, cols))

TupleVector(len::Int, cols::AbstractVector...) =
    TupleVector(len, map(Column, cols))

TupleVector(cols::AbstractVector...) =
    TupleVector(map(Column, cols))

function RecordVector(len::Int, fcols::Vector{Tuple{Symbol,Column}})
    nfields = length(fcols)
    fieldnames = Vector{Symbol}(nfields)
    fieldcols = Vector{Column}(nfields)
    fieldtypes = Vector{Type}(nfields)
    seen = Set{Symbol}()
    for (k, (field, col)) in enumerate(fcols)
        if field == Symbol("")
            field = Symbol("_", k)
        end
        if field in seen
            field = Symbol(field, "′")
        end
        fieldnames[k] = field
        fieldcols[k] = col
        fieldtypes[k] = datatype(col)
        push!(seen, field)
    end
    R = recordtype(fieldnames){fieldtypes...}
    cols = (fieldcols...)
    return TupleVector{R,typeof(cols)}(len, cols)
end

RecordVector(fcols::Vector{Tuple{Symbol,Column}}) =
    RecordVector(!isempty(fcols) ? length(fcols[1][2]) : 0, fcols)

RecordVector(len::Int; kwargs...) =
    RecordVector(len, collect(Tuple{Symbol,Column}, kwargs))

RecordVector(; kwargs...) =
    RecordVector(collect(Tuple{Symbol,Column}, kwargs))

@inline columns(tv::TupleVector) = tv.cols

@inline size(tv::TupleVector) = (tv.len,)

@inline length(tv::TupleVector) = tv.len

@inline getindex{R}(tv::TupleVector{R}, i::Int) =
    getdata(R, tv.cols, i)

@inline getindex{R,C}(tv::TupleVector{R,C}, idxs::AbstractVector{Int}) =
    TupleVector{R,C}(length(idxs), map(col -> col[idx], tv.cols))

# Composite vector with a hidden signature.

immutable DataVector <: AbstractCompositeVector{Any}
    meta::Any
    len::Int
    cols::ColumnVector

    function DataVector(meta, len::Int, cols::ColumnVector)
        @boundscheck for col in cols
            length(col) == len || error("unexpected column length")
        end
        new(meta, len, cols)
    end
end

DataVector(len::Int, cols::ColumnVector) =
    DataVector(Any, len, cols)

DataVector(cols::ColumnVector) =
    DataVector(Any, guesslength(cols), cols)

DataVector(len::Int, cols::AbstractVector...) =
    DataVector(len, collect(Column, cols))

DataVector(cols::AbstractVector...) =
    DataVector(collect(Column, cols))

convert(::Type{TupleVector}, dv::DataVector) =
    let cols = (dv.cols...), R = guesstype(dv.meta, dv.cols)
        TupleVector{R,typeof(cols)}(dv.len, cols)
    end

@inline columns(dv::DataVector) =
    dv.cols

@inline eltype(dv::DataVector) =
    guesstype(dv.meta, dv.cols)

@inline size(dv::DataVector) = (dv.len,)

@inline length(dv::DataVector) = dv.len

getindex(dv::DataVector, i::Int) =
    getdata(dv.meta, dv.cols, i)

getindex(dv::DataVector, idxs::AbstractVector{Int}) =
    DataVector(dv.meta, length(idxs), Column[col[idxs] for col in dv.cols])

Base.mapfoldl(f, op, v0, dv::DataVector) =
    mapfoldl(f, op, v0, convert(TupleVector, dv))

Base.mapfoldl(f, op, dv::DataVector) =
    mapfoldl(f, op, convert(TupleVector, dv))

Base.mapfoldr(f, op, v0, dv::DataVector) =
    mapfoldr(f, op, v0, convert(TupleVector, dv))

Base.mapfoldr(f, op, dv::DataVector) =
    mapfoldr(f, op, convert(TupleVector, dv))

Base.mapreduce(f, op, v0, dv::DataVector) =
    mapreduce(f, op, v0, convert(TupleVector, dv))

Base.mapreduce(f, op, dv::DataVector) =
    mapreduce(f, op, convert(TupleVector, dv))

Base.foreach(f, dv::DataVector) =
    foreach(f, convert(TupleVector, dv))

Base.map(f, dv::DataVector) =
    map(f, convert(TupleVector, dv))

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    LoggerSerializationError(logger::Logger)

Provides a helpful error message when a logger cannot be serialized.
"""
struct LoggerSerializationError <: Exception
    logger::Logger
end

function Base.showerror(io::IO, e::LoggerSerializationError)
    print(io, string(
        "$(typeof(e)): $(e.logger) was unable to be serialized. ",
        "Perhaps try not using shared loggers between processes?",
    ))
end

"""
    EscalationError(msg)

An error type for log records that have been escalated to errors.
"""
struct EscalationError <: Exception
    msg
end


text:
fname = "winequality.csv"

using DelimitedFiles
data = DelimitedFiles.readdlm(fname, ';')
# 1600x12 Array{Any,2}:
#    "fixed acidity"   "volatile acidity"  …    "alcohol"   "quality"
#   7.4               0.7                      9.4         5.0
#   7.8               0.88                     9.8         5.0
#   7.8               0.76                     9.8         5.0
#  11.2               0.28                     9.8         6.0

# variant:
data3 = DelimitedFiles.readdlm(fname, ';', Float64, '\n', header=true)
#> ([7.4 0.7 … 9.4 5.0; 7.8 0.88 … 9.8 5.0; … ; 5.9 0.645 … 10.2 5.0;
# 6.0 0.31 … 11.0 6.0],
# AbstractString["fixed acidity" "volatile acidity" … "alcohol" "quality"])

# extracting data:
row3 = data[3, :] # measurements for wine xyz
# 12-element Array{Any,1}: 7.8  0.88  0  2.6  0.098 25 67  0.9968  3.2  0.68  9.8  5
col3 = data[ :, 3]  # measurements of citric acid
# 1600-element Array{Any,1}:
# "citric acid"
# 0.0 0.0 0.04 0.56 0.0 0.0 ⋮ 0.08 0.08 0.1 0.13 0.12 0.47 ...
x = data[:, 2:4]
# 1600×3 Array{Any,2}:
#   "volatile acidity"   "citric acid"   "residual sugar"
#  0.7                  0               1.9
#  0.88                 0               2.6
#  0.76                 0.04            2.3
y = data[70:75, 2:4]
# 6x3 Array{Any,2}:
 # 0.32   0.57  2.0
 # 0.705  0.05  1.9
 # 0.63   0.08  1.9
 # 0.67   0.23  2.1
 # 0.69   0.22  1.9
 # 0.675  0.26  2.1
 z = [data[:,3] data[:,6] data[:,11]]
 # 1600x3 Array{Any,2}:
 #  "citric acid"    "free sulfur dioxide"    "alcohol"
 # 0.0             11.0                      9.4
 # 0.0             25.0                      9.8
 # 0.04            15.0                      9.8
 # 0.56            17.0                      9.8
 # ⋮
 # 0.08            32.0                     10.5
 # 0.1             39.0                     11.2
 # 0.13            29.0                     11.0
 # 0.12            32.0                     10.2
 # 0.47            18.0                     11.0

 # without headers:
 z = [data[2:end,3] data[2:end,6] data[2:end,11]]

 mutable struct Wine
    fixed_acidity::Array{Float64}
    volatile_acidity::Array{Float64}
    citric_acid::Array{Float64}
    #...define the other columns here
    quality::Array{Float64}
end
# wine1 = Wine(7.4, 0.7, 0, ..., 5)
# wine1 = Wine(data[1, :]...) # ... = splat-operator

data = z
# writing to csv file:
writedlm("partial.dat", data, ';')


text:
type masks_NaN_2dim <: masks_Forward_Backward

  Dim1 :: Array{Int64,1}
  Dim2 :: Array{Int64,1}

  mask :: Array{Bool}
  maskp :: Array{Bool}

  rD1 :: Array{Int64,1}
  mask2 :: Array{Bool,2}
end

#= ForwardD2(xar, Dim1 ) = ForwardD2(xar, [Dim1...;]); =#
function ForwardD2(xar :: Array{Float64}, Dim1 :: Array{Int64})
  mask = isnan(xar)
  D1   = Dim1
  D2   = filter(x->!in(x,D1),[1:ndims(xar);])
  xp   = permutedims(xar,[D1,D2;])
  maskp = isnan(xp)

  x2 = reshape(xp, reduce(*,size(xar,D1...)),reduce(*,size(xar,D2...)))
  rD1 = find(sum(isnan(x2),2).==0)
  mask2 = isnan(x2)
  M    = masks_NaN_2dim(D1,D2,mask,maskp,rD1,mask2)
  return (M , x2[rD1,:])
end

function BackwardD2( r2::Array{Float64,2} , M :: masks_NaN_2dim;)
  x2 = fill(NaN,size(M.mask2))
  x2[M.rD1,:] = r2
  Xp = reshape(x2,size(M.mask,[M.Dim1, M.Dim2]...))
  #= Xar = permutedims(Xp, indexin([1:ndims(Xp)],[M.Dim1,M.Dim2])) =#
  #= Xar = permutedims(Xp, indexin([1:ndims(Xp)],[M.Dim1,M.Dim2])) =#
  Xar = ipermutedims(Xp, [M.Dim1,M.Dim2])
  return Xar
end

function LBackwardD2( r2::Array{Float64,1} , M :: masks_NaN_2dim;)
  x2 = fill(NaN,size(M.mask2,1))
  x2[M.rD1] = r2
  Xp = reshape(x2,size(M.mask,[M.Dim1]...))
  return Xp
end

function LBackwardD2( r2::Array{Float64,2} , M :: masks_NaN_2dim;)
  x2 = fill(NaN,size(M.mask2,1),size(r2,2))
  x2[M.rD1,:] = r2

  Xp = reshape(x2,[size(M.mask,M.Dim1...)...,size(r2,2);]...)
  return Xp
end
export ForwardD2,BackwardD2,LBackwardD2;

function DFevaluate(EE :: _EOF_Functor,xarr :: Array{Float64},
                    Dim1 :: Array{Int64,1} ;nsv  :: Int64 = 3)
  (M,real1D) = ForwardD2(xarr, Dim1)
  (S,V,D,ratio) = Fevaluate(EE, real1D,nsv = nsv)
  S2 = LBackwardD2(S,M)
  return (S2,V,D,ratio)
end

function DFevaluate(SS :: _SVD_Functor,
                     xarr :: Array{Float64}, yarr :: Array{Float64},
                     Dim1x :: Array{Int64,1}, Dim1y :: Array{Int64,1} = Dim1x
                     ;nsv :: Int64 = 3)

  (Mx, rx1) = ForwardD2(xarr,Dim1x)
  (My, ry1) = ForwardD2(yarr,Dim1y)
  (Sx, Sy, V, ratio, PC1x, PC1y) = Fevaluate(SS, rx1,ry1, nsv= nsv)
  Sx = LBackwardD2(Sx,Mx)
  Sy = LBackwardD2(Sy,My)
  Pctx = reshape(PC1x,size(Mx.mask,Mx.Dim2...))
  Pcty = reshape(PC1y,size(My.mask,My.Dim2...))
  return (Sx,Sy,ratio,V,Pctx,Pcty)
end

function DFevaluate( MM :: _Mean_Functor , xarr :: Array{Float64}, Dim1 :: Array{Int64,1})
  xa = copy(xarr)
  ref = fill(1.0,size(xa)); ref[isnan(xa)] = 0.0
  xa[isnan(xa)] = 0.0
  return mean(xa,Dim1) ./ mean(ref, Dim1)
end

export masks_NaN_Forward_2dim
export Fevaluate
export DFevaluate


text:
# # DrWatson Workflow Tutorial

# ```@setup workflow
# cd(@__DIR__)
# ```

# *Disclaimer: DrWatson assumes basic knowledge of how Julia's
# [project manager](https://pkgdocs.julialang.org/v1/environments/) works.*

# This example page demonstrates how DrWatson's functions help a typical scientific
# workflow, as illustrated below:

# ![](workflow.png)

# Blue text comes from DrWatson. Of course, not all of DrWatson's functionality
# will be highlighted in this tutorial nor is shown in the above figure!

# ## 1. Setup the project

# So, let's start a new scientific project. You want your project to be contained
# in a folder. So let's create a new project, located at current working directory
using DrWatson
initialize_project("DrWatson Example"; authors="Datseris", force=true)

# Alright now we have a project set up. The project has a default reasonable structure,
# as illustrated in the [Default Project Setup](@ref) page:
#
# ```@setup workflow
# using DrWatson
# struct ShowFile
#     file::String
# end
# function Base.show(io::IO, ::MIME"text/plain", f::ShowFile)
#     write(io, read(f.file))
# end
# ```
# ```@example workflow
# ShowFile(dirname(pathof(DrWatson))*"/defaults/project_structure.txt") # hide
# ```

# For example, folders exist for data, plots, scripts, source code, etc.
# Three files are noteworthy:
# * Project.toml: Defines project
# * Manifest.toml: Contains exact list of project dependencies
# * .git (hidden folder): Contains reversible and searchable history of the project

# The scientific project we have created is also a [Julia project environment](https://docs.julialang.org/en/v1/manual/code-loading/#Environments-1).
# This means that it has its own dedicated dependencies and versions of dependencies.
# This project is now active by default so we can start adding packages
# that we will be using in the project. We'll add the following for demonstrating
using Pkg
Pkg.add(["Statistics", "JLD2"])

# ## 2. Write some scripts

# We start by writing some script for our project that will do some dummy calculations.
# Let's create `scripts/example.jl` in our project. All following code is supposed to
# exist in that file.

# ```@setup workflow
# cd(joinpath(@__DIR__, "DrWatson Example"))
# ```

# Now, with DrWatson every script (typically) starts with the following two lines:
# ```@setup workflow
# quickactivate("DrWatson Example", "DrWatson Example")
# ```

# ```julia
# using DrWatson
# @quickactivate "DrWatson Example" # <- project name
# ```
# This command does something simple: it searches the folder of the script, and its
# parent folders, until it finds a Project.toml. It activates that project, but
# if the project name doesn't match the given name (here `"DrWatson Example"`)
# it throws an error. Let's see the project we activated:
projectname()

# This is **extremely useful** for two reasons. First, it is guaranteed
# that our scripts run within the context of the project and thus use the correct
# package versions.

# Second, DrWatson provides the powerful function [`projectdir`](@ref) and its derivatives
# like `datadir, plotsdir, srcdir`, etc.
projectdir()

# `projectdir` will **always** return the path to the currently active project. It doesn't
# matter where its called from, or where the active project actually is.
# So, by using DrWatson, you don't care anymore where your current script is, you only
# care about the target directory.

datadir()

#

datadir("sims", "electron_gas")

# Giving arguments to `projectdir` and derivatives joins paths.

# ## 3. Prepare simulations

# Let's say we write a simple simulation function, that creates some data
function fakesim(a, b, v, method = "linear")
    if method == "linear"
        r = @. a + b * v
    elseif method == "cubic"
        r = @. a*b*v^3
    end
    y = sqrt(b)
    return r, y
end

# and we create some parameters in our scripts and run the simulation
a, b = 2, 3
v = rand(5)
method = "linear"
r, y = fakesim(a, b, v, method)

# Okay, that is fine, but it is typically the case that in scientific context
# some simulations are done for several different combinations of parameters.
# It is convenient to group all parameters in a dictionary, with the keys
# being the parameters. Depending on the package you will use to actually save the data,
# the key type should be either `String` or `Symbol`. Here we will be using
# JLD2.jl and therefore the key type will be `String`.
params = @strdict a b v method

# Above we used the ultra-cool [`@strdict`](@ref) macro, which creates a
# dictionary from existing variables!

# Now, for every simulation we want to do, we would create such a container.
# We can use the [`dict_list`](@ref) to ease up the process of preparing several
# of these parameter containers

allparams = Dict(
    "a" => [1, 2], # it is inside vector. It is expanded.
    "b" => [3, 4],
    "v" => [rand(5)],     # single element inside vector; no expansion
    "method" => "linear", # not in vector = not expanded, even if naturally iterable
)

dicts = dict_list(allparams)

# using `dict_list` is great, because it has a very clear design on how to expand
# containers, while not caring whether the parameter values are iterable or not.
# In short **everything in a `Vector` is expanded once** (`Vector`s of length 1
# are not expanded naturally). See [`dict_list`](@ref) for more details.

# The resulting dictionaries are then typically given into a `main` or `makesim`
# function that actually does the simulation given some input parameters.

# ## 4. Run and save
# Alright, we now have to actually save the results, so we first define:

function makesim(d::Dict)
    @unpack a, b, v, method = d
    r, y = fakesim(a, b, v, method)
    fulld = copy(d)
    fulld["r"] = r
    fulld["y"] = y
    return fulld
end

# and then we can save the results by once again leveraging [`projectdir`](@ref)
# ```julia
# for (i, d) in enumerate(dicts)
#     f = makesim(d)
#     wsave(datadir("simulations", "sim_$(i).jld2"), f)
# end
# ```

# *(`wsave` is a function from DrWatson, that ensures that the directory you try to
# save the data exists. It then calls `FileIO.save`)*

# Here each simulation was named according to a number.
# But this is not how we do it in science... We typically want the input parameters
# to be part of the file name. E.g. here we would want the file name to be something like
# `a=2_b=3_method=linear.jld2`. It would be also nice that such a naming scheme would
# apply to arbitrary input parameters so that we don't have to manually write
# `a=$(a)_b=$(b)_method=$(method)` and change this code every time we change
# a parameter name...

# Enter [`savename`](@ref):

savename(params)

# `savename` takes as an input pretty much *any* Julia composite container with key-value
# pairs and transforms it into such a name. We can even do

savename(dicts[1], "jld2")

# `savename` is flexible and smart. As you noticed, even though the vector `v` with
# 5 numbers is part of the input, it wasn't included in the name (on purpose).
# See the [`savename`](@ref) documentation for more.

# We now transform our make+save loop into

for (i, d) in enumerate(dicts)
    f = makesim(d)
    wsave(datadir("simulations", savename(d, "jld2")), f)
end

readdir(datadir("simulations"))

# That is cool, but we can do better. In fact, **much better**.
# Remember that the project initialized by DrWatson is a git repository.
# So, now we quickly go into git and commit the script we have created and all changes
# (not shown here).

# Then we make+save again, but now instead of `wsave` we use [`@tagsave`](@ref):

# ```@setup workflow
# for (i, d) in enumerate(dicts)
#     f = makesim(d)
#     @tagsave(datadir("simulations", savename(d, "jld2")), f; gitpath = "../..")
# end
# ```

# ```julia
# for (i, d) in enumerate(dicts)
#     f = makesim(d)
#     @tagsave(datadir("simulations", savename(d, "jld2")), f)
# end
# ```

# and let's load the first simulation

firstsim = readdir(datadir("simulations"))[1]

wload(datadir("simulations", firstsim))

# So what happened is that `tagsave` automatically added git-related information
# into the file we saved (the field `:gitcommit`), enabling reproducibility!

# It gets even better! Because [`@tagsave`](@ref) is a macro, it deduced automatically
# where the script that called `@tagsave` was located. It even includes the exact line of code
# that called the `@tagsave` command. This information is in the `:script` field of the
# saved data!

# ## 5. Analyze results

# Cool, now we can start analyzing some simulations. The actual analysis is your job,
# but DrWatson can help you get started with the [`collect_results`](@ref) function.
# Notice that you need to be `using DataFrames` to access the function!
Pkg.add(["DataFrames"])
using DataFrames

df = collect_results(datadir("simulations"))

# Some things to note:
# * the returned object is a `DataFrame` for further analysis.
# * the input to `collect_result` is a folder, **not** a dataframe!
#   The function does the loading and combining for you.
# * If you create new simulations, you can iteratively (or all from scratch) add them to
#   this dataframe.
# * If you create new simulations that have **new parameters**, that don't exist in the
#   simulations already saved, that's no problem. `collect_results` will appropriately
#   and automatically add `missing` to all parameter values that don't exist in previews
#   and/or current simulations. This is demonstrated explicitly in the
#   [Adapting to new data/parameters](@ref) real world example, so it is not repeated here.
# * Similarly with e.g. `savename`, `collect_results` is a flexible function. It has
#   several configuration options.

# Great! so now we are doing some analysis and we want to save some results...
# It is very often the case in science that the same analysis may be done again, and again,
# and some times even with the same parameters... And poor scientitists sometimes forget
# to change the name of the output file, and end up overwritting previous work! Devastating!

# To avoid such scenarios, we can use the function [`safesave`](@ref), e.g.

analysis = 42

safesave(datadir("ana", "linear.jld2"), @strdict analysis)

# If a file `linear.jld2` exists in that folder, it is not overwritten. Instead, it is
# renamed to `linear#1.jld2`, and a new `linear.jld2` file is made!

# ## 6. Share your project

# This is already discussed in the [Reproducibility](@ref) section of the docs
# so there is no reason to copy/paste everything here.
# What is demonstrated there is that it is truly trivial to share your project
# with a colleague, and this project is guaranteed to work for them!

# ---

# And that's it! We hope that DrWatson will take some stress out of the absurdly stressfull
# scientific life!


text:
using Test

@testset "ReinforcementLearningZoo.jl" begin
    # include("cfr/cfr.jl")
end


text:
## Excel

@info "XLSX.jl is available - including functionality to read / write Excel (xlsx) files"

using .XLSX

function read_table(::TableIOInterface.ExcelFormat, filename:: AbstractString, sheetname:: AbstractString; kwargs...)
    xf = XLSX.readxlsx(filename)
    try
        sheet = xf[sheetname]
        return XLSX.eachtablerow(sheet)
    finally
        close(xf)
    end
end

function read_table(t::TableIOInterface.ExcelFormat, filename:: AbstractString; kwargs...)
    xf = XLSX.readxlsx(filename)
    try
        table_list = _list_tables(t, xf)
        length(table_list) > 1 && @warn "File contains more than one table, the alphabetically first one is taken"
        sheet = xf[first(table_list)]
        return DataFrame(XLSX.eachtablerow(sheet); copycols=false) # this would be no valid Table.jl output if not converted to DataFrame
    finally
        close(xf)
    end
end

function write_table!(::TableIOInterface.ExcelFormat, filename:: AbstractString, tablename:: AbstractString, table:: DataFrame; kwargs...)
    _checktable(table)
    _checktablename(tablename)
    XLSX.writetable(filename, table; overwrite=true, sheetname=tablename, kwargs...)
    nothing
end

# XLSX supports only DataFrames, not arbitrary Tables.jl inputs. For export, the table is converted to a DataFrame first.
write_table!(::TableIOInterface.ExcelFormat, filename:: AbstractString, tablename:: AbstractString, table; kwargs...) = write_table!(TableIOInterface.ExcelFormat(), filename, tablename, DataFrame(table); kwargs...)

function list_tables(t::TableIOInterface.ExcelFormat, filename:: AbstractString)
    xf = XLSX.readxlsx(filename)
    try
        return _list_tables(t, xf)
    finally
        close(xf)
    end
end

_list_tables(::TableIOInterface.ExcelFormat, xf) = XLSX.sheetnames(xf) |> sort


text:
######################################################
#
# MSM Plot
#
# Date:   2014-12-02
# Author: Hessel Tuinhof
#
# Simulate and plot a Markov-Switching Multifractal (MSM)
#
# Required packages: 	PyPlot
# Required functions:	MSM.jl
#
######################################################

# load packages & functions
using PyPlot
require("../src/msm/MSM.jl")

# define colors
rouge = (0.79, 0.0, 0.1)

# configure parameters
num_sim = 1
sim_length = 5000
kbar = 8
m0 = 1.4
sigma = 0.5
gamma_kbar = 0.95
b = 3

# set seed
srand(314159265358979323846)

# plot MSM (returns & log price)
figure()
msm = MSM(num_sim, sim_length, m0, kbar, sigma, gamma_kbar, b)
plot(msm, linewidth = 0.5, color = rouge)
xlim(xmax = sim_length)
figure()
plot(cumsum(msm), linewidth = 1, color = rouge)
xlim(xmax = sim_length)

max_stars_repo_name: hesseltuinhof/masterthesis

max_issues_repo_name: hesseltuinhof/masterthesis

max_forks_repo_name: hesseltuinhof/masterthesis


text:
mutable struct Rule
    color::String
    containedBags::Dict{String,Int}
    Rule() = new("",Dict())
end

function parseLine(line)
    rule = Rule()
    containedBags = Dict()
    sections = split(line, "contain")
    colorsection = split(sections[1], ' ')
    rule.color = colorsection[1]*colorsection[2]
    rulesections = split(sections[2], ',')
    if rulesections[1] != " no other bags."
        for rulesection in rulesections
            words = split(rulesection, ' ')
            amount = parse(Int,words[2])
            rulecolor = words[3]*words[4]
            containedBags[rulecolor] = amount
        end
    end
    rule.containedBags = containedBags
    return rule

end

allRules = []
open("data.txt") do file
    for line in eachline(file)
        push!(allRules, parseLine(line))
    end
end

function countGoldenBags()
    goldRelatedColors = Set()
    currentColors = ["shinygold"]
    while length(currentColors) != 0
        newColors::Array{String} = []
        for currColor in currentColors
            for rule in allRules
                if haskey(rule.containedBags, currColor)
                    push!(newColors, rule.color)
                end
            end
        end
        println(newColors)
        for newColor in newColors
            push!(goldRelatedColors, newColor)
        end
        currentColors = newColors
    end
    println(length(goldRelatedColors))
end

countGoldenBags()


text:
@field Links = FieldTraits.ComposedDict()

"""
Composable, that allows one to link fields and register callbacks to field changes
"""
macro reactivecomposed(expr)
    composed_type(expr, [Links], ReactiveComposable)
end

immutable _NoUpdate end
const NoUpdate = _NoUpdate()

@propagate_inbounds function setindex!{F <: Field}(ct::ReactiveComposable, value, field::Type{F}, ::_NoUpdate)
    _setindex!(ct, value, field)
end
@propagate_inbounds function setindex!{F <: Field}(ct::ReactiveComposable, value, field::Type{F})
    _setindex!(ct, value, field)
    links = ct[Links]
    if haskey(links, field)
        link = links[field]
        for (func, fields, args) in link
            func(map(f-> ct[f], fields)..., args...)
        end
    end
    value
end

"""
Links a field with another in to Composable types.
After this operation, setindex!(a, val) will result in
setindex!(b, val) being executed.
link!(Scale, a => b)
"""
function link!{F <: Field}(::Type{F}, pair::Pair{<: ReactiveComposable, <: Any})
    a, b = pair
    on(a, F) do val
        b[F] = val
    end
end

function on(F, object::ReactiveComposable, head, tail...)
    links = object[Links]
    args = (head, tail...)
    _fields = []; _args = []
    field, state = first(args), start(args)
    nomorefields = false
    for elem in args
        if !nomorefields && isa(elem, Type) && elem <: Field
            push!(_fields, elem)
        else
            nomorefields = true
            push!(_args, elem)
        end
    end
    args = (_args...,)
    fields = (_fields...,)
    for field in fields
        fieldlinks = get!(links, field, [])
        # adds a callback to the field
        push!(fieldlinks, (F, fields, args))
    end
end

function convertfor(::Type{Links}, x::Partial, val)
    # never copy links
    default(Links, x)
end


text:
# XXX: additional backedges for JET analysis
# currently invalidation of JET cache just relies on invalidation of native code cache, i.e.
# it happens only when there is a backedge from erroneous frame, which is not always true
# e.g. native interpreter doesn't add backedge when the return type is `Any`, but invalidation
# may still change the result of JET analysis
# so we may need to add additional backedges for frames from which some errors are reported

# we just overload `cache_result!` and so we don't need to set up our own cache
CC.code_cache(interp::JETInterpreter) = code_cache(interp.native)

# TODO better to use `overload_cache_result!!` hack ?
function CC.cache_result!(interp::JETInterpreter, result::InferenceResult, valid_worlds::WorldRange)
    linfo = result.linfo

    # check if the existing linfo metadata is also sufficient to describe the current inference result
    # to decide if it is worth caching this
    already_inferred = CC.already_inferred_quick_test(interp, linfo) # always false
    if !already_inferred && CC.haskey(WorldView(code_cache(interp), valid_worlds), linfo)
        already_inferred = true
    end

    # TODO: also don't store inferred code if we've previously decided to interpret this function
    if !already_inferred
        inferred_result = CC.transform_result_for_cache(interp, linfo, result.src)
        CC.setindex!(code_cache(interp), CodeInstance(result, inferred_result, valid_worlds), linfo)
    end
    unlock_mi_inference(interp, linfo)

    @static :backedges ∉ fieldnames(MethodInstance) && return nothing

    # this function is called from `_typeinf(interp::AbstractInterpreter, frame::InferenceState)`
    # and so at this point `linfo` is not registered in `ANALYZED_LINFOS`
    # (unless inference happens multiple times for this frame)
    if linfo ∉ ANALYZED_LINFOS
        add_jet_callback!(linfo)
    end

    return nothing
end

function add_jet_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_jet_cache!]
    else
        if !any(cb->cb!==invalidate_jet_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_jet_cache!)
        end
    end
end

function invalidate_jet_cache!(replaced, max_world, depth = 0)
    replaced ∉ ANALYZED_LINFOS && return
    delete!(ANALYZED_LINFOS, replaced)
    delete!(JET_GLOBAL_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            invalidate_jet_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end


text:
function spre(O::AbstractArray)
    N = size(O, 1)
    A = spdiagm( ones(ComplexF64, N) )
    B = O
    return kron(A, B)
end

function spost(O::AbstractArray)
    N = size(O, 1)
    A = sparse( transpose(O) )
    B = spdiagm( ones(ComplexF64, N) )
    return kron(A, B)
end

function sprepost(A, B::AbstractArray)
    return sparse( spre(A) * spost(B) )
end

function lindblad_dissipator(O::AbstractArray)
    Od_O = adjoint(O) * O
    return sprepost(O, adjoint(O)) - 0.5 * spre(Od_O) - 0.5 * spost(Od_O)
end

function destroy(N::Number)
    return spdiagm(1 => Array{ComplexF64}(sqrt.(1:N - 1)))
end

function eye(N::Number)
    return spdiagm(ones(ComplexF64, N))
end

function fock(N::Number, pos::Number)
    array = zeros(N)
    array[pos + 1] = 1
    return Array{ComplexF64}( array )
end

function projection(N::Number, i::Number, j::Number, shift = false)
    if shift
        return sparse( fock(N, i - 1) * adjoint(fock(N, j - 1)) )
    else
        return sparse( fock(N, i) * adjoint(fock(N, j)) )
    end
end

# function exp_gpu(M)
#     M_tmp = copy(M)
#     E_gpu, U_gpu = CUDA.CUSOLVER.heevd!('V','U', M_tmp)
#     return U_gpu * Diagonal(exp.(E_gpu)) * adjoint(U_gpu)
# end

function sinm(O::AbstractArray)
    M = collect(O)
    return sparse( 0.5im * (exp(1im * M) - exp(-1im * M)) )
end

function cosm(O::AbstractArray)
    M = collect(O)
    return sparse( 0.5 * (exp(1im * M) + exp(-1im * M)) )
end

function expect(op::AbstractArray, state::AbstractArray)
    if ishermitian(op)
        return real(adjoint(state) * op * state)
    else
        return adjoint(state) * op * state
    end
end


text:
include("euler/euler.jl")
using .Primes: get_primes

compute(n::Int)::Int = sum(get_primes(n))

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
######################
### Mutable states ###
######################

mutable struct StanHMCAdaptorState
    i           :: Int
    window_size :: Int
    next_window :: Int
end

################
### Adaptors ###
################

# TODO: currently only StanHMCAdaptor has the filed `n_adapts`. maybe we could unify all
# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.
struct StanHMCAdaptor{M<:AbstractPreconditioner, Tssa<:StepSizeAdaptor} <: AbstractAdaptor
    n_adapts    :: Int
    pc          :: M
    ssa         :: Tssa
    init_buffer :: Int
    term_buffer :: Int
    state       :: StanHMCAdaptorState
end
Base.show(io::IO, a::StanHMCAdaptor) =
    print(io, "StanHMCAdaptor(n_adapts=$(a.n_adapts), pc=$(a.pc), ssa=$(a.ssa), init_buffer=$(a.init_buffer), term_buffer=$(a.term_buffer))")

function StanHMCAdaptor(
    n_adapts::Int,
    pc::M,
    ssa::StepSizeAdaptor,
    init_buffer::Int=75,
    term_buffer::Int=50,
    window_size::Int=25
) where {M<:AbstractPreconditioner}
    next_window = init_buffer + window_size - 1
    return StanHMCAdaptor(n_adapts, pc, ssa, init_buffer, term_buffer, StanHMCAdaptorState(0, window_size, next_window))
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
function is_in_window(tp::StanHMCAdaptor)
    return (tp.state.i >= tp.init_buffer) &&
           (tp.state.i < tp.n_adapts - tp.term_buffer) &&
           (tp.state.i != tp.n_adapts)
end

function is_window_end(tp::StanHMCAdaptor)
    return (tp.state.i == tp.state.next_window) &&
           (tp.state.i != tp.n_adapts)
end

function compute_next_window!(tp::StanHMCAdaptor)
    if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
        tp.state.window_size *= 2
        tp.state.next_window = tp.state.i + tp.state.window_size
        if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
            next_window_boundary = tp.state.next_window + 2 * tp.state.window_size
            if (next_window_boundary >= tp.n_adapts - tp.term_buffer)
                tp.state.next_window = tp.n_adapts - tp.term_buffer - 1
            end
        end
    end
end

getM⁻¹(adaptor::StanHMCAdaptor) = getM⁻¹(adaptor.pc)
getϵ(adaptor::StanHMCAdaptor)   = getϵ(adaptor.ssa)
finalize!(adaptor::StanHMCAdaptor) = finalize!(adaptor.ssa)
# reset!(adaptor::StanHMCAdaptor) # Not used.

function adapt!(tp::StanHMCAdaptor, θ::AbstractVector{<:Real}, α::AbstractFloat)
    tp.state.i += 1

    adapt!(tp.ssa, θ, α)

    resize!(tp.pc, θ) # Resize pre-conditioner if necessary.

    # Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
    if is_in_window(tp)
        # We accumlate stats from θ online and only trigger the update of M⁻¹ in the end of window.
        is_update_M⁻¹ = is_window_end(tp)
        adapt!(tp.pc, θ, α, is_update_M⁻¹)
    end

    if is_window_end(tp)
        reset!(tp.ssa)
        reset!(tp.pc)
        # If window ends, compute next window
        compute_next_window!(tp)
    end
end


text:
module Pattern

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

abstract Pattern

type PVar <: Pattern
    name::Symbol
end

type TypeAssert <: Pattern
    T
end

type Product <: Pattern
    factors::Vector{Pattern}
end

type TuplePat <: Pattern
    ps::Vector{Pattern}
end

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

abstract Node

# asserts that all args are egal
type EgalNode <: Node
    args::Vector{Node}
end

# asserts that valueof(v) is of type valueof(T)
type IsaNode <: Node
    v::Node
    T::Node
end

# asserts that valueof(t) is a tuple(ps...)
type TupleNode <: Node
    t::Node
    ps::Vector{Node}
end

type Atom <: Node
    value
end

type Var <: Node
    name::Symbol
end

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

abstract Value

type Atom{T} <: Value
    value::T
end
type Var <: Value
    name::Symbol
end

abstract Guard

# asserts that all args are egal
type Egal <: Guard
    args::Vector{Value}
end

# asserts that valueof(v)::valueof(T)
type Isa <: Guard
    v::Value
    T::Value
end

# asserts that valueof(t) is a tuple(ps...)
type TupleGuard <: Guard
    t::Value
    ps::Vector{Value}
end

end # module


text:
gen_op(CAT)

c  = @select CAT where Close;
cy = collapse(c, year, last)

# daily
Cd   = percentchange(c).values
GMMd = GMM(3, Cd, method=:kmeans, nInit=50, nIter=10)
MMd  = MixtureModel(GMMd)

#yearly
Cy   = percentchange(cy).values
GMMy = GMM(3, Cy, method=:kmeans, nInit=50, nIter=10)
MMy  = MixtureModel(GMMy)

function rday(n::Int=13089)
    rand(MMd, n)
end

function ryear(n::Int=51)
    rand(MMy, n)
end


text:
# kernel
X = [1 1; 2 2;  1 -1]
exp_K = [1.0 0.367879  0.135335; 0.367879 1.0 0.00673795; 0.135335 0.00673795 1.0]
@test_approx_eq_eps rbf(X,X) exp_K 1e-6

# test mercer map
X = [1 1; 2 2;  1 -1]
K_centered = [0.444513 -0.144742 -0.299771; -0.144742 0.530245 -0.385503; -0.299771 -0.385503 0.685274]
V = [0.768623 -0.275472; -0.622878 -0.527911; -0.145746 0.803383]
d = [0.786544; 1.02048]
X_mapped = [0.604556 -0.281114; -0.489921 -0.538723; -0.114635 0.819837]

map = MercerMap(X, :rbf)

@test_approx_eq_eps map.K K_centered 1e-6
@test_approx_eq_eps map.d d 1e-6
@test_approx_eq_eps map.V V 1e-6
@test_approx_eq_eps RegERMs.apply(map) X_mapped 1e-6

X = [1 1; 0 0;  -1 -1]
XT = rand(100,2)
V = [-0.7071067811865475; 0.0; 0.7071067811865477]
d = [2.0]
X_mapped = [-1.414213562373095; 0.0; 1.4142135623730954]

map = MercerMap(X, :linear)
@test_approx_eq_eps map.K linear(X,X) 1e-6 # kernel is already centered
@test_approx_eq_eps map.d d 1e-6
@test_approx_eq_eps map.V V 1e-6
@test_approx_eq_eps RegERMs.apply(map, X) X_mapped 1e-6
RegERMs.apply(map, XT)


text:
# 2D Heat Diffusion
# LBM
#
# Generalize version on complicated mesh, with 1D indexing of lattices.
#
# Hongyang Zhou, hyzhou@umich.edu 05/01/2019

using PyPlot
include("Parameters.jl")
include("lattice.jl")
include("BC.jl")

using .Parameters, .BC, .lattice

function solveHeatDiffusion2D!(param::Param,lat::Lattice,mesh::Mesh)

# Set aliases
ω,Δt,TotalTime = param.ω,param.Δt,param.TotalTime
BCtype = param.BCtype
BCDirichlet,BCRiemann  = Parameters.BCDirichlet, Parameters.BCRiemann
Tleft,Tright,Ttop,Tbot = param.Tleft, param.Tright, param.Ttop, param.Tbot

nw,w,nN  = lat.nw, lat.w, lat.nN
T   = zeros(1,nN)
feq = Array{Float64,2}(undef,nw,nN)
# There will be an issue if I allocate with undef due to values on the boundary.
ftemp = zeros(nw,nN)
t   = 0.0     # Current simulation time

# Advance
while t < TotalTime

   # Calculate equilibrium distribution
   for iw = 1:nw
      @. feq[iw,:] = w[iw]*T[1,:]
   end

   # Collision
   for iw = 1:nw
      @. lat.f[iw,:] = (1.0-ω)*lat.f[iw,:] + ω.*feq[iw,:]
   end

   # Streaming
   for i=1:nN, iw = 1:nw
      iN = lat.IJ[iw,i]
      (iN > nN || iN < 1) && continue
      ftemp[iw,iN] = lat.f[iw,i]
   end
   #copy!(lat.f,ftemp) # This won't work in Julia 1.0!
   lat.f .= ftemp

   # Remaining BCs
   setBCs!(lat,param,mesh)

   # Calculate macroscopic T
   T .= sum(lat.f, dims=1)

   t += Δt
end

end

function init()

param,mesh = setParameters()

Δx,Δt,LT = param.Δx,param.Δt,param.LatticeType

if LT == 1
   lat = LatticeD2Q4(Δx,Δt,mesh)
elseif LT == 2
   lat = LatticeD2Q5(Δx,Δt,mesh)
elseif LT == 3
   lat = LatticeD2Q9(Δx,Δt,mesh)
end

return param,lat,mesh

end

function plotSol(lat::Lattice,mesh::Mesh)
# Visualization
#x = range(0,stop=lat.nI-1,length=lat.nI)
#y = range(0,stop=lat.nJ-1,length=lat.nJ)
T = sum(lat.f, dims=1)

tricontourf(mesh.coord[1,:],mesh.coord[2,:],T[1,:]); colorbar()
title("2D Heat Diffusion in a Plate")
xlabel("x")
ylabel("y")
grid(true)

# 1D cut plot
#plot(0:nx-1,T[1,:,50],linestyle="-",marker="o",label="LBM")
#xlim(0,30)
end

function main()
   param,lat,mesh = init()

   solveHeatDiffusion2D!(param,lat,mesh)

   if param.DoPlot
      plotSol(lat,mesh)
   end

end


text:
using Plots
using NeRCA
using BenchmarkTools

# pyplot()

dom_positions = [NeRCA.Position(0, 0, z)
                 for z in range(120, length=18, stop=800)]

track = NeRCA.Track([0, 1, 0], [0, 500, 500], 0)

plot(dom_positions, track)
plot!(dom_positions, NeRCA.Track([0,0,-1], [0,0,800]))
plot!(dom_positions, NeRCA.Track([0,-1,0], [0,1000,400]))

function eventplot(dom_positions, track::NeRCA.Track)
        plot3d([p.x for p in dom_positions],
               [p.y for p in dom_positions],
               [p.z for p in dom_positions],
               marker=2)
        plot3d!(
            [track.pos.x, track.pos.x + track.dir.x * 1000],
            [track.pos.y, track.pos.y + track.dir.y * 1000],
            [track.pos.z, track.pos.z + track.dir.z * 1000],
        )
end

eventplot(dom_positions, track)


text:
import .SIMD: Vec

@inline Vec(g::Base.Generator{<:LengthRange}) = Vec(Tuple(g))


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

JLLWrappers.@generate_main_file_header("RDKit")
JLLWrappers.@generate_main_file("RDKit", UUID("03d1d220-30e6-562a-9e1a-3062d7b56d75"))
end  # module RDKit_jll


text:
module DAGFlowColGen

using JuMP
using CPLEX
using Gurobi
using ..MyGraphModule
using ..DagFlowInstanceModule

mutable struct RestrictedMaster
    LP::JuMP.Model
    consRef
    outflowRef
    lambda
end

function setupMaster(instance::DagFlowInstance, dummyCost::Float64)
    LP = Model(solver = CplexSolver(CPX_PARAM_SIMDISPLAY=0))
    #LP = Model(solver = GurobiSolver())

    nSets = length(instance.nodePartition)
    outFlowLimits = instance.outFlowLimits
    # Defining the variables. One per customer, corresponding to a route just containing that customer.
    @variable(LP, lambda[1:nSets] >= 0)
    # Setting the objective
    @objective(LP, Min, sum(dummyCost*lambda[i] for i=1:nSets))
    # Each partition should be "covered"
    @constraint(LP, consRef[i=1:nSets], lambda[i] == 1)
    # out flow constraints
    @constraint(LP, outflowRef[i=1:length(outFlowLimits)], 0 <= outFlowLimits[i].limit)
    RMAS = RestrictedMaster(LP, consRef, outflowRef, lambda)
    return RMAS
end

function solvePricing(instance::DagFlowInstance, dualsPartition, dualsFlowLimit, maxNPaths)
    n,m = MyGraphModule.size(instance.g)
    # pricing problem:
    nodeBonus = zeros(Float64,n)
    for i = 1:length(instance.nodePartition)
        for node = instance.nodePartition[i]
            nodeBonus[node] = dualsPartition[i]
        end
    end
    for i = 1:length(instance.outFlowLimits)
        nodeBonus[instance.outFlowLimits[i].node] += dualsFlowLimit[i]
    end
    #shortestPath(g::MyGraph, nodeOrder::Vector{Int32}, nodeBonus::Vector{Float64}, source::Int32, sink::Int32, maxNPaths::Int64, eps::Float64)
    shortestPath(instance.g, instance.nodeOrder, nodeBonus, instance.sourceNode, instance.sinkNode, maxNPaths, -0.00001)
end

function addPathsToRMAS(instance::DagFlowInstance, RMAS::RestrictedMaster, paths)

    #try to print out the nodes in the path using nodeInfoDict from instance.
    #=for i=1:length(paths)
        path = paths[i]
        println("Path $i: ")
        for node in path
            println("$node: $(instance.nodeInfoDict[node])")
        end
    end=#
    for i=1:length(paths)
        nVisits = zeros(Int64, length(instance.nodePartition))
        outFlowVisits = zeros(Int64, length(instance.outFlowLimits))
        path = paths[i]
        origCost = 0
        prevNode = -1
        for node in path
            if instance.node2Partion[node]  > 0
                nVisits[instance.node2Partion[node]] += 1
            end
            if prevNode != -1
                origCost += MyGraphModule.cost(instance.g, prevNode, node)
                #println("c($prevNode,$node) = ", MyGraphModule.cost(instance.g, prevNode, node))
            end
            if prevNode in instance.setOutFlowNodes
                for j=1:length(instance.outFlowLimits)
                    if prevNode == instance.outFlowLimits[j].node
                        outFlowVisits[j] += 1
                    end
                end
            end
            prevNode = node
        end
        #println("nVisits: $nVisits")

        touchedConstraints = ConstraintRef[]
        vals = Float64[]
        for j=1:length(nVisits)
            if nVisits[j] > 0
                push!(touchedConstraints, RMAS.consRef[j])
                push!(vals,nVisits[j])
            end
        end
        for j=1:length(instance.outFlowLimits)
            if outFlowVisits[j] > 0
                push!(touchedConstraints, RMAS.outflowRef[j])
                push!(vals,outFlowVisits[j])
            end
        end

        @variable(
            RMAS.LP,                              # Model to be modified
            xNew >= 0,                            # New variable to be added
            objective=origCost,                   # cost coefficient of new varaible in the objective
            inconstraints=touchedConstraints ,    # constraints to be modified
            coefficients=vals                     # the coefficients of the variable in those constraints
        )
        push!(RMAS.lambda, xNew) # Pushing the new variable in the array of new variables
    end
end

function colgen(instance::DagFlowInstance)
    startTime = time()
    RMAS = setupMaster(instance, 10000.0)

    iter = 0
    done = false
    timeLP = 0
    timePP = 0
    timeAddCol = 0
    colAdded = 0
    while  !done #&& iter < 100
        timeRMSStart = time()
        solve(RMAS.LP)
        timeLP += time()-timeRMSStart

        println("$iter: RMAS objective value: ", getobjectivevalue(RMAS.LP))
        #writeLP(RMAS.LP, "master.lp", genericnames=true)
        dualsPartition = getdual(RMAS.consRef)
        dualsFlowLimit = getdual(RMAS.outflowRef)
        #println("dualsPartition: $dualsPartition")
        #println("dualsFlowLimit: $dualsFlowLimit")
        timePPStart = time()
        paths, costs = solvePricing(instance, dualsPartition, dualsFlowLimit, 5)
        timePP += time()-timePPStart
        if length(paths) > 0
            timeAddColStart = time()
            addPathsToRMAS(instance, RMAS, paths)
            timeAddCol += time()-timeAddColStart
            colAdded += length(paths)
        else
            done = true
        end
        iter +=1
        println("iter = $iter, done = $done")
    end
    totTime = time()-startTime
    println("Total time: $totTime, RMAS LP time: $timeLP, PP time: $timePP, add columns time: $timeAddCol, columns added: $colAdded")

end

end


text:
# This file was generated, do not modify it. # hide
@show round(minimum(curve.measurements), sigdigits=3)
@show round(maximum(curve.measurements), sigdigits=3)