text:
struct adVAE
	encoder
	generator
	transformer
end

function adVAE(
		;idim::Int=1,
		zdim::Int=1,
		hdim::Int=64,
		activation = "relu",
		nlayers::Int=3,
		init_seed=nothing,
		kwargs...
		)
	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing

	encoder = Flux.Chain(
		build_mlp(idim, hdim, hdim, nlayers-1, activation=activation),
		ConditionalDists.SplitLayer(hdim, [zdim,zdim], [identity,softplus])
	)
	transformer = Flux.Chain(
		build_mlp(zdim*2, hdim, hdim, 2, activation=activation),
		ConditionalDists.SplitLayer(hdim, [zdim,zdim], [identity,softplus])
	)
	generator = build_mlp(zdim, hdim, idim, nlayers, activation=activation, lastlayer="linear")

	(init_seed !== nothing) ? Random.seed!() : nothing
	return adVAE(encoder, generator, transformer)
end

function Conv_adVAE(
		;isize::Int=1,
		in_ch::Int=1,
		zdim::Int=64,
		nf::Int=16,
		extra_layers::Int=0,
		depth::Int=3,
		activation="relu",
		init_seed=nothing,
		kwargs...
		)

	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing
	encoder = Flux.Chain(
		ConvEncoder(isize, in_ch, zdim, nf, extra_layers),
		x->reshape(x, (zdim,size(x)[end])), # equivalent to flatten
		ConditionalDists.SplitLayer(zdim, [zdim, zdim], [identity, softplus])
	)
	transformer = Flux.Chain(
		build_mlp(zdim*2, zdim, zdim, depth, activation=activation),
		ConditionalDists.SplitLayer(zdim, [zdim,zdim], [identity,softplus])
	)
	generator = Flux.Chain(
		Flux.Dense(zdim, zdim, eval(:($(Symbol(activation)))) ),
		x->reshape(x, (1,1,zdim,size(x)[end])), # inverse to flatten
		ConvDecoder(isize, zdim, in_ch, nf, extra_layers)
	)
	(init_seed !== nothing) ? Random.seed!() : nothing
	return adVAE(encoder, generator, transformer)
end

Flux.@functor adVAE

function (advae::adVAE)(x)
	μ, Σ = advae.encoder(x)
	z = μ + Σ * randn(Float32)
	xᵣ = advae.generator(z)
end

kl_divergence(μ, Σ) = - Flux.mean(0.5f0 * sum(1f0 .+ log.(Σ.^2) - μ.^2  - Σ.^2, dims=1))
kl_divergence(μ₁, Σ₁, μ₂, Σ₂) = Flux.mean(sum(log.(Σ₂) - log.(Σ₁) + (Σ₁.^2 + (μ₁ - μ₂).^2) ./ (2*Σ₂.^2) .- 0.5f0, dims=1))

function loss(advae::adVAE, x; γ=1e-3, λ=1e-2, mx=1, mz=1)
	μ, Σ = advae.encoder(x)
	μₜ, Σₜ = advae.transformer(cat(μ, Σ, dims=1))

	z = μ + Σ * randn(Float32)
	zₜ = μₜ + Σₜ * randn(Float32)

	xᵣ = advae.generator(z)
	xₜᵣ = advae.generator(zₜ)

	μᵣ, Σᵣ = advae.encoder(xᵣ)
	μₜᵣ, Σₜᵣ = advae.encoder(xₜᵣ)

	# 𝓛 for generator => 𝓛 = 𝓛_z + 𝓛_zₜ
	𝓛_z = Flux.Losses.mse(x, xᵣ) .+ γ * kl_divergence(μᵣ, Σᵣ)
	𝓛_zₜ = max(0f0, (mx - Flux.Losses.mse(xᵣ, xₜᵣ))) + γ * max(0f0, (mz - kl_divergence(μₜᵣ, Σₜᵣ)))
	𝓛 = 𝓛_z + 𝓛_zₜ
	# 𝓛ₜ for transformer
	𝓛ₜ = kl_divergence(μ, Σ, μₜ, Σₜ)
	# 𝓛ₑ for encoder
	𝓛ₑ = Flux.Losses.mse(x, xᵣ) .+ λ * kl_divergence(μ, Σ)
		+ γ * max(0f0, (mz - kl_divergence(μᵣ, Σᵣ)))
		+ γ * max(0f0, (mz - kl_divergence(μₜᵣ, Σₜᵣ)))

	return 𝓛 + λ*𝓛ₜ, 𝓛ₑ, 𝓛, 𝓛ₜ
end

function anomaly_score(advae::adVAE, real_input; L=100, dims=3, batch_size=64, to_testmode::Bool=true)
	real_input = Flux.Data.DataLoader(real_input, batchsize=batch_size)
	(to_testmode == true) ? Flux.testmode!(advae) : nothing
	advae = advae |> gpu
	output = Array{Float32}([])
	for x in real_input
        x = x |> gpu
		X = Array{Float32}(undef, size(x)[end], 1)
		μ, Σ = advae.encoder(x)
		for l=1:L
			z = μ + Σ * randn(Float32)
			xᵣ = advae.generator(z)
			rec_loss = vec(Flux.mse(x, xᵣ, agg=x->Flux.mean(x, dims=dims))) |> cpu # loss per batch
			X = cat(X, rec_loss, dims=2)
		end
		output = cat(output, Flux.mean(X, dims=2), dims=1)
	end
	(to_testmode == true) ? Flux.testmode!(advae, false) : nothing
	return output
end

function StatsBase.fit!(advae::adVAE, data, params)
	# prepare batches & loaders
	train_loader, val_loader = prepare_dataloaders(data, params)
	# training info logger
	history = MVHistory()
	# prepare for early stopping
	best_adVAE = deepcopy(advae)
	patience = params.patience
	best_val_loss = 1e10
	val_batches = length(val_loader)

	# ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) = Optimiser(ADAM(η, β), WeightDecay(decay))
	opt_step1 = haskey(params, :decay) ? ADAMW(params.lr, (0.9, 0.999), params.decay) : ADAM(params.lr)
	opt_step2 = haskey(params, :decay) ? ADAMW(params.lr, (0.9, 0.999), params.decay) : ADAM(params.lr)

	ps_step1 = Flux.params(advae.generator, advae.transformer)
	ps_step2 = Flux.params(advae.encoder)

	progress = Progress(length(train_loader))
	for (iter, X) in enumerate(train_loader)
		# training step 1
		loss_1, back = Flux.pullback(ps_step1) do
			loss(advae|>gpu , getobs(X)|>gpu, γ=params.gamma, λ=params.lambda, mx=params.mx, mz=params.mz)
		end
		grad_step1 = back((1f0, 0f0, 0f0, 0f0))
		Flux.Optimise.update!(opt_step1, ps_step1, grad_step1)
		# training step 2
		loss_2, back = Flux.pullback(ps_step2) do
			loss(advae|>gpu , getobs(X)|>gpu, γ=params.gamma, λ=params.lambda, mx=params.mx, mz=params.mz)
		end
		grad_step2 = back((0f0, 1f0, 0f0, 0f0))
		Flux.Optimise.update!(opt_step1, ps_step2, grad_step2)

		push!(history, :loss_step_1, iter, loss_1[1])
		push!(history, :encoder_loss, iter, loss_2[2])
		push!(history, :generator_loss, iter, loss_1[3])
		push!(history, :transformer_loss, iter, loss_1[4])

		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:loss_step_1, loss_1[1]),
			(:loss_step_2, loss_2[2])
			])
		# TODO: check and prepare early stopping
		if mod(iter, params.check_every) == 0
			total_val_loss = 0
			Flux.testmode!(advae)
			for X_val in val_loader
				X_val = X_val |> gpu
				(𝓛, 𝓛ₑ, _, _) = loss(advae|>gpu , X_val, γ=params.gamma, λ=params.lambda, mx=params.mx, mz=params.mz)
				total_val_loss += 𝓛 + 𝓛ₑ
			end
			Flux.testmode!(advae, false)
			push!(history, :validation_loss, iter, total_val_loss/val_batches)

			if total_val_loss < best_val_loss
				best_val_loss = total_val_loss
				patience = params.patience
				best_adVAE = deepcopy(advae)
			else
				patience -= 1
				if patience == 0
					@info "Stopped training after $(iter) iterations"
					global iters = iter - params.check_every*params.patience
					break
				end
			end
		end
		if iter == params.iters
			global iters = params.iters
		end
	end
	return history, best_adVAE, sum(map(p->length(p), Flux.params(advae))), iters
end


text:
using SymbolicUtils: Rewriters

const KEEP = typemin(Int)

include("compat/bareiss.jl")

function alias_elimination(sys)
    sys = initialize_system_structure(sys; quick_cancel=true)
    s = structure(sys)
    is_linear_equations, eadj, cadj = find_linear_equations(sys)

    v_eliminated, v_types, n_null_vars, degenerate_equations, linear_equations = alias_eliminate_graph(
        s, is_linear_equations, eadj, cadj
    )

    s = structure(sys)
    @unpack fullvars, graph = s

    n_reduced_states = length(v_eliminated) - n_null_vars
    subs = OrderedDict()
    if n_reduced_states > 0
        for (i, v) in enumerate(@view v_eliminated[n_null_vars+1:end])
            subs[fullvars[v]] = iszeroterm(v_types, v) ? 0.0 :
                                isalias(v_types, v) ? fullvars[alias(v_types, v)] :
                                -fullvars[negalias(v_types, v)]
        end
    end

    dels = Set{Int}()
    eqs = copy(equations(sys))
    for (ei, e) in enumerate(linear_equations)
        vs = 𝑠neighbors(graph, e)
        if isempty(vs)
            push!(dels, e)
        else
            rhs = 0
            for vj in eachindex(vs)
                var = fullvars[vs[vj]]
                rhs += cadj[ei][vj] * var
            end
            eqs[e] = 0 ~ rhs
        end
    end
    dels = sort(collect(dels))
    deleteat!(eqs, dels)

    dict = Dict(subs)
    for (ieq, eq) in enumerate(eqs)
        eqs[ieq] = eq.lhs ~ fixpoint_sub(eq.rhs, dict)
    end

    newstates = []
    sts = states(sys)
    for j in eachindex(fullvars)
        if isirreducible(v_types, j)
            isdervar(s, j) || push!(newstates, fullvars[j])
        end
    end

    @set! sys.eqs = eqs
    @set! sys.states = newstates
    @set! sys.observed = [observed(sys); [lhs ~ rhs for (lhs, rhs) in pairs(subs)]]
    @set! sys.structure = nothing
    return sys
end

"""
    SparseMatrixCLIL{T, Ti}

The SparseMatrixCLIL represents a sparse matrix in two distinct ways:

1. As a sparse (in both row and column) n x m matrix
2. As a row-dense, column-sparse k x m matrix

The data structure keeps a permutation between the row order of the two representations.
Swapping the rows in one does not affect the other.

On construction, the second representation is equivalent to the first with fully-sparse
rows removed, though this may cease being true as row permutations are being applied
to the matrix.

The default structure of the `SparseMatrixCLIL` type is the second structure, while
the first is available via the thin `AsSubMatrix` wrapper.
"""
struct SparseMatrixCLIL{T, Ti<:Integer} <: AbstractSparseMatrix{T, Ti}
    nparentrows::Int
    ncols::Int
    nzrows::Vector{Ti}
    row_cols::Vector{Vector{Ti}}
    row_vals::Vector{Vector{T}}
end
Base.size(S::SparseMatrixCLIL) = (length(S.nzrows), S.ncols)
Base.copy(S::SparseMatrixCLIL{T, Ti}) where {T, Ti} =
    SparseMatrixCLIL(S.nparentrows, S.ncols, copy(S.nzrows), copy(S.row_cols), copy(S.row_vals))
function swaprows!(S::SparseMatrixCLIL, i, j)
    swap!(S.nzrows, i, j)
    swap!(S.row_cols, i, j)
    swap!(S.row_vals, i, j)
end

