text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots

# Condiciones iniciales del problema
u0 = Float32[2.0; 0.0]
datasize = 30
tspan = (0.0f0, 1.5f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

# Se define y resuelve el problema con parámetros reales
function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))

# Se define una red neuronal de con dos capas densas/lineales.
# a ésta se le agrega la asunción de que el sistema transforma
# x hacia x^3. Esto como conocimiento previo que mejore el tiempo.
dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, tanh),
                  FastDense(50, 2))

# Se define el problema junto con su método de resolución y tiempo al cual parar.
prob_neuralode = NeuralODE(dudt2, tspan, Tsit5(), saveat = tsteps)

# Funciones que guardarán la predicción de la red neuronal y su valor de error.
function predict_neuralode(p)
  Array(prob_neuralode(u0, p))
end

function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end

# Utilizado para graficar a tiempo real las predicciones realizadas.
callback = function (p, l, pred; doplot = true)
  display(l)
  # plot current prediction against data
  plt = scatter(tsteps, ode_data[1,:], label = "data")
  scatter!(plt, tsteps, pred[1,:], label = "prediction")
  if doplot
    display(plot(plt))
  end
  return false
end

# Se entrena la red neuronal utilizando la data de la resolución con parámetros reales.
# Ésto se hace en dos pasos, primero con el optimizador ADAM para acercarnos rápidamente
# al mínimo, y una vez cerca utilizamos Limited-memory BFGS para un ajuste más eficiente.
result_neuralode = DiffEqFlux.sciml_train(loss_neuralode, prob_neuralode.p,
                                          ADAM(0.05), cb = callback,
                                          maxiters = 300)

result_neuralode2 = DiffEqFlux.sciml_train(loss_neuralode,
                                           result_neuralode.minimizer,
                                           LBFGS(),
                                           cb = callback,
                                           allow_f_increases = false)

max_stars_repo_name: anathompson/taller-julia-macaotech

max_issues_repo_name: anathompson/taller-julia-macaotech

max_forks_repo_name: anathompson/taller-julia-macaotech


text:
# testing a set of points
# 1d case
using Cubature

include("model.jl")
include("plotting/plot_distribution_single.jl")
include("root_finding/zero_finding.jl")

function BTG_main(training_data, testing_data, kernel, nonlinfun)
    # training points
    x_train = training_data[1]
    z_train = training_data[2]
    n_train = size(x_train,1)

    # Test data
    x_test = testing_data[1]
    z_true = testing_data[2]
    n_test = size(x_test,1)

    for i in 1:n_test
        # select where to predict
        x_test_temp = x_test[i, :]
        z_true_temp = z_true[i]
        println("Current i = $i \n x_pred = $x_test_temp")
        flush(stdout)
        # single point prediction
        start = time()
        distribution_temp, alpha = model(x_test_temp, training_data, kernel, nonlinfun)
        elapsed = time() - start
        # print("Alpha = $alpha")
        print("time taken for single point prediction: $elapsed \n")

        # check pdf integral == 1
        # get the current pdf
        pdf_temp = distribution_temp.pdf
        # get the upper bound for integral
        b = distribution_temp.upperbd
        int_temp = hquadrature(pdf_temp, 1e-6, b)[1]
        print("Current integral interval: [0, $b], int(pdf) = $int_temp \n")

        # plot
        PyPlot.clf()
        plot_distribution_single(distribution_temp, x_test_temp, z_true_temp)
        PyPlot.savefig("figures/test_$i.pdf")
        println("------------------------")
    end
end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using MyPkg
using QuadGK
using Test

integral, err = quadgk(MyPkg.g,0,1,rtol=1e-5)

@test Simpson(MyPkg.g,0,1) == integral
@test Simpson(MyPkg.g,0,1) == 1
@test Simpson(MyPkg.g,1,2) == 3


text:
using Lux, Random, Test

include("../utils.jl")

rng = Random.default_rng()
Random.seed!(rng, 0)

@testset "Dropout" begin
    layer = Dropout(0.5f0)
    println(layer)
    ps, st = Lux.setup(rng, layer)
    x = randn(Float32, 5, 2)

    x_, st_ = layer(x, ps, st)
    x__, st__ = layer(x, ps, st)
    x___, st___ = layer(x_, ps, st_)

    @test st_.rng != st.rng
    @test st_.rng == st__.rng
    @test x_ == x__
    @test x_ != x___

    run_JET_tests(layer, x, ps, st)
    test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                  rtol=1.0f-3)

    st = Lux.testmode(st)

    @test first(layer(x, ps, st)) == x
end

@testset "VariationalHiddenDropout" begin
    layer = VariationalHiddenDropout(0.5f0)
    println(layer)
    ps, st = Lux.setup(rng, layer)
    x = randn(Float32, 5, 2)

    x_, st_ = layer(x, ps, st)
    x__, st__ = layer(x, ps, st)
    x___, st___ = layer(x_, ps, st_)

    @test st_.rng != st.rng
    @test st_.rng == st__.rng
    @test st_.mask == st__.mask
    @test x_ == x__
    @test x_ != x___

    run_JET_tests(layer, x, ps, st)
    run_JET_tests(layer, x, ps, st_)
    test_gradient_correctness_fdm(x -> sum(layer(x, ps, st)[1]), x; atol=1.0f-3,
                                  rtol=1.0f-3)
    test_gradient_correctness_fdm(x -> sum(layer(x, ps, st_)[1]), x; atol=1.0f-3,
                                  rtol=1.0f-3)

    st__ = Lux.update_state(st_, :update_mask, Val(true))
    x___, st___ = layer(x, ps, st__)

    @test st___.mask != st__.mask
    @test x___ != x_

    run_JET_tests(layer, x, ps, st__)
    test_gradient_correctness_fdm(x -> sum(layer(x, ps, st__)[1]), x; atol=1.0f-3,
                                  rtol=1.0f-3)
end


text:
using DiffEqGPU, CuArrays, OrdinaryDiffEq, Test

function lorenz(du,u,p,t)
 @inbounds begin
     du[1] = p[1]*(u[2]-u[1])
     du[2] = u[1]*(p[2]-u[3]) - u[2]
     du[3] = u[1]*u[2] - p[3]*u[3]
 end
 nothing
end

CuArrays.allowscalar(false)
u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = (10.0f0,28.0f0,8/3f0)
prob = ODEProblem(lorenz,u0,tspan,p)
const pre_p = [rand(Float32,3) for i in 1:10]
prob_func = (prob,i,repeat) -> remake(prob,p=pre_p[i].*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func)