function bareiss_update_virtual_colswap_mtk!(zero!, M::SparseMatrixCLIL, k, swapto, pivot, last_pivot; pivot_equal_optimization=true)
    # for ei in nzrows(>= k)
    eadj = M.row_cols
    old_cadj = M.row_vals
    vpivot = swapto[2]

    ## N.B.: Micro-optimization
    #
    # For rows that do not have an entry in the eliminated column, all this
    # update does is multiply the row in question by `pivot/last_pivot` (the
    # result of which is guaranteed to be integer by general properties of the
    # bareiss algorithm, even if `pivot/last_pivot` is not).
    #
    # Thus, when `pivot == last pivot`, we can skip the update for any rows that
    # do not have an entry in the eliminated column (because we'd simply be
    # multiplying by 1).
    #
    # As an additional MTK-specific enhancement, we further allow the case
    # when the absolute values are equal, i.e. effectively multiplying the row
    # by `-1`. To ensure this is legal, we need to show two things.
    # 1. The multiplication does not change the answer and
    # 2. The multiplication does not affect the fraction-freeness of the Bareiss
    #    algorithm.
    #
    # For point 1, remember that we're working on a system of linear equations,
    # so it is always legal for us to multiply any row by a sclar without changing
    # the underlying system of equations.
    #
    # For point 2, note that the factorization we're now computing is the same
    # as if we had multiplied the corresponding row (accounting for row swaps)
    # in the original matrix by `last_pivot/pivot`, ensuring that the matrix
    # itself has integral entries when `last_pivot/pivot` is integral (here we
    # have -1, which counts). We could use the more general integrality
    # condition, but that would in turn disturb the growth bounds on the
    # factorization matrix entries that the bareiss algorithm guarantees. To be
    # conservative, we leave it at this, as this captures the most important
    # case for MTK (where most pivots are `1` or `-1`).
    pivot_equal = pivot_equal_optimization && abs(pivot) == abs(last_pivot)

    for ei in k+1:size(M, 1)
        # elimate `v`
        coeff = 0
        ivars = eadj[ei]
        vj = findfirst(isequal(vpivot), ivars)
        if vj !== nothing
            coeff = old_cadj[ei][vj]
            deleteat!(old_cadj[ei], vj)
            deleteat!(eadj[ei], vj)
        elseif pivot_equal
            continue
        end

        # the pivot row
        kvars = eadj[k]
        kcoeffs = old_cadj[k]
        # the elimination target
        ivars = eadj[ei]
        icoeffs = old_cadj[ei]

        tmp_incidence = similar(eadj[ei], 0)
        tmp_coeffs = similar(old_cadj[ei], 0)
        vars = union(ivars, kvars)

        for v in vars
            v == vpivot && continue
            ck = getcoeff(kvars, kcoeffs, v)
            ci = getcoeff(ivars, icoeffs, v)
            ci = (pivot*ci - coeff*ck) ÷ last_pivot
            if ci !== 0
                push!(tmp_incidence, v)
                push!(tmp_coeffs, ci)
            end
        end

        eadj[ei] = tmp_incidence
        old_cadj[ei] = tmp_coeffs
    end

    # Swap pivots to the front of the coefficient list
    # TODO: This prevents the coefficient list from being sorted, making
    # the rest of the algorithm much more expensive
    pivot_idx = findfirst(==(vpivot), eadj[k])
    deleteat!(eadj[k], pivot_idx)
    deleteat!(old_cadj[k], pivot_idx)
    pushfirst!(eadj[k], vpivot)
    pushfirst!(old_cadj[k], pivot)
end

function bareiss_update_virtual_colswap_mtk!(zero!, M::AbstractMatrix, k, swapto, pivot, last_pivot; pivot_equal_optimization=true)
    if pivot_equal_optimization
        error("MTK pivot micro-optimization not implemented for `$(typeof(M))`.
            Turn off the optimization for debugging or use a different matrix type.")
    end
    bareiss_update_virtual_colswap!(zero!, M, k, swapto, pivot, last_pivot)
end

struct AsSubMatrix{T, Ti<:Integer} <: AbstractSparseMatrix{T, Ti}
    M::SparseMatrixCLIL{T, Ti}
end
Base.size(S::AsSubMatrix) = (S.M.nparentrows, S.M.ncols)

function Base.getindex(S::SparseMatrixCLIL{T}, i1, i2) where {T}
    checkbounds(S, i1, i2)

    nncol = findfirst(==(i2), S.row_cols[i1])
    isnothing(nncol) && return zero(T)

    return S.row_vals[i1][nncol]
end

function Base.getindex(S::AsSubMatrix{T}, i1, i2) where {T}
    checkbounds(S, i1, i2)
    S = S.M

    nnrow = findfirst(==(i1), S.nzrows)
    isnothing(nnrow) && return zero(T)

    nncol = findfirst(==(i2), S.row_cols[nnrow])
    isnothing(nncol) && return zero(T)

    return S.row_vals[nnrow][nncol]
end

"""
$(SIGNATURES)

Find the first linear variable such that `𝑠neighbors(adj, i)[j]` is true given
the `constraint`.
"""
@inline function find_first_linear_variable(
        M::SparseMatrixCLIL,
        range,
        mask,
        constraint,
    )
    eadj = M.row_cols
    for i in range
        vertices = eadj[i]
        if constraint(length(vertices))
            for (j, v) in enumerate(vertices)
                (mask === nothing || mask[v]) && return (CartesianIndex(i, v), M.row_vals[i][j])
            end
        end
    end
    return nothing
end

@inline function find_first_linear_variable(
        M::AbstractMatrix,
        range,
        mask,
        constraint,
    )
    for i in range
        row = @view M[i, :]
        if constraint(count(!iszero, row))
            for (v, val) in enumerate(row)
                iszero(val) && continue
                if mask === nothing || mask[v]
                    return CartesianIndex(i, v), val
                end
            end
        end
    end
    return nothing
end

function find_masked_pivot(variables, M, k)
    r = find_first_linear_variable(M, k:size(M,1), variables, isequal(1))
    r !== nothing && return r
    r = find_first_linear_variable(M, k:size(M,1), variables, isequal(2))
    r !== nothing && return r
    r = find_first_linear_variable(M, k:size(M,1), variables, _->true)
    return r
end

function alias_eliminate_graph(s::SystemStructure, is_linear_equations, eadj, cadj)
    @unpack graph, varassoc = s
    invvarassoc = inverse_mapping(varassoc)

    old_cadj = map(copy, cadj)

    is_not_potential_state = iszero.(varassoc)
    is_linear_variables = copy(is_not_potential_state)
    for i in 𝑠vertices(graph); is_linear_equations[i] && continue
        for j in 𝑠neighbors(graph, i)
            is_linear_variables[j] = false
        end
    end
    solvable_variables = findall(is_linear_variables)

    linear_equations = findall(is_linear_equations)

    mm = SparseMatrixCLIL(nsrcs(graph),
                          ndsts(graph),
                          linear_equations, eadj, old_cadj)

    function do_bareiss!(M, cadj=nothing)
        rank1 = rank2 = nothing
        function find_pivot(M, k)
            if rank1 === nothing
                r = find_masked_pivot(is_linear_variables, M, k)
                r !== nothing && return r
                rank1 = k - 1
            end
            if rank2 === nothing
                r = find_masked_pivot(is_not_potential_state, M, k)
                r !== nothing && return r
                rank2 = k - 1
            end
            return find_masked_pivot(nothing, M, k)
        end
        function myswaprows!(M, i, j)
            cadj !== nothing && swap!(cadj, i, j)
            swaprows!(M, i, j)
        end
        bareiss_ops = ((M,i,j)->nothing, myswaprows!, bareiss_update_virtual_colswap_mtk!, bareiss_zero!)
        rank3 = bareiss!(M, bareiss_ops; find_pivot)
        rank1 = something(rank1, rank3)
        rank2 = something(rank2, rank3)
        (rank1, rank2, rank3)
    end

    # mm2 = Array(copy(mm))
    # @show do_bareiss!(mm2)
    # display(mm2)

    (rank1, rank2, rank3) = do_bareiss!(mm, cadj)

    v_solved = [eadj[i][1] for i in 1:rank1]
    v_eliminated = setdiff(solvable_variables, v_solved)
    n_null_vars = length(v_eliminated)

    v_types = fill(KEEP, ndsts(graph))
    for v in v_eliminated
        v_types[v] = 0
    end

    # kind of like the backward substitution
    for ei in reverse(1:rank2)
        locally_structure_simplify!(
                                    (eadj[ei], old_cadj[ei]),
                                    invvarassoc, v_eliminated, v_types
                                   )
    end

    reduced = false
    for ei in 1:rank2
        if length(cadj[ei]) >= length(old_cadj[ei])
            cadj[ei] = old_cadj[ei]
        else
            # MEMORY ALIAS of a vector
            eadj[ei] = 𝑠neighbors(graph, linear_equations[ei])
            reduced |= locally_structure_simplify!(
                                                   (eadj[ei], cadj[ei]),
                                                   invvarassoc, v_eliminated, v_types
                                                  )
        end
    end

    while reduced
        for ei in 1:rank2
            if !isempty(eadj[ei])
                reduced |= locally_structure_simplify!(
                                                       (eadj[ei], cadj[ei]),
                                                       invvarassoc, v_eliminated, v_types
                                                      )
                reduced && break # go back to the begining of equations
            end
        end
    end

    for ei in rank2+1:length(linear_equations)
        cadj[ei] = old_cadj[ei]
    end

    for (ei, e) in enumerate(linear_equations)
        graph.fadjlist[e] = eadj[ei]
    end

    degenerate_equations = rank3 < length(linear_equations) ? linear_equations[rank3+1:end] : Int[]
    return v_eliminated, v_types, n_null_vars, degenerate_equations, linear_equations
end

iszeroterm(v_types, v) = v_types[v] == 0
isirreducible(v_types, v) = v_types[v] == KEEP
isalias(v_types, v) = v_types[v] > 0 && !isirreducible(v_types, v)
alias(v_types, v) = v_types[v]
negalias(v_types, v) = -v_types[v]

function locally_structure_simplify!(
        (vars, coeffs),
        invvarassoc, v_eliminated, v_types
       )
    while length(vars) > 1 && any(!isequal(KEEP), (v_types[v] for v in @view vars[2:end]))
        for vj in 2:length(vars)
            v = vars[vj]
            if isirreducible(v_types, v)
                continue
            elseif iszeroterm(v_types, v)
                deleteat!(vars, vj)
                deleteat!(coeffs, vj)
                break
            else
                coeff = coeffs[vj]
                if isalias(v_types, v)
                    v = alias(v_types, v)
                else
                    v = negalias(v_types, v)
                    coeff = -coeff
                end

                has_v = false
                for vi in 2:length(vars)
                    (vi !== vj && vars[vi] == v) || continue
                    has_v = true
                    c = (coeffs[vi] += coeff)
                    if c == 0
                        if vi < vj
                            deleteat!(vars, [vi, vj])
                            deleteat!(coeffs, [vi, vj])
                        else
                            deleteat!(vars, [vj, vi])
                            deleteat!(coeffs, [vj, vi])
                        end
                    end
                    break
                end # for vi

                if has_v
                    break
                else
                    vars[vj] = v
                    coeffs[vj] = coeff
                end # if
            end # else
        end # for
    end # while

    v = first(vars)

    # Do not attempt to eliminate derivatives
    invvarassoc[v] != 0 && return false

    if length(vars) == 1
        push!(v_eliminated, v)
        v_types[v] = 0
        empty!(vars); empty!(coeffs)
        return true
    elseif length(vars) == 2 && abs(coeffs[1]) == abs(coeffs[2])
        if (coeffs[1] > 0 && coeffs[2] < 0) || (coeffs[1] < 0 && coeffs[2] > 0)
            # positive alias
            push!(v_eliminated, v)
            v_types[v] = vars[2]
        else
            # negative alias
            push!(v_eliminated, v)
            v_types[v] = -vars[2]
        end
        empty!(vars); empty!(coeffs)
        return true
    end
    return false
end

swap!(v, i, j) = v[i], v[j] = v[j], v[i]

function getcoeff(vars, coeffs, var)
    for (vj, v) in enumerate(vars)
        v == var && return coeffs[vj]
    end
    return 0
end

function inverse_mapping(assoc)
    invassoc = zeros(Int, length(assoc))
    for (i, v) in enumerate(assoc)
        v <= 0 && continue
        invassoc[v] = i
    end
    return invassoc
end

"""
$(SIGNATURES)

Use Kahn's algorithm to topologically sort observed equations.

Example:
```julia
julia> @variables t x(t) y(t) z(t) k(t)
(t, x(t), y(t), z(t), k(t))

julia> eqs = [
           x ~ y + z
           z ~ 2
           y ~ 2z + k
       ];

julia> ModelingToolkit.topsort_equations(eqs, [x, y, z, k])
3-element Vector{Equation}:
 Equation(z(t), 2)
 Equation(y(t), k(t) + 2z(t))
 Equation(x(t), y(t) + z(t))
```
"""
function topsort_equations(eqs, states; check=true)
    graph, assigns = observed2graph(eqs, states)
    neqs = length(eqs)
    degrees = zeros(Int, neqs)

    for 𝑠eq in 1:length(eqs); var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            # 𝑠eq => 𝑑eq
            degrees[𝑑eq] += 1
        end
    end

    q = Queue{Int}(neqs)
    for (i, d) in enumerate(degrees)
        d == 0 && enqueue!(q, i)
    end

    idx = 0
    ordered_eqs = similar(eqs, 0); sizehint!(ordered_eqs, neqs)
    while !isempty(q)
        𝑠eq = dequeue!(q)
        idx+=1
        push!(ordered_eqs, eqs[𝑠eq])
        var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            degree = degrees[𝑑eq] = degrees[𝑑eq] - 1
            degree == 0 && enqueue!(q, 𝑑eq)
        end
    end

    (check && idx != neqs) && throw(ArgumentError("The equations have at least one cycle."))

    return ordered_eqs
end

function observed2graph(eqs, states)
    graph = BipartiteGraph(length(eqs), length(states))
    v2j = Dict(states .=> 1:length(states))

    # `assigns: eq -> var`, `eq` defines `var`
    assigns = similar(eqs, Int)

    for (i, eq) in enumerate(eqs)
        lhs_j = get(v2j, eq.lhs, nothing)
        lhs_j === nothing && throw(ArgumentError("The lhs $(eq.lhs) of $eq, doesn't appear in states."))
        assigns[i] = lhs_j
        vs = vars(eq.rhs)
        for v in vs
            j = get(v2j, v, nothing)
            j !== nothing && add_edge!(graph, i, j)
        end
    end

    return graph, assigns
end

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(eqs, dict)
    sub = Base.Fix2(fixpoint_sub, dict)
    map(eq->eq.lhs ~ sub(eq.rhs), eqs)
end

ext: jl


text:
# WIP

@fastmath function ltc_syn_state!(du,u,p,t, n_in, n_neurons, w_sens, w_syns)
  v, nt_sens, nt_syns = ltc_syn_state_unpack_u(n_in, n_neurons, u)
  dv, dnt_sens, dnt_syns = ltc_syn_state_unpack_u(n_in, n_neurons, du)

  stim, C, G_L, E_L,
  A_ntr_sens, vh_ntr_sens, s_ntr_sens, τ_accum_ntc_sens, τ_elim_ntc_sens, r_ra_sens, nth_ra_sens, G_sens, E_sens,
  A_ntr_syns, vh_ntr_syns, s_ntr_syns, τ_accum_ntc_syns, τ_elim_ntc_syns, r_ra_syns, nth_ra_syns, G_syns, E_syns = ltc_syn_state_unpack_p(n_in,n_neurons,p)

  I_senss = zeros(eltype(stim), n_in, n_neurons)
  I_synss = zeros(eltype(stim), n_neurons, n_neurons)

  @inbounds @simd for d in 1:n_neurons
    for s in 1:n_in
      I_senss[s,d] = w_sens[s,d] * G_sens[s,d] * sensa[s,d] * (v[d] - E_sens[s,d])
      α_sensa[s,d] = α_sensa_A[s,d] * exp((stim[s]-Nt_vh[s,d])/Nt_k[s,d])
      β_sensa[s,d] = τ_accum[s,d] * ((stim[s] -τ_elim[s,d]) / (exp((stim[s] -τ_elim[s,d])/β_sensa_k[s,d]) - 1f0))
    end
    for s in 1:n_neurons
      I_synss[s,d] = w_syns[s,d] * G_syns[s,d] * synsa[s,d] * (v[d] - E_syns[s,d])
      α_synsa[s,d] = α_synsa_A[s,d] * exp((v[s]-α_synsa_vh[s,d])/α_synsa_k[s,d])
      β_synsa[s,d] = β_synsa_A[s,d] * ((v[s] -β_synsa_vh[s,d]) / (exp((v[s] -β_synsa_vh[s,d])/β_synsa_k[s,d]) - 1f0))
    end
  end
  I_sens = reshape(sum(I_senss, dims=1), :)
  I_syns = reshape(sum(I_synss, dims=1), :)

  I_tot = (G_L .* (v .- E_L)) .+ I_sens .+ I_syns

  ntr_sens = A_ntr_sens / (1 + exp(-s_ntr_sens * (v - vh_ntr_sens)))

  @.. dv = -I_tot / C
  @.. dnt_sens = ntr_sens / τ_accum_sens - nt_sens / τ_elim_sens
  @.. dnt_syns = ntr_syns / τ_accum_syns - nt_syns / τ_elim_syns

  nothing
end

function ltc_syn_state_init_u0p(n_in, n_neurons; T=Float32)
  stim   = fill(T(0), n_in)
  #w_sens = zeros(T, n_in, n_neurons)
  #w_syns = zeros(T, n_neurons, n_neurons)
  C      = rand_uniform(T,  1.0,     1.002,   n_neurons)
  G_L    = rand_uniform(T,  0.0001,  0.1,     n_neurons)
  E_L    = rand_uniform(T, -0.3,     0.3,     n_neurons) #.± 0.1f0
  G_sens = rand_uniform(T,  0.0001,  0.1,     n_in, n_neurons)
  α_sensa_A  = rand_uniform(T, 0.1, 1,     n_in, n_neurons)
  Nt_vh = rand_uniform(T, 0.2, 0.4, n_in, n_neurons)
  Nt_k  = rand_uniform(T, 0.1, 0.9,     n_in, n_neurons)
  τ_accum  = rand_uniform(T, 0.1, 1,     n_in, n_neurons)
  τ_elim = rand_uniform(T, 0.2, 0.4, n_in, n_neurons)
  β_sensa_k  = rand_uniform(T, 0.1, 0.9,     n_in, n_neurons)
  E_sens = rand_uniform(T, -0.5,     0.5,     n_in, n_neurons)
  G_syns = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)
  α_synsa_A  = rand_uniform(T, 0.1, 1,     n_neurons, n_neurons)
  α_synsa_vh = rand_uniform(T, 0.2, 0.4, n_neurons, n_neurons)
  α_synsa_k  = rand_uniform(T, 0.1, 0.9,     n_neurons, n_neurons)
  β_synsa_A  = rand_uniform(T, 0.1, 1,     n_neurons, n_neurons)
  β_synsa_vh = rand_uniform(T, 0.2, 0.4, n_neurons, n_neurons)
  β_synsa_k  = rand_uniform(T, 0.1, 0.9,     n_neurons, n_neurons)
  E_syns = rand_uniform(T, -0.5,     0.5,     n_neurons, n_neurons)

  v = rand_uniform(T, 0.01, 0.1, n_neurons) #.± 0.1f0
  sensa = rand_uniform(T, 0.001, 0.01, n_in, n_neurons)
  synsa = rand_uniform(T, 0.001, 0.01, n_neurons, n_neurons)

  u0 = vcat(v, vec(sensa), vec(synsa))
  p = vcat(stim,
    C, G_L, E_L,
    vec(G_sens), vec(α_sensa_A),vec(Nt_vh),vec(Nt_k),vec(τ_accum),vec(τ_elim),vec(β_sensa_k), vec(E_sens),
    vec(G_syns), vec(α_synsa_A),vec(α_synsa_vh),vec(α_synsa_k),vec(β_synsa_A),vec(β_synsa_vh),vec(β_synsa_k), vec(E_syns),
  )
  return u0, p
end

function ltc_syn_state_bounds(n_in, n_neurons; T=Float32)
  lb = vcat(
    [1.0 for _ in 1:n_neurons],               # C
    [0 for _ in 1:n_neurons],                 # G_L
    [-1 for _ in 1:n_neurons],                # E_L
    [0 for _ in 1:n_in*n_neurons],            # G_sens
    [0.01 for _ in 1:n_in*n_neurons],            # α_sensa_A
    [0 for _ in 1:n_in*n_neurons],            # Nt_vh
    [0.01 for _ in 1:n_in*n_neurons],            # Nt_k
    [0.01 for _ in 1:n_in*n_neurons],            # τ_accum
    [0 for _ in 1:n_in*n_neurons],            # τ_elim
    [0.01 for _ in 1:n_in*n_neurons],            # β_sensa_k
    [-1 for _ in 1:n_in*n_neurons],           # E_sens
    [0 for _ in 1:n_neurons*n_neurons],       # G_syns
    [0.01 for _ in 1:n_neurons*n_neurons],            # α_synsa_A
    [0 for _ in 1:n_neurons*n_neurons],            # α_synsa_vh
    [0.01 for _ in 1:n_neurons*n_neurons],            # α_synsa_k
    [0.01 for _ in 1:n_neurons*n_neurons],            # β_synsa_A
    [0 for _ in 1:n_neurons*n_neurons],            # β_synsa_vh
    [0.01 for _ in 1:n_neurons*n_neurons],            # β_synsa_k
    [-1 for _ in 1:n_neurons*n_neurons],      # E_syns

    [-0.4 for _ in 1:n_neurons],             # v
    [0 for _ in 1:n_in*n_neurons],            # sensa
    [0 for _ in 1:n_neurons*n_neurons])             # synsa

  ub = vcat(
    [2.0 for _ in 1:n_neurons],
    [0.4 for _ in 1:n_neurons],
    [1 for _ in 1:n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [1.1 for _ in 1:n_in*n_neurons],            # α_sensa_A
    [1 for _ in 1:n_in*n_neurons],            # Nt_vh
    [1 for _ in 1:n_in*n_neurons],            # Nt_k
    [1 for _ in 1:n_in*n_neurons],            # τ_accum
    [1 for _ in 1:n_in*n_neurons],            # τ_elim
    [1 for _ in 1:n_in*n_neurons],            # β_sensa_k
    [1 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [1.1 for _ in 1:n_neurons*n_neurons],            # α_synsa_A
    [1 for _ in 1:n_neurons*n_neurons],            # α_synsa_vh
    [1 for _ in 1:n_neurons*n_neurons],            # α_synsa_k
    [1 for _ in 1:n_neurons*n_neurons],            # β_synsa_A
    [1 for _ in 1:n_neurons*n_neurons],            # β_synsa_vh
    [1 for _ in 1:n_neurons*n_neurons],            # β_synsa_k
    [1 for _ in 1:n_neurons*n_neurons],

    [0.2 for _ in 1:n_neurons],
    [0.1 for _ in 1:n_in*n_neurons],            # sensa
    [0.1 for _ in 1:n_neurons*n_neurons],

    ) # v)

  T.(lb), T.(ub)
end

function ltc_syn_state_unpack_u(n_in,n_neurons,u)
  vl = n_neurons
  sensal = n_in*n_neurons
  synsal = n_neurons*n_neurons

  s = 1
  v = @view u[s:s+vl-1]
  s += vl
  sensa = @view u[reshape(s:s+sensal-1, n_in, n_neurons)]
  s += sensal
  synsa = @view u[reshape(s:s+synsal-1, n_neurons, n_neurons)]
  return v, sensa, synsa
end

function ltc_syn_state_unpack_p(n_in,n_neurons,p)
  stiml = n_in
  Cl = n_neurons
  G_Ll = n_neurons
  E_Ll = n_neurons
  sensl = n_in*n_neurons
  synsl = n_neurons*n_neurons
  G_sensl = sensl
  α_sensa_Al = sensl
  Nt_vhl = sensl
  Nt_kl = sensl
  τ_accuml = sensl
  τ_eliml = sensl
  β_sensa_kl = sensl
  E_sensl = sensl
  G_synsl = synsl
  α_synsa_Al = synsl
  α_synsa_vhl = synsl
  α_synsa_kl = synsl
  β_synsa_Al = synsl
  β_synsa_vhl = synsl
  β_synsa_kl = synsl
  E_synsl = synsl

  s = 1
  stim = @view p[s:s+stiml-1]
  s += stiml
  C = @view p[s:s+Cl-1]
  s += Cl
  G_L = @view p[s:s+G_Ll-1]
  s += G_Ll
  E_L = @view p[s:s+E_Ll-1]
  s += E_Ll
  G_sens = @view p[reshape(s:s+G_sensl-1, n_in, n_neurons)]
  s += G_sensl
  α_sensa_A = @view p[reshape(s:s+α_sensa_Al-1, n_in, n_neurons)]
  s += α_sensa_Al
  Nt_vh = @view p[reshape(s:s+Nt_vhl-1, n_in, n_neurons)]
  s += Nt_vhl
  Nt_k = @view p[reshape(s:s+Nt_kl-1, n_in, n_neurons)]
  s += Nt_kl
  τ_accum = @view p[reshape(s:s+τ_accuml-1, n_in, n_neurons)]
  s += τ_accuml
  τ_elim = @view p[reshape(s:s+τ_eliml-1, n_in, n_neurons)]
  s += τ_eliml
  β_sensa_k = @view p[reshape(s:s+β_sensa_kl-1, n_in, n_neurons)]
  s += β_sensa_kl
  E_sens = @view p[reshape(s:s+E_sensl-1, n_in, n_neurons)]
  s += E_sensl

  G_syns = @view p[reshape(s:s+G_synsl-1, n_neurons, n_neurons)]
  s += G_synsl
  α_synsa_A = @view p[reshape(s:s+α_synsa_Al-1, n_neurons, n_neurons)]
  s += α_sensa_Al
  α_synsa_vh = @view p[reshape(s:s+α_synsa_vhl-1, n_neurons, n_neurons)]
  s += α_synsa_vhl
  α_synsa_k = @view p[reshape(s:s+α_synsa_kl-1, n_neurons, n_neurons)]
  s += α_synsa_kl
  β_synsa_A = @view p[reshape(s:s+β_synsa_Al-1, n_neurons, n_neurons)]
  s += β_synsa_Al
  β_synsa_vh = @view p[reshape(s:s+β_synsa_vhl-1, n_neurons, n_neurons)]
  s += β_synsa_vhl
  β_synsa_k = @view p[reshape(s:s+β_synsa_kl-1, n_neurons, n_neurons)]
  s += β_synsa_kl
  E_syns = @view p[reshape(s:s+E_synsl-1, n_neurons, n_neurons)]
  s += E_synsl
  return stim, C, G_L, E_L, G_sens, α_sensa_A, Nt_vh, Nt_k, τ_accum, τ_elim, β_sensa_k, E_sens, G_syns, α_synsa_A, α_synsa_vh, α_synsa_k, β_synsa_A, β_synsa_vh, β_synsa_k, E_syns
end

function LTCSynState(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),
  #wiring=nothing,
  mtkize=true, gen_jac=false, kwargs...)

  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]

  s_in = wiring.s_in
  n_neurons = wiring.n_total

  lb, ub = ltc_syn_state_bounds(s_in, n_neurons; T)
  u0, p = ltc_syn_state_init_u0p(s_in, n_neurons; T)

  dudt!(du,u,p,t) = ltc_syn_state!(du,u,p,t, s_in, n_neurons, w_sens, w_syns)

  rnncell = BCTRNNCell(wiring, solver, sensealg, dudt!, u0, tspan, p, lb, ub; mtkize, gen_jac, kwargs...)
  MyRecur(rnncell)
end