#Performance check with nvvp
# CUDAnative.CUDAdrv.@profile
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@time sol = solve(monteprob,ROCK4(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@time solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,
                                                 batch_size=5,saveat=1.0f0)
@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@time solve(monteprob,Tsit5(),EnsembleCPUArray(),trajectories=10,saveat=1.0f0)
@time solve(monteprob,Tsit5(),EnsembleThreads(), trajectories=10,saveat=1.0f0)
@time solve(monteprob,Tsit5(),EnsembleSerial(),  trajectories=10,saveat=1.0f0)

#=
solve(monteprob,TRBDF2(),EnsembleCPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
solve(monteprob,TRBDF2(),EnsembleGPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
@test_broken solve(monteprob,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=2,saveat=1.0f0)
=#

function lorenz_jac(J,u,p,t)
 @inbounds begin
     σ = p[1]
     ρ = p[2]
     β = p[3]
     x = u[1]
     y = u[2]
     z = u[3]
     J[1,1] = -σ
     J[2,1] = ρ - z
     J[3,1] = y
     J[1,2] = σ
     J[2,2] = -1
     J[3,2] = x
     J[1,3] = 0
     J[2,3] = -x
     J[3,3] = -β
 end
 nothing
end

function lorenz_tgrad(J,u,p,t)
 nothing
end

func = ODEFunction(lorenz,jac=lorenz_jac,tgrad=lorenz_tgrad)
prob_jac = ODEProblem(func,u0,tspan,p)
monteprob_jac = EnsembleProblem(prob_jac, prob_func = prob_func)

@time solve(monteprob_jac,Rodas5(linsolve=LinSolveGPUSplitFactorize()),EnsembleCPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,Rodas5(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleCPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)
@time solve(monteprob_jac,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),EnsembleGPUArray(),dt=0.1,trajectories=10,saveat=1.0f0)

condition = function (u,t,integrator)
    @inbounds u[1] > 5
end

affect! = function (integrator)
    @inbounds integrator.u[1] = -4
end

callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=DiscreteCallback(condition,affect!,save_positions=(false,false)))
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)

c_condition = function (u,t,integrator)
    @inbounds u[1] - 3
end

c_affect! = function (integrator)
    @inbounds integrator.u[1] += 20
end

callback_prob = ODEProblem(lorenz,u0,tspan,p,callback=ContinuousCallback(c_condition,c_affect!,save_positions=(false,false)))
callback_monteprob = EnsembleProblem(callback_prob, prob_func = prob_func)
@test_broken solve(callback_monteprob,Tsit5(),EnsembleGPUArray(),trajectories=2,saveat=1.0f0).retcode == :Success

#=
using OrdinaryDiffEq, LinearAlgebra, ParameterizedFunctions
LinearAlgebra.BLAS.set_num_threads(1)
rober = @ode_def begin
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
=#

function rober_f(internal_var___du,internal_var___u,internal_var___p,t)
    @inbounds begin
        internal_var___du[1] = -(internal_var___p[1]) * internal_var___u[1] + internal_var___p[3] * internal_var___u[2] * internal_var___u[3]
        internal_var___du[2] = (internal_var___p[1] * internal_var___u[1] - internal_var___p[2] * internal_var___u[2] ^ 2) - internal_var___p[3] * internal_var___u[2] * internal_var___u[3]
        internal_var___du[3] = internal_var___p[2] * internal_var___u[2] ^ 2
    end
    nothing
end

function rober_jac(internal_var___J,internal_var___u,internal_var___p,t)
    @inbounds begin
        internal_var___J[1, 1] = -(internal_var___p[1])
        internal_var___J[1, 2] = internal_var___p[3] * internal_var___u[3]
        internal_var___J[1, 3] = internal_var___p[3] * internal_var___u[2]
        internal_var___J[2, 1] = internal_var___p[1] * 1
        internal_var___J[2, 2] = -2 * internal_var___p[2] * internal_var___u[2] - internal_var___p[3] * internal_var___u[3]
        internal_var___J[2, 3] = -(internal_var___p[3]) * internal_var___u[2]
        internal_var___J[3, 1] = 0 * 1
        internal_var___J[3, 2] = 2 * internal_var___p[2] * internal_var___u[2]
        internal_var___J[3, 3] = 0 * 1
    end
    nothing
end

rober_prob = ODEProblem(ODEFunction(rober_f,jac=rober_jac),Float32[1.0,0.0,0.0],(0.0,1f5),(0.04,3f7,1f4))
sol = solve(rober_prob,Rodas5(),abstol=1f-8,reltol=1f-8)
rober_monteprob = EnsembleProblem(rober_prob, prob_func = prob_func)

@time sol = solve(rober_monteprob,TRBDF2(linsolve=LinSolveGPUSplitFactorize()),
                  EnsembleGPUArray(),trajectories=10,saveat=1.0f0,abstol=1f-8,reltol=1f-8)
@time sol = solve(monteprob,TRBDF2(),EnsembleThreads(),trajectories=10,
                  abstol=1e-8,reltol=1e-8,saveat=1.0f0)


text:
register_size(::Type{T}) where {T} = REGISTER_SIZE
register_size(::Type{T}) where {T<:Union{Signed,Unsigned}} = SIMD_INTEGER_REGISTER_SIZE

intlog2(N::I) where {I <: Integer} = (8sizeof(I) - one(I) - leading_zeros(N)) % I
intlog2(::Type{T}) where {T} = intlog2(sizeof(T))
ispow2(x::Integer) = (x & (x - 1)) == zero(x)
nextpow2(W) = vshl(one(W), vsub(8sizeof(W), leading_zeros(vsub(W, one(W)))))
prevpow2(W) = vshl(one(W), vsub(vsub((8sizeof(W)) % UInt, one(UInt)), leading_zeros(W) % UInt))
# prevpow2(W) = (one(W) << ((((8sizeof(W)) % UInt) - one(UInt)) - (leading_zeros(W) % UInt)))
prevpow2(W::Signed) = prevpow2(W % Unsigned) % Signed

# @generated function pick_vector_width(::Type{T}) where {T<:NativeTypes}
#     max(1, register_size(T) >>> intlog2(T))
# end
function pick_vector_width_shift(::Type{T}) where {T<:NativeTypes}
    # W = pick_vector_width(T)
    Wshift = intlog2(register_size(T)) - intlog2(T)
    1 << Wshift, Wshift
end
function pick_vector_width_shift(N::Integer, ::Type{T}) where {T<:NativeTypes}
    Wshift_N = VectorizationBase.intlog2(2N - 1)
    Wshift_st = intlog2(register_size(T)) - VectorizationBase.intlog2(sizeof(T))
    Wshift = min(Wshift_N, Wshift_st)
    W = 1 << Wshift
    W, Wshift
end

# For the sake of convenient mask support, we allow 8 so that the mask can be a full byte
# max_vector_width(::Type{T}) where {T} = max(8, pick_vector_width(T))

# @inline pick_vector_width(::Type{T1}, ::Type{T2}) where {T1,T2} = min(pick_vector_width(T1), pick_vector_width(T2))
# @inline pick_vector_width(::Type{T1}, ::Type{T2}, ::Type{T3}, args::Vararg{Any,K}) where {T1,T2,T3,K} = min(pick_vector_width(T1), pick_vector_width(T2, T3, args...))

function pick_vector_width_shift_from_size(N::Int, size_T::Int)
    Wshift_N = VectorizationBase.intlog2(2N - 1)
    Wshift_st = intlog2(REGISTER_SIZE) - VectorizationBase.intlog2(size_T)
    Wshift = min(Wshift_N, Wshift_st)
    W = 1 << Wshift
    W, Wshift
end

# @inline function pick_vector_width(N::Integer, args::Vararg{Any,K}) where {K}
#     min(nextpow2(N), pick_vector_width(args...))
# end
# @inline function pick_vector_width_shift(N::Integer, args::Vararg{Any,K}) where {K}
#     W = pick_vector_width(N, args...)
#     W, intlog2(W)
# end

# @generated pick_vector_width(::Union{Val{N},StaticInt{N}}, ::Type{T}) where {N,T} = pick_vector_width(N,T)
# @generated function pick_vector_width(::Union{Val{N},StaticInt{N}}, ::Type{T}, ::Type{T2}, args::Vararg{Any,K}) where {N,T,T2,K}
#     pvw = pick_vector_width(N, T)
#     quote
#         $(Expr(:meta,:inline))
#         min($pvw, pick_vector_width_val(T2, args...))
#     end
# end
# pick_vector_width(::Union{Val{N},StaticInt{N}}, ::Type{T}, args::Vararg{Any,K}) where {N,T,K} = min(@show(nextpow2(N)), @show(pick_vector_width(T, args...)))
# pick_vector_width(::Union{Val{N},StaticInt{N}}, arg, args::Vararg{Any,K}) where {N,K} = min(nextpow2(N), pick_vector_width(arg, args...))
# pick_vector_width(::Val{N}, args::Vararg{Any,K}) where {N,K} = ((Nmax,a,W) = @show((N, args, pick_vector_width(N, args...))); W)
# @inline pick_vector_width_val(::Type{T}) where {T} = StaticInt{pick_vector_width(T)}()
# @inline pick_vector_width_val(::Union{Val{N},StaticInt{N}}, ::Type{T}) where {N,T} = StaticInt{pick_vector_width(Val(N), T)}()
# @inline pick_vector_width_val(::Type{T1}, ::Type{T2}, args::Vararg{Any,K}) where {T1,T2,K} = StaticInt{pick_vector_width(T1,T2,args...)}()
# @generated function pick_vector_width_val(::Union{Val{N},StaticInt{N}}, ::Type{T1}, ::Type{T2}, args::Vararg{Any,K}) where {T1,T2,N,K}
#     Wstart = pick_vector_width(N, T1, T2)
#     iszero(K) && return Expr(:call, Expr(:curly, :StaticInt, Wstart))
#     Expr(:block, Expr(:call, Expr(:curly, :StaticInt, Expr(:call, :min, Wstart, :(pick_vector_width(args...))))))

#     # W = min(StaticInt{N}(), min(pick_vector_width_val(T1), pick_vector_width_val(T2, args...)))
#     # StaticInt{W}()
#     # StaticInt{pick_vector_width(Val{N}(),T1,T2,args...)}()
# end

function _pick_vector_width(vargs...)
    min_W = 1
    max_W = REGISTER_SIZE
    for v ∈ vargs
        T = v.parameters[1]
        if T === Bit
            min_W = 8
        elseif (SIMD_INTEGER_REGISTER_SIZE != REGISTER_SIZE) && T <: Integer # only check subtype if it matters
            max_W = min(max_W, SIMD_INTEGER_REGISTER_SIZE ÷ sizeof(T))
        else
            max_W = min(max_W, REGISTER_SIZE ÷ sizeof(T))
        end
    end
    W = max(min_W, max_W)
end
@generated function pick_vector_width_val(vargs...)
    W = _pick_vector_width(vargs...)
    Expr(:call, Expr(:curly, :StaticInt, W))
end
adjust_W(N, W) = min(nextpow2(N), W)
@generated function pick_vector_width_val(::Union{Val{N},StaticInt{N}}, vargs...) where {N}
    W = adjust_W(N, _pick_vector_width(vargs...))
    Expr(:call, Expr(:curly, :StaticInt, W))
end
pick_vector_width(::Union{StaticInt{N},Val{N}}, args...) where {N} = Int(pick_vector_width_val(StaticInt{N}(), args...))

pick_vector_width(::Type{T}) where {T} = Int(pick_vector_width_val(T))
pick_vector_width(N::Integer, T) = min(nextpow2(N), pick_vector_width(T))

function int_type_symbol(W)
    bits = 8*(SIMD_INTEGER_REGISTER_SIZE ÷ W)
    if bits ≤ 8
        :Int8
    elseif bits ≤ 16
        :Int16
    elseif bits ≤ 32
        :Int32
    else # even if Int === Int32? Or should this be `Int`?
        :Int64
    end
end
@generated int_type(::Union{Val{W},StaticInt{W}}) where {W} = int_type_symbol(W)

pick_vector(::Type{T}) where {T} = _pick_vector(pick_vector_width_val(T), T)
_pick_vector(::StaticInt{W}, ::Type{T}) where {W,T} = Vec{W,T}
@generated pick_vector(::Val{N}, ::Type{T}) where {N, T} =  Expr(:curly, :Vec, pick_vector_width(N, T), T)
pick_vector(N::Int, ::Type{T}) where {T} = pick_vector(Val(N), T)

@inline MM(::Union{Val{W},StaticInt{W}}) where {W} = MM{W}(0)
@inline MM(::Union{Val{W},StaticInt{W}}, i) where {W} = MM{W}(i)
# @inline MM{W}(a::LazyMul) where {W} = MM{W}(data(a))
@inline gep(ptr::Ptr, i::MM) = gep(ptr, i.i)

@inline staticm1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vsub(i.i, one(I)))
@inline staticp1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vadd(i.i, one(I)))
@inline vadd(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vadd(i.i, j))
@inline vadd(i::Integer, j::MM{W,X}) where {W,X} = MM{W,X}(vadd(i, j.i))
@inline vadd(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd(i.i, j))
@inline vadd(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd(i, j.i))
@inline vadd(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
@inline vadd(::StaticInt{0}, j::MM{W,X}) where {W,X} = j
# @inline vadd(i::MM{W,X}, j::MM{W,S}) where {W,X,S} = MM{W}(vadd(i.i, j.i), StaticInt{X}() + StaticInt{S}())
@inline vsub(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vsub(i.i, j))
@inline vsub(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub(i.i, j))
@inline vsub(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
# @inline vmul_no_promote(i, j) = vmul(i,j)
# @inline vmul_no_promote(i, j::MM{W,X}) where {W,X} = MM{W,X}(vmul(i, j.i))
# @inline vmuladdnp(a, b, c) = vadd(vmul(a,b), c)
# @inline vmuladdnp(a, b::MM{W,X}, c) where {W,X} = vadd(MM{W,X}(vmul(a,b.i)), c)
# @inline vmuladdnp(a, b::MM{W,X}, c::MM{W,X}) where {W,X} = vadd(vmul(a,b), c)
# @inline vmuladdnp(a, b::MM{W,X}, c::Vec) where {W,X} = vadd(vmul(a,b), c)
# @inline vmuladdnp(a, b::MM{W,X}, c::_Vec) where {W,X} = vadd(vmul(a,b), c)

@inline Base.:(+)(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vadd(i.i, j))
@inline Base.:(+)(i::Integer, j::MM{W,X}) where {W,X} = MM{W,X}(vadd(i, j.i))
@inline Base.:(+)(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd(i.i, j))
@inline Base.:(+)(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd(i, j.i))
# @inline Base.:(+)(i::MM{W}, j::MM{W}) where {W} = MM{W}(i.i + j.i)
@inline Base.:(-)(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vsub(i.i, j))
# @inline Base.:(-)(i::Integer, j::MM{W}) where {W} = MM{W}(i - j.i)
@inline Base.:(-)(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub(i.i, j))
# @inline Base.:(-)(i::MM, ::StaticInt{0}) = i
@inline Base.:(-)(i::MM) = i * StaticInt{-1}()
# @inline Base.:(+)(i::MM, ::StaticInt{0}) = i
# @inline Base.:(+)(::StaticInt{0}, i::MM) = i
@inline Base.:(*)(::StaticInt{M}, i::MM{W,X}) where {M,W,X} = MM{W}(i.i * StaticInt{M}(), StaticInt{X}() * StaticInt{M}())
@inline Base.:(*)(i::MM{W,X}, ::StaticInt{M}) where {M,W,X} = MM{W}(i.i * StaticInt{M}(), StaticInt{X}() * StaticInt{M}())
# @inline Base.:(-)(::StaticInt{i}, j::MM{W}) where {W,i} = MM{W}(i - j.i)
# @inline Base.:(-)(i::MM{W}, j::MM{W}) where {W} = MM{W}(i.i - j.i)
# @inline Base.:(*)(i::MM{W}, j) where {W} = MM{W}(i.i * j)
# @inline Base.:(*)(i, j::MM{W}) where {W} = MM{W}(i * j.i)
# @inline Base.:(*)(i::MM{W}, ::StaticInt{j}) where {W,j} = MM{W}(i.i * j)
# @inline Base.:(*)(::StaticInt{i}, j::MM{W}) where {W,i} = MM{W}(i * j.i)
# @inline Base.:(*)(i::MM{W}, j::MM{W}) where {W} = MM{W}(i.i * j.i)

# @inline scalar_less(i, j) = i < j
# @inline scalar_less(i::MM, j::Integer) = i.i < j
# @inline scalar_less(i::Integer, j::MM) = i < j.i
# @inline scalar_less(i::MM, ::StaticInt{j}) where {j} = i.i < j
# @inline scalar_less(::StaticInt{i}, j::MM) where {i} = i < j.i
# @inline scalar_less(i::MM, j::MM) = i.i < j.i
# @inline scalar_greater(i, j) = i > j
# @inline scalar_greater(i::MM, j::Integer) = i.i > j
# @inline scalar_greater(i::Integer, j::MM) = i > j.i
# @inline scalar_greater(i::MM, ::StaticInt{j}) where {j} = i.i > j
# @inline scalar_greater(::StaticInt{i}, j::MM) where {i} = i > j.i
# @inline scalar_greater(i::MM, j::MM) = i.i > j.i
# @inline scalar_equal(i, j) = i == j
# @inline scalar_equal(i::MM, j::Integer) = i.i == j
# @inline scalar_equal(i::Integer, j::MM) = i == j.i
# @inline scalar_equal(i::MM, ::StaticInt{j}) where {j} = i.i == j
# @inline scalar_equal(::StaticInt{i}, j::MM) where {i} = i == j.i
# @inline scalar_equal(i::MM, j::MM) = i.i == j.i
# @inline scalar_notequal(i, j) = i != j
# @inline scalar_notequal(i::MM, j::Integer) = i.i != j
# @inline scalar_notequal(i::Integer, j::MM) = i != j.i
# @inline scalar_notequal(i::MM, ::StaticInt{j}) where {j} = i.i != j
# @inline scalar_notequal(::StaticInt{i}, j::MM) where {i} = i != j.i
# @inline scalar_notequal(i::MM, j::MM) = i.i != j.i

@inline Base.:(==)(::AbstractIrrational, ::MM{W,<:Integer}) where {W} = zero(Mask{W})
@inline Base.:(==)(x::AbstractIrrational, i::MM{W}) where {W} = x == Vec(i)
@inline Base.:(==)(::MM{W,<:Integer}, ::AbstractIrrational) where {W} = zero(Mask{W})
@inline Base.:(==)(i::MM{W}, x::AbstractIrrational) where {W} = Vec(i) == x

@generated function Base.promote_rule(::Type{MM{W,X,I}}, ::Type{T2}) where {W,X,I,T2<:NativeTypes}
    if register_size(T2) ≥ sizeof(T2) * W
        return :(Vec{$W,$T2})
    elseif T2 <: Signed
        return :(Vec{$W,$(int_type_symbol(W))})
    elseif T2 <: Unsigned
        return :(Vec{$W,unsigned($(int_type_symbol(W)))})
    else
        return :(Vec{$W,$T2})
    end
end

# @inline _vload(ptr::Ptr{T}, i::Integer) where {T} = vload(ptr + vmul(sizeof(T), i))
# @inline _vload(ptr::Ptr, v::Vec{<:Any,<:Integer}) = vload(ptr, v.data)
# @inline _vload(ptr::Ptr, v::Vec{<:Any,<:Integer}) = vload(ptr, v)
# @inline _vload(ptr::Ptr{T}, i::MM{W}) where {W,T} = vload(Val{W}(), ptr + vmul(sizeof(T), i.i))
# @inline _vload(ptr::AbstractPointer, i) = _vload(ptr.ptr, offset(ptr, i))
# @inline vload(ptr::AbstractPointer{T}, i::Tuple) where {T} = _vload(ptr.ptr, offset(ptr, i))


text:
#from synapseclient.entity import Entity, Project, Folder, File, split_entity_namespaces, is_container
#from synapseclient.exceptions import *
import SynapseClient: Entity
using SynapseClient.entity

@testset "Entity" begin
    # Test the basics of creating and accessing properties on an entity
    for i in 0:1
        e = Entity(name="Test object", description="I hope this works",
                   annotations = Dict("foo"=>123, "nerds"=>["chris","jen","janey"], "annotations"=>"How confusing!"),
                   properties  = Dict("annotations"=>"/repo/v1/entity/syn1234/annotations",
                                      "md5"=>"cdef636522577fc8fb2de4d95875b27c",
                                      "parentId"=>"syn1234"),
                   concreteType="org.sagebionetworks.repo.model.Data")

        # Should be able to create an Entity from an Entity
        if i == 1
            e = create(Entity,e)
        end
        @test e.parentId == "syn1234"

        @test e.properties["parentId"] =="syn1234"

        @test e.foo == 123

        @test e.annotations["foo"] == 123

        @test hasattr(e, "parentId")
        @test hasattr(e, "foo")
        @test hasattr(e, "qwerqwer") == false

        # Annotations is a bit funny, because there is a property call
        # "annotations", which will be masked by a member of the object
        # called "annotations". Because annotations are open-ended, we
        # might even have an annotations called "annotations", which gets
        # really confusing.
        @test typeof(e.annotations) <: AbstractDict

        @test e.properties["annotations"] == "/repo/v1/entity/syn1234/annotations"
        @test e.annotations["annotations"] == "How confusing!"

        @test e.nerds == ["chris","jen","janey"]
        @test all(k->hasattr(e,k), [:name, :description, :foo, :nerds, :annotations, :md5, :parentId])

        # Test modifying properties
        e.description = "Working, so far"
        @test e.description == "Working, so far"
        e.description = "Wiz-bang flapdoodle"
        @test e.description == "Wiz-bang flapdoodle"

        # Test modifying annotations
        e.foo = 999
        @test e.annotations["foo"] == 999
        e.foo = 12345
        @test e.annotations["foo"] == 12345

        # Test creating a new annotation
        e.bar = 888
        @test e.annotations["bar"] == 888
        e.bat = 7788
        @test e.annotations["bat"] == 7788

        # # Test replacing annotations object
        e.annotations = Dict("splat"=>"a totally new set of annotations", "foo"=>456)
        @test e.foo == 456

        @test typeof(e.annotations) <: AbstractDict

        @test e.annotations["foo"] == 456

        @test e.properties["annotations"] == "/repo/v1/entity/syn1234/annotations"

        ## test unicode properties
        e.train = "時刻表には記載されない　月への列車が来ると聞いて"
        e.band = "Motörhead"
        e.lunch = "すし"

        println(e)

    end
end
# def test_subclassing():
#     """Test ability to subclass and add a member variable"""

#     ## define a subclass of Entity to make sure subclassing and creating
#     ## a new member variable works
#     class FoobarEntity(Entity):
#         def __init__(self, x):
#             self.__dict__['x'] = x

#     foobar = FoobarEntity(123)
#     assert foobar.x == 123
#     assert 'x' in foobar.__dict__
#     assert foobar.__dict__['x'] == 123
#     foobar.id = 'syn999'
#     assert foobar.properties['id'] == 'syn999'
#     foobar.n00b = 'henry'
#     assert foobar.annotations['n00b'] == 'henry'

@testset "entity_creation" begin
    props = Dict(
        "id"=>"syn123456",
        "concreteType"=>"org.sagebionetworks.repo.model.Folder",
        "parentId"=>"syn445566",
        "name"=>"Testing123"
    )
    annos = Dict("testing"=>123)
    folder = create(Entity, props, annos)

    @test folder.concreteType == "org.sagebionetworks.repo.model.Folder"
    @test typeof(folder) == Folder
    @test folder.name == "Testing123"
    @test folder.testing == 123

    ## In case of unknown concreteType, fall back on generic Entity object
    props = Dict(
        "id"=>"syn123456",
        "concreteType"=>"org.sagebionetworks.repo.model.DoesntExist",
        "parentId"=>"syn445566",
        "name"=>"Whatsits"
    )
    whatsits = create(Entity,props)

    @test whatsits.concreteType == "org.sagebionetworks.repo.model.DoesntExist"
    @test typeof(whatsits) == Entity
end

@testset "parent_id_required" begin
    xkcd1 = File("http://xkcd.com/1343/", name="XKCD: Manuals", parent="syn1000001", synapseStore=false)
    @test xkcd1.parentId == "syn1000001"

    xkcd2 = File("http://xkcd.com/1343/", name="XKCD: Manuals", parentId="syn1000002", synapseStore=false)
    @test xkcd2.parentId == "syn1000002"

    @test_pythrows SynapseMalformedEntityError File("http://xkcd.com/1343/", name="XKCD: Manuals", synapseStore=false)
end

@testset "entity_constructors" begin
    project = Project("TestProject", id="syn1001", foo="bar")
    @test project.name == "TestProject"
    @test project.foo == "bar"

    folder = Folder("MyFolder", parent=project, foo="bat", id="syn1002")
    @test folder.name == "MyFolder"
    @test folder.foo == "bat"
    @test folder.parentId == "syn1001"

    a_file = File("/path/to/fabulous_things.zzz", parent=folder, foo="biz", contentType="application/cattywampus")
    #@test a_file.name == "fabulous_things.zzz"
    @test a_file.concreteType == "org.sagebionetworks.repo.model.FileEntity"
    @test a_file.path == "/path/to/fabulous_things.zzz"
    @test a_file.foo == "biz"
    @test a_file.parentId == "syn1002"
    @test a_file.contentType == "application/cattywampus"
    @test "contentType" in keys(a_file._file_handle)
end

@testset "property_keys" begin
    @test "parentId" in PyFile._property_keys
    @test "versionNumber" in PyFile._property_keys
    @test "dataFileHandleId" in PyFile._property_keys
end

@testset "keys" begin
    f = File("foo.xyz", parent="syn1234", foo="bar")

    # iter_keys = collect(keys(f)) # TODO: implement keys(e::AbstractEntity)

    # @test "parentId" in iter_keys == true
    # @test "name" in iter_keys == true
    # @test "foo" in iter_keys == true
    # @test "concreteType" in iter_keys == true
end

@testset "attrs" begin
    f = File("foo.xyz", parent="syn1234", foo="bar")
    @test hasattr(f, "parentId")
    @test hasattr(f, "foo")
    @test hasattr(f, "path")
end

@testset "split_entity_namespaces" begin
    # """Test split_entity_namespaces"""

    e = Dict("concreteType"=>"org.sagebionetworks.repo.model.Folder",
         "name"=>"Henry",
         "color"=>"blue",
         "foo"=>1234,
         "parentId"=>"syn1234")
    (properties,annotations,local_state) = split_entity_namespaces(e)

    @test Set(keys(properties)) == Set(["concreteType", "name", "parentId"])
    @test properties["name"] == "Henry"
    @test Set(keys(annotations)) == Set(["color", "foo"])
    @test annotations["foo"] == 1234
    @test length(local_state) == 0

    e = Dict("concreteType"=>"org.sagebionetworks.repo.model.FileEntity",
         "name"=>"Henry",
         "color"=>"blue",
         "foo"=>1234,
         "parentId"=>"syn1234",
         "dataFileHandleId"=>54321,
         "cacheDir"=>"/foo/bar/bat",
         "files"=>["foo.xyz"],
         "path"=>"/foo/bar/bat/foo.xyz")
    (properties,annotations,local_state) = split_entity_namespaces(e)

    @test Set(keys(properties)) == Set(["concreteType", "name", "parentId", "dataFileHandleId"])
    @test properties["name"] == "Henry"
    @test properties["dataFileHandleId"] == 54321
    @test Set(keys(annotations)) == Set(["color", "foo"])
    @test annotations["foo"] == 1234
    @test Set(keys(local_state)) == Set(["cacheDir", "files", "path"])
    @test local_state["cacheDir"] == "/foo/bar/bat"

    f = create(Entity,properties,annotations,local_state)
    @test f.properties["dataFileHandleId"] == 54321
    @test f.properties["name"] == "Henry"
    @test f.annotations["foo"] == 1234
    @test f.__dict__["cacheDir"] == "/foo/bar/bat"
    @test f.__dict__["path"] == "/foo/bar/bat/foo.xyz"
end

@testset "concrete_type" begin
    f1 = File("http://en.wikipedia.org/wiki/File:Nettlebed_cave.jpg", name="Nettlebed Cave", parent="syn1234567", synapseStore=false)
    @test f1.concreteType == "org.sagebionetworks.repo.model.FileEntity"
end

@testset "is_container" begin
    ## result from a Synapse entity annotation query
    ## Note: prefix may be capitalized or not, depending on the from clause of the query
    result = Dict("entity.versionNumber"=>1,
              "entity.nodeType"=>"project",
              "entity.concreteType"=>["org.sagebionetworks.repo.model.Project"],
              "entity.createdOn"=>1451512703905,
              "entity.id"=>"syn5570912",
              "entity.name"=>"blah")
    @test is_container(result)

    result = Dict("Entity.nodeType"=>"project",
              "Entity.id"=>"syn5570912",
              "Entity.name"=>"blah")
    @test is_container(result)

    result = Dict("entity.concreteType"=>["org.sagebionetworks.repo.model.Folder"],
              "entity.id"=>"syn5570914",
              "entity.name"=>"flapdoodle")
    @test is_container(result)

    result = Dict("File.concreteType"=>["org.sagebionetworks.repo.model.FileEntity"],
              "File.id"=>"syn5570914",
              "File.name"=>"flapdoodle")
    @test is_container(result) == false

    @test is_container(Folder("Stuff", parentId="syn12345"))
    @test is_container(Project("My Project", parentId="syn12345"))
    @test is_container(File("asdf.png", parentId="syn12345")) == false
end


text:
@testset "successive_f_tol" begin
    alg = GradientDescent(
        alphaguess = LineSearches.InitialStatic(),
        linesearch = LineSearches.Static(),
    )
    opt = Optim.Options(
        iterations = 10,
        successive_f_tol = 5,
        f_tol = 3,
        g_tol = -1,
    )
    result = Optim.optimize(
        sum,
        (y, _) -> fill!(y, 1),
        [0.0, 0.0],
        alg,
        opt,
    )
    @test result.iterations == opt.successive_f_tol + 1
end


text:
using Revise
using Plots
using Basins
using DynamicalSystems

# struct MagPendulum{T<:AbstractFloat}
#     magnets::Vector{SVector{2, T}}
# end
#
# function (m::MagPendulum)(u, p, t)
#     x, y, vx, vy = u
#     γ, d, α, ω = p
#     dx, dy = vx, vy
#     dvx, dvy = @. -ω^2*(x, y) - α*(vx, vy)
#     for ma in m.magnets
#         δx, δy = (x - ma[1]), (y - ma[2])
#         D = sqrt(δx^2 + δy^2 + d^2)
#         dvx -= γ*(x - ma[1])/D^5
#         dvy -= γ*(y - ma[2])/D^5
#     end
#     return SVector(dx, dy, dvx, dvy)
# end
#
# function mag_pendulum(u = [sincos(0.12553*2π)..., 0, 0];
#     γ = 1.0, d = 0.3, α = 0.2, ω = 0.5, N = 3)
#     m = MagPendulum([SVector(cos(2π*i/N), sin(2π*i/N)) for i in 1:N])
#     p = [γ, d, α, ω]
#     ds = ContinuousDynamicalSystem(m, u, p)
# end

# https://cdn.ima.org.uk/wp/wp-content/uploads/2020/03/Chaos-in-the-Magnetic-Pendulum-from-MT-April-2020.pdf
#ds = mag_pendulum(γ=1, d=0.5, α=0.175, ω=1., N=4)

ds = Systems.magnetic_pendulum(γ=1, d=0.3, α=0.2, ω=0.5, N=4)
integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9)
xg=range(-2,2,length=500)
yg=range(-2,2,length=500)
@time bsn = Basins.basins_general(xg, yg, integ; dt=1., idxs=1:2)

# Uncertainty exponent for these parameter and grid
bd = box_counting_dim(xg, yg, bsn)
α = 2 - bd

plot(xg,yg,bsn.basin',seriestype=:heatmap)


text:
"""
Main module of TinyModia.

* Developer: Hilding Elmqvist, Mogram AB
* First version: December 2020
* License: MIT (expat)

"""
module TinyModia

# Defalut switch settings
logStatistics = false
logExecution = false
logCalculations = false

useNewCodeGeneration = false

using Reexport

export instantiateModel, @instantiateModel, assert, stringifyDefinition

using Base.Meta: isexpr
using OrderedCollections: OrderedDict

using ModiaBase.Symbolic
using ModiaBase.Simplify
using ModiaBase.BLTandPantelidesUtilities
using ModiaBase.BLTandPantelides
using ModiaBase.Differentiate
using ModiaBase
export ModiaBase

@reexport using Unitful
using  Measurements
import MonteCarloMeasurements
using JSON
#using Profile
using TimerOutputs
using InteractiveUtils

global to = TimerOutput()

export stripUnit
"""
    stripUnit(v)

Convert the unit of `v` to the preferred units (default are the SI units),
and then strip the unit. For details see `upreferred` and `preferunits` in
[Unitful](https://painterqubits.github.io/Unitful.jl/stable/conversion/)

The function is defined as: `stripUnit(v) = ustrip(upreferred.(v))`.
"""
stripUnit(v) = ustrip.(upreferred.(v))

include("ModelCollections.jl")
include("EvaluateParameters.jl")
include("EventHandler.jl")
include("CodeGeneration.jl")
# include("GenerateGetDerivatives.jl")
include("Synchronous.jl")
include("SimulateAndPlot.jl")
include("ReverseDiffInterface.jl")
include("PathPlanning.jl")

# include("IncidencePlot.jl")
# using .IncidencePlot
# Base.Experimental.@optlevel 0

const drawIncidence = false

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory

const Version = "0.8.0-dev"
const Date = "2021-07-27"

#println(" \n\nWelcome to Modia - Dynamic MODeling and Simulation in julIA")
print(" \n\nWelcome to ")
printstyled("Tiny", color=:light_black)
print("Mod")
printstyled("ia", bold=true, color=:red)
print(" - ")
printstyled("Dynamic ", color=:light_black)
print("Mod")
printstyled("eling and Simulation with Jul", color=:light_black)
printstyled("ia", bold=true, color=:red)

println()
println("Version $Version ($Date)")

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

function assert(condition, message)
    if ! condition
        println("ASSERTION: ", message)
    end
end

function printArray(array, heading; order=1:length(array), log=false, number=true, indent=false, spacing=" ")
    if length(array) > 0 && log
        println(heading)
        n = 0
        for v in array
            n += 1
            if number
                print(lpad(order[n], 4), ": ")
            end
            if indent
                print(("|"*spacing)^(n-1))
            end
            println(replace(replace(string(v), "\n" => " "), "  " => " "))
        end
    end
end

function performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log=false)
    nUnknowns = length(unknowns) - length(states)
    nEquations = length(equations)

    if nEquations > 0 && nUnknowns != nEquations
        printstyled("The number of unknowns ($nUnknowns) and equations ($nEquations) is not equal!\n", bold=true, color=:red)
        printArray(parameters, "Parameters:", log=true)
        printArray(states, "Potential states:", log=true)
        printArray(setdiff(unknowns,states), "Unknowns:", log=true)
        printArray(equations, "Equations:", log=true)
    end

    ok = nUnknowns == nEquations

    # Check consistency
    EG::Array{Array{Int,1},1} = [G; buildExtendedSystem(Avar)]
    assignEG = matching(EG, length(Avar))
#        @assert all(assignEG .> 0) "  The DAE is structurally singular\n  $assignEG."
    assigned = [a > 0 for a in assignEG]
    if ! (all(assigned)) || ! ok
        if nUnknowns > nEquations
            missingEquations = nUnknowns - nEquations
            printstyled("\nThe DAE has no unique solution, because $missingEquations equation(s) missing!\n", bold=true, color=:red)
        elseif nEquations > nUnknowns
            tooManyEquations = nEquations - nUnknowns
            printstyled("\nThe DAE has no unique solution, because $tooManyEquations equation(s) too many\n", bold=true, color=:red)
        else
            printstyled("\nThe DAE has no unique solution, because it is structurally inconsistent\n" *
                        "  (the number of equations is identical to the number of unknowns)\n", bold=true, color=:red)
        end

        singular = true
        hequ = []
        # Create dummy equations for printing: integrate(x, der(x)) = 0
        for i in 1:length(Avar)
            a = Avar[i]
            if a > 0
                n = unknowns[i]
                der = unknowns[a]
                push!(hequ, "integrate($n, $der) = 0")
            end
        end
        equationsEG = vcat(equations, hequ)
        println("\nAssigned equations after consistency check:")
        printAssignedEquations(equationsEG, unknowns, 1:length(EG), assignEG, Avar, fill(0, length(EG)))
        unassigned = []
        for i in eachindex(assignEG)
            if assignEG[i] == 0
                push!(unassigned, unknowns[i])
            end
        end
        printArray(unassigned, "Not assigned variables after consistency check:", log=true)
        if nUnknowns > nEquations
            missingEquations = nUnknowns - nEquations
            error("Aborting, because $missingEquations equation(s) missing!")
        elseif nEquations > nUnknowns
            tooManyEquations = nEquations - nUnknowns
            error("Aborting, because $tooManyEquations equation(s) too many!")
        else
            error("Aborting, because DAE is structurally inconsistent")
        end
    elseif log
        println("  The DAE is structurally nonsingular.")
    end
end

function assignAndBLT(equations, unknowns, parameters, Avar, G, states, logDetails, log=false, logTiming=false)
    vActive = [a == 0 for a in Avar]

    if drawIncidence
        showIncidence(equations, unknowns, G, vActive, [], [], [], [[e] for e in 1:length(G)])
    end

    if false # logAssignmentBeforePantelides
        @show G length(Avar) vActive
        assign = matching(G, length(Avar), vActive)
        @show assign

        println("\nAssigned equations:")
        Bequ = fill(0, length(G))
        printAssignedEquations(equations, unknowns, 1:length(G), assign, Avar, Bequ)

        printUnassigned(equations, unknowns, assign, Avar, Bequ, vActive)
    end

    if logTiming
        println("\nCheck consistency of equations by matching extended equation set")
        @timeit to "performConsistencyCheck" performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log)
    else
        performConsistencyCheck(G, Avar, vActive, parameters, unknowns, states, equations, log)
    end

    if logTiming
        println("\nIndex reduction (Pantelides)")
        @timeit to "pantelides!" assign, Avar, Bequ = pantelides!(G, length(Avar), Avar)
    else
        assign, Avar, Bequ = pantelides!(G, length(Avar), Avar)
    end
    if logDetails
        @show assign Avar Bequ
    end

    moreVariables = length(Avar) > length(unknowns)
    unknowns = vcat(unknowns, fill(:x, length(Avar) - length(unknowns)))
    for i in 1:length(Avar)
        if Avar[i] > 0
            unknowns[Avar[i]] = derivative(unknowns[i], keys(parameters))
        end
    end

    moreEquations = length(Bequ) > length(equations)
    equations = vcat(equations, fill(:(a=b), length(Bequ) - length(equations)))
    for i in 1:length(Bequ)
        if Bequ[i] > 0
            if log
                println("Differentiating: ", equations[i])
            end
            equations[Bequ[i]] = simplify(derivative(equations[i], keys(parameters)))
            if log
                println("  Differentiated equation: ", equations[Bequ[i]])
            end
        end
    end

    if moreVariables && log
        printArray(unknowns, "Unknowns after index reduction:", log=logDetails)
    end
    if moreEquations && log
        printArray(equations, "Equations after index reduction:", log=logDetails)
    end

    HG = Array{Array{Int64,1},1}()
    originalEquationIndex = []
    highestDerivativeEquations = []
    for i in 1:length(G)
        if Bequ[i] == 0
            push!(HG, G[i])
            push!(highestDerivativeEquations, equations[i])
            push!(originalEquationIndex, i)
        end
    end

    vActive = [a == 0 for a in Avar]
    if logTiming
        println("Assign")
        @timeit to "matching" assign = matching(HG, length(Avar), vActive)
    else
        assign = matching(HG, length(Avar), vActive)
    end
    if logTiming
        println("BLT")
        @timeit to "BLT" bltComponents = BLT(HG, assign)
    else
       bltComponents = BLT(HG, assign)
    end
    if logDetails
        @show HG
        @show bltComponents
    end

    if log
        println("\nSorted highest derivative equations:")
        printSortedEquations(highestDerivativeEquations, unknowns, bltComponents, assign, Avar, Bequ)
    end

    if drawIncidence
        assignedVar, unAssigned = invertAssign(assign)
        showIncidence(equations, unknowns, HG, vActive, [], assign, assignedVar, bltComponents, sortVariables=true)
    end

    for b in bltComponents
        if length(b) > 1 && logDetails
            println("\nTearing of BLT block")
            @show b
            es::Array{Int64,1} = b
            assignedVar, unAssigned = invertAssign(assign)
            vs = assignedVar[b]
            td = TearingSetup(HG, length(unknowns))
            (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, (e,v) -> v in HG[e], es, vs)
            @show vTear eResidue vSolved eSolved
            printArray(unknowns[vTear], "torn variables:", log=logDetails)
            printArray(equations[eResidue], "residue equations:", log=logDetails)
            printArray(unknowns[vSolved], "solved variables:", log=logDetails)
            printArray(equations[eSolved], "solved equations:", log=logDetails)
        end
    end

    blt = Vector{Int}[]  # TODO: Change type
    for c in bltComponents
        push!(blt, originalEquationIndex[c])
    end
    for i in eachindex(assign)
        if assign[i] != 0
            assign[i] = originalEquationIndex[assign[i]]
        end
    end
    return (unknowns, equations, G, Avar, Bequ, assign, blt, parameters)