text:
@block JuliusKrumbiegel ["legend"] begin
    using MakieLayout
    using AbstractPlotting: px
    @cell "Basic legend" [legend] begin

        using AbstractPlotting: px

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        xs = 0:0.5:10
        ys = sin.(xs)
        lin = lines!(ax, xs, ys, color = :blue)
        sca = scatter!(ax, xs, ys, color = :red, markersize = 15px)

        leg = LLegend(scene, [lin, sca, [lin, sca]], ["a line", "some dots", "both together"])
        layout[1, 2] = leg

        scene
    end

    @cell "Inline legend" [legend] begin
        scene, layout = layoutscene()

        ax = layout[1, 1] = LAxis(scene)

        ls = [lines!(ax, 1:10, (1:10) .* i, color = rand(RGBf0)) for i in 1:5]

        leg = layout[1, 1] = LLegend(scene, ls, ["line $i" for i in 1:5];
            width = Auto(false), height = Auto(false), halign = :left, valign = :top,
            margin = (10, 10, 10, 10))

        scene
    end

    @cell "Inline legend postioning" [legend] begin

        haligns = [:left, :right, :center]
        valigns = [:top, :bottom, :center]

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        xs = 0:0.1:10
        lins = [lines!(ax, xs, sin.(xs .* i), color = color)
            for (i, color) in zip(1:3, [:red, :blue, :green])]

        legends = [LLegend(
                scene, lins, ["Line $i" for i in 1:3],
                "$ha & $va",
                width = Auto(false),
                margin = (10, 10, 10, 10),
                halign = ha, valign = va, orientation = :horizontal
            ) for (j, ha, va) in zip(1:3, haligns, valigns)]

        for leg in legends
            layout[1, 1] = leg
        end

        scene
    end

    @cell "Multi-column legend" [legend] begin
        using AbstractPlotting: px

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        xs = 0:0.1:10
        lins = [lines!(ax, xs, sin.(xs .+ 3v), color = RGBf0(v, 0, 1-v)) for v in 0:0.1:1]

        leg = LLegend(scene, lins, string.(1:length(lins)), ncols = 3)

        layout[1, 2] = leg

        scene
    end

    @cell "Manual legend" [legend] begin
        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        elem_1 = [LineElement(color = :red, linestyle = nothing),
                  MarkerElement(color = :blue, marker = 'x', strokecolor = :black)]

        elem_2 = [PolyElement(color = :red, strokecolor = :blue),
                  LineElement(color = :black, linestyle = :dash)]

        elem_3 = LineElement(color = :green, linestyle = nothing,
                linepoints = Point2f0[(0, 0), (0, 1), (1, 0), (1, 1)])

        elem_4 = MarkerElement(color = :blue, marker = 'π',
                strokecolor = :transparent,
                markerpoints = Point2f0[(0.2, 0.2), (0.5, 0.8), (0.8, 0.2)])

        elem_5 = PolyElement(color = :green, strokecolor = :black,
                polypoints = Point2f0[(0, 0), (1, 0), (0, 1)])

        leg = layout[1, 2] = LLegend(scene,
            [elem_1, elem_2, elem_3, elem_4, elem_5],
            ["Line & Marker", "Poly & Line", "Line", "Marker", "Poly"])

        scene
    end

    @cell "Horizontal legend" [legend] begin
        using AbstractPlotting: px

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        xs = 0:0.5:10
        ys = sin.(xs)
        lin = lines!(ax, xs, ys, color = :blue)
        sca = scatter!(ax, xs, ys, color = :red, markersize = 15px)

        leg = LLegend(scene, [lin, sca, lin], ["a line", "some dots", "line again"])
        layout[1, 2] = leg

        leg_horizontal = LLegend(scene, [lin, sca, lin], ["a line", "some dots", "line again"],
        orientation = :horizontal, width = Auto(false), height = Auto(true))
        layout[2, 1] = leg_horizontal
        scene
    end

    @cell "Multi-bank legend" [legend] begin

        # You can control the number of banks with the `nbanks` attribute.
        # Banks are columns when in vertical mode,
        # and rows when in horizontal mode.

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        xs = 0:0.1:10
        lins = [lines!(ax, xs, sin.(xs .+ 3v), color = RGBf0(v, 0, 1-v)) for v in 0:0.1:1]

        leg = LLegend(scene, lins, string.(1:length(lins)), nbanks = 3)
        layout[1, 2] = leg

        scene

    end

    @cell "Multi-group legends" [legend] begin

        scene, layout = layoutscene(resolution = (1400, 900))

        ax = layout[1, 1] = LAxis(scene)

        markersizes = [5, 10, 15, 20]
        colors = [:red, :green, :blue, :orange]

        for ms in markersizes, color in colors
            scatter!(ax, randn(5, 2), markersize = ms * px, color = color)
        end

        group_size = [MarkerElement(marker = :circle, color = :black, strokecolor = :transparent,
            markersize = ms * px) for ms in markersizes]

        group_color = [PolyElement(color = color, strokecolor = :transparent)
            for color in colors]

        legends = [LLegend(scene,
            [group_size, group_color],
            [string.(markersizes), string.(colors)],
            ["Size", "Color"]) for _ in 1:6]

        layout[1, 2:4] = legends[1:3]
        layout[2:4, 1] = legends[4:6]

        for l in legends[4:6]
            l.orientation = :horizontal
            l.height = Auto(true)
            l.width = Auto(false)
        end

        legends[2].titleposition = :left
        legends[5].titleposition = :left

        legends[3].nbanks = 2
        legends[5].nbanks = 2
        legends[6].nbanks = 2

        scene
    end
end


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

name = "PDFHighlights"
version = v"0.1.1"