end

function setAvar(unknowns)
    variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))
    Avar = fill(0, length(unknowns))
    states = []
    derivatives = []
    for i in 1:length(unknowns)
        v = unknowns[i]
        if isexpr(v, :call)
            push!(derivatives, v)
            x = v.args[2]
            ix = variablesIndices[x]
            Avar[ix] = i
            push!(states, unknowns[ix])
        end
    end
    return Avar, states, derivatives
end

mutable struct ModelStructure
    parameters::OrderedDict{Any,Any}
    mappedParameters::OrderedDict{Any,Any}
    init::OrderedDict{Any,Any}
    start::OrderedDict{Any,Any}
    variables::OrderedDict{Any,Any}
    flows::OrderedDict{Any,Any}
    inputs::OrderedDict{Any,Any}
    outputs::OrderedDict{Any,Any}
#    components
#    extends
    equations::Array{Expr,1}
end

ModelStructure() = ModelStructure(OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), Expr[])

function printDict(label, d)
    if length(d) > 0
        print(label)
        for (k,v) in d
            print(k)
            if v != nothing
                print(" => ", v)
            end
            print(", ")
        end
        println()
    end
end

function printModelStructure(m; label="", log=false)
    if log
        println("ModelStructure ", label)
        printDict("  parameters:  ", m.parameters)
        printDict("  mappedParameters:  ", m.mappedParameters)
        printDict("  init:        ", m.init)
        printDict("  start:       ", m.start)
        printDict("  variables:   ", m.variables)
        printDict("  flows:       ", m.flows)
        printDict("  inputs:      ", m.inputs)
        printDict("  outputs:     ", m.outputs)
        println("  equations: ")
        for e in m.equations
            println("    ", e)
        end
    end