# Collection of sources required to complete build
sources = [
	GitSource("https://github.com/paveloom-j/PDFHighlights.jl.git", "327d2f73aafbda615101dc3e7b6a296de491dee8"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/PDFHighlights.jl/deps/

CC="gcc"
LIB_NAME="PDFHighlightsWrapper"
OBJECTS=(get_author_title get_lines_comments_pages)

OBJECTS_O=()
for index in ${!OBJECTS[@]}; do
    OBJECTS_O[${index}]=${OBJECTS[${index}]}.o
done

INCLUDE_FLAGS="`pkg-config --cflags poppler-glib`"
LINK_FLAGS="-lgio-2.0 `pkg-config --libs poppler-glib`"

if [[ "${target}" == *-darwin* ]]; then
    OBJ_FLAGS="-std=c99 -O3 -fPIC"
    LIB_FLAGS="-shared"
elif [[ "${target}" == *-mingw* ]]; then
    OBJ_FLAGS="-std=c99 -O3"
    LIB_FLAGS="-shared"
    LINK_FLAGS="${LINK_FLAGS} -Wl,--out-implib,${prefix}/lib/${LIB_NAME}.dll.a"
else
    OBJ_FLAGS="-std=c99 -O3 -fPIC"
    LIB_FLAGS="-shared -Wl,--no-undefined"
fi

for file in ${OBJECTS[@]}; do
    ${CC} ${OBJ_FLAGS} -c ${file}.c -o ${file}.o ${INCLUDE_FLAGS}
done

${CC} ${LIB_FLAGS} -o ${libdir}/${LIB_NAME}.${dlext} ${OBJECTS_O[@]} ${LINK_FLAGS}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Poppler_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct DataFile{T}
    filename::String
    savename::String
    units::Vector{String}
    legend_units::String
    idx::Int
end

function DataFile(filename::String, ext, delim, rename=true)
    filename = preprocess(filename, rename)
    savename = joinpath(dirname(filename), basename(filename) * ext)
    units = extract_units(filename, delim)

    idx = 1
    if occursin("CV", filename)
        legend_units = "mV/s"
        idx = 2
        T = Val{:CV}
    elseif occursin("C&D", filename)
        legend_units = "A"
        T = Val{Symbol("C&D")}
    elseif occursin("EIS", filename)
        legend_units = "mA/cm^2"
        T = Val{:EIS}
    else
        legend_units = ""
        T = Val{:unknown}
    end

    DataFile{T}(filename, savename, units, legend_units, idx)
end

function preprocess(filename, rename)
    !rename && return filename
    if occursin("C&D", filename)
        if !endswith(filename, "_C") && !endswith(filename, "_D")
            new_filename = filename * "_D"
            mv(filename, new_filename)
            @info "Renamed $filename to $new_filename"
            filename = new_filename
        end
    end

    return filename
end

function extract_units(filename, delim)
    firstline = readline(filename)
    parts = split(firstline, delim)
    units = String[]

    for p in parts
        m = match(r" \((?<unit>\w)\)", p)
        if !isnothing(m)
            push!(units, m[:unit])
        end
    end

    return units
end


text:
function _vcov(x::AbstractVector, y::AbstractVector, corrected::Bool, μᵪ::Number, μᵧ::Number)
    # Calculate covariance
    σᵪᵧ = zero(promote_type(typeof(μᵪ), typeof(μᵧ), Int))
    @turbo for i ∈ indices((x,y))
            δᵪ = x[i] - μᵪ
            δᵧ = y[i] - μᵧ
            σᵪᵧ += δᵪ * δᵧ
    end
    σᵪᵧ = σᵪᵧ / (length(x)-corrected)
    return σᵪᵧ
end

function _vtcov(x::AbstractVector, y::AbstractVector, corrected::Bool, μᵪ::Number, μᵧ::Number)
    # Calculate covariance
    σᵪᵧ = zero(promote_type(typeof(μᵪ), typeof(μᵧ), Int))
    @tturbo for i ∈ indices((x,y))
            δᵪ = x[i] - μᵪ
            δᵧ = y[i] - μᵧ
            σᵪᵧ += δᵪ * δᵧ
    end
    σᵪᵧ = σᵪᵧ / (length(x)-corrected)
    return σᵪᵧ
end

"""
```julia
vcov(x::AbstractVector, y::AbstractVector; corrected::Bool=true, multithreaded=:auto)
```
Compute the covariance between the vectors `x` and `y`.
As `Statistics.cov`, but vectorized and (optionally) multithreaded.

If `corrected` is `true` as is the default, _Bessel's correction_ will be applied,
such that the sum is scaled by `n-1` rather than `n`, where `n = length(x)`.
"""
function vcov(x::AbstractVector, y::AbstractVector; corrected::Bool=true, multithreaded=:auto)
    # Check lengths
    nᵪ = length(x)
    nᵧ = length(y)
    @assert nᵪ == nᵧ

    if (multithreaded==:auto && length(x) > 4095) || multithreaded==true
        μᵪ = _vtmean(x,:)
        μᵧ = _vtmean(y,:)
        σᵪᵧ = _vtcov(x, y, corrected, μᵪ, μᵧ)
    else
        μᵪ = _vmean(x,:)
        μᵧ = _vmean(y,:)
        σᵪᵧ = _vcov(x, y, corrected, μᵪ, μᵧ)
    end
    return σᵪᵧ
end

"""
```julia
vcov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true, multithreaded=:auto)
```
Compute the covariance matrix of the matrix `X`, along dimension `dims`.
As `Statistics.cov`, but vectorized and (optionally) multithreaded.

If `corrected` is `true` as is the default, _Bessel's correction_ will be applied,
such that the sum is scaled by `n-1` rather than `n`, where `n = length(x)`.
"""
function vcov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true, multithreaded=:auto)
    if (multithreaded===:auto && length(A) > 4095) || multithreaded===true
        _vtcov(X, dims, corrected)
    else
        _vcov(X, dims, corrected)
    end
end
export vcov

function _vcov(X, dims, corrected)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Σ = similar(X, Tₒ, (m, m))
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means for each column
        μ = ntuple(m) do d
            vmean(view(X,:,d))
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vcov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    elseif dims == 2
        # Precalculate means for each row
        μ = ntuple(m) do d
            vmean(view(X,d,:))
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vcov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    else
        throw("Dimension not in range")
    end
    return Σ
end

function _vtcov(X, dims, corrected)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Σ = similar(X, Tₒ, (m, m))
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means for each column
        μ = ntuple(m) do d
            _vtmean(view(X,:,d),:)
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vtcov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    elseif dims == 2
        # Precalculate means for each row
        μ = ntuple(m) do d
            _vtmean(view(X,d,:),:)
        end
        # Fill covariance matrix symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vtcov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Σ[i,j] = Σ[j,i] = σᵢⱼ
            end
        end
    else
        throw("Dimension not in range")
    end
    return Σ
end

"""
```julia
vcor(x::AbstractVector, y::AbstractVector, multithreaded=:auto)
```
Compute the (Pearson's product-moment) correlation between the vectors `x` and `y`.
As `Statistics.cor`, but vectorized and (optionally) multithreaded.

Equivalent to `cov(x,y) / (std(x) * std(y))`.
"""
function vcor(x::AbstractVector, y::AbstractVector; corrected::Bool=true, multithreaded=:auto)
    # Check lengths
    nᵪ = length(x)
    nᵧ = length(y)
    @assert nᵪ == nᵧ

    if (multithreaded==:auto && length(x) > 4095) || multithreaded==true
        μᵪ = _vtmean(x,:)
        μᵧ = _vtmean(y,:)
        σᵪ = _vtstd(μᵪ, corrected, x, :)
        σᵧ = _vtstd(μᵧ, corrected, y, :)
        σᵪᵧ = _vtcov(x, y, corrected, μᵪ, μᵧ)
    else
        μᵪ = _vmean(x,:)
        μᵧ = _vmean(y,:)
        σᵪ = _vstd(μᵪ, corrected, x, :)
        σᵧ = _vstd(μᵧ, corrected, y, :)
        σᵪᵧ = _vcov(x, y, corrected, μᵪ, μᵧ)
    end
    ρᵪᵧ = σᵪᵧ / (σᵪ * σᵧ)
    return ρᵪᵧ
end

"""
```julia
vcor(X::AbstractMatrix; dims::Int=1, multithreaded=:auto)
```
Compute the (Pearson's product-moment) correlation matrix of the matrix `X`,
along dimension `dims`. As `Statistics.cor`, but vectorized and (optionally)
multithreaded.
"""
function vcor(X::AbstractMatrix; dims::Int=1, corrected::Bool=true, multithreaded=:auto)
    if (multithreaded===:auto && length(A) > 4095) || multithreaded===true
        _vtcor(X, dims, corrected)
    else
        _vcor(X, dims, corrected)
    end
end
export vcor

function _vcor(X, dims, corrected)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Ρ = similar(X, Tₒ, (m, m))
    # Diagonal must be unity
    @inbounds for i = 1:m
        Ρ[i,i] = one(Tₒ)
    end
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            _vmean(view(X,:,d),:)
        end
        σ = ntuple(m) do d
            _vstd(μ[d], corrected, view(X,:,d), :)
        end
        # Fill off-diagonals symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vcov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    elseif dims == 2
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            _vmean(view(X,d,:),:)
        end
        σ = ntuple(m) do d
            _vstd(μ[d], corrected, view(X,d,:), :)
        end
        @inbounds for i = 1:m
            for j = 1:i-1
                σᵢⱼ = _vcov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    else
        throw("Dimension not in range")
    end
    return Ρ
end

function _vtcor(X, dims, corrected)
    Tₒ = Base.promote_op(/, eltype(X), Int)
    n = size(X, dims)
    m = size(X, mod(dims,2)+1)
    Ρ = similar(X, Tₒ, (m, m))
    # Diagonal must be unity
    @inbounds for i = 1:m
        Ρ[i,i] = one(Tₒ)
    end
    # Only two dimensions are possible, so handle each manually
    if dims == 1
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            _vtmean(view(X,:,d),:)
        end
        σ = ntuple(m) do d
            _vtstd(μ[d], corrected, view(X,:,d), :)
        end
        # Fill off-diagonals symmetrically
        @inbounds for i = 1:m
            for j = 1:i
                σᵢⱼ = _vtcov(view(X,:,i), view(X,:,j), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    elseif dims == 2
        # Precalculate means and standard deviations
        μ = ntuple(m) do d
            _vtmean(view(X,d,:),:)
        end
        σ = ntuple(m) do d
            _vtstd(μ[d], corrected, view(X,d,:), :)
        end
        @inbounds for i = 1:m
            for j = 1:i-1
                σᵢⱼ = _vtcov(view(X,i,:), view(X,j,:), corrected, μ[i], μ[j])
                Ρ[i,j] = Ρ[j,i] = σᵢⱼ / (σ[i] * σ[j])
            end
        end
    else
        throw("Dimension not in range")
    end
    return Ρ
end


text:
@defstruct ConcatLayer Layer (
  name :: AbstractString = "concat",
  (dim :: Int = 3, dim >= 1),
  (bottoms :: Vector{Symbol} = [], length(bottoms) >= 2),
  (tops :: Vector{Symbol} = [], length(tops) == 1)
)
@characterize_layer(ConcatLayer,
  can_do_bp  => true
)

type ConcatLayerState <: LayerState
  layer      :: ConcatLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::ConcatLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  # make sure all input blobs has the same tensor-dim
  for i = 2:length(inputs)
    @assert ndims(inputs[i]) == ndims(inputs[1])
  end
  tensor_dim = ndims(inputs[1])
  @assert layer.dim <= tensor_dim

  dims = map(size, inputs)
  dim_total = 0
  for i = 1:tensor_dim
    if i == layer.dim
      dim_total = sum(map(x -> x[i], dims))
    else
      # make sure the non-concating dims are equal
      for j = 2:length(dims)
        @assert dims[j][i] == dims[1][i]
      end
    end
  end

  my_dim = [dims[1]...]
  my_dim[layer.dim] = dim_total
  my_dim = tuple(my_dim...)
  data_type = eltype(inputs[1])

  blobs = Blob[make_blob(backend, data_type, my_dim)]
  blobs_diff = map(diffs) do blob
    if isa(blob, NullBlob)
      NullBlob()
    else
      make_blob(backend, data_type, my_dim)
    end
  end

  return ConcatLayerState(layer, blobs, blobs_diff)
end

function shutdown(backend::Backend, state::ConcatLayerState)
  map(destroy, state.blobs)
  map(destroy, state.blobs_diff)
end

function forward(backend::CPUBackend, state::ConcatLayerState, inputs::Vector{Blob})
  output = state.blobs[1]
  idx = map(x -> 1:x, [size(output)...])
  idx_start = 1
  for i = 1:length(inputs)
    idx_end = idx_start + size(inputs[i], state.layer.dim) - 1
    idx[state.layer.dim] = idx_start:idx_end
    setindex!(output.data, inputs[i].data, idx...)
    idx_start = idx_end + 1
  end
end

function backward(backend::CPUBackend, state::ConcatLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  idx = map(x -> 1:x, [size(state.blobs_diff[1])...])
  idx_start = 1
  for i = 1:length(diffs)
    idx_end = idx_start + size(inputs[i], state.layer.dim) - 1
    idx[state.layer.dim] = idx_start:idx_end
    if !isa(diffs[i], NullBlob)
      diffs[i].data[:] = getindex(state.blobs_diff[1].data, idx...)
    end
    idx_start = idx_end + 1
  end
end


text:
function calc_tderivative!(integrator, cache, dtd1, repeat_step)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack du2,fsalfirst,dT,tf,linsolve_tmp = cache

    # Time derivative
    if !repeat_step # skip calculation if step is repeated
      if DiffEqBase.has_tgrad(f)
        f.tgrad(dT, uprev, p, t)
      else
        tf.uprev = uprev
        tf.p = p
        derivative!(dT, tf, t, du2, integrator, cache.grad_config)
      end
    end

    f(fsalfirst, uprev, p, t)
    integrator.destats.nf += 1
    @.. linsolve_tmp = fsalfirst + dtd1*dT
  end
end

function calc_tderivative(integrator, cache)
  @unpack t,dt,uprev,u,f,p = integrator

  # Time derivative
  if DiffEqBase.has_tgrad(f)
    dT = f.tgrad(uprev, p, t)
  else
    tf = cache.tf
    tf.u = uprev
    tf.p = p
    dT = derivative(tf, t, integrator)
  end
  dT
end

"""
    calc_J(integrator,cache,is_compos)

Interface for calculating the jacobian.

For constant caches, a new jacobian object is returned whereas for mutable
caches `cache.J` is updated. In both cases, if `integrator.f` has a custom
jacobian update function, then it will be called for the update. Otherwise,
either ForwardDiff or finite difference will be used depending on the
`jac_config` of the cache.
"""
function calc_J(integrator, cache::OrdinaryDiffEqConstantCache, is_compos)
  @unpack t,dt,uprev,u,f,p = integrator
  if DiffEqBase.has_jac(f)
    J = f.jac(uprev, p, t)
  else
    J = jacobian(cache.uf,uprev,integrator)
  end
  integrator.destats.njacs += 1
  is_compos && (integrator.eigen_est = opnorm(J, Inf))
  return J
end

"""
    calc_J!(integrator,cache,is_compos)

Interface for calculating the jacobian.

For constant caches, a new jacobian object is returned whereas for mutable
caches `cache.J` is updated. In both cases, if `integrator.f` has a custom
jacobian update function, then it will be called for the update. Otherwise,
either ForwardDiff or finite difference will be used depending on the
`jac_config` of the cache.
"""
function calc_J!(integrator, cache::OrdinaryDiffEqMutableCache, is_compos)
  @unpack t,dt,uprev,u,f,p = integrator
  J = cache.J
  if DiffEqBase.has_jac(f)
    f.jac(J, uprev, p, t)
  else
    @unpack du1,uf,jac_config = cache
    uf.t = t
    uf.p = p
    jacobian!(J, uf, uprev, du1, integrator, jac_config)
  end
  integrator.destats.njacs += 1
  is_compos && (integrator.eigen_est = opnorm(J, Inf))
end

"""
    WOperator(mass_matrix,gamma,J[;transform=false])

A linear operator that represents the W matrix of an ODEProblem, defined as

```math
W = MM - \\gamma J
```

or, if `transform=true`:

```math
W = \\frac{1}{\\gamma}MM - J
```

where `MM` is the mass matrix (a regular `AbstractMatrix` or a `UniformScaling`),
`γ` is a real number proportional to the time step, and `J` is the Jacobian
operator (must be a `AbstractDiffEqLinearOperator`). A `WOperator` can also be
constructed using a `*DEFunction` directly as

    WOperator(f,gamma[;transform=false])

`f` needs to have a jacobian and `jac_prototype`, but the prototype does not need
to be a diffeq operator --- it will automatically be converted to one.

`WOperator` supports lazy `*` and `mul!` operations, the latter utilizing an
internal cache (can be specified in the constructor; default to regular `Vector`).
It supports all of `AbstractDiffEqLinearOperator`'s interface.
"""
mutable struct WOperator{T,
  MType <: Union{UniformScaling,AbstractMatrix},
  GType <: Real,
  JType <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  } <: DiffEqBase.AbstractDiffEqLinearOperator{T}
  mass_matrix::MType
  gamma::GType
  J::JType
  transform::Bool       # true => W = mm/gamma - J; false => W = mm - gamma*J
  inplace::Bool
  _func_cache           # cache used in `mul!`
  _concrete_form         # non-lazy form (matrix/number) of the operator
  WOperator(mass_matrix, gamma, J, inplace; transform=false) = new{eltype(J),typeof(mass_matrix),
    typeof(gamma),typeof(J)}(mass_matrix,gamma,J,transform,inplace,nothing,nothing)
end
function WOperator(f::DiffEqBase.AbstractODEFunction, gamma, inplace; transform=false)
  @assert DiffEqBase.has_jac(f) "f needs to have an associated jacobian"
  if isa(f, Union{SplitFunction, DynamicalODEFunction})
    error("WOperator does not support $(typeof(f)) yet")
  end
  # Convert mass matrix, if needed
  mass_matrix = f.mass_matrix
  if !isa(mass_matrix, Union{AbstractMatrix,UniformScaling})
    mass_matrix = convert(AbstractMatrix, mass_matrix)
  end
  # Convert jacobian, if needed
  J = deepcopy(f.jac_prototype)
  if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
    J = DiffEqArrayOperator(J; update_func=f.jac)
  end
  return WOperator(mass_matrix, gamma, J, inplace; transform=transform)
end

set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W)
DiffEqBase.update_coefficients!(W::WOperator,u,p,t) = (update_coefficients!(W.J,u,p,t); W)
function Base.convert(::Type{AbstractMatrix}, W::WOperator)
  if W._concrete_form === nothing || !W.inplace
    # Allocating
    if W.transform
      W._concrete_form = -W.mass_matrix / W.gamma + convert(AbstractMatrix,W.J)
    else
      W._concrete_form = -W.mass_matrix + W.gamma * convert(AbstractMatrix,W.J)
    end
  else
    # Non-allocating
    _W = W._concrete_form
    J = convert(AbstractMatrix,W.J)
    if W.transform
      if _W isa Diagonal # axpby doesn't specialize on Diagonal matrix
        @inbounds for i in axes(W._concrete_form, 1)
          _W[i, i] = J[i, i] - inv(W.gamma) * W.mass_matrix[i, i]
        end
      else
        copyto!(_W, W.mass_matrix)
        axpby!(one(W.gamma), J, -inv(W.gamma), _W)
      end
    else
      if _W isa Diagonal # axpby doesn't specialize on Diagonal matrix
        @inbounds for i in axes(W._concrete_form, 1)
          _W[i, i] = W.gamma*J[i, i] - W.mass_matrix[i, i]
        end
      else
        copyto!(_W, W.mass_matrix)
        axpby!(W.gamma, J, -one(W.gamma), W._concrete_form)
      end
    end
  end
  W._concrete_form
end
function Base.convert(::Type{Number}, W::WOperator)
  if W.transform
    W._concrete_form = -W.mass_matrix / W.gamma + convert(Number,W.J)
  else
    W._concrete_form = -W.mass_matrix + W.gamma * convert(Number,W.J)
  end
  W._concrete_form
end
Base.size(W::WOperator, args...) = size(W.J, args...)
function Base.getindex(W::WOperator, i::Int)
  if W.transform
    -W.mass_matrix[i] / W.gamma + W.J[i]
  else
    -W.mass_matrix[i] + W.gamma * W.J[i]
  end
end
function Base.getindex(W::WOperator, I::Vararg{Int,N}) where {N}
  if W.transform
    -W.mass_matrix[I...] / W.gamma + W.J[I...]
  else
    -W.mass_matrix[I...] + W.gamma * W.J[I...]
  end
end
function Base.:*(W::WOperator, x::Union{AbstractVecOrMat,Number})
  if W.transform
    (W.mass_matrix*x) / -W.gamma + W.J*x
  else
    -W.mass_matrix*x + W.gamma * (W.J*x)
  end
end
function Base.:\(W::WOperator, x::Union{AbstractVecOrMat,Number})
  if size(W) == () # scalar operator
    convert(Number,W) \ x
  else
    convert(AbstractMatrix,W) \ x
  end
end

function LinearAlgebra.mul!(Y::AbstractVecOrMat, W::WOperator, B::AbstractVecOrMat)
  if W._func_cache === nothing
    # Allocate cache only if needed
    W._func_cache = Vector{eltype(W)}(undef, size(Y, 1))
  end
  if W.transform
    # Compute mass_matrix * B
    if isa(W.mass_matrix, UniformScaling)
      a = -W.mass_matrix.λ / W.gamma
      @.. Y = a * B
    else
      mul!(Y, W.mass_matrix, B)
      lmul!(-1/W.gamma, Y)
    end
    # Compute J * B and add
    mul!(W._func_cache, W.J, B)
    Y .+= W._func_cache
  else
    # Compute mass_matrix * B
    if isa(W.mass_matrix, UniformScaling)
      @.. Y = W.mass_matrix.λ * B
    else
      mul!(Y, W.mass_matrix, B)
    end
    # Compute J * B
    mul!(W._func_cache, W.J, B)
    # Add result
    axpby!(W.gamma, W._func_cache, -one(W.gamma), Y)
  end
end

function do_newJ(integrator, alg::T, cache, repeat_step)::Bool where T # any changes here need to be reflected in FIRK
  repeat_step && return false
  !integrator.opts.adaptive && return true
  !alg_can_repeat_jac(alg) && return true
  isnewton = T <: NewtonAlgorithm
  isnewton && (T <: RadauIIA5 ? ( nlstatus = cache.status ) : ( nlstatus = DiffEqBase.get_status(cache.nlsolver) ))
  nlsolvefail(nlstatus) && return true
  # reuse J when there is fast convergence
  fastconvergence = nlstatus === FastConvergence
  return !fastconvergence
end

function do_newW(integrator, nlsolver::T, new_jac, W_dt)::Bool where T # any changes here need to be reflected in FIRK
  integrator.iter <= 1 && return true
  new_jac && return true
  # reuse W when the change in stepsize is small enough
  dt = integrator.dt
  new_W_dt_cutoff = T <: NLSolver ? nlsolver.cache.new_W_dt_cutoff : #= FIRK =# nlsolver.new_W_dt_cutoff
  smallstepchange = (dt/W_dt-one(dt)) <= new_W_dt_cutoff
  return !smallstepchange
end

@noinline _throwWJerror(W, J) = throw(DimensionMismatch("W: $(axes(W)), J: $(axes(J))"))
@noinline _throwWMerror(W, mass_matrix) = throw(DimensionMismatch("W: $(axes(W)), mass matrix: $(axes(mass_matrix))"))

@inline function jacobian2W!(W::AbstractMatrix, mass_matrix::MT, dtgamma::Number, J::AbstractMatrix, W_transform::Bool)::Nothing where MT
  # check size and dimension
  iijj = axes(W)
  @boundscheck (iijj === axes(J) && length(iijj) === 2) || _throwWJerror(W, J)
  mass_matrix isa UniformScaling || @boundscheck axes(mass_matrix) === axes(W) || _throwWMerror(W, mass_matrix)
  @inbounds if W_transform
    invdtgamma = inv(dtgamma)
    if MT <: UniformScaling
      copyto!(W, J)
      idxs = diagind(W)
      λ = -mass_matrix.λ
      @.. @view(W[idxs]) = muladd(λ, invdtgamma, @view(J[idxs]))
    else
      @.. W = muladd(-mass_matrix, invdtgamma, J)
    end
  else
    if MT <: UniformScaling
      idxs = diagind(W)
      @.. W = dtgamma*J
      λ = -mass_matrix.λ
      @.. @view(W[idxs]) = @view(W[idxs]) + λ
    else
      @.. W = muladd(dtgamma, J, -mass_matrix)
    end
  end
  return nothing
end

function calc_W!(integrator, cache::OrdinaryDiffEqMutableCache, dtgamma, repeat_step, W_transform=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack J,W = cache
  alg = unwrap_alg(integrator, true)
  mass_matrix = integrator.f.mass_matrix
  is_compos = integrator.alg isa CompositeAlgorithm
  isnewton = alg isa NewtonAlgorithm

  # fast pass
  # we only want to factorize the linear operator once
  new_jac = true
  new_W = true
  if (f isa ODEFunction && islinear(f.f)) || (integrator.alg isa SplitAlgorithms && f isa SplitFunction && islinear(f.f1.f))
    new_jac = false
    @goto J2W # Jump to W calculation directly, because we already have J
  end

  # check if we need to update J or W
  W_dt = isnewton ? cache.nlsolver.cache.W_dt : dt # TODO: RosW
  new_jac = isnewton ? do_newJ(integrator, alg, cache, repeat_step) : true
  new_W = isnewton ? do_newW(integrator, cache.nlsolver, new_jac, W_dt) : true

  # calculate W
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    isnewton || DiffEqBase.update_coefficients!(W,uprev,p,t) # we will call `update_coefficients!` in NLNewton
    @label J2W
    W.transform = W_transform; set_gamma!(W, dtgamma)
  else # concrete W using jacobian from `calc_J!`
    new_jac && calc_J!(integrator, cache, is_compos)
    new_W && jacobian2W!(W, mass_matrix, dtgamma, J, W_transform)
  end
  if isnewton
    set_new_W!(cache.nlsolver, new_W) && DiffEqBase.set_W_dt!(cache.nlsolver, dt)
  end
  new_W && (integrator.destats.nw += 1)
  return nothing
end

function calc_W!(integrator, cache::OrdinaryDiffEqConstantCache, dtgamma, repeat_step, W_transform=false)
  @unpack t,uprev,p,f = integrator
  @unpack uf = cache
  mass_matrix = integrator.f.mass_matrix
  isarray = typeof(uprev) <: AbstractArray
  # calculate W
  uf.t = t
  is_compos = typeof(integrator.alg) <: CompositeAlgorithm
  if (f isa ODEFunction && islinear(f.f)) || (f isa SplitFunction && islinear(f.f1.f))
    J = f.f1.f
    W = WOperator(mass_matrix, dtgamma, J, false; transform=W_transform)
  elseif DiffEqBase.has_jac(f)
    J = f.jac(uprev, p, t)
    if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
      J = DiffEqArrayOperator(J)
    end
    W = WOperator(mass_matrix, dtgamma, J, false; transform=W_transform)
    integrator.destats.nw += 1
  else
    integrator.destats.nw += 1
    J = calc_J(integrator, cache, is_compos)
    W_full = W_transform ? -mass_matrix*inv(dtgamma) + J :
                           -mass_matrix + dtgamma*J
    W = W_full isa Number ? W_full : lu(W_full)
  end
  is_compos && (integrator.eigen_est = isarray ? opnorm(J, Inf) : abs(J))
  W
end

function calc_rosenbrock_differentiation!(integrator, cache, dtd1, dtgamma, repeat_step, W_transform)
  calc_tderivative!(integrator, cache, dtd1, repeat_step)
  calc_W!(integrator, cache, dtgamma, repeat_step, W_transform)
  return nothing
end

# update W matrix (only used in Newton method)
update_W!(integrator, cache, dt, repeat_step) =
  update_W!(cache.nlsolver, integrator, cache, dt, repeat_step)

function update_W!(nlsolver::NLSolver, integrator, cache::OrdinaryDiffEqMutableCache, dt, repeat_step)
  if isnewton(nlsolver)
    calc_W!(integrator, cache, dt, repeat_step, true)
  end
  nothing
end

function update_W!(nlsolver::NLSolver, integrator, cache::OrdinaryDiffEqConstantCache, dt, repeat_step)
  if isnewton(nlsolver)
    DiffEqBase.set_W!(nlsolver, calc_W!(integrator, cache, dt, repeat_step, true))
  end
  nothing
end

iip_get_uf(alg::OrdinaryDiffEqAlgorithm,nf,t,p) = DiffEqDiffTools.UJacobianWrapper(nf,t,p)
oop_get_uf(alg::OrdinaryDiffEqAlgorithm,nf,t,p) = DiffEqDiffTools.UDerivativeWrapper(nf,t,p)


text:
#@static Sys.iswindows() ? (const thelib = "C:/programs/compa_libs/webview/compiled/VC12_64/webview_w64") : (const thelib = "libgmt")  # Name of the shared lib.
@static Sys.iswindows() ? (Sys.WORD_SIZE == 64 ? (const thelib = "webview_w64") : (const thelib = "webview_w32")) : (const thelib = "webview")  # Name of the shared lib.

const Void = Cvoid
const NULL = C_NULL

mutable struct RECT
	left::Int32
	top::Int32
	right::Int32
	bottom::Int32
end

mutable struct webview_priv
	hwnd::Ptr{Cvoid}
#	IOleObject **browser;
	browser::Ptr{Ptr{Cvoid}};
	is_fullscreen::Int32
	saved_style::UInt32
	saved_ex_style::UInt32
	saved_rect::RECT
end

mutable struct webview
	url::Ptr{UInt8};
	title::Ptr{UInt8};
	width::Int32
	height::Int32
	resizable::Int32
	debug::Int32
#	external_invoke_cb::webview_external_invoke_cb_t
	external_invoke_cb::Ptr{Cvoid}
	priv::webview_priv
	userdata::Ptr{Cvoid}
end

function webview(title::Ptr{UInt8}, url::Ptr{UInt8}, width::Integer=800, height::Integer=600, resizable::Integer=1)
	ccall((:webview, thelib), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint, Cint),
	                                 title, url, width, height, resizable)
end
function webview(title::String, url::String, width::Integer=800, height::Integer=600, resizable::Integer=1)
	ccall((:webview, thelib), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint, Cint, Cint),
	                                 title, url, width, height, resizable)
end

function webview_init(wv::Ptr{webview})
	ccall((:webview_init, thelib), Cint, (Ptr{webview},), wv)
end

# This one is used so far only to detect if we have Webview lib present or not.
function webview_print_log(t::String="")
	ccall((:webview_print_log, thelib), Cvoid, (Ptr{UInt8},), t)
end


text:
using RecipesBase
using Printf

@recipe function plot(r::Roc{T}; convex_hull=false, traditional=false) where T
    ## defaults
    title --> "ROC"
    xguide --> "Pfa"
    yguide --> (traditional ? "Phit" : "Pmiss")
    legend --> false
    if convex_hull
        x = r.pfa[r.ch]
        y = r.pmiss[r.ch]
    else
        x = r.pfa
        y = r.pmiss
    end
    if traditional
        y = 1 .- y
    end
    (x, y)
end

@userplot DetPlot

@recipe function detplot(d::DetPlot; convex_hull=false)
    if (length(d.args) != 1 || ! (typeof(d.args[1]) <: Roc))
        error("Argument should be of type Roc, got $(length(d.args)) $(typeof(d.args[1]))", )
    end
    r, = d.args
    ##default
    title --> "DET plot"
    xguide --> "Pfa (%)"
    yguide --> "Pmiss (%)"
    legend --> false
    xlims --> (qnorm(0.001), qnorm(0.5))
    ylims --> (qnorm(0.001), qnorm(0.5))
    linewidth --> 2
    ## Not trivial in Plots.jl to control ticks and ticklabels
    ticks --> qnorm.([0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 40] ./ 100)
    formatter --> x -> begin
       p = 100 * pnorm(x)
       if p < 1
           Printf.@sprintf("%3.1f", p)
       else
           Printf.@sprintf("%2.0f", p)
       end
   end
    ## compulsory
    aspect_ratio := :equal

    if convex_hull
        x = qnorm.(r.pfa[r.ch])
        y = qnorm.(r.pmiss[r.ch])
    else
        x = qnorm.(r.pfa)
        y = qnorm.(r.pmiss)
    end
    (x, y)
end

@userplot ApePlot

@recipe function apeplot(ape::ApePlot; xmin = -7, xmax = 7)
    if !(length(ape.args) == 1 && typeof(ape.args[1]) <: Roc)
        error("Expecting one argument of type Roc, got $(typeof(ape.args[1]))")
    end
    r, = ape.args
    title --> "Applied Probability of Error"
    xguide --> "prior log odds"
    yguide --> "Bayes' error rate"
    label --> ["Bayes' Error" "Minimum" "Trivial"]
    linecolor --> [:red :green :black]
    lo = collect(xmin:0.01:xmax)
    be = ber(r, lo)
    mbe = minber(r, lo)
    defbe = 1 ./ ROCAnalysis.normfactor(lo)
    ylims --> (0, 1.1maximum(be))
    x = lo
    y = hcat(be, mbe, defbe)
    (x,y)
end

@userplot NbePlot

@recipe function nbeplot(nbe::NbePlot; xmin = -10, xmax = 5)
    if !(length(nbe.args) == 1 && typeof(nbe.args[1]) <: Roc)
        error("Expecting one argument of type Roc, got $(typeof(nbe.args[1]))")
    end
    r, = nbe.args
    title --> "Normalized Bayes' Error"
    xguide --> "prior log odds"
    yguide --> "normalized Bayes' error"
    label --> ["Bayes' Error" "Minimum BE" "False Alarms" "Misses" "Minimum FA" "Minimum miss"]
    legend --> :left
    linestyle --> [:solid :solid :dash :dot :dash :dot]
    linecolor --> [:red :green :red :red :green :green]
    lo = collect(xmin:0.01:xmax)
    norm = ROCAnalysis.normfactor(lo)
    bfa, bmiss = [be .* norm for be in ROCAnalysis.ber_famiss(r, lo)]
    nbe = bfa + bmiss
    mbfa, mbmiss = [mbe .* norm for mbe in ROCAnalysis.minber_famiss(r, lo)]
    mnbe = minber(r, lo) .* norm
    x = lo
    y = hcat(nbe, mnbe, bfa, bmiss, mbfa, mbmiss)
    (x, y)
end

@userplot LlrPlot

@recipe function llrplot(llr::LlrPlot)
    if !(length(llr.args) == 1 && typeof(llr.args[1]) <: Roc)
        error("Expecting one argument of type Roc, got $(typeof(nbe.args[1]))")
    end
    r, = llr.args
    title --> "LLR plot"
    xguide --> "score"
    yguide --> "log likelihood ratio"
    legend --> false
    llr_extrema = extrema(filter(x -> !isinf(x), r.llr))
    mi = max(llr_extrema[1], minimum(r.θ)) ## max of mins
    ma = min(llr_extrema[2], maximum(r.θ)) ## min of maxs
    xlims --> (mi, ma)
    ylims --> (mi, ma)
    (r.θ, r.llr)
end


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using BasicPOMCP
using POMCPOW
using ParticleFilters
using CPUTime
using DataFrames
using ARDESPOT
using Distributions
using DiscreteValueIteration
using SubHunt
using QMDP

@everywhere begin
    using POMDPs
    using ParticleFilters
    using POMDPToolbox
    using ContinuousPOMDPTreeSearchExperiments
end

pomdp = SubHuntPOMDP()
vs = ValueIterationSolver()
if !isdefined(:vp) || vp.mdp != pomdp
    vp = solve(vs, pomdp, verbose=true)
end
qs = QMDPSolver()
qp = QMDP.create_policy(qs, pomdp)
qp.alphas[:] = vp.qmat

function gen_sims(x::Vector{Float64}, n, k, seed)
    c = x[1]
    @assert c >= 0.0
    k_obs = x[2]
    @assert k_obs >= 1.0
    inv_alpha_obs = x[3]
    @assert inv_alpha_obs >= 0.1

    sims = []

    for i in 1:n
        rng = MersenneTwister(0)
        ro = ValueIterationSolver()

        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(c),
                               final_criterion=MaxTries(),
                               max_depth=20,
                               max_time=1.0,
                               enable_action_pw=false,
                               k_observation=k_obs,
                               alpha_observation=1/inv_alpha_obs,
                               estimate_value=FOValue(vp),
                               check_repeat_obs=false,
                               default_action=ReportWhenUsed(qp),
                               rng=rng
                              )

        planner = solve(solver, pomdp)
        filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(10_000),
                                           0.05,
                                           MersenneTwister(i+10000*k))

        srand(planner, i+40000*k)
        sim = Sim(deepcopy(pomdp),
                  planner,
                  filter,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=100,
                  metadata=Dict(:i=>i, :k=>k)
                 )
        push!(sims, sim)
    end

    return sims
end

# start_mean = [100.0, 2.0, 10.0]
# start_cov = diagm([100.0^2, 10.0^2, 20.0^2])
start_mean = [17.0, 6.0, 50.0]
start_cov = diagm([20.0^2, 10.0^2, 50.0^2])
d = MvNormal(start_mean, start_cov)
K = 120 # 60 # number of parameter samples
n = 100 # 100 # number of evaluation simulations
m = 15  # 15 # number of elite samples
max_iters = 100

for i in 1:max_iters
    sims = []
    params = Vector{Vector{Float64}}(K)
    print("creating $K simulation sets")
    for k in 1:K
        p = rand(d)
        p[1] = max(0.0, p[1])
        p[2] = max(1.0, p[2])
        p[3] = max(0.1, p[3])
        params[k] = p
        k_sims = gen_sims(p, n, k, i)
        print(".")
        append!(sims, k_sims)
    end
    println()
    results = run_parallel(sims)
    # results = run(sims)
    combined = by(results, :k) do df
        DataFrame(mean_reward=mean(df[:reward]))
    end
    @show mean(combined[:mean_reward])
    order = sortperm(combined[:mean_reward])
    elite = params[combined[:k][order[K-m:end]]]
    elite_matrix = Matrix{Float64}(length(start_mean), m)
    for k in 1:m
        elite_matrix[:,k] = elite[k]
    end
    try
        d = fit(typeof(d), elite_matrix)
    catch ex
        if ex isa Base.LinAlg.PosDefException
            println("pos def exception")
            d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
        else
            rethrow(ex)
        end
    end
    println("iteration $i")
    @show mean(d)
    @show det(cov(d))
    @show ev = eigvals(cov(d))
    for j in 1:length(ev)
        @show eigvecs(cov(d))[:,j]
    end