end

prependDict(dict, prefix) = OrderedDict([prepend(k, prefix) => prepend(v, prefix) for (k,v) in dict])

function mergeModelStructures(parent::ModelStructure, child::ModelStructure, prefix)
    merge!(parent.parameters, child.parameters)
    parent.mappedParameters[prefix] = child.mappedParameters

    merge!(parent.init, child.init)
    parent.mappedParameters[prefix] = child.mappedParameters

    merge!(parent.start, child.start)
    parent.mappedParameters[prefix] = child.mappedParameters

    merge!(parent.variables, child.variables)
    merge!(parent.flows, child.flows)

    push!(parent.equations, prepend(child.equations, prefix)...)
end

function replaceLinearIntegerEquations(Gint, eInt, GcInt, unknowns, equations)
    for i = 1:length(Gint)
        e = Gint[i]
        if length(e) > 0
            # Construct equation
            rhs = 0
            for j = 1:length(e)
                v  = e[j]
                vc = GcInt[i][j]
                rhs = add(rhs, mult(vc, unknowns[v]))
            end
            equ = :(0 = $rhs)
            equations[eInt[i]] = equ
        else
            equations[eInt[i]] = :(0 = 0)
        end
    end
end

function performAliasReduction(unknowns, equations, Avar, logDetails, log)

    @timeit to "enumerate(unknowns)" variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))

    linearVariables = []
    nonlinearVariables = []
    linearEquations = Vector{Int}()
    nonlinearEquations = []
    G = Vector{Int}[] # Array{Array{Int64,1},1}()
    Gint = Array{Array{Int64,1},1}()
    GcInt = Array{Array{Int64,1},1}()

    @timeit to "build graph information" for i in 1:length(equations)
        e = equations[i]
        nameIncidence, coefficients, rest, linear = getCoefficients(e)
        incidence = [variablesIndices[n] for n in nameIncidence if n in keys(variablesIndices)]
        unique!(incidence)
        push!(G, incidence)

        if linear && all([c == 1 || c == -1 for c in coefficients]) && rest == 0 && all(n in keys(variablesIndices) for n in nameIncidence)
            push!(linearVariables, nameIncidence...)
            push!(linearEquations, i)
            push!(Gint, incidence)
            push!(GcInt, coefficients)
        else
            push!(nonlinearVariables, nameIncidence...)
            push!(nonlinearEquations, e)
        end
    end

    if logDetails
        @show G
        @show Avar
    end

    @timeit to "unique!(linearVariables)" unique!(linearVariables)
    @timeit to "unique!(nonlinearVariables)" unique!(nonlinearVariables)
    @timeit to "setdiff" vSolveInLinearEquations = setdiff(linearVariables, nonlinearVariables)
    if logDetails
        printArray(linearVariables, "linearVariables:", number=false)
        printArray(equations[linearEquations], "linearEquations:", number=false)