end


text:
using PQEcolaPoint
using Test

# -- Default test ------------------------------------------------------ #
function default_pqecolapoint_test()
    return true
end
# ------------------------------------------------------------------------------- #

@testset "default_test_set" begin
    @test default_pqecolapoint_test() == true
end


text:
include("all_different/bipartite.jl")
include("all_different/scc.jl")

function init_constraint_struct(::AllDifferentSetInternal, internals)
    AllDifferentConstraint(
        internals,
        Int[], # pval_mapping will be filled later
        Int[], # vertex_mapping => later
        Int[], # vertex_mapping_bw => later
        Int[], # di_ei => later
        Int[], # di_ej => later
        MatchingInit(),
        SCCInit(),
        Int[],
    )
end

"""
    init_constraint!(com::CS.CoM, constraint::AllDifferentConstraint, fct::MOI.VectorOfVariables, set::AllDifferentSetInternal)

Initialize the AllDifferentConstraint by filling matching_init
"""
function init_constraint!(
    com::CS.CoM,
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal
)
    pvals = constraint.pvals
    nindices = length(constraint.indices)

    min_pvals, max_pvals = extrema(pvals)
    len_range = max_pvals - min_pvals + 1

    # we need the + length(pvals) for scc edges from the new vertex
    num_edges = length(pvals) * (nindices + 1)

    constraint.pval_mapping = zeros(Int, length(pvals))
    constraint.vertex_mapping = zeros(Int, len_range)
    constraint.vertex_mapping_bw = zeros(Int, nindices + length(pvals))
    constraint.di_ei = zeros(Int, num_edges)
    constraint.di_ej = zeros(Int, num_edges)

    # fill matching_init
    m = nindices
    n = len_range

    constraint.matching_init = MatchingInit(
        0,
        zeros(Int, m),
        zeros(Int, n),
        zeros(Int, m + 1),
        zeros(Int, m + n),
        zeros(Int, m + n),
        zeros(Int, m + n),
        zeros(Bool, m),
        zeros(Bool, n),
    )

    constraint.scc_init = SCCInit(
        zeros(Int, m + n + 2),
        zeros(Int, m + n + 1),
        zeros(Int, m + n + 1),
        zeros(Bool, m + n + 1),
        zeros(Int, m + n + 1),
    )

    # check if lp model exists and then add an equality constraint for better bound computation
    com.lp_model === nothing && return true # return feasibility

    lp_backend = backend(com.lp_model)
    lp_vidx = create_lp_variable!(com.lp_model, com.lp_x)
    # create == constraint with sum of all variables equal the newly created variable
    sats =
        [MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(vidx)) for vidx in constraint.indices]
    push!(sats, MOI.ScalarAffineTerm(-1.0, MOI.VariableIndex(lp_vidx)))
    saf = MOI.ScalarAffineFunction(sats, 0.0)
    MOI.add_constraint(lp_backend, saf, MOI.EqualTo(0.0))

    constraint.bound_rhs = [BoundRhsVariable(lp_vidx, typemin(Int), typemax(Int))]

    # if constraints are part of the all different constraint
    # the all different constraint can be split more parts to get better bounds
    # i.e https://github.com/Wikunia/ConstraintSolver.jl/issues/114
    for sc_idx in constraint.sub_constraint_idxs
        lp_vidx = create_lp_variable!(com.lp_model, com.lp_x)
        # create == constraint with sum of all variables equal the newly created variable
        sats = [
            MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(vidx))
            for vidx in com.constraints[sc_idx].indices
        ]
        push!(sats, MOI.ScalarAffineTerm(-1.0, MOI.VariableIndex(lp_vidx)))
        saf = MOI.ScalarAffineFunction(sats, 0.0)
        MOI.add_constraint(lp_backend, saf, MOI.EqualTo(0.0))
        push!(constraint.bound_rhs, BoundRhsVariable(lp_vidx, typemin(Int), typemax(Int)))
    end
    return true # still feasible
end

"""
    get_alldifferent_extrema(sorted_min, sorted_max, len)

`sorted_max` should be desc and `sorted_min` asc
Return the minimum and maximum sum using `len` values of sorted_min/sorted_max while satisfying the all different constraint
"""
function get_alldifferent_extrema(sorted_min, sorted_max, len)
    max_sum = sorted_max[1]
    last_val = max_sum
    for i in 2:len
        if sorted_max[i] >= last_val
            last_val -= 1
        else
            last_val = sorted_max[i]
        end
        max_sum += last_val
    end

    min_sum = sorted_min[1]
    last_val = min_sum
    for i in 2:len
        if sorted_min[i] <= last_val
            last_val += 1
        else
            last_val = sorted_min[i]
        end
        min_sum += last_val
    end

    return min_sum, max_sum
end

function get_sorted_extrema(
    com,
    constraint::AllDifferentConstraint,
    vidx::Int,
    lb::Int,
    ub::Int,
)
    max_vals = zeros(Int, length(constraint.indices))
    min_vals = zeros(Int, length(constraint.indices))
    search_space = com.search_space
    for i in 1:length(constraint.indices)
        v_idx = constraint.indices[i]
        if v_idx == vidx
            max_vals[i] = ub
            min_vals[i] = lb
        else
            min_vals[i] = search_space[v_idx].min
            max_vals[i] = search_space[v_idx].max
        end
    end

    # sort the max_vals desc and obtain bound by enforcing all different
    sort!(max_vals; rev = true)
    # sort the min_vals asc and obtain bound by enforcing all different
    sort!(min_vals)
    return min_vals, max_vals
end

"""
    update_best_bound_constraint!(com::CS.CoM,
        constraint::AllDifferentConstraint,
        fct::MOI.VectorOfVariables,
        set::AllDifferentSetInternal,
        vidx::Int,
        lb::Int,
        ub::Int
    )

Update the bound constraint associated with this constraint. This means that the `bound_rhs` bounds will be changed according to
the possible values the all different constraint allows.
i.e if we have 4 variables all between 1 and 10 the maximum sum is 10+9+8+7 and the minimum sum is 1+2+3+4
Additionally one of the variables can be bounded using `vidx`, `lb` and `ub`
"""
function update_best_bound_constraint!(
    com::CS.CoM,
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal,
    vidx::Int,
    lb::Int,
    ub::Int,
)
    isempty(constraint.bound_rhs) && return
    search_space = com.search_space

    # compute bounds
    # get the maximum/minimum value for each variable
    min_vals, max_vals = get_sorted_extrema(com, constraint, vidx, lb, ub)

    min_sum, max_sum =
        get_alldifferent_extrema(min_vals, max_vals, length(constraint.indices))

    constraint.bound_rhs[1].lb = min_sum
    constraint.bound_rhs[1].ub = max_sum

    i = 1
    for sc_idx in constraint.sub_constraint_idxs
        i += 1
        sub_constraint = com.constraints[sc_idx]
        min_sum, max_sum =
            get_alldifferent_extrema(min_vals, max_vals, length(sub_constraint.indices))
        constraint.bound_rhs[i].lb = min_sum
        constraint.bound_rhs[i].ub = max_sum
    end
end

"""
    prune_constraint!(com::CS.CoM, constraint::AllDifferentConstraint, fct::MOI.VectorOfVariables, set::AllDifferentSetInternal; logs = true)

Reduce the number of possibilities given the `AllDifferentConstraint`.
Return whether still feasible and throws a warning if infeasible and `logs` is set to `true`
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal;
    logs = true,
)
    indices = constraint.indices
    pvals = constraint.pvals
    nindices = length(indices)

    search_space = com.search_space
    fixed_vals, unfixed_indices = fixed_vs_unfixed(search_space, indices)

    fixed_vals_set = Set(fixed_vals)
    # check if one value is used more than once
    if length(fixed_vals_set) < length(fixed_vals)
        logs && @warn "The problem is infeasible"
        return false
    end

    bfixed = true
    current_fixed_vals = fixed_vals

    while bfixed
        new_fixed_vals = Int[]
        bfixed = false
        for i in 1:length(unfixed_indices)
            pi = unfixed_indices[i]
            vidx = indices[pi]
            @views c_search_space = search_space[vidx]
            if !CS.isfixed(c_search_space)
                for pv in current_fixed_vals
                    if has(c_search_space, pv)
                        if !rm!(com, c_search_space, pv)
                            logs && @warn "The problem is infeasible"
                            return false
                        end

                        if nvalues(c_search_space) == 1
                            only_value = CS.value(c_search_space)
                            push!(fixed_vals, only_value)
                            push!(new_fixed_vals, only_value)
                            bfixed = true
                            break
                        end
                    end
                end
            end
        end
        current_fixed_vals = new_fixed_vals
    end

    if length(fixed_vals) == nindices
        return true
    end

    min_pvals, max_pvals = extrema(pvals)
    len_range = max_pvals - min_pvals + 1

    # find maximum_matching for infeasible check and Berge's lemma
    # building graph
    # ei = indices, ej = possible values
    # i.e ei=[1,2,1] ej = [1,2,2] => 1->[1,2], 2->[2]
    min_pvals_m1 = min_pvals - 1

    pval_mapping = constraint.pval_mapping
    vertex_mapping = constraint.vertex_mapping
    vertex_mapping_bw = constraint.vertex_mapping_bw

    vc = 0
    for i in indices
        vc += 1
        vertex_mapping_bw[vc] = i
    end
    pvc = 0
    for pv in pvals
        vc += 1
        pvc += 1
        pval_mapping[pvc] = pv
        vertex_mapping[pv - min_pvals_m1] = vc
        vertex_mapping_bw[vc] = pv
    end
    num_nodes = vc

    # count the number of edges
    num_edges = 0
    @inbounds for i in indices
        num_edges += nvalues(search_space[i])
    end

    di_ei = constraint.di_ei
    di_ej = constraint.di_ej
    di_ei .= 0
    di_ej .= 0

    # add edge from each index to the possible values
    edge_counter = 0
    vc = 0
    @inbounds for i in indices
        vc += 1
        if isfixed(search_space[i])
            edge_counter += 1
            di_ei[edge_counter] = vc
            di_ej[edge_counter] =
                vertex_mapping[CS.value(search_space[i]) - min_pvals_m1] - nindices
        else
            for pv in view_values(search_space[i])
                edge_counter += 1
                di_ei[edge_counter] = vc
                di_ej[edge_counter] = vertex_mapping[pv - min_pvals_m1] - nindices
            end
        end
    end

    matching_init = constraint.matching_init
    matching_init.l_in_len = num_edges
    # there have to be at least as many values `len_range` as variables `vc`
    vc > len_range && return false
    maximum_matching = bipartite_cardinality_matching(
        di_ei,
        di_ej,
        vc,
        len_range;
        l_sorted = true,
        matching_init = matching_init,
    )
    if maximum_matching.weight != nindices
        logs && @warn "Infeasible (No maximum matching was found)"
        return false
    end

    # directed edges for strongly connected components
    # Important: di_ei must be sorted asc for strongly connected components
    vc = 0
    edge_counter = 0
    @inbounds for i in indices
        vc += 1
        if isfixed(search_space[i])
            edge_counter += 1
            di_ei[edge_counter] = vc
            di_ej[edge_counter] = vertex_mapping[CS.value(search_space[i]) - min_pvals_m1]
        else
            edge_counter += 1
            pv = pval_mapping[maximum_matching.match[vc]]
            di_ei[edge_counter] = vc
            di_ej[edge_counter] = vertex_mapping[pv - min_pvals_m1]
        end
    end

    new_vertex = num_nodes + 1
    used_in_maximum_matching = zeros(Bool, len_range)
    @inbounds for pv in pvals
        vc = 0
        for i in indices
            vc += 1
            if has(search_space[i], pv)
                if pv != pval_mapping[maximum_matching.match[vc]]
                    edge_counter += 1
                    di_ei[edge_counter] = vertex_mapping[pv - min_pvals_m1]
                    di_ej[edge_counter] = vc
                elseif length(pvals) > nindices # if we have more values than indices
                    edge_counter += 1
                    di_ei[edge_counter] = vertex_mapping[pv - min_pvals_m1]
                    di_ej[edge_counter] = new_vertex
                    used_in_maximum_matching[pv - min_pvals_m1] = true
                end
            end
        end
    end

    # INFO for both: length(pvals) > nindices
    # add extra node which has all values as input which are used in the maximum matching
    # and the ones in maximum matching as inputs see: http://www.minicp.org (Part 6)
    # the direction is opposite to that of minicp

    # if we have more values than indices
    if length(pvals) > nindices
        for pv in pvals
            # value not in maximum matching
            if !used_in_maximum_matching[pv - min_pvals_m1]
                edge_counter += 1
                di_ei[edge_counter] = new_vertex
                di_ej[edge_counter] = vertex_mapping[pv - min_pvals_m1]
            end
        end
    end

    # Important di_ei must be sorted asc !!!
    sccs_map = scc(di_ei, di_ej, constraint.scc_init)

    # remove the left over edges from the search space
    vmb = vertex_mapping_bw
    for (src, dst) in zip(di_ei, di_ej)
        # only zeros coming afterwards
        if src == 0
            break
        end
        # edges to the extra node don't count
        if src > num_nodes || dst > num_nodes
            continue
        end

        # if in same strong component -> part of a cycle -> part of a maximum matching
        if sccs_map[src] == sccs_map[dst]
            continue
        end
        #  remove edges in maximum matching and then edges which are part of a cycle
        if src <= nindices &&
           dst == vertex_mapping[pval_mapping[maximum_matching.match[src]] - min_pvals_m1]
            continue
        end

        cind = vmb[dst]
        if !rm!(com, search_space[cind], vmb[src])
            logs && @warn "The problem is infeasible"
            return false
        end
    end

    return true
end

"""
    still_feasible(com::CoM, constraint::AllDifferentConstraint, fct::MOI.VectorOfVariables, set::AllDifferentSetInternal, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal,
    vidx::Int,
    value::Int,
)
    indices = constraint.indices
    for i in 1:length(indices)
        if indices[i] == vidx
            continue
        end
        if issetto(com.search_space[indices[i]], value)
            return false
        end
    end
    return true
end

function is_constraint_solved(
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal,
    values::Vector{Int},
)
    return allunique(values)
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::AllDifferentConstraint,
        fct::MOI.VectorOfVariables,
        set::AllDifferentSetInternal
    )

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::AllDifferentConstraint,
    fct::MOI.VectorOfVariables,
    set::AllDifferentSetInternal,
)
    return !allunique(
        CS.value(var) for var in com.search_space[constraint.indices] if isfixed(var)
    )
end


text:
############################################################################################
"""
$(SIGNATURES)
Check gradient computations of different backends against 'objective'.

# Examples
```julia
```

"""
function check_gradients(
    _rng::Random.AbstractRNG,
    objective::Objective,
    ADlibraries = [:ForwardDiff, :ReverseDiff, :Zygote],
    θᵤ = randn(_rng, length(objective)),
    difftune = map(backend -> DiffObjective(objective, AutomaticDiffTune(backend, objective)), ADlibraries);
    printoutput = true
)
## Compute Gradients
    ℓobjectiveresults = ℓGradientResult[]
    for iter in eachindex(difftune)
        push!(ℓobjectiveresults, log_density_and_gradient(difftune[iter], θᵤ))
        if printoutput
            println(ADlibraries[iter], " gradient call succesfull.")
        end
    end
## Check differences
    ℓobjective_diff = map(
        iter -> ℓobjectiveresults[1].ℓθᵤ - ℓobjectiveresults[iter].ℓθᵤ, eachindex(ℓobjectiveresults)
    )
    for iter in eachindex(ℓobjectiveresults)
        if printoutput
            println("Log objective result difference of ", ADlibraries[1], " against ", ADlibraries[iter], ": ", ℓobjective_diff[iter])
        end
    end

    ℓobjective_gradient_diff = map(
        iter -> sum(abs.(ℓobjectiveresults[1].∇ℓθᵤ .- ℓobjectiveresults[iter].∇ℓθᵤ)), eachindex(ℓobjectiveresults)
    )
    for iter in eachindex(ℓobjective_gradient_diff)
        if printoutput
            println("Log objective gradient difference of ", ADlibraries[1], " against ", ADlibraries[iter], ": ", ℓobjective_gradient_diff[iter])
        end
    end
## Compare against base Forward and ReverseDiff
    grad_fd = ForwardDiff.gradient(objective, θᵤ)
    grad_rd = ReverseDiff.gradient(objective, θᵤ)
    fdrd_diff = map(result ->
        (sum(abs.(result.∇ℓθᵤ .- grad_fd)), sum(abs.(result.∇ℓθᵤ .- grad_rd))), ℓobjectiveresults
    )
    for iter in eachindex(fdrd_diff)
        if printoutput
            println("Log objective gradient difference of ", ADlibraries[iter], " against Forward/Reverse call: ", fdrd_diff[iter])
        end
    end
## Return differences
    return (
        names = ADlibraries,
        difftune = difftune,
        ℓobjectiveresults = ℓobjectiveresults,
        ℓobjective_diff = ℓobjective_diff,
        ℓobjective_gradient_diff = ℓobjective_gradient_diff,
        Forward_Reverse_diff = fdrd_diff,
    )
end

############################################################################################
# Export
export check_gradients


text:
mutable struct Timer
    starttime_ns::typeof(Base.time_ns())
    paused_elapsed_ns::typeof(Base.time_ns())
    Timer() = new(0,0)
end

function start!(timer::Timer)
    timer.starttime_ns = (Base.time_ns)()
    return nothing
end
started(timer::Timer) = (timer.starttime_ns ≠ 0)

""" Return seconds since timer was started or 0 if not yet started. """
function elapsed(timer::Timer)
    local elapsedtime_ns = (Base.time_ns)() - timer.starttime_ns
    return started(timer) * float(elapsedtime_ns) / 1000000000
end

function pause!(timer::Timer)
    timer.paused_elapsed_ns = (Base.time_ns)() - timer.starttime_ns
    return nothing
end
function unpause!(timer::Timer)
    timer.starttime_ns = (Base.time_ns)()
    timer.starttime_ns -= timer.paused_elapsed_ns;
    return nothing
end

t = Timer()
start!(t)
elapsed(t)
pause!(t)
unpause!(t)
elapsed(t)

max_stars_repo_name: NHDaly/pongClone

max_issues_repo_name: NHDaly/pongClone

max_forks_repo_name: NHDaly/pongClone


text:
using ControlSystems, OrdinaryDiffEq, NLopt, BlackBoxOptim, ForwardDiff
p0          = [0.2,0.8,1] # Initial guess
K(kp,ki,kd) = pid(kp=kp, ki=ki, kd=kd)
K(p)        = K(p...)

# Define process model
ζ  = 0.1
ω  = 1.; ω² = ω^2
P  = tf(ω²,[1, 2ζ*ω, ω²])*tf(1,[1,1])

Ω  = exp10.(LinRange(-1,2,150))  # Frequency vector to eval constraints
h  = 0.1 # Sample time for time-domain evaluation
Tf = 60.  # Time horizon
t  = 0:h:Tf-h

Ms = 1.4 # Maximum allowed magnitude of sensitivity function
Mt = 1.4 # Maximum allowed magnitude of complimentary sensitivity function

p  = copy(p0)

function timedomain(p)
    C     = K(p[1], p[2], p[3])
    S     = 1/(1+P*C) # Sensitivity fun
    PS    = ss(P*S)   # TF from load disturbance to output
    s     = Simulator(PS, (t,x) -> [1]) # Sim. unit step load disturbance
    ty    = eltype(p) # So that all inputs to solve have same numerical type (ForwardDiff.Dual)
    x0    = zeros(PS.nx) .|> ty
    tspan = (ty(0.),ty(Tf))
    sol   = solve(s, x0, tspan, Tsit5())
    y     = PS.C*sol(t) # y = C*x
    C,y
end

function freqdomain(p)
    C     = K(p[1], p[2], p[3])
    S     = 1/(1+P*C) # Sensitivity fun
    T     = tf(1.) - S# Comp. Sensitivity fun
    Sw    = vec(bode(S,Ω)[1]) # Freq. domain constraints
    Tw    = vec(bode(T,Ω)[1]) # Freq. domain constraints
    Sw,Tw
end

# Evaluates and plots candidate controller
function evalsol(p::Vector)
    C,y = timedomain(p)
    Sw,Tw = freqdomain(p)
    plot(t,y', layout=2, show=false)
    plot!(Ω, [Sw Tw] , lab=["Sw" "Tw"], subplot=2, xscale=:log10, yscale=:log10, show=false)
    plot!([Ω[1],Ω[end]], [Ms,Ms], c = :black, l=:dash, subplot=2, show=false, lab="Ms")
    plot!([Ω[1],Ω[end]], [Mt,Mt], c = :purple, l=:dash, subplot=2, lab="Mt")
    gui()
    false
end
evalsol(res::BlackBoxOptim.OptimizationResults) = evalsol(best_candidate(res))

function constraintfun(p)
    Sw,Tw = freqdomain(p)
    [maximum(Sw)-Ms; maximum(Tw)-Mt]
end

function costfun(p)
    C,y = timedomain(p)
    mean(abs,y) # ~ Integrated absolute error IAE
end

function runopt(p, costfun, constraintfun;
    f_tol = 1e-5,
    x_tol = 1e-3,
    c_tol = 1e-8,
    f_cfg = ForwardDiff.GradientConfig(costfun, p),
    g_cfg = ForwardDiff.JacobianConfig(constraintfun, p),
    lb = zeros(length(p)))

    c1 = constraintfun(p)
    np = length(p)
    nc = length(c1)

    function f(p::Vector, grad::Vector)
        if length(grad) > 0
            grad .= ForwardDiff.gradient(costfun,p,f_cfg)
        end
        costfun(p)
    end

    function c(result, p::Vector, grad)
        if length(grad) > 0
            grad .= ForwardDiff.jacobian(constraintfun,p,g_cfg)'
        end
        result .= constraintfun(p)
    end

    opt = Opt(:LD_SLSQP, np)
    lower_bounds!(opt, lb)
    xtol_rel!(opt, x_tol)
    ftol_rel!(opt, f_tol)

    min_objective!(opt, f)
    inequality_constraint!(opt, c, c_tol*ones(nc))
    minf,minx,ret = NLopt.optimize(opt, p)
end

f_cfg = ForwardDiff.GradientConfig(costfun, p)
g_cfg = ForwardDiff.JacobianConfig(constraintfun, p)
@time minf,minx,ret = runopt(1p0, costfun, constraintfun, x_tol=1e-6, c_tol=1e-12, f_cfg=f_cfg, g_cfg=g_cfg)
evalsol(minx)

# # Optimize costfun using derivative-free method
# res1 = compare_optimizers(costfun; SearchRange = (0.,2.), NumDimensions = length(p0), MaxTime = 20.0)
# res2 = bboptimize(costfun, NumDimensions = length(p0), MaxTime = 20.0)
# evalsol(res2)
# p = best_candidate(res2)