#        printArray(nonlinearVariables, "nonlinearVariables:", number=false)
#        printArray(nonlinearEquations, "nonlinearEquations:", number=false)
        @show linearEquations vSolveInLinearEquations Gint GcInt
    end

    GCopy = deepcopy(G)
    AvarCopy = deepcopy(Avar)
    @timeit to "simplifyLinearIntegerEquations!" (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(GCopy, linearEquations, GcInt, AvarCopy)

    if logDetails
        @show vEliminated vProperty nvArbitrary redundantEquations
    end
    if log
        printSimplifiedLinearIntegerEquations(GCopy, linearEquations, GcInt, vEliminated, vProperty, nvArbitrary, redundantEquations, v->string(unknowns[v]), printTest=false)
    end

    substitutions = OrderedDict()
    @timeit to "build substitutions" if length(vEliminated) - nvArbitrary > 0
        for v in vEliminated[nvArbitrary+1:end]
            if isZero(vProperty, v)
                substitutions[unknowns[v]] = 0.0
            elseif isAlias(vProperty, v)
                substitutions[unknowns[v]] = unknowns[alias(vProperty,v)]
            else
                substitutions[unknowns[v]] = :(- $(unknowns[negAlias(vProperty,v)]))
            end
        end
    end

    @timeit to "replaceLinearIntegerEquations" replaceLinearIntegerEquations(GCopy[linearEquations], linearEquations, GcInt, unknowns, equations)
#    printArray(equations, "new equations:", log=log)

    reducedEquations = []
    @timeit to "substitute" for ei in 1:length(equations)
        e = equations[ei]
        if e != :(0 = 0)
            push!(reducedEquations, substitute(substitutions, e))
        end
    end

    reducedUnknowns = []
    for vi in 1:length(unknowns)
        if isNotEliminated(vProperty, vi)
            push!(reducedUnknowns, unknowns[vi])
        end
    end

    reducedVariablesIndices = OrderedDict(key => k for (k, key) in enumerate(reducedUnknowns))
    reducedG = Vector{Int}[] # Array{Array{Int64,1},1}()
    @timeit to "build reducedG" for i in 1:length(reducedEquations)
        e = reducedEquations[i]
        nameIncidence, coefficients, rest, linear = getCoefficients(e)
        incidence = [reducedVariablesIndices[n] for n in nameIncidence if n in keys(reducedVariablesIndices)]
        unique!(incidence)
        push!(reducedG, incidence)
    end

    reducedAvar, reducedStates, reducedDerivatives = setAvar(reducedUnknowns)
    return reducedEquations, reducedUnknowns, reducedAvar, reducedG, reducedStates, vEliminated, vProperty
end

function stateSelectionAndCodeGeneration(modStructure, name, modelModule, FloatType, init, start, inputs, outputs, vEliminated, vProperty, unknownsWithEliminated, mappedParameters;
    unitless=false, logStateSelection=false, logCode=false, logExecution=false, logCalculations=false, logTiming=false, evaluateParameters=false)
    (unknowns, equations, G, Avar, Bequ, assign, blt, parameters) = modStructure

    function getSolvedEquationAST(e, v)
        (solution, solved) = solveEquation(equations[e], unknowns[v])
        unknown = solution.args[1]
        expr = solution.args[2]
###        solution = :(_variables[$(v+1)] = $unknown = $expr)
        equ = string(equations[e])
        @assert solved "Equation not solved for $(unknowns[v]): $(equations[e])"
#        sol = string(solution)
#        solution = prepend(makeDerivativeVar(solution, components), :m)
        solution = substituteForEvents(solution)
        if false
    #        solution = :(begin println("Executing: ", $sol); $solution end)
            if sol != equ
                solution = :(try $solution; catch e; println("\n\nFailure executing: ", $sol); println("Original equation: ", $equ, "\n\n"); rethrow(e) end)
            else
                solution = :(try $solution; catch e; println("\n\nFailure executing: ", $sol, "\n\n"); rethrow(e) end)
            end
    #        solution = :(try $solution; catch e; println("Failure executing: ", $sol); println(sprint(showerror,e)); rethrow(e) end)
    #        solution = :(try $solution; catch e; println("Failure executing: ", $sol); printstyled(stderr,"ERROR: ", bold=true, color=:red);
    #        printstyled(stderr,sprint(showerror,e), color=:light_red); println(stderr); end)
        end
        solution = makeDerVar(solution, parameters, inputs, evaluateParameters)
        if logCalculations
            var = string(unknowns[v])
            solutionString = string(solution)
            return :(println("Calculating: ", $solutionString); $solution; println("  Result: ", $var, " = ", upreferred.($(solution.args[1]))))
        else
            return solution
        end
    end

    function getResidualEquationAST(e, residualName)
        eq = equations[e] # prepend(makeDerivativeVar(equations[e], components), :m)
        eqs = sub(eq.args[2], eq.args[1])
        resid = makeDerVar(:(ustrip($eqs)), parameters, inputs, evaluateParameters)
        residual = :($residualName = $resid)
        residString = string(resid)
        if logCalculations
            return :(println("Calculating residual: ", $residString); $residualName = $resid; println("  Residual: ", $residualName) )
#            return makeDerVar(:(dump($(makeDerVar(eq.args[2]))); dump($(makeDerVar(eq.args[1]))); $residual; println($residualName, " = ", upreferred.(($(eq.args[2]) - $(eq.args[1]))))))
        else
            return residual
        end
    end

    function isLinearEquation(e_original, v_original)
        equ = equations[e_original]
        var = unknowns[v_original]
        return isLinear(equ, var)
    end

    function isSolvableEquation(e_original, v_original)
        equ = equations[e_original]
        var = unknowns[v_original]
        (solution, solved) = solveEquation(equ, var)
        return solved
    end

    hasParticles(value) = typeof(value) <: MonteCarloMeasurements.StaticParticles ||
                          typeof(value) <: MonteCarloMeasurements.Particles

    function var_unit(v)
        var = unknowns[v]
        if var in keys(init)
            value = eval(init[var])
        elseif var in keys(start)
            value = eval(start[var])
        else
            value = 0.0
        end
        if hasParticles(value)  # Units not yet support for particles
            return ""
        end
        # if length(value) == 1
        if ! (typeof(value) <: Array)
            un = unit(value)
        else
            un = unit.(value)   # un = [unit(v) for v in value]  # unit.(value) does not work for MonteCarloMeasurements
            @assert all([un[i] == un[1] for i in 2:length(un)]) "The unit of all elements of state vector must be equal: $var::$(value)"
            un = un[1]
        end
        return replace(string(un), " " => "*") # Fix since Unitful removes * in unit strings
    end

    function var_has_start(v_original)
        return unknowns[v_original] in keys(init) || unknowns[v_original] in keys(start)
    end

    function var_fixed(v_original)
        return unknowns[v_original] in keys(init)
    end

    function var_length(v_original)
        v = unknowns[v_original]
        if v in keys(init)
            value = eval(init[v])
            len = hasParticles(value) ? 1 : length(value)
        elseif v in keys(start)
            value = eval(start[v])
            len = hasParticles(value) ? 1 : length(value)
        else
            len = 1
        end
        return len
    end

    function eq_equation(e)
        return replace(replace(string(equations[e]), "\n" => " "), "  " => " ")
    end

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

    solvedAST = []
    Gsolvable = copy(G)
    juliaVariables  = [Symbol(u) for u in unknowns]
    stringVariables = [String(Symbol(v)) for v in unknowns]
    allEquations    = equations
    #    @show stringVariables equations G blt assign Avar Bequ Gsolvable

    stateSelectionFunctions = StateSelectionFunctions(
        var_name               = v -> stringVariables[v],
        var_julia_name         = v -> juliaVariables[v],
        var_unit               = var_unit,
        var_has_start          = var_has_start,
        var_fixed              = var_fixed,
        var_nominal            = v_original -> NaN,
        var_unbounded          = v_original -> false,
        var_length             = var_length,
        equation               = eq_equation,
        isSolvableEquation     = isSolvableEquation,
        isLinearEquation       = isLinearEquation,
        getSolvedEquationAST   = getSolvedEquationAST,
        getResidualEquationAST = getResidualEquationAST,
        showMessage            = (message;severity=0,from="???",details="",variables=Int[],equations=Int[]) ->
                                    ModiaBase.showMessage(message, severity, from, details,
                                                          stringVariables[variables],
                                                          string.(allEquations[equations]))
        )

    if length(equations) == 0
        return nothing
    end

    if logTiming
        println("Get sorted and solved AST")
        @timeit to "getSortedAndSolvedAST" (success,AST,equationInfo) = getSortedAndSolvedAST(
            G, Array{Array{Int64,1},1}(blt), assign, Avar, Bequ, stateSelectionFunctions;
            unitless, log = logStateSelection, modelName = name)
    else
        (success,AST,equationInfo) = getSortedAndSolvedAST(
            G, Array{Array{Int64,1},1}(blt), assign, Avar, Bequ, stateSelectionFunctions;
            unitless, log = logStateSelection, modelName = name)
    end

    if ! success
        error("Aborting")
    end

    if useNewCodeGeneration
        AST, newFunctions = prepareGenerateGetDerivatives(AST, modelModule, functionSize, juliaVariables)
    else
        newFunctions = Vector{Expr}()
    end

#    @show AST
#    @show equationInfo

    selectedStates = [v.x_name_julia for v in equationInfo.x_info]

    startValues = []
    for v in selectedStates
        v = Meta.parse(string(v))
        if v in keys(init)
            value = eval(init[v])
#            value = [0.0 0.0]  # TODO: Fix
#            @show value
            if hasParticles(value)
                push!(startValues, value)
            else
                push!(startValues, value...)
            end
        elseif v in keys(start)
            value = eval(start[v])
            if hasParticles(value)
                push!(startValues, value)
            else
                push!(startValues, value...)
            end
        else
            push!(startValues, 0.0)
        end
    end
    if logCode
        println("startValues = ", startValues)
    end

    vSolvedWithInit = equationInfo.vSolvedWithFixedTrue
    vSolvedWithInitValuesAndUnit = OrderedDict{String,Any}()
    for name in vSolvedWithInit
        nameAsExpr = Meta.parse(name)
        if haskey(init, nameAsExpr)
            vSolvedWithInitValuesAndUnit[name] = eval( init[nameAsExpr] )
        else
            @warn "Internal issue of TinyModia: $name is assumed to have an init-value, but it is not found."
        end
    end

    # Generate code
    nCrossingFunctions, nAfter, nClocks, nSamples, previousVars, preVars, holdVars = getEventCounters()
    previousVars = Symbol.(previousVars)
    preVars = Symbol.(preVars)
    holdVars = Symbol.(holdVars)

    # Variables added to result
    extraResults = vcat(:time, setdiff([Symbol(u) for u in unknowns],
                                        Symbol[Symbol(xi_info.x_name_julia)     for xi_info in equationInfo.x_info],
                                        Symbol[Symbol(xi_info.der_x_name_julia) for xi_info in equationInfo.x_info]))

    if true # logTiming
#        println("Generate code")
        if useNewCodeGeneration
            @timeit to "generate_getDerivativesNew!" code = generate_getDerivativesNew!(AST, newFunctions, modelModule, equationInfo, [:(_p)], vcat(:time, [Symbol(u) for u in unknowns]), previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
        else
            @timeit to "generate_getDerivatives!" code = generate_getDerivatives!(AST, equationInfo, [:(_p)], extraResults, previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
        end
    else
#        code = generate_getDerivatives!(AST, equationInfo, Symbol.(keys(parameters)), extraResults, :getDerivatives, hasUnits = !unitless)
        code = generate_getDerivatives!(AST, newFunctions, modelModule, equationInfo, [:(_p)], extraResults, previousVars, preVars, holdVars, :getDerivatives, hasUnits = !unitless)
    end
    if logCode
        @show mappedParameters
        showCodeWithoutComments(code)
    end

    # Compile code

#    generatedFunction = @RuntimeGeneratedFunction(modelModule, code)
    #getDerivatives = Core.eval(modelModule, code)

    if logTiming
        println("eval code")
        @time @timeit to "eval(code)" Core.eval(modelModule, code)
    else
        Core.eval(modelModule, code)
    end
    getDerivatives = modelModule.getDerivatives

    # If generatedFunction is not packed inside a function, DifferentialEquations.jl crashes
#    getDerivatives(derx,x,m,time) = generatedFunction(derx, x, m, time)

    convertedStartValues = convert(Vector{FloatType}, [ustrip(v) for v in startValues])  # ustrip.(value) does not work for MonteCarloMeasurements
#    @show mappedParameters

#    println("Build SimulationModel")

    model = @timeit to "build SimulationModel" SimulationModel{FloatType, OrderedDict{Symbol,Any}}(modelModule, name, getDerivatives, equationInfo, convertedStartValues, previousVars, preVars, holdVars,
#                                         parameters, vcat(:time, [Symbol(u) for u in unknowns]);
                                         OrderedDict(:(_p) => mappedParameters ), extraResults;
                                         vSolvedWithInitValuesAndUnit, vEliminated, vProperty,
                                         var_name = (v)->string(unknownsWithEliminated[v]),
                                         nz=nCrossingFunctions, nAfter=nAfter,  unitless=unitless)

    # Execute code
    if logExecution
        println("\nExecute getDerivatives")
    #    @show startValues
        derx = deepcopy(convertedStartValues) # To get the same type as for x (deepcopy is needed for MonteCarloMeasurements)
        println("First executions of getDerivatives")
        @timeit to "execute getDerivatives" try
            TimeType = timeType(model)
            @time Base.invokelatest(getDerivatives, derx, convertedStartValues, model, convert(TimeType, 0.0))
            @time Base.invokelatest(getDerivatives, derx, convertedStartValues, model, convert(TimeType, 0.0))
#            @show derx
        catch e
            error("Failed: ", e)
            return nothing
        end
    end
    return model
end

"""
    modelInstance = @instantiateModel(model; FloatType = Float64, aliasReduction=true, unitless=false,
        log=false, logModel=false, logDetails=false, logStateSelection=false, logCode=false,
        logExecution=logExecution, logCalculations=logCalculations, logTiming=false)

Instantiates a model, i.e. performs structural and symbolic transformations and generates a function for calculation of derivatives suitable for simulation.

* `model`: model (declarations and equations)
* `FloatType`: Variable type for floating point numbers, for example: Float64, Measurements{Float64}, StaticParticles{Float64,100}, Particles{Float64,2000}
* `aliasReduction`: Perform alias elimination and remove singularities
* `unitless`: Remove units (useful while debugging models and needed for MonteCarloMeasurements)
* `log`: Log the different phases of translation
* `logModel`: Log the variables and equations of the model
* `logDetails`: Log internal data during the different phases of translation
* `logStateSelection`: Log details during state selection
* `logCode`: Log the generated code
* `logExecution`: Log the execution of the generated code (useful for timing compilation)
* `logCalculations`: Log the calculations of the generated code (useful for finding unit bugs)
* `logTiming`: Log timing of different phases
* `return modelInstance prepared for simulation`
"""
macro instantiateModel(model, kwargs...)
    modelName = string(model)
    code = :( instantiateModel($model; modelName=$modelName, modelModule=@__MODULE__, source=@__FILE__, $(kwargs...) ) )
    return esc(code)
end

"""
See documentation of macro @instatiateModel
"""
function instantiateModel(model; modelName="", modelModule=nothing, source=nothing, FloatType = Float64, aliasReduction=true, unitless=false,
    log=false, logModel=false, logDetails=false, logStateSelection=false, logCode=false,
    logExecution=logExecution, logCalculations=logCalculations, logTiming=false, evaluateParameters=false)
    try
        println("\nInstantiating model $modelModule.$modelName")
        resetEventCounters()
        global to = TimerOutput()

        modelStructure = ModelStructure()

        if isexpr(model, :quote)
            model = eval(model) # model defined with macro
        end

        if typeof(model) <: NamedTuple || typeof(model) <: Dict || typeof(model) <: OrderedDict
            if logModel
                @showModel(model)
            end
            if false
                println("drawModel")
                jsonDiagram = drawModel(modelName, model)
                println()
#                println("jsonModel")
#                JSON.print(jsonDiagram, 2)
            end

            if logTiming
                println("Flatten")
                @timeit to "flatten" flattenModelTuple!(model, modelStructure, modelName, to; unitless, log)
            else
                flattenModelTuple!(model, modelStructure, modelName, to; unitless, log)
            end
            pars = OrderedDict{Symbol, Any}([(Symbol(k),v) for (k,v) in modelStructure.parameters])
            if length(modelStructure.equations) > 0
                flatModel = Model() | pars | Model(equations = [removeBlock(e) for e in modelStructure.equations])
            else
                flatModel = Model() | pars
            end
            printModelStructure(modelStructure, log=false)
            name = modelName
        else
            @show model
            error("Invalid model format")
        end

        allVariables = Incidence[]
        if logTiming
            println("Find incidence")
            @timeit to "findIncidence!" findIncidence!(modelStructure.equations, allVariables)
        else
            findIncidence!(modelStructure.equations, allVariables)
        end
        unique!(allVariables)

        unknowns = setdiff(allVariables, keys(modelStructure.parameters), keys(modelStructure.inputs), [:time, :instantiatedModel, :_leq_mode, :_x])
        Avar, states, derivatives = setAvar(unknowns)
        vActive = [a == 0 for a in Avar]

        if logStatistics
            println("Number of states:    ", length(states))
            if length(unknowns)-length(states) != length(modelStructure.equations)
                println("Number of unknowns:  ", length(unknowns)-length(states))
            end
            println("Number of equations: ", length(modelStructure.equations))
        end
        printArray(modelStructure.parameters, "Parameters:", log=log)
        printArray(modelStructure.inputs, "Inputs:", log=log)
        printArray(modelStructure.outputs, "Outputs:", log=log)
        printArray(states, "Potential states:", log=log)
        printArray(setdiff(unknowns, states), "Unknowns:", log=log)
        printArray(modelStructure.equations, "Equations:", log=log)
        if logDetails
            @show modelStructure.parameters modelStructure.mappedParameters modelStructure.init modelStructure.start modelStructure.variables modelStructure.flows modelStructure.inputs modelStructure.outputs
        end

        unknownsWithEliminated = unknowns
        if aliasReduction
            if log || logTiming
                println("Perform alias elimination and remove singularities")
            end
            if logTiming
                @timeit to "performAliasReduction" equations, unknowns, Avar, G, states, vEliminated, vProperty = performAliasReduction(unknowns, modelStructure.equations, Avar, logDetails, log)
                println("Number of reduced unknowns:  ", length(unknowns)-length(states))
                println("Number of reduced equations: ", length(equations))
            else
                equations, unknowns, Avar, G, states, vEliminated, vProperty = performAliasReduction(unknowns, modelStructure.equations, Avar, logDetails, log)
            end
            printArray(states, "States:", log=log)
            printArray(setdiff(unknowns, states), "Unknowns after alias reduction:", log=log)
            printArray(equations, "Equations after alias reduction:", log=log)
        else
            equations = modelStructure.equations
            G = Vector{Int}[] # Array{Array{Int64,1},1}()
            variablesIndices = OrderedDict(key => k for (k, key) in enumerate(unknowns))

            println("Build incidence matrix")
            unknownsSet = Set(unknowns)
            @timeit to "build incidence matrix" for i in 1:length(equations)
                e = equations[i]
                nameIncidence, coefficients, rest, linear = getCoefficients(e)

                incidence = [] # [variablesIndices[n] for n in nameIncidence if n in unknownsSet]
                for n in nameIncidence
                    if n in keys(variablesIndices)
                        push!(incidence, variablesIndices[n])
                    end
                end
                unique!(incidence)
                push!(G, incidence)
            end
            vEliminated = Int[]
            vProperty   = Int[]
        end

        modStructure = assignAndBLT(equations, unknowns, modelStructure.parameters, Avar, G, states, logDetails, log, logTiming)

        inst = stateSelectionAndCodeGeneration(modStructure, name, modelModule, FloatType, modelStructure.init, modelStructure.start, modelStructure.inputs, modelStructure.outputs,
            vEliminated, vProperty, unknownsWithEliminated, modelStructure.mappedParameters;
            unitless, logStateSelection, logCode, logExecution, logCalculations, logTiming, evaluateParameters)

        if logTiming
            show(to, compact=true)
            println()
        end

        inst #, flatModel

    catch e
        if isa(e, ErrorException)
            println()
            printstyled("Model error: ", bold=true, color=:red)
            printstyled(e.msg, "\n", bold=true, color=:red)
            printstyled("Aborting instantiateModel for $modelName in $modelModule\n", bold=true, color=:red)
            println()
#            Base.rethrow()
        else
            Base.rethrow()
        end
        show(to, compact=true)
        println()
    end

end

end


text:
export
    Problem,
    change_integration

"""$(TYPEDEF) Trajectory Optimization Problem.
Contains the full definition of a trajectory optimization problem, including:
* dynamics model (`Model`)
* objective (`Objective`)
* constraints (`ConstraintSet`)
* initial and final states
* Primal variables (state and control trajectories)
* Discretization information: knot points (`N`), time step (`dt`), and total time (`tf`)

# Constructors:
```julia
Problem(model, obj, constraints, x0, xf, Z, N, tf) # defaults to RK3 integration
Problem{Q}(model, obj, constraints, x0, xf, Z, N, tf) where Q<:QuadratureRule
Problem(model, obj, xf, tf; x0, constraints, N, X0, U0, dt, integration)
Problem{Q}(prob::Problem)  # change integration
```
where `Z` is a trajectory (Vector of `KnotPoint`s)

# Arguments
* `model`: Dynamics model. Can be either `Discrete` or `Continuous`
* `obj`: Objective
* `X0`: Initial state trajectory. If omitted it will be initialized with NaNs, to be later overwritten by the solver.
* `U0`: Initial control trajectory. If omitted it will be initialized with zeros.
* `x0`: Initial state. Defaults to zeros.
* `xf`: Final state. Defaults to zeros.
* `dt`: Time step
* `tf`: Final time. Set to zero to specify a time penalized problem.
* `N`: Number of knot points. Defaults to 51, unless specified by `dt` and `tf`.
* `integration`: One of the defined integration types to discretize the continuous dynamics model.
Both `X0` and `U0` can be either a `Matrix` or a `Vector{Vector}`, but must be the same.
At least 2 of `dt`, `tf`, and `N` need to be specified (or just 1 of `dt` and `tf`).
"""
struct Problem{Q<:QuadratureRule,T<:AbstractFloat}
    model::AbstractModel
    obj::AbstractObjective
    constraints::ConstraintSet{T}
    x0::SVector
    xf::SVector
    Z::Traj
    N::Int
    t0::T
    tf::T
    function Problem{Q}(model::AbstractModel, obj::AbstractObjective,
            constraints::ConstraintSet,
            x0::SVector, xf::SVector,
            Z::Traj, N::Int, t0::T, tf::T) where {Q,T}
        n,m = size(model)
        @assert length(x0) == length(xf) == n
        @assert length(Z) == N
        @assert tf > t0
        new{Q,T}(model, obj, constraints, x0, xf, Z, N, t0, tf)
    end
end

"Use RK3 as default integration"
Problem(model, obj, constraints, x0, xf, Z, N, t0, tf) =
    Problem{RK3}(model, obj, constraints, x0, xf, Z, N, t0, tf)

function Problem(model::L, obj::O, xf::AbstractVector, tf;
        constraints=ConstraintSet(size(model)...,length(obj)),
        t0=zero(tf),
        x0=zero(xf), N::Int=length(obj),
        X0=[x0*NaN for k = 1:N],
        U0=[@SVector zeros(size(model)[2]) for k = 1:N-1],
        dt=fill((tf-t0)/(N-1),N),
        integration=DEFAULT_Q) where {L,O}
    n,m = size(model)
    if dt isa Real
        dt = fill(dt,N)
    end
    if X0 isa AbstractMatrix
        X0 = [X0[:,k] for k = 1:size(X0,2)]
    end
    if U0 isa AbstractMatrix
        U0 = [U0[:,k] for k = 1:size(U0,2)]
    end
    t = range(t0, tf, length=N)
    Z = Traj(X0,U0,dt,t)

    Problem{integration}(model, obj, constraints, SVector{n}(x0), SVector{n}(xf),
        Z, N, t0, tf)
end

"$(TYPEDSIGNATURES)
Get number of states, controls, and knot points"
Base.size(prob::Problem) = size(prob.model)..., prob.N

"""```julia
integration(::Problem)
integration(::DynamicsConstraint)
```
Get the integration rule"""
integration(prob::Problem{Q}) where Q = Q

"```julia
controls(::Problem)
controls(::AbstractSolver)
controls(::Traj)
```
Get the control trajectory
"
controls(prob::Problem) = controls(prob.Z)

"```julia
states(::Problem)
states(::AbstractSolver)
states(::Traj)
```
Get the state trajectory
"
states(prob::Problem) = states(prob.Z)

"```julia
initial_trajectory!(::Problem, Z)
initial_trajectory!(::AbstractSolver, Z)
```
Copy the trajectory "
function initial_trajectory!(prob::Problem, Z::Traj)
    for k = 1:prob.N
        prob.Z[k].z = Z[k].z
    end
end

"```julia
initial_states!(::Union{Problem,AbstractSolver}, X0::Vector{<:AbstractVector})
initial_states!(::Union{Problem,AbstractSolver}, X0::AbstractMatrix)
```
Copy the state trajectory "
function initial_states!(prob::Problem, X0::Vector{<:AbstractVector})
    set_states!(prob.Z, X0)
end

function initial_states!(prob::Problem, X0::AbstractMatrix)
    X0 = [X0[:,k] for k = 1:size(X0,2)]
    set_states!(prob.Z, X0)
end

"```julia
initial_controls!(::Union{Problem,AbstractSolver}, U0::Vector{<:AbstractVector})
initial_controls!(::Union{Problem,AbstractSolver}, U0::AbstractMatrx)
```
Copy the control trajectory "
function initial_controls!(prob::Problem, U0::Vector{<:AbstractVector})
    set_controls!(prob.Z, U0)
end

function initial_controls!(prob::Problem, u0::AbstractVector{<:Real})
    U0 = [copy(u0) for k = 1:prob.N]
    initial_controls!(prob, U0)
end

"```julia
cost(::Problem)
cost(::AbstractSolver)
```
Compute the cost for the current trajectory"
@inline cost(prob::Problem) = cost(prob.obj, prob.Z)

"Copy the problem"
function copy(prob::Problem{Q}) where Q
    Problem{Q}(prob.model, copy(prob.obj), copy(prob.constraints), prob.x0, prob.xf,
        copy(prob.Z), prob.N, prob.t0, prob.tf)
end

function max_violation(prob::Problem, Z::Traj=prob.Z)
    conSet = get_constraints(prob)
    evaluate!(conSet, Z)
    max_violation!(conSet)
    return maximum(conSet.c_max)
end

num_constraints(prob::Problem) = get_constraints(prob).p

@inline get_constraints(prob::Problem) = prob.constraints

"```julia
change_integration(prob::Problem, Q<:QuadratureRule)
```
Change dynamics integration for the problem"
change_integration(prob::Problem, ::Type{Q}) where Q<:QuadratureRule =
    Problem{Q}(prob)

function Problem{Q}(p::Problem) where Q
    Problem{Q}(p.model, p.obj, p.constraints, p.x0, p.xf, p.Z, p.N, p.t0, p.tf)
end

@inline rollout!(prob::Problem) = rollout!(prob.model, prob.Z, prob.x0)

function Problem(p::Problem; model=p.model, obj=p.obj, constraints=p.constraints,
    x0=p.x0, xf=p.xf, t0=p.t0, tf=p.tf)
    Problem(model, obj, constraints, x0, xf, p.Z, p.N, t0, tf)
end


text:
using DataStructures
using Random
using BSON
using DataDeps

export fit!, predict

function pos_tagger_datadep_register()
    register(DataDep("POS Perceptron Tagger Weights",
        """
        The trained weights for the average Perceptron Tagger on Part of Speech Tagging task.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/pretrainedMod.bson.zip",
        "52519cb3aea5d8f74368faedea831471e5df34567de4748d15decea7424743d3",
        post_fetch_method = function(fn)
            unpack(fn)
            rm("__MACOSX", recursive=true)
            file = readdir()[1]
            mv(file, "POSWeights.bson")
        end
    ))
end

"""
This file contains the Average Perceptron model and Perceptron Tagger
which was original implemented by Matthew Honnibal.

The model learns by basic perceptron algorithm
but after all iterations weights are being averaged

AVERAGE PERCEPTRON MODEL

This struct contains the actual Average Perceptron Model
"""
mutable struct AveragePerceptron
    classes :: Set
    weights :: Dict
    _totals :: DefaultDict
    _tstamps :: DefaultDict
    i :: Int64
    START :: Array

    function AveragePerceptron()
        new(Set(), Dict(), DefaultDict(0), DefaultDict(0), 1, ["-START-", "-START2-"])
    end
end

"""
Predicting the class using current weights by doing Dot-product of features
and weights and return the scores
"""
function predict(self::AveragePerceptron, features)
    scores = DefaultDict(0.0)
    for (feature, value) in features
        if feature ∉ keys(self.weights)
            continue
        end
        weights = self.weights[feature]
        for (label, weight) in weights
            scores[label] += value * weight
        end
    end
    function custmax(scores)
        a = [scores[class] for class in self.classes]
        zipped = collect(zip(a, self.classes))
        currmax = zipped[1]
        for i=2:(length(zipped))
            currmax = max(currmax, zipped[i])
        end
        return currmax[2]
    end
    return custmax(scores)
end

"""
Applying the perceptron learning algorithm
Increment the truth weights and decrementing the guess weights,
if the guess is wrong
"""
function update(self::AveragePerceptron, truth, guess, features)
    function upd_feat(c, f, w, v)
        param = (f, c)
        n_iters_at_this_weight = self.i - self._tstamps[param]
        self._totals[param] += n_iters_at_this_weight * w
        self.weights[f][c] = w + v
        self._tstamps[param] = self.i
    end

    self.i += 1
    if truth === guess
        return nothing
    end
    for (f, value) in features
        if f in keys(self.weights)
            weights = self.weights[f]
        else
            self.weights[f] = Dict()
            weights = Dict()
        end
        upd_feat(truth, f, get(weights, truth, 0.0), 1.0)
        upd_feat(guess, f, get(weights, guess, 0.0), -1.0)
    end
    return nothing
end

"""
Averaging the weights over all time stamps
"""
function average_weights(self::AveragePerceptron)
    function newRound(fl, in)
        temp = fl*(10^in)
        return (float(round(temp))/(10^in))
    end
    for (feature , weights) in self.weights
        new_feat_weights = Dict()
        for (clas, weight) in weights
            param = (feature, clas)
            total = self._totals[param]
            total += (self.i - self._tstamps[param])*weight
            averaged = newRound(total/float(self.i-1), 3)
            if averaged != nothing
                new_feat_weights[clas] = averaged
            end
        self.weights[feature] = new_feat_weights
        end
    end
    return nothing
end

"""
# PERCEPTRON TAGGER

This struct contains the POS tagger "PerceptronTagger" which uses model in "AveragePerceptron"
In this training can be done and weights can be saved
Or a pretrain weights can be used (which are trained on same features)
and train more or can be used to predict

## To train:

```julia
julia> tagger = PerceptronTagger(false)

julia> fit!(tagger, [[("today","NN"),("is","VBZ"),("good","JJ"),("day","NN")]])
```

## To load pretrain model:

```julia
julia> tagger = PerceptronTagger(true)
```

## To predict tag:

```julia
julia> predict(tagger, ["today", "is"])
```
"""
mutable struct PerceptronTagger
    model :: AveragePerceptron
    tagdict :: Dict
    classes :: Set
    START :: Array
    END :: Array
    _sentences

        PerceptronTagger() = new(AveragePerceptron(), Dict(), Set(), ["-START-", "-START2-"], ["-END-", "-END2-"], [])
end

function PerceptronTagger(load::Bool)
    self = PerceptronTagger()

    # If load is true then a pretrain model will be import from location
    if load
        location = joinpath(datadep"POS Perceptron Tagger Weights", "POSWeights.bson")
        pretrained = BSON.load(location)
        self.model.weights = pretrained[:weights]
        self.tagdict = pretrained[:tagdict]
        self.classes = self.model.classes = Set(pretrained[:classes])
        println("loaded successfully")
    end

    return self
end

"""
makes a dictionary for single-tag words
params : sentences - an array of tuples which contains word and correspinding tag
"""
function makeTagDict(self::PerceptronTagger, sentences)
    counts = DefaultDict(()->DefaultDict(0))
    for sentence in sentences
        append!(self._sentences, sentences)
        for (word, tag) in sentence
            counts[word][tag] += 1
            push!(self.classes, tag)
        end
    end
    freq_thresh = 20
    ambiguity_thresh = 0.97
    for (word, tag_freqs) in counts
        mode, tag = findmax(collect(values(tag_freqs))); tag = collect(keys(tag_freqs))[tag]
        n = sum(values(tag_freqs))
        if (n >= freq_thresh) && ((mode/n) >= ambiguity_thresh)
            self.tagdict[word] = tag
        end
    end
end

"""
This function is used to normalize the given word
params : word - String
"""
function normalize(word)
    word = string(word)
    if occursin("-", word) && (word[1] != "-")
        return "!HYPHEN"
    elseif occursin(r"^[\d]{4}$", word)
        return "!YEAR"
    elseif occursin(r"^[\d]$", string(word[1]))
        return "!DIGITS"
    else
        return lowercase(word)
    end
end

"""
Converting the token into a feature representation, implemented as Dict
If the features change, a new model should be trained

# Arguments:

- `i` - index of word(or token) in sentence
- `word` - token
- `context` - array of tokens with starting and ending specifiers
- `prev` == "-START-" prev2 == "-START2-" - Start specifiers
"""
function getFeatures(self::PerceptronTagger, i, word, context, prev, prev2)
    function add(sep, name, args...)
        str = name
        for arg in args
            str *= sep * arg
        end
        if str in keys(features)
            features[str] += 1
        else
            features[str] = 1
        end
        return nothing
    end
    i += length(self.START)
    features = OrderedDict()
    add(" ", "bias")
    if length(word) >= 3
        add(" ", "i suffix", word[end-2:end])
    else
        add(" ", "i suffix", word)
    end
    add(" ", "i pref1", word[1])
    add(" ", "i-1 tag", prev)
    add(" ", "i-2 tag", prev2)
    add(" ", "i tag+i-2 tag", prev, prev2)
    add(" ", "i word", context[i])
    add(" ", "i-1 tag+i word", prev, context[i])
    add(" ", "i-1 word", context[i-1])
    if length(context[i-1]) >= 3
        add(" ", "i-1 suffix", context[i-1][end-2:end])
    else
        add(" ", "i-1 suffix", context[i-1])
    end
    add(" ", "i-2 word", context[i-2])
    add(" ", "i+1 word", context[i+1])
    if length(context[i+1]) >= 3
        add(" ", "i+1 suffix", context[i+1][end-2:end])
    else
        add(" ", "i+1 suffix", context[i+1])
    end
    add(" ", "i+2 word", context[i+2])
    return features
end

"""
    predict(::PerceptronTagger, tokens)
    predict(::PerceptronTagger, sentence)

Used for predicting the tags for given sentence or array of tokens
"""
function predict(self::PerceptronTagger, tokens::Vector{String})
    prev, prev2 = self.START
    output = []

    context = vcat(self.START, [normalize(word) for word in tokens], self.END)
    for (i, word) in enumerate(tokens)
        tag = get(self.tagdict, word, nothing)
        if tag === nothing
            features = getFeatures(self, i, word, context, prev, prev2)
            tag = predict(self.model, features)
        end
        push!(output, (word, tag))
        prev2 = prev
        prev = tag
    end
    return output
end

function (tagger::PerceptronTagger)(input)
    predict(tagger, input)
end

predict(tagger::PerceptronTagger, sentence::String) =
        predict(tagger, tokenize(Languages.English(), sentence))
predict(tagger::PerceptronTagger, sd::StringDocument) =
        predict(tagger, text(sd))
predict(tagger::PerceptronTagger, fd::FileDocument) =
        predict(tagger, text(fd))
predict(tagger::PerceptronTagger, td::TokenDocument) =
        predict(tagger, tokens(td))
function predict(tagger::PerceptronTagger, ngd::NGramDocument)
    @warn "POS tagging for NGramDocument not available."
end

"""
    fit!(::PerceptronTagger, sentences::Vector{Vector{Tuple{String, String}}}, save_loc::String, nr_iter::Integer)

Used for training a new model or can be used for training
an existing model by using pretrained weigths and classes

Contains main training loop for number of epochs.
After training weights, tagdict and classes are stored in the specified location.

# Arguments:
- `::PerceptronTagger` : Input PerceptronTagger model
- `sentences::Vector{Vector{Tuple{String, String}}}` : Array of the all token seqeunces with target POS tag
- `save_loc::String` : To specify the saving location
- `nr_iter::Integer` : Total number of training iterations for given sentences(or number of epochs)
"""
function fit!(self::PerceptronTagger, sentences::Vector{Vector{Tuple{String, String}}}, save_loc::String, nr_iter::Integer)
    self._sentences = []
    makeTagDict(self, sentences)
    self.model.classes = self.classes
    for iter=1:nr_iter
        c = 0; n = 0;
        for sentence in self._sentences
            words, tags = [x[1] for x in sentence], [x[2] for x in sentence]
            prev, prev2 = self.START
            context = vcat(self.START, [normalize(w) for w in words], self.END)
            for (i, word) in enumerate(words)
                guess = get(self.tagdict, word, nothing)
                if guess == nothing
                    feats = getFeatures(self, i, word, context, prev, prev2)
                    guess = predict(self.model, feats)
                    update(self.model, tags[i], guess, feats)
                end
                prev2 = prev
                prev = guess
                c += (guess == tags[i])
                n += 1
            end
        end
        shuffle(self._sentences)
        println("iteration : $iter")
    end
    self._sentences = nothing
    average_weights(self.model)

    if save_loc != ""
        bson(save_loc, weights = self.model.weights, tagdict = self.tagdict, classes = collect(self.classes))
    end
end

fit!(self::PerceptronTagger, sentences::Vector{Vector{Tuple{String, String}}}, nr_iter::Integer) = fit!(self::PerceptronTagger, sentences, "", nr_iter)
fit!(self::PerceptronTagger, sentences::Vector{Vector{Tuple{String, String}}}, save_loc::String) = fit!(self::PerceptronTagger, sentences, save_loc, 5)
fit!(self::PerceptronTagger, sentences::Vector{Vector{Tuple{String, String}}}) = fit!(self::PerceptronTagger, sentences, "", 5)


text:
#=

# Error handling

Error handling is achieved via a Rust-like mechanism through [ResultTypes.jl](https://github.com/iamed2/ResultTypes.jl). All Vulkan API functions that return a Vulkan [result code](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkResult.html) are wrapped into a `ResultTypes.Result`, which holds either the result or a [`VulkanError`](@ref) if a non-zero status code is encountered. Note that `ResultTypes.Result` is distinct from `Vulkan.Result` which is the wrapped version of `VkResult`. Errors can be manually handled with the following pattern

=#

using Vulkan

res = create_instance(["VK_LAYER_KHRONOS_validation"], [])
if iserror(res)
    err = unwrap_error(res)
    if err.code == ERROR_INCOMPATIBLE_DRIVER
        error("""
              No driver compatible with the requested API version could be found.
              Please make sure that a driver supporting Vulkan is installed, and
              that it is up to date with the requested version.
              """)
    elseif err.code == ERROR_LAYER_NOT_PRESENT
        error("Validation layers not available.")
    else
        throw(err)
    end
else # get the instance
    unwrap(res)
end

# Note that calling `unwrap` directly on the result will throw any contained `VulkanError` if there is one. So, if you just want to throw an exception when encountering an error, you can do

unwrap(create_instance([], []))

# Because it may be tedious to unwrap everything by hand and explicitly set the create info structures, [convenience constructors](@ref expose-create-info-args) are defined for handle types so that you can just do

Instance([], [])

#=

However, note that exceptions are thrown whenever the result is an error with this shorter approach.

Furthermore, all functions that may return non-success (but non-error) codes return a `Vulkan.Result` type along with any other returned value, since the return code may still be of interest.

For more details on the `ResultTypes.Result` type and how to handle it, please consult the [ResultTypes.jl documentation](https://iamed2.github.io/ResultTypes.jl/stable/).

=#

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using QuickStructs
using FactCheck

facts("CircularBuffer") do

  cb = CircularBuffer(Int, 5)
  @fact length(cb) => 0
  @fact capacity(cb) => 5
  @fact_throws first(cb)
  @fact isfull(cb) => false

  push!(cb, 1)
  @fact length(cb) => 1
  @fact capacity(cb) => 5
  @fact isfull(cb) => false

  append!(cb, 2:8)
  @fact length(cb) => capacity(cb)
  @fact isfull(cb) => true
  @fact toarray(cb) => Int[4,5,6,7,8]
  @fact cb[2] => 5
  @fact_throws cb[6]
  @fact_throws cb[3:6]
  @fact cb[3:4] => Int[6,7]
  @fact cb[[1,5]] => Int[4,8]

end

facts("RollingWindow") do

  w = RollingWindow(Int, 0:2:4)
  @fact length(w) => 0
  @fact capacity(w) => 3
  @fact_throws first(w)
  @fact isfull(w) => false

  push!(w, 1)
  @fact length(w) => 0
  @fact capacity(w) => 3
  @fact isfull(w) => false

  append!(w, 2:8)
  @fact length(w) => capacity(w)
  @fact isfull(w) => true
  @fact toarray(w) => Int[8,6,4]
  @fact w[2] => 6
  @fact_throws w[4]
  @fact_throws w[3:6]
  @fact w[2:3] => Int[6,4]
  @fact w[[1,3]] => Int[8,4]

end

FactCheck.exitstatus()


text:
package java.lang;

public class IllegalArgumentException extends RuntimeException {
    public IllegalArgumentException() { }
    public IllegalArgumentException(String s) { super(s); }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase