text:
export get_independent_variables
export get_harmonic_equations
export is_rearranged
export slow_flow, slow_flow!
export _equations_without_brackets

show(eom::HarmonicEquation) = show_fields(eom)

"""
    harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")

Expand each variable of `diff_eom` using the harmonics assigned to it with `time` as the time variable.
For each harmonic of each variable, an instance of `HarmonicVariable` (describing a pair of variables (u,v)) is automatically created and named.

`coordinates` allows for using different coordinate systems (e.g. 'polars') - CURRENTLY INACTIVE
"""
harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian") = harmonic_ansatz(eom, time, Dict([[var, coordinates] for var in get_variables(eom)]))

function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num, trafo_types::Dict{Num,String})
    !is_harmonic(diff_eom, time) && error("The differential equation is not harmonic in ", time, " !")
    old_eqs = [diff_eom.equations[var] for var in get_variables(diff_eom)]
    eqs, vars = deepcopy(old_eqs), []
    rules = Dict()

    idx_counter = 1 # keep count to label new variables
    for var in get_variables(diff_eom)
        to_substitute = 0 # holds all the subtitution rules
        for ω in diff_eom.harmonics[var] # 2 new variables are created for each ω of each variable
            these_names = coordinate_names[trafo_types[var]] # stores two strings denoting the variables as either "a", "ϕ", "u", "v"
            unique_symbols = [  these_names[1] * string(idx_counter) , these_names[2] * string(idx_counter) ]
            repl_rule, new_vars = _rotate_variable(var, ω, time, trafo_types[var], new_symbols=unique_symbols) # create the new variables
            to_substitute += repl_rule
            append!(vars, [new_vars])
            idx_counter += 1
        end
        rules[var] = to_substitute
    end
    eqs = substitute_all(eqs, rules)
    HarmonicEquation(eqs, Vector{HarmonicVariable}(vars), diff_eom)
end

function slow_flow!(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)
    eom.equations = expand_derivatives.(eom.equations) # expand all the derivatives

    # fast_time => slow_time for derivatives up to degree-1
    replace0 = [var => substitute_all(var, fast_time => slow_time) for var in get_variables(eom)]
    replace_degrees = [Differential(fast_time)^deg => Differential(slow_time)^deg for deg in 1:degree-1]
    replace = flatten([replace0, replace_degrees])

    # degree derivatives are removed
    drop = [d(var, fast_time, degree) => 0 for var in get_variables(eom)]

    eom.equations = substitute_all(substitute_all(eom.equations, drop), replace)
    eom.variables = substitute_all(eom.variables, replace)
end

"""
    slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)

Removes all derivatives w.r.t `fast_time` (and their products) in `eom` of power `degree`.
In the remaining derivatives, `fast_time` is replaced by `slow_time`.
"""
function slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)::HarmonicEquation
    new_eq = deepcopy(eom)
    slow_flow!(new_eq, fast_time=fast_time, slow_time=slow_time, degree=degree)
    new_eq
end

#   Drop powers of `var` of degree >= `deg` from the equation set in `eom`.
function drop_powers(eom::HarmonicEquation, var, deg::Integer)
    new_eom = deepcopy(eom)
    new_eom.equations = drop_powers(eom.equations, var, deg)
    new_eom
end

"Rearrange an equation system such that the field equations is equal to the vector specified in new_lhs"
function rearrange!(eom::HarmonicEquation, new_rhs::Vector{Num})
    soln = Symbolics.solve_for(eom.equations, new_rhs, simplify=false,check=true)
    eom.equations = soln .~ new_rhs
    return nothing
end

function rearrange(eom::HarmonicEquation, new_rhs::Vector{Num})
    new_eom = deepcopy(eom)
    rearrange!(new_eom, new_rhs)
    return new_eom
end

"""
$(TYPEDSIGNATURES)
Check if `eom` is rearranged to the standard form, such that the derivatives of the variables are on one side.
"""
function is_rearranged(eom::HarmonicEquation)
    tvar = get_independent_variables(eom)[1]
    isequal(getfield.(eom.equations, :rhs), d(get_variables(eom), tvar))
end

"""
$(TYPEDSIGNATURES)
Rearrange `eom` to the standard form, such that the derivatives of the variables are on one side.
"""
function rearrange_standard(eom::HarmonicEquation)
    tvar = get_independent_variables(eom)[1]
    dvars =  d(get_variables(eom), tvar)
    rearrange(eom, dvars)
end

"""
$(TYPEDSIGNATURES)
Get the internal symbols of the independent variables of `eom`.
"""
function get_variables(eom::HarmonicEquation)
    return flatten(get_variables.(eom.variables))
end

"Get the parameters (not time nor variables) of a HarmonicEquation"
function _parameters(eom::HarmonicEquation)
    all_symbols = flatten([cat(get_variables(eq.lhs), get_variables(eq.rhs), dims=1) for eq in eom.equations])
    # subtract the set of independent variables (i.e., time) from all free symbols
    setdiff(all_symbols, get_variables(eom), get_independent_variables(eom))
end

###
# Extending Symbolics.jl's simplify and substitute
###

"Apply `rules` to both `equations` and `variables` field of `eom`"
function substitute_all!(eom::HarmonicEquation, rules::Dict)
    eom.equations = replace(eom.equations, rules)
    eom.variables = substitute_all(eom.variables, rules)
end

"Simplify the equations in HarmonicEquation."
function simplify!(eom::HarmonicEquation)
    eom.equations = [simplify(eq) for eq in eom.equations]
end

"""
$(TYPEDSIGNATURES)
Return the independent variables (typically time) of `eom`.
"""
function get_independent_variables(eom::HarmonicEquation)::Vector{Num}
    dynamic_vars = flatten(getfield.(eom.variables, Symbol("symbols")))
    return flatten(unique([SymbolicUtils.arguments(var.val) for var in dynamic_vars]))
end

"""
$(TYPEDSIGNATURES)
Extract the Fourier components of `eom` corresponding to the harmonics specified in `eom.variables`.
For each harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients).
`time` does not appear in the resulting equations anymore.

Underlying assumption: all time-dependences are harmonic.
"""
function fourier_transform(eom::HarmonicEquation, time::Num)
    new_eom = deepcopy(eom)
    fourier_transform!(new_eom, time)
    return new_eom
end

function fourier_transform!(eom::HarmonicEquation, time::Num)
    avg_eqs = Vector{Equation}(undef, 2*length(eom.variables))

    natural_variables = flatten([[x,x] for x in getfield.(eom.variables, :natural_variable)])
    equation_assignments = flatten([findall(x -> isequal(x, var), collect(keys(eom.natural_equation.equations))) for var in natural_variables])

    for i in 1:length(avg_eqs)
        j = Int(ceil(i/2))
        ω = eom.variables[j].ω
        eq = eom.equations[equation_assignments[i]]

        if isodd(i)
            avg_eqs[i] = fourier_cos_term(eq, ω, time)
        else
            avg_eqs[i] = fourier_sin_term(eq, ω, time)
        end
    end

    eom.equations = avg_eqs
end

"""
    get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothing, slow_time=nothing)

Apply the harmonic ansatz, followed by the slow-flow, Fourier transform and dropping
higher-order derivatives to obtain
a set of ODEs (the harmonic equations) governing the harmonics of `diff_eom`.

The harmonics evolve in `slow_time`, the oscillating terms themselves in `fast_time`.
If no input is used, a variable T is defined for `slow_time` and `fast_time` is taken as the independent variable
of `diff_eom`.

By default, all products of order > 1 of `slow_time`-derivatives are dropped,
which means the equations are linear in the time-derivatives.

# Example
```julia-repl
julia> @variables t, x(t), ω0, ω, F;

# enter the simple harmonic oscillator
julia> diff_eom = DifferentialEquation( d(x,t,2) + ω0^2 * x ~ F *cos(ω*t), x);

# expand x in the harmonic ω
julia> add_harmonic!(diff_eom, x, ω);

# get equations for the harmonics evolving in the slow time T
julia> harmonic_eom = get_harmonic_equations(diff_eom)

A set of 2 harmonic equations
Variables: u1(T), v1(T)
Parameters: ω0, ω, F

Harmonic ansatz:
x(t) = u1*cos(ωt) + v1*sin(ωt)

Harmonic equations:

(ω0^2)*u1(T) + (2//1)*ω*Differential(T)(v1(T)) - (ω^2)*u1(T) ~ F

(ω0^2)*v1(T) - (ω^2)*v1(T) - (2//1)*ω*Differential(T)(u1(T)) ~ 0
```

"""
function get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothing, slow_time=nothing)

    slow_time = isnothing(slow_time) ? (@variables T; T) : slow_time
    fast_time = isnothing(fast_time) ? get_independent_variables(diff_eom)[1] : fast_time

    all(isempty.(values(diff_eom.harmonics))) && error("No harmonics specified!")
    eom = harmonic_ansatz(diff_eom, fast_time); # substitute trig functions into the differential equation
    eom = slow_flow(eom, fast_time=fast_time, slow_time=slow_time); # drop 2nd order time derivatives
    fourier_transform!(eom, fast_time); # perform averaging over the frequencies originally specified in dEOM
    ft_eom_simplified = drop_powers(eom, d(get_variables(eom), slow_time), 2); # drop higher powers of the first-order derivatives
    return ft_eom_simplified
end

"Rearrange `eq` to have zero on the right-hand-side."
_set_zero_rhs(eq::Equation) = eq.lhs - eq.rhs ~ 0
_set_zero_rhs(eqs::Vector{Equation}) = [_set_zero_rhs(eq) for eq in eqs]

"Returns the equation system in `eom`, dropping all argument brackets (i.e., u(T) becomes u)."
function _equations_without_brackets(eom::HarmonicEquation)
    variable_rules = [var => declare_variable(var_name(var)) for var in get_variables(eom)]
    equations_lhs = Num.(getfield.(eom.equations, :lhs)  - getfield.(eom.equations, :rhs))
    substitute_all(equations_lhs, variable_rules)
end

function drop_powers(eom::HarmonicEquation, terms, deg::Int)
    new_eom = deepcopy(eom)
    new_eom.equations = drop_powers(eom.equations, terms, deg)
    new_eom
end


text:
"""
    blockcat(n, A)
    blockcat(n, A::AbstractMatrix{SparseMatrixCSC{Ti,Tv}}) where {Ti,Tv}

Concatenates the n-by-n matrix of matrices A into a global matrix.
If eltype(A) = SparseMatrixCSC, blockcat returns a global sparse matrix.
Otherwise, blockcat returns a global dense array of type Array{eltype(first(A))}.
"""
function blockcat(n, A)
    A0 = first(A)
    nA = size(A0, 1)
    B  = Array{eltype(A0),2}(undef, n * nA, n * nA)
    for a in 1:n, j in 1:n
        B[(a-1)*nA .+ (1:nA), (j-1)*nA .+ (1:nA)] .= A[a,j]
    end
    return B
end

function blockcat(n, A::AbstractMatrix{SparseMatrixCSC{Ti,Tv}}) where {Ti,Tv}
    A0 = first(A)
    nA = size(A0, 1)
    B  = spzeros(n * nA, n * nA)
    for ii in 1:n, jj in 1:n
        for ijv in zip(findnz(A[ii,jj])...)
            i,j,Aij = ijv
            B[(ii-1)*nA + i, (jj-1)*nA + j] = Aij
        end
    end
    return B
end

#####
##### Jacobian functions
#####

function scale_factor(hadamard_product_type::Symbol)
    if hadamard_product_type == :skew
        return -1.0I
    else
        return 1.0I
    end
end

"""
    hadamard_jacobian(A_list, dF::Fxn, U, Fargs...;
                      skip_index=(i,j)->false) where {N,Fxn}
    hadamard_jacobian(A_list, hadamard_product_type, dF::Fxn, U,
                      Fargs...; skip_index=(i,j)->false) where {N,Fxn}

Computes Jacobian of the flux differencing term ∑_i sum(Ai.*Fi). Outputs array whose
entries are blocks of the Jacobian matrix corresponding to components of flux vectors.

Inputs:
- `A_list` = tuple of operators
- `hadamard_product_type` = `:skew`, `:sym`. Specifies if `Ai.*Fi` is skew or symmetric.
- `dF` = Jacobian of the flux function `F(uL,uR)` with respect to `uR`
- `U` = solution at which to evaluate the Jacobian
- `Fargs` = extra args for `dF(uL,uR)`
- (optional) `skip_index(i,j)` = optional function to skip computation of (i,j)th entry
"""
function hadamard_jacobian(A_list,dF::Fxn,U,Fargs...;skip_index=(i,j)->false) where {N,Fxn}

    # make symmetric and skew parts, call jacobian on each part
    Asym_list  = (A->.5*(A+A')).(A_list)
    Askew_list = (A->.5*(A-A')).(A_list)
    return ((x,y)->x .+ y).(
           hadamard_jacobian(Asym_list,:sym,dF,U,Fargs...;skip_index=skip_index),
           hadamard_jacobian(Askew_list,:skew,dF,U,Fargs...;skip_index=skip_index)
           )
end

function hadamard_jacobian(A_list,hadamard_product_type,dF::Fxn,U,Fargs...;
                           skip_index=(i,j)->false) where {Fxn}
    Nfields = length(U)
    A = SMatrix{Nfields,Nfields}([zero(first(A_list)) for i=1:Nfields, j=1:Nfields])
    hadamard_jacobian!(A, A_list, hadamard_product_type, dF, U, Fargs...; skip_index=skip_index)
    return A
end

"""
    hadamard_jacobian(A_list, dF::Fxn, U, ::Val{Nfields}, Fargs...;
                      skip_index=(i,j)->false) where {N,Nfields, Fxn}
    hadamard_jacobian(A_list, hadamard_product_type, dF::Fxn, U, ::Val{Nfields},
                      Fargs...; skip_index=(i,j)->false) where {N,Fxn}

Version with static information `Nfields` as a value-type argument.

# Example
```julia
julia> hadamard_jacobian(A_list, dF, U, Val(4)) # Nfields = 4
```
"""
function hadamard_jacobian(A_list,dF::Fxn,U,::Val{Nfields},Fargs...;
                           skip_index=(i,j)->false) where {N,Nfields,Fxn}

    # make symmetric and skew parts, call jacobian on each part
    Asym_list  = (A->.5*(A+A')).(A_list)
    Askew_list = (A->.5*(A-A')).(A_list)
    return ((x,y)->x .+ y).(
           hadamard_jacobian(Asym_list,:sym,dF,U,Val(Nfields),Fargs...;skip_index=skip_index),
           hadamard_jacobian(Askew_list,:skew,dF,U,Val(Nfields),Fargs...;skip_index=skip_index)
           )
end

function hadamard_jacobian(A_list,hadamard_product_type,dF::Fxn,U,::Val{Nfields},Fargs...;
                           skip_index=(i,j)->false) where {Nfields,Fxn}
    A = SMatrix{Nfields,Nfields}([zero(first(A_list)) for i=1:Nfields, j=1:Nfields])
    hadamard_jacobian!(A, A_list, hadamard_product_type, dF, U, Fargs...; skip_index=skip_index)
    return A
end

"""
    hadamard_jacobian!(A::SMatrix{N,N},A_list,hadamard_product_type::Symbol,
                       dF::Fxn,U,Fargs...; skip_index=(i,j)->false) where {N,Nd,Fxn}

Mutating version of [`hadamard_jacobian`](@ref). `A` = matrix for storing Jacobian
output, with each entry storing a block of the Jacobian.
"""
function hadamard_jacobian!(A::SMatrix{Nfields,Nfields},A_list,hadamard_product_type::Symbol,
                            dF::Fxn,U,Fargs...; skip_index=(i,j)->false) where {Fxn,Nfields}

    rows,cols = axes(first(A_list))

    # accumulator for sum(Q.*dF,1) over jth column
    Atype = eltype(first(A_list))
    dFaccum = MMatrix{Nfields,Nfields}(zeros(Atype,Nfields,Nfields))

    # loop over cols + non-zero ids
    for j in cols
       Uj = getindex.(U,j)

       fill!(dFaccum,zero(Atype))
       for i in rows
           if skip_index(i,j)==false
               Ui = getindex.(U,i)
               A_ij_list = getindex.(A_list,i,j)
               dFij = dF(Ui,Uj,getindex.(Fargs,i)...,getindex.(Fargs,j)...)

               @inbounds for n = 1:Nfields, m=1:Nfields
                   dFijQ = sum(map((x,y)->x.*y,getindex.(dFij,m,n),A_ij_list)) # sum result for multiple operators
                   A[m,n][i,j] += dFijQ
                   dFaccum[m,n] += dFijQ # accumulate column sums on-the-fly
               end
           end
       end

       # add diagonal entry for each block
       @inbounds for n=1:Nfields, m=1:Nfields
           A[m,n][j,j] += scale_factor(hadamard_product_type)*dFaccum[m,n]
       end
    end
end

function hadamard_jacobian!(A::SMatrix{N,N},A_list::NTuple{D,SparseMatrixCSC},
                            hadamard_product_type::Symbol,dF::Fxn,U,Fargs...;
                            skip_index=(i,j)->false) where {N,D,Fxn}
    for (i,Ai) in enumerate(A_list)
        dF_i = (x->getindex(x,i))∘dF
        hadamard_jacobian!(A,Ai,hadamard_product_type,dF_i,U,Fargs...)
    end
end

function hadamard_jacobian!(A::SMatrix{Nfields,Nfields},Amat::SparseMatrixCSC,
                            hadamard_product_type::Symbol,dF::Fxn,U,Fargs...) where {Nfields,Fxn}

    rows = rowvals(Amat)
    vals = nonzeros(Amat)

    # accumulator for sum(Q.*dF,1) over jth column
    Atype = eltype(Amat)
    dFaccum = MMatrix{Nfields,Nfields}(zeros(Atype,Nfields,Nfields))

    # loop over cols + non-zero ids
    for j = 1:size(Amat,2)
       Uj = getindex.(U,j)
       fill!(dFaccum,zero(Atype))
       for row_id in nzrange(Amat,j)
           i = rows[row_id]
           A_ij = vals[row_id]

           Ui = getindex.(U,i)
           dFij = dF(Ui,Uj,getindex.(Fargs,i)...,getindex.(Fargs,j)...)

           @inbounds for n = 1:Nfields, m=1:Nfields
               dFijQ = dFij[m,n].*A_ij # sum result for single op
               A[m,n][i,j] += dFijQ
               dFaccum[m,n] += dFijQ # accumulate column sums on-the-fly
           end
       end

       # add diagonal entry for each block
       @inbounds for n=1:Nfields, m=1:Nfields
           A[m,n][j,j] += scale_factor(hadamard_product_type)*dFaccum[m,n]
       end
    end
end

ext: jl


text:
using BinDeps
@BinDeps.setup

libbnet = library_dependency("libbnet")

libdir = BinDeps.libdir(libbnet)
srcdir = joinpath(BinDeps.srcdir(libbnet), "binary_networks")
provides(Sources, URI("https://raw.githubusercontent.com/afternone/CommunityDetection.jl/master/deps/binary_networks.tar.gz"), libbnet)
provides(BuildProcess,
    (@build_steps begin
        GetSources(libbnet)
        CreateDirectory(libdir)
        @build_steps begin
            ChangeDirectory(srcdir)
            FileRule(joinpath(libdir, "libbnet.so"), @build_steps begin
                `make`
                `cp libbnet.so $libdir`
            end)
        end
    end), libbnet)

@BinDeps.install Dict(:libbnet => :libbnet)


text:
using Squeezing
using Compat.Test

@testset "squeezing sum" begin
    @test @inferred @squeezing(sum(ones(3, 2), 1)) == [3.0, 3.0]
end

@testset "squeezing reducedim" begin
    A = reshape(Vector(1:16), (4,4))
    @test @inferred @squeezing(reducedim(max, A, 1)) == collect(4:4:16)
    @test @inferred @squeezing(reducedim(max, A, 2)) == collect(13:16)
end


text:
# The following was inspired by AtariAlgos.jl
# Modifications by Dennis Wilson @d9w and Erwan Lecarpentier @erwanlecarpentier

using ArcadeLearningEnvironment
using Colors
using ImageCore
using ImageTransformations
using Interpolations
export
    Game,
    close!,
    reset!,
    draw,
    get_inputs,
    get_state_buffer,
    get_observation_buffer,
    get_state!,
    get_observation!,
    get_state,
    get_screen,
    get_rgb,
    get_grayscale,
    get_state_ref,
    get_ram,
    save_screen_png

struct Game
    ale::ALEPtr
    width::Int
    height::Int
    actions::Array{Int32}
end

function Game(romfile::String, seed::Int64; kwargs...)
    ale = ALE_new()
    setInt(ale, "random_seed", Cint(seed))
    setInt(ale, "repeat_action_probability", Cint(0))
    # setBool(ale, "color_averaging", true)
    # setInt(ale, "frame_skip", Int32(1)) # 1 means no frame skip
    setFloat(ale, "repeat_action_probability", Float32(0.0))
    kwargs_dict = Dict(kwargs)
    if haskey(kwargs_dict, :lck)  # Thread safe
        lock(kwargs_dict[:lck]) do
            loadROM(ale, romfile)
        end
        # unlock(kwargs_dict[:lck])
    else
        loadROM(ale, romfile)
    end
    if haskey(kwargs_dict, :state_ref)  # Initial state
        ArcadeLearningEnvironment.restoreSystemState(ale, kwargs_dict[:state_ref])
    end
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    actions = getMinimalActionSet(ale)
    #=
    for _ in 1:3 # 1st actions are noop
        act(ale, 0)
    end
    =#
    Game(ale, w, h, actions)
end

function reset!(game::Game)
    reset_game(game.ale)
end

function get_state_ref(romfile::String, seed::Int64)
    ale = ALE_new()
    setInt(ale, "random_seed", Cint(seed))
    loadROM(ale, romfile)
    state_ref = ArcadeLearningEnvironment.cloneSystemState(ale)
    ALE_del(ale)
    state_ref
end

function close!(game::Game)
    ALE_del(game.ale)
end

function draw(game::Game)
    rawscreen = getScreenRGB(game.ale)
    colorview(RGB, Float64.(reshape(rawscreen/256.,
                                    (3, game.width, game.height))))';
end

function get_inputs(game::Game)
    screen = getScreen(game.ale)/(0xff*1.0)
    screen = reshape(screen, (game.width, game.height))'
    # imresize(screen, (42, 32))/256.
    screen
end

function scale(x::Array{UInt8,2}; factor::Int64=2)
    outsz = ceil.(Int, size(x)./factor)
    tilesz = (factor,factor)
    out = Array{UInt8,2}(undef, outsz)
    R = TileIterator(axes(x), tilesz)
    i = 1
    for tileaxs in R
       out[i] = maximum(view(x, tileaxs...))
       i += 1
    end
    out
end

function get_state_buffer(game::Game, grayscale::Bool)
    size = grayscale ? game.width * game.height : game.width * game.height * 3
    Vector{UInt8}(undef, size)
end

function get_observation_buffer(game::Game, grayscale::Bool, downscale::Bool)
    n_obs = grayscale ? 1 : 3
    if downscale
        # (l+1)÷2 for odd l, and l÷2 + 1 for even l
        w = isodd(game.width) ? (game.width + 1) ÷ 2 : game.width ÷ 2 + 1
        h = isodd(game.height) ? (game.height + 1) ÷ 2 : game.height ÷ 2 + 1
        obs_size = (w, h)
    else
        obs_size = (game.width, game.height)
    end
    [Array{UInt8, 2}(undef, obs_size) for _ in 1:n_obs]
end

function get_state!(
    s::Vector{UInt8},
    game::Game,
    grayscale::Bool
)
    if grayscale
        ArcadeLearningEnvironment.getScreenGrayscale!(game.ale, s)
    else
        ArcadeLearningEnvironment.getScreenRGB!(game.ale, s)
    end
end

"""
grayscale downscale - @btime
1 1 - 41.143 μs (16 allocations: 276.27 KiB)
1 0 - 77.752 ns (3 allocations: 192 bytes)
0 0 - 54.647 μs (7 allocations: 98.97 KiB)
0 1 - 158.778 μs (43 allocations: 828.58 KiB)
"""
function get_observation!(
    o::Vector{Matrix{UInt8}},
    s::Vector{UInt8},
    game::Game,
    grayscale::Bool,
    downscale::Bool
)
    if grayscale
        if downscale
            o .= [convert(Matrix{UInt8},
                floor.(imresize(reshape(s, (game.width, game.height)),
                ratio=0.5, method=BSpline(Linear()))))]
            #    convert(Matrix{UInt8},
            #    imresize(reshape(s, (game.width, game.height)),
            #    ratio=0.5, method=BSpline(Constant())))]
        else
            o .= [reshape(s, (game.width, game.height))]
        end
    else
        if downscale
            o .= [
                convert(Matrix{UInt8},
                    floor.(imresize(
                        reshape(
                            @view(s[i:3:length(s)]),
                            (game.width, game.height)
                        ),
                        ratio=0.5,
                        method=BSpline(Linear())
                    ))
                ) for i in 1:3
            ]
        else
            o .= [
                convert(
                    Array{UInt8,2},
                    reshape(
                        @view(s[i:3:length(s)]),
                        (game.width, game.height)
                    )
                ) for i in 1:3
            ]
        end
    end
end

function get_state(
    game::Game,
    grayscale::Bool,
    downscale::Bool
)
    s = grayscale ? get_grayscale(game) : get_rgb(game)
    if downscale
        @inbounds for i in eachindex(s)
            # s[i] = scale(s[i])
            #s[i] = convert(Matrix{UInt8}, imresize(s[i], ratio=0.5, method=BSpline(Constant())))
            s[i] = convert(Matrix{UInt8}, floor.(imresize(s[i], ratio=0.5, method=BSpline(Linear()))))
            #s[i] = convert(Matrix{UInt8}, floor.(restrict(s[i]))) # Memory issues ?
        end
    end
    s
end

function save_screen_png(game::Game, filename::String)
    saveScreenPNG(game.ale, filename)
end

function get_screen(game::Game)
    rawscreen = getScreen(game.ale)
    convert(Array{UInt8,2}, reshape(@view(rawscreen[1:length(rawscreen)]), (game.width, game.height)))
end

function get_rgb(game::Game)
    rawscreen = getScreenRGB(game.ale)
    [convert(Array{UInt8,2}, reshape(@view(rawscreen[i:3:length(rawscreen)]), (game.width, game.height))) for i in 1:3]
end

function get_grayscale(game::Game)
    gs = getScreenGrayscale(game.ale)
    gs = reshape(gs, (game.width, game.height))
    [gs]
end

function get_ram(game::Game)
    getRAM(game.ale) ./ typemax(UInt8)
end


text:
############################################################################
### Types
############################################################################

struct ImplyEmpty{V, V2} <: SemImply
    start_val::V
    identifier::V2
end

############################################################################
### Constructors
############################################################################

function ImplyEmpty(;
        specification = nothing,
        start_val = start_fabin3,
        kwargs...)

        if !isa(start_val, Vector)
            if specification isa RAMMatrices
                ram_matrices = specification
                identifier = StructuralEquationModels.identifier(ram_matrices)
            elseif specification isa ParameterTable
                ram_matrices = RAMMatrices!(specification)
                identifier = StructuralEquationModels.identifier(ram_matrices)
            else
                throw(ErrorException("The RAM constructor does not know how to handle your specification object.
                \n Please specify your model as either a ParameterTable or RAMMatrices."))
            end
            start_val = start_val(;ram_matrices = ram_matrices, specification = specification, kwargs...)
        end

        return ImplyEmpty(start_val, identifier)
end

############################################################################
### functors
############################################################################

function (imply::ImplyEmpty)(par, F, G, H, model) end

############################################################################
### Recommended methods
############################################################################

identifier(imply::ImplyEmpty) = imply.identifier
n_par(imply::ImplyEmpty) = imply.n_par

update_observed(imply::ImplyEmpty, observed::SemObs; kwargs...) = imply

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::ImplyEmpty)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
# ---
# title: 20. Valid Parentheses
# id: problem20
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String, Stack
# link: <https://leetcode.com/problems/valid-parentheses/description/>
# hidden: true
# ---
#
# Given a string `s` containing just the characters `'('`, `')'`, `'{'`, `'}'`,
# `'['` and `']'`, determine if the input string is valid.
#
# An input string is valid if:
#
#   1. Open brackets must be closed by the same type of brackets.
#   2. Open brackets must be closed in the correct order.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "()"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: s = "()[]{}"
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: s = "(]"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: s = "([)]"
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: s = "{[]}"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of parentheses only `'()[]{}'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    GaN Potential

See 10.1088/1361-648X/ab6cbe
"""
mutable struct GaN <: Potential
    lj_Ga_Ga::LennardJones
    lj_N_N::LennardJones
    bm_Ga_N::BornMayer
    c::Coulomb
    no_Ga::Int64
    no_N::Int64
end

function GaN(params::Dict)
    # Creates the GaN model
    lj_Ga_Ga = LennardJones(params["ε_Ga_Ga"], params["σ_Ga_Ga"])
    lj_N_N = LennardJones(params["ε_N_N"], params["σ_N_N"])
    bm_Ga_N = BornMayer(params["A_Ga_N"], params["ρ_Ga_N"])
    c = Coulomb(params["q_Ga"], params["q_N"], params["ε0"])
    gan = GaN(lj_Ga_Ga, lj_N_N, bm_Ga_N, c, params["no_Ga"], params["no_N"])
    return gan
end

"""
    potential_energy(p::GaN, r::Position, args...)

Calculates GaN potential energy.
"""
function potential_energy(p::GaN, r::Position, args...)
    i = args[1]; j = args[2]
    if i <= p.no_Ga && j <= p.no_N # Ga-Ga interaction
        return potential_energy(p.c, r) + potential_energy(p.lj_Ga_Ga, r)
    elseif i > p.no_Ga && j > p.no_N # N-N interaction
        return potential_energy(p.c, r) + potential_energy(p.lj_N_N, r)
    else # Ga-N or N-Ga interaction
        return potential_energy(p.c, r) + potential_energy(p.bm_Ga_N, r)
    end

end


text:
using Documenter, ReactiveMP

makedocs(
    modules  = [ ReactiveMP ],
    clean    = true,
    sitename = "ReactiveMP.jl",
    pages    = [
        "Home"            => "index.md",
        "Getting started" => "getting-started.md",
        "Extending"       => "extending.md",
        "Distributions"   => "distributions.md",
        "API"             => [
            "Node" => "api/node.md"
        ]
        # "TODO"         => "todo.md",
        # "Contributing" => "contributing.md",
        # "Utils"        => "utils.md"
    ],
    format   = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/biaslab/ReactiveMP.jl.git"
    )
end


text:
using RobotOS
@rosimport std_msgs.msg: Float64MultiArray
@rosimport geometry_msgs.msg: Point
rostypegen()
using std_msgs.msg       # subscribe to topic
using geometry_msgs.msg

#TODO: use proper data type(s) here
#       to receive message from ROS (and actually from Gazebo)
function callback(msg::Float64MultiArray)
   println(msg)
end
# https://github.com/jdlangs/RobotOS.jl

# TODO: here is an example of pushing data in fixed rate
#       (just like there is a timer)
function loop(pub_obj)
   loop_rate = Rate(1)
   npt = Float64MultiArray()
   npt.data = Array{Float64}(6)
   npt.data[6] = 0.0
   while ! is_shutdown()
     npt.data[1] += 0.2
     npt.data[2] = 0.0
     npt.data[3] = 0.0
     npt.data[4] = 0.0
     npt.data[5] = 0.0
     npt.data[6] += 0.2  # NOTE this might need to be set to 0

     # send information to ROS (then Gazebo)
     publish(pub_obj, npt)
     rossleep(loop_rate)
   end
end

function main()
   #init_node("rosjl_example") # what is this?

   pub = Publisher{Float64MultiArray}("gazebo_client/test", queue_size=10)
   sub = Subscriber{Float64MultiArray}("gazebo_client/sendscan", callback, (), queue_size=10)

   loop(pub)
end

if ! isinteractive()
   main()
end

# set up roscore
#@pyimport roslaunch

# setup for the Python script
#using PyCall
#@pyimport imp
#Node=imp.load_source("Node","pythonNode.py");

include("utils.jl") # functions

function __init__()
  addprocs(3);

  # start up a roscore ?? or is it already running
  #roslaunch.roscore

  # run the ptthon script Echo.py to start ROS node to cummunicate between julia and Gazebo
  node=remotecall(Node[:main],2,());
  println("ROS node Loaded!");

  # run the julia script handler.jl to get point cloud information from Gazebo
  data=remotecall(main,3,());
  println("Ready to Collect Data From Gazebo!");

  # start Gazebo
  gazebo=remotecall(run,4,(`gazebo velodyne.world`));
  println("Gazebo LiDAR model running!");

end

export data


text:
# Assess sensitivity of mescososm model to different alpha values,

# Define differential equation with first order kinetics
# and a power law residence time distribution
@sam_model("firstOrder", "powerLaw")

# Specify time span, background, initial condition, and parameters
tspan = (0.0, 1800) # model run time
h(p,t) = 310.0 # background specific conductance
u0 = [325.0] # concentration in channel after spike
τ₋ = 0.1
τ₊ = 900.0
k = 0.0
s = 27.0
sHZ = s/3
sCH = s - sHZ

# Solve model for two different alpha values
α = 1.5
qPerUnitHZ = qCalc_powerLaw(1, τ₋, τ₊, α)
p = [τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH]
prob2 = remake(prob, p = p)
sol1 = solve(prob2, MethodOfSteps(RK4()), saveat = 1)

α = 1.7
qPerUnitHZ = qCalc_powerLaw(1, τ₋, τ₊, α)
p = [τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH]
prob4 = remake(prob, p = p)
sol2 = solve(prob4, MethodOfSteps(RK4()), saveat = 1)

# Plot model output
plot(sol1, label = "alpha = 1.5", ylabel = "Nitrate Concentration (ppb N)", lw = 3, legend = :topleft)
plot!(sol2, label = "alpha = 1.7", xlabel = "Seconds", lw = 3)

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
module CIAOAlgorithms

const RealOrComplex{R} = Union{R,Complex{R}}
const Maybe{T} = Union{T, Nothing}

# utulities
include("utilities/indexingUtilities.jl")

# algorithms
include("algorithms/Finito/Finito.jl")
include("algorithms/ProShI/ProShI.jl")
include("algorithms/SVRG/SVRG.jl")
include("algorithms/SAGA_SAG/SAGA.jl")

end # module

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
using Documenter
using Quon
using DocThemeIndigo

indigo = DocThemeIndigo.install(Quon)

makedocs(;
    modules = [Quon],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical="https://QuantumBFS.github.io/Quon.jl",
        assets=String[indigo],
    ),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.com/QuantumBFS/Quon.jl",
    sitename = "Quon.jl",
)

deploydocs(; repo = "github.com/QuantumBFS/Quon.jl")


text:
using Flux
using Flux: @functor
using NamedTupleTools

@views function gate(x::AbstractArray{<:Real,3})
   split = size(x, 2) ÷ 2
   tanh.(x[:, 1:split, :]) .* σ.(x[:, (split+1):end, :])
end

struct ResBlock{D, R, S, C}
   dilconv  :: D
   resconv  :: R
   skipconv :: S
   condconv :: C
end

struct LastBlock{D, S, C}
   dilconv  :: D
   skipconv :: S
   condconv :: C
end

@functor ResBlock
@functor LastBlock

function ResBlock(dilation::Integer, nch::NamedTuple{(:res, :skip, :cond),<:NTuple{3,Integer}})
   dilconv  = Conv((2,), nch.res => 2nch.res; dilation = dilation, pad = (dilation, 0)) |> gpu
   resconv  = Conv((1,), nch.res => nch.res)  |> gpu
   skipconv = Conv((1,), nch.res => nch.skip) |> gpu
   condconv = Conv((1,), nch.cond => 2nch.res, tanh) |> gpu
   ResBlock(dilconv, resconv, skipconv, condconv)
end

function LastBlock(dilation::Integer, nch::NamedTuple{(:res, :skip, :cond),<:NTuple{3,Integer}})
   dilconv  = Conv((2,), nch.res => 2nch.res; dilation = dilation, pad = (dilation, 0)) |> gpu
   skipconv = Conv((1,), nch.res => nch.skip) |> gpu
   condconv = Conv((1,), nch.cond => 2nch.res, tanh) |> gpu
   LastBlock(dilconv, skipconv, condconv)
end

function Base.show(io::IO, block::ResBlock)
   dc, rc, sc, cc = block.dilconv, block.resconv, block.skipconv, block.condconv
   dilation = first(dc.dilation)
   _, res, skip = size(sc.weight)
   cond = size(cc.weight, 2)
   print(io,
      """
      ResBlock($dilation, (res=$res, skip=$skip, cond=$cond))
      │  Dilated     $dc
      │  Residual    $rc
      │  Skip        $sc
      └─ Conditional $cc
      """
   )
end

function Base.show(io::IO, block::LastBlock)
   dc, sc, cc = block.dilconv, block.skipconv, block.condconv
   dilation = first(dc.dilation)
   _, res, skip = size(sc.weight)
   cond = size(cc.weight, 2)
   print(io,
      """
      LastBlock($dilation, (res=$res, skip=$skip, cond=$cond))
      │  Dilated     $dc
      │  Skip        $sc
      └─ Conditional $cc
      """
   )
end

"""                             residual connection
        ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
        ┃                              ┏━━ tanh ━━┓                       ┃
input ━━┻━━> 2×1 dilatedconv ━━━ + ━━━ ❖          × ━━┳━━ 1×1 resconv ━━━ + ━━━> output
                                 ┃     ┗━━━ σ ━━━━┛   ┃
conditional input ━━━> 1×1 condconv             1×1 skipconv
                                                      ┃
Σskip ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━> + ━━━━━━━━> Σskip
"""
function (m::ResBlock)((Σskip, x)::NTuple{2,DenseArray{<:Real,3}})
   y = gate(m.dilconv(x))
   return (Σskip + m.skipconv(y), x + m.resconv(y))
end

function (m::ResBlock)((Σskip, x, h)::NTuple{3,DenseArray{<:Real,3}})
   y = gate(m.dilconv(x) + m.condconv(h))
   return (Σskip + m.skipconv(y), x + m.resconv(y), h)
end

function (m::LastBlock)((Σskip, x)::NTuple{2,DenseArray{<:Real,3}})
   y = gate(m.dilconv(x))
   return Σskip + m.skipconv(y)
end

function (m::LastBlock)((Σskip, x, h)::NTuple{3,DenseArray{<:Real,3}})
   y = gate(m.dilconv(x) + m.condconv(h))
   return Σskip + m.skipconv(y)
end

struct WaveNet{I, R, O}
   inblock   :: I
   resblocks :: R
   outblock  :: O
   depth     :: Int
end

@functor WaveNet

function WaveNet(nlayers::Integer = 30,
      nch::NamedTuple{(:res, :skip, :cond, :out),<:NTuple{4,Integer}} = (res=64, skip=256, cond=0, out=256); depth::Integer = 10
   )
   nch′ = delete(nch, :out)
   inblock = Conv((2,), 1 => nch.res, tanh; pad=(1,0)) |> gpu

   resblocks = Chain(
      (ResBlock(2^(i % depth), nch′) for i ∈ 0:(nlayers-2))...,
      LastBlock(2^((nlayers-1) % depth), nch′)
   ) |> gpu

   outblock = Chain(
      Conv((1,), nch.skip => nch.out, leakyrelu),
      Conv((1,), nch.out => nch.out),
      softmax
   ) |> gpu
   WaveNet(inblock, resblocks, outblock, depth)
end

function Base.show(io::IO, m::WaveNet)
   nlayers = length(m.resblocks)
   depth = m.depth
   res = length(m.inblock.bias)
   cond = size(first(m.resblocks).condconv.weight, 2)
   _, skip, out = size(first(m.outblock).weight)
   print(io, "WaveNet($nlayers, (res=$res, skip=$skip, cond=$cond, out=$out); depth=$depth)")
end

function (m::WaveNet)(x::DenseArray{<:Real,3})
   y = m.inblock(x)
   len, _, batch = size(y)
   skip₀ = gpu(zeros(Float32, len, nch.skip, batch))
   Σskip = m.resblocks((skip₀, y))
   return m.outblock(leakyrelu.(Σskip))
end

function (m::WaveNet)(x::DenseArray{<:Real,3}, h::DenseArray{<:Real,3})
   y = m.inblock(x)
   len, _, batch = size(y)
   skip₀ = zeros(Float32, len, nch.skip, batch)
   Σskip = m.resblocks((skip₀, y, h))
   return m.outblock(leakyrelu.(Σskip))
end


text:
function ScatterMatrixPlot(olddf::DataFrame;colorido=[],filepath::AbstractString="scattermatrix",mime::AbstractString="svg",xwidth=0cm,ywidth=0cm,legenda::Bool=false)
  pl1=ScatterMatrix1(olddf, colorido, legenda)
  pl2=ScatterMatrix2(olddf, colorido, legenda)

  filepath1=string(filepath,"_1")
  filepath2=string(filepath,"_2")

   if(xwidth==0cm)
	xwidth=((size(olddf)[2]+1)*5)cm
   end
   if(ywidth==0cm)
	ywidth=((size(olddf)[2]+1)*5)cm
   end

    if( (mime=="svg"||mime=="SVG")==true)
        filepath1=string(filepath1,".svg")
        filepath2=string(filepath2,".svg")

        image1=draw(SVG(filepath1, xwidth, ywidth), pl1)
        image2=draw(SVG(filepath2, xwidth, ywidth), pl2)
    elseif ( (mime=="png"||mime=="PNG")==true)
        filepath1=string(filepath1,".png")
        filepath2=string(filepath2,".png")

        image1=draw(PNG(filepath1, xwidth, ywidth), pl1)
        image2=draw(PNG(filepath2, xwidth, ywidth), pl2)
    elseif( (mime=="js.svg"||mime=="JS.SVG")==true)
        filepath1=string(filepath1,".js.svg")
        filepath2=string(filepath2,".js.svg")

        image1=draw(SVG(filepath1, xwidth, ywidth), pl1)
        image2=draw(SVG(filepath2, xwidth, ywidth), pl2)
    else
        filepath1=string(filepath1,".svg")
        filepath2=string(filepath2,".svg")

        image1=draw(SVG(filepath1, xwidth, ywidth), pl1)
        image2=draw(SVG(filepath2, xwidth, ywidth), pl2)
    end

    open_imagefile(filepath1)
    open_imagefile(filepath2)

end

max_stars_repo_name: wakakusa/ScatterMatrixPlots

max_issues_repo_name: wakakusa/ScatterMatrixPlots

max_forks_repo_name: wakakusa/ScatterMatrixPlots


text:
#julia4 --machinefile=$PBS_NODEFILE -- run100.jl 1 2 110
#run on all available cores using:
#addprocs(readcsv(ENV["PBS_NODEFILE"], ASCIIString)[:, 1])
@everywhere args = @fetchfrom 1 ARGS
import LatentGaussianMixtureModel
using DelimitedFiles
#include(joinpath(Pkg.dir("LatentGaussianMixtureModel"), "biccompare/runbic.jl"))
include(joinpath(dirname(pathof(LatentGaussianMixtureModel)), "..", "biccompare/runbic.jl"))
@everywhere C=parse(Int,args[2])
teststat = pmap((b) -> Brun(b, C), 1:parse(Int, args[1]))
writedlm("bic$(C).csv", teststat, ',')
run(`sed -i "s/\[//g" bic$(C).csv`)
run(`sed -i "s/\]//g" bic$(C).csv`)


text:
"""
Boundary bbox

- st_bbox(xmin, ymin, xmax, ymax)
- st_bbox(lon, lat)
- st_bbox(ga::GeoArray)
- st_bbox!(ga::GeoArray, b::bbox) : reset bbox
- st_bbox(ncfile::String)
"""

function st_bbox(lon, lat)
    cellsize_x = abs(lon[2] - lon[1])
    cellsize_y = abs(lat[2] - lat[1])
    bbox(minimum(lon) - cellsize_x/2, minimum(lat) - cellsize_y/2,
        maximum(lon) + cellsize_x/2, maximum(lat) + cellsize_y/2)
end

function st_bbox(ga::GeoArray)
    ax, ay = ga.f(SVector(0,0))
    bx, by = ga.f(SVector(size(ga)[1:2]))
    # (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by))
    bbox(min(ax, bx), min(ay, by), max(ax, bx), max(ay, by))
end

function st_bbox(ncfile::String)
    lat = NetCDF.ncread(ncfile, "lat")
    lon = NetCDF.ncread(ncfile, "lon")
    st_bbox(lon, lat)
end

function bbox_to_affine(size::Tuple{Integer, Integer}, b::bbox)
    AffineMap(
        SMatrix{2,2}((b.xmax - b.xmin) / size[1], 0, 0, -(b.ymax - b.ymin)/size[2]),
        SVector(b.xmin, b.ymax)
        )
end

"""Set geotransform of `GeoArray` by specifying a bounding box.
Note that this only can result in a non-rotated or skewed `GeoArray`."""
function st_bbox!(ga::GeoArray, b::bbox)
    ga.f = bbox_to_affine(size(ga)[1:2], b)
    ga
end

export bbox, st_bbox, st_bbox!


text:
function ms()

    global MS_START
    if ! @isdefined MS_START
        MS_START = time()
    end
    ms = round(Int, 1000*(time() - MS_START))
    return (ms/1000)

end

function ms_reset()

    global MS_START = time()

end

macro ms(cmd)

    cmdstr = string(cmd)
    # information about the file that called this macro
    if isinteractive()
        # interactive mode => no file to call
        fileinfo = "(mode repl)"
    else
        # retrieve the relative name of the file that called this macro
        fname = string(__source__.file)
        fname = fname[length("$APPDIR")+2 : end]
        # retrieve the line number where the macro is called
        fline = string(__source__.line)
        fileinfo = "$(fname):$fline"
    end

    quote
        local t0 = time()
        print("@ms START ", round(ms(), digits=3), "s ", $fileinfo, ":", $(cmdstr))
        println(" ... ")
        local val = $(esc(cmd))
        local t1 = time()
        print("@ms END ", round(ms(), digits=3), "s ", $fileinfo, ":", $(cmdstr))
        println(" in ", round(t1-t0, digits=3), "s")
        val
    end

end

# start the chronometer
ms()

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
# This file is programmatically generated, do not edit it!
#   Generated from `base/utf8ienc.dtx` from the LaTeX 2e source tree.
#   Available under the LaTeX Project Public License (LPPL).
#   File commit hash: e6ffcaf0
const u2latex = Dict(
    ' ' => "\\nobreakspace",
    '¡' => "\\textexclamdown",
    '¢' => "\\textcent",
    '£' => "\\textsterling",
    '¤' => "\\textcurrency",
    '¥' => "\\textyen",
    '¦' => "\\textbrokenbar",
    '§' => "\\textsection",
    '¨' => "\\textasciidieresis",
    '©' => "\\textcopyright",
    'ª' => "\\textordfeminine",
    '«' => "\\guillemotleft",
    '¬' => "\\textlnot",
    '\uad' => "\\-",
    '®' => "\\textregistered",
    '¯' => "\\textasciimacron",
    '°' => "\\textdegree",
    '±' => "\\textpm",
    '²' => "\\texttwosuperior",
    '³' => "\\textthreesuperior",
    '´' => "\\textasciiacute",
    'µ' => "\\textmu",
    '¶' => "\\textparagraph",
    '·' => "\\textperiodcentered",
    '¸' => "\\c\\ ",
    '¹' => "\\textonesuperior",
    'º' => "\\textordmasculine",
    '»' => "\\guillemotright",
    '¼' => "\\textonequarter",
    '½' => "\\textonehalf",
    '¾' => "\\textthreequarters",
    '¿' => "\\textquestiondown",
    'À' => "\\`A",
    'Á' => "\\'A",
    'Â' => "\\^A",
    'Ã' => "\\~A",
    'Ä' => "\\\"A",
    'Å' => "\\r A",
    'Æ' => "\\AE",
    'Ç' => "\\c C",
    'È' => "\\`E",
    'É' => "\\'E",
    'Ê' => "\\^E",
    'Ë' => "\\\"E",
    'Ì' => "\\`I",
    'Í' => "\\'I",
    'Î' => "\\^I",
    'Ï' => "\\\"I",
    'Ð' => "\\DH",
    'Ñ' => "\\~N",
    'Ò' => "\\`O",
    'Ó' => "\\'O",
    'Ô' => "\\^O",
    'Õ' => "\\~O",
    'Ö' => "\\\"O",
    '×' => "\\texttimes",
    'Ø' => "\\O",
    'Ù' => "\\`U",
    'Ú' => "\\'U",
    'Û' => "\\^U",
    'Ü' => "\\\"U",
    'Ý' => "\\'Y",
    'Þ' => "\\TH",
    'ß' => "\\ss",
    'à' => "\\`a",
    'á' => "\\'a",
    'â' => "\\^a",
    'ã' => "\\~a",
    'ä' => "\\\"a",
    'å' => "\\r a",
    'æ' => "\\ae",
    'ç' => "\\c c",
    'è' => "\\`e",
    'é' => "\\'e",
    'ê' => "\\^e",
    'ë' => "\\\"e",
    'ì' => "\\`\\i",
    'í' => "\\'\\i",
    'î' => "\\^\\i",
    'ï' => "\\\"\\i",
    'ð' => "\\dh",
    'ñ' => "\\~n",
    'ò' => "\\`o",
    'ó' => "\\'o",
    'ô' => "\\^o",
    'õ' => "\\~o",
    'ö' => "\\\"o",
    '÷' => "\\textdiv",
    'ø' => "\\o",
    'ù' => "\\`u",
    'ú' => "\\'u",
    'û' => "\\^u",
    'ü' => "\\\"u",
    'ý' => "\\'y",
    'þ' => "\\th",
    'ÿ' => "\\\"y",
    'Ā' => "\\=A",
    'ā' => "\\=a",
    'Ă' => "\\u A",
    'ă' => "\\u a",
    'Ą' => "\\k A",
    'ą' => "\\k a",
    'Ć' => "\\'C",
    'ć' => "\\'c",
    'Ĉ' => "\\^C",
    'ĉ' => "\\^c",
    'Ċ' => "\\.C",
    'ċ' => "\\.c",
    'Č' => "\\v C",
    'č' => "\\v c",
    'Ď' => "\\v D",
    'ď' => "\\v d",
    'Đ' => "\\DJ",
    'đ' => "\\dj",
    'Ē' => "\\=E",
    'ē' => "\\=e",
    'Ĕ' => "\\u E",
    'ĕ' => "\\u e",
    'Ė' => "\\.E",
    'ė' => "\\.e",
    'Ę' => "\\k E",
    'ę' => "\\k e",
    'Ě' => "\\v E",
    'ě' => "\\v e",
    'Ĝ' => "\\^G",
    'ĝ' => "\\^g",
    'Ğ' => "\\u G",
    'ğ' => "\\u g",
    'Ġ' => "\\.G",
    'ġ' => "\\.g",
    'Ģ' => "\\c G",
    'ģ' => "\\c g",
    'Ĥ' => "\\^H",
    'ĥ' => "\\^h",
    'Ĩ' => "\\~I",
    'ĩ' => "\\~\\i",
    'Ī' => "\\=I",
    'ī' => "\\=\\i",
    'Ĭ' => "\\u I",
    'ĭ' => "\\u\\i",
    'Į' => "\\k I",
    'į' => "\\k\\i",
    'İ' => "\\.I",
    'ı' => "\\i",
    'Ĳ' => "\\IJ",
    'ĳ' => "\\ij",
    'Ĵ' => "\\^J",
    'ĵ' => "\\^\\j",
    'Ķ' => "\\c K",
    'ķ' => "\\c k",
    'Ĺ' => "\\'L",
    'ĺ' => "\\'l",
    'Ļ' => "\\c L",
    'ļ' => "\\c l",
    'Ľ' => "\\v L",
    'ľ' => "\\v l",
    'Ł' => "\\L",
    'ł' => "\\l",
    'Ń' => "\\'N",
    'ń' => "\\'n",
    'Ņ' => "\\c N",
    'ņ' => "\\c n",
    'Ň' => "\\v N",
    'ň' => "\\v n",
    'Ŋ' => "\\NG",
    'ŋ' => "\\ng",
    'Ō' => "\\=O",
    'ō' => "\\=o",
    'Ŏ' => "\\u O",
    'ŏ' => "\\u o",
    'Ő' => "\\H O",
    'ő' => "\\H o",
    'Œ' => "\\OE",
    'œ' => "\\oe",
    'Ŕ' => "\\'R",
    'ŕ' => "\\'r",
    'Ŗ' => "\\c R",
    'ŗ' => "\\c r",
    'Ř' => "\\v R",
    'ř' => "\\v r",
    'Ś' => "\\'S",
    'ś' => "\\'s",
    'Ŝ' => "\\^S",
    'ŝ' => "\\^s",
    'Ş' => "\\c S",
    'ş' => "\\c s",
    'Š' => "\\v S",
    'š' => "\\v s",
    'Ţ' => "\\c T",
    'ţ' => "\\c t",
    'Ť' => "\\v T",
    'ť' => "\\v t",
    'Ũ' => "\\~U",
    'ũ' => "\\~u",
    'Ū' => "\\=U",
    'ū' => "\\=u",
    'Ŭ' => "\\u U",
    'ŭ' => "\\u u",
    'Ů' => "\\r U",
    'ů' => "\\r u",
    'Ű' => "\\H U",
    'ű' => "\\H u",
    'Ų' => "\\k U",
    'ų' => "\\k u",
    'Ŵ' => "\\^W",
    'ŵ' => "\\^w",
    'Ŷ' => "\\^Y",
    'ŷ' => "\\^y",
    'Ÿ' => "\\\"Y",
    'Ź' => "\\'Z",
    'ź' => "\\'z",
    'Ż' => "\\.Z",
    'ż' => "\\.z",
    'Ž' => "\\v Z",
    'ž' => "\\v z",
    'ƒ' => "\\textflorin",
    'Ǎ' => "\\v A",
    'ǎ' => "\\v a",
    'Ǐ' => "\\v I",
    'ǐ' => "\\v \\i",
    'Ǒ' => "\\v O",
    'ǒ' => "\\v o",
    'Ǔ' => "\\v U",
    'ǔ' => "\\v u",
    'Ǣ' => "\\=\\AE",
    'ǣ' => "\\=\\ae",
    'Ǧ' => "\\v G",
    'ǧ' => "\\v g",
    'Ǩ' => "\\v K",
    'ǩ' => "\\v k",
    'Ǫ' => "\\k O",
    'ǫ' => "\\k o",
    'ǰ' => "\\v\\j",
    'Ǵ' => "\\'G",
    'ǵ' => "\\'g",
    'Ș' => "\\textcommabelow S",
    'ș' => "\\textcommabelow s",
    'Ț' => "\\textcommabelow T",
    'ț' => "\\textcommabelow t",
    'Ȳ' => "\\=Y",
    'ȳ' => "\\=y",
    'ȷ' => "\\j",
    'ˆ' => "\\textasciicircum",
    'ˇ' => "\\textasciicaron",
    '˜' => "\\textasciitilde",
    '˘' => "\\textasciibreve",
    '˝' => "\\textacutedbl",
    'Ѐ' => "\\`\\CYRE",
    'Ё' => "\\CYRYO",
    'Ђ' => "\\CYRDJE",
    'Ѓ' => "\\'\\CYRG",
    'Є' => "\\CYRIE",
    'Ѕ' => "\\CYRDZE",
    'І' => "\\CYRII",
    'Ї' => "\\CYRYI",
    'Ј' => "\\CYRJE",
    'Љ' => "\\CYRLJE",
    'Њ' => "\\CYRNJE",
    'Ћ' => "\\CYRTSHE",
    'Ќ' => "\\'\\CYRK",
    'Ѝ' => "\\`\\CYRI",
    'Ў' => "\\CYRUSHRT",
    'Џ' => "\\CYRDZHE",
    'А' => "\\CYRA",
    'Б' => "\\CYRB",
    'В' => "\\CYRV",
    'Г' => "\\CYRG",
    'Д' => "\\CYRD",
    'Е' => "\\CYRE",
    'Ж' => "\\CYRZH",
    'З' => "\\CYRZ",
    'И' => "\\CYRI",
    'Й' => "\\CYRISHRT",
    'К' => "\\CYRK",
    'Л' => "\\CYRL",
    'М' => "\\CYRM",
    'Н' => "\\CYRN",
    'О' => "\\CYRO",
    'П' => "\\CYRP",
    'Р' => "\\CYRR",
    'С' => "\\CYRS",
    'Т' => "\\CYRT",
    'У' => "\\CYRU",
    'Ф' => "\\CYRF",
    'Х' => "\\CYRH",
    'Ц' => "\\CYRC",
    'Ч' => "\\CYRCH",
    'Ш' => "\\CYRSH",
    'Щ' => "\\CYRSHCH",
    'Ъ' => "\\CYRHRDSN",
    'Ы' => "\\CYRERY",
    'Ь' => "\\CYRSFTSN",
    'Э' => "\\CYREREV",
    'Ю' => "\\CYRYU",
    'Я' => "\\CYRYA",
    'а' => "\\cyra",
    'б' => "\\cyrb",
    'в' => "\\cyrv",
    'г' => "\\cyrg",
    'д' => "\\cyrd",
    'е' => "\\cyre",
    'ж' => "\\cyrzh",
    'з' => "\\cyrz",
    'и' => "\\cyri",
    'й' => "\\cyrishrt",
    'к' => "\\cyrk",
    'л' => "\\cyrl",
    'м' => "\\cyrm",
    'н' => "\\cyrn",
    'о' => "\\cyro",
    'п' => "\\cyrp",
    'р' => "\\cyrr",
    'с' => "\\cyrs",
    'т' => "\\cyrt",
    'у' => "\\cyru",
    'ф' => "\\cyrf",
    'х' => "\\cyrh",
    'ц' => "\\cyrc",
    'ч' => "\\cyrch",
    'ш' => "\\cyrsh",
    'щ' => "\\cyrshch",
    'ъ' => "\\cyrhrdsn",
    'ы' => "\\cyrery",
    'ь' => "\\cyrsftsn",
    'э' => "\\cyrerev",
    'ю' => "\\cyryu",
    'я' => "\\cyrya",
    'ѐ' => "\\`\\cyre",
    'ё' => "\\cyryo",
    'ђ' => "\\cyrdje",
    'ѓ' => "\\'\\cyrg",
    'є' => "\\cyrie",
    'ѕ' => "\\cyrdze",
    'і' => "\\cyrii",
    'ї' => "\\cyryi",
    'ј' => "\\cyrje",
    'љ' => "\\cyrlje",
    'њ' => "\\cyrnje",
    'ћ' => "\\cyrtshe",
    'ќ' => "\\'\\cyrk",
    'ѝ' => "\\`\\cyri",
    'ў' => "\\cyrushrt",
    'џ' => "\\cyrdzhe",
    'Ѣ' => "\\CYRYAT",
    'ѣ' => "\\cyryat",
    'Ѫ' => "\\CYRBYUS",
    'ѫ' => "\\cyrbyus",
    'Ѳ' => "\\CYRFITA",
    'ѳ' => "\\cyrfita",
    'Ѵ' => "\\CYRIZH",
    'ѵ' => "\\cyrizh",
    'Ѷ' => "\\C\\CYRIZH",
    'ѷ' => "\\C\\cyrizh",
    'Ҍ' => "\\CYRSEMISFTSN",
    'ҍ' => "\\cyrsemisftsn",
    'Ҏ' => "\\CYRRTICK",
    'ҏ' => "\\cyrrtick",
    'Ґ' => "\\CYRGUP",
    'ґ' => "\\cyrgup",
    'Ғ' => "\\CYRGHCRS",
    'ғ' => "\\cyrghcrs",
    'Ҕ' => "\\CYRGHK",
    'ҕ' => "\\cyrghk",
    'Җ' => "\\CYRZHDSC",
    'җ' => "\\cyrzhdsc",
    'Ҙ' => "\\CYRZDSC",
    'ҙ' => "\\cyrzdsc",
    'Қ' => "\\CYRKDSC",
    'қ' => "\\cyrkdsc",
    'Ҝ' => "\\CYRKVCRS",
    'ҝ' => "\\cyrkvcrs",
    'Ҟ' => "\\CYRKHCRS",
    'ҟ' => "\\cyrkhcrs",
    'Ҡ' => "\\CYRKBEAK",
    'ҡ' => "\\cyrkbeak",
    'Ң' => "\\CYRNDSC",
    'ң' => "\\cyrndsc",
    'Ҥ' => "\\CYRNG",
    'ҥ' => "\\cyrng",
    'Ҧ' => "\\CYRPHK",
    'ҧ' => "\\cyrphk",
    'Ҩ' => "\\CYRABHHA",
    'ҩ' => "\\cyrabhha",
    'Ҫ' => "\\CYRSDSC",
    'ҫ' => "\\cyrsdsc",
    'Ҭ' => "\\CYRTDSC",
    'ҭ' => "\\cyrtdsc",
    'Ү' => "\\CYRY",
    'ү' => "\\cyry",
    'Ұ' => "\\CYRYHCRS",
    'ұ' => "\\cyryhcrs",
    'Ҳ' => "\\CYRHDSC",
    'ҳ' => "\\cyrhdsc",
    'Ҵ' => "\\CYRTETSE",
    'ҵ' => "\\cyrtetse",
    'Ҷ' => "\\CYRCHRDSC",
    'ҷ' => "\\cyrchrdsc",
    'Ҹ' => "\\CYRCHVCRS",
    'ҹ' => "\\cyrchvcrs",
    'Һ' => "\\CYRSHHA",
    'һ' => "\\cyrshha",
    'Ҽ' => "\\CYRABHCH",
    'ҽ' => "\\cyrabhch",
    'Ҿ' => "\\CYRABHCHDSC",
    'ҿ' => "\\cyrabhchdsc",
    'Ӏ' => "\\CYRpalochka",
    'Ӂ' => "\\U\\CYRZH",
    'ӂ' => "\\U\\cyrzh",
    'Ӄ' => "\\CYRKHK",
    'ӄ' => "\\cyrkhk",
    'Ӆ' => "\\CYRLDSC",
    'ӆ' => "\\cyrldsc",
    'Ӈ' => "\\CYRNHK",
    'ӈ' => "\\cyrnhk",
    'Ӌ' => "\\CYRCHLDSC",
    'ӌ' => "\\cyrchldsc",
    'Ӎ' => "\\CYRMDSC",
    'ӎ' => "\\cyrmdsc",
    'Ӑ' => "\\U\\CYRA",
    'ӑ' => "\\U\\cyra",
    'Ӓ' => "\\\"\\CYRA",
    'ӓ' => "\\\"\\cyra",
    'Ӕ' => "\\CYRAE",
    'ӕ' => "\\cyrae",
    'Ӗ' => "\\U\\CYRE",
    'ӗ' => "\\U\\cyre",
    'Ә' => "\\CYRSCHWA",
    'ә' => "\\cyrschwa",
    'Ӛ' => "\\\"\\CYRSCHWA",
    'ӛ' => "\\\"\\cyrschwa",
    'Ӝ' => "\\\"\\CYRZH",
    'ӝ' => "\\\"\\cyrzh",
    'Ӟ' => "\\\"\\CYRZ",
    'ӟ' => "\\\"\\cyrz",
    'Ӡ' => "\\CYRABHDZE",
    'ӡ' => "\\cyrabhdze",
    'Ӣ' => "\\=\\CYRI",
    'ӣ' => "\\=\\cyri",
    'Ӥ' => "\\\"\\CYRI",
    'ӥ' => "\\\"\\cyri",
    'Ӧ' => "\\\"\\CYRO",
    'ӧ' => "\\\"\\cyro",
    'Ө' => "\\CYROTLD",
    'ө' => "\\cyrotld",
    'Ӭ' => "\\\"\\CYREREV",
    'ӭ' => "\\\"\\cyrerev",
    'Ӯ' => "\\=\\CYRU",
    'ӯ' => "\\=\\cyru",
    'Ӱ' => "\\\"\\CYRU",
    'ӱ' => "\\\"\\cyru",
    'Ӳ' => "\\H\\CYRU",
    'ӳ' => "\\H\\cyru",
    'Ӵ' => "\\\"\\CYRCH",
    'ӵ' => "\\\"\\cyrch",
    'Ӷ' => "\\CYRGDSC",
    'ӷ' => "\\cyrgdsc",
    'Ӹ' => "\\\"\\CYRERY",
    'ӹ' => "\\\"\\cyrery",
    'Ӻ' => "\\CYRGDSCHCRS",
    'ӻ' => "\\cyrgdschcrs",
    'Ӽ' => "\\CYRHHK",
    'ӽ' => "\\cyrhhk",
    'Ӿ' => "\\CYRHHCRS",
    'ӿ' => "\\cyrhhcrs",
    '฿' => "\\textbaht",
    'Ḃ' => "\\.B",
    'ḃ' => "\\.b",
    'ẞ' => "\\SS",
    '\u200c' => "\\textcompwordmark",
    '‐' => "-",
    '‒' => "\\textendash",
    '–' => "\\textendash",
    '—' => "\\textemdash",
    '―' => "\\textemdash",
    '‖' => "\\textbardbl",
    '‘' => "\\textquoteleft",
    '’' => "\\textquoteright",
    '‚' => "\\quotesinglbase",
    '“' => "\\textquotedblleft",
    '”' => "\\textquotedblright",
    '„' => "\\quotedblbase",
    '†' => "\\textdagger",
    '‡' => "\\textdaggerdbl",
    '•' => "\\textbullet",
    '…' => "\\textellipsis",
    '‰' => "\\textperthousand",
    '‱' => "\\textpertenthousand",
    '‹' => "\\guilsinglleft",
    '›' => "\\guilsinglright",
    '※' => "\\textreferencemark",
    '‽' => "\\textinterrobang",
    '⁄' => "\\textfractionsolidus",
    '⁎' => "\\textasteriskcentered",
    '⁒' => "\\textdiscount",
    '₡' => "\\textcolonmonetary",
    '₤' => "\\textlira",
    '₦' => "\\textnaira",
    '₩' => "\\textwon",
    '₫' => "\\textdong",
    '€' => "\\texteuro",
    '₱' => "\\textpeso",
    '℃' => "\\textcelsius",
    '№' => "\\textnumero",
    '℗' => "\\textcircledP",
    '℞' => "\\textrecipe",
    '℠' => "\\textservicemark",
    '™' => "\\texttrademark",
    'Ω' => "\\textohm",
    '℧' => "\\textmho",
    '℮' => "\\textestimated",
    '←' => "\\textleftarrow",
    '↑' => "\\textuparrow",
    '→' => "\\textrightarrow",
    '↓' => "\\textdownarrow",
    '〈' => "\\textlangle",
    '〉' => "\\textrangle",
    '␢' => "\\textblank",
    '␣' => "\\textvisiblespace",
    '◦' => "\\textopenbullet",
    '◯' => "\\textbigcircle",
    '♪' => "\\textmusicalnote",
    'Ḡ' => "\\=G",
    'ḡ' => "\\=g",
    '\ufeff' => "\\ifhmode\\nobreak\\fi",
    '£' => "\\textsterling",
    '§' => "\\textsection",
    '©' => "\\textcopyright",
    '±' => "\\textpm",
    '¶' => "\\textparagraph",
    '†' => "\\textdagger",
    '‡' => "\\textdaggerdbl",
    '…' => "\\textellipsis",
)


text:
abstract type BuildTargets end
struct JuliaTarget <: BuildTargets end
struct StanTarget <: BuildTargets end
struct CTarget <: BuildTargets end
struct MATLABTarget <: BuildTargets end

abstract type ParallelForm end
struct SerialForm <: ParallelForm end
struct MultithreadedForm <: ParallelForm end
struct DistributedForm <: ParallelForm end
struct DaggerForm <: ParallelForm end

"""
`build_function`

Generates a numerically-usable function from a ModelingToolkit `Expression`.
If the `Expression` is an `Operation`, the generated function is a function
with a scalar output, otherwise if it's an `AbstractArray{Operation}`, the output
is two functions, one for out-of-place AbstractArray output and a second which
is a mutating function. The outputted functions match the given argument order,
i.e., f(u,p,args...) for the out-of-place and scalar functions and
`f!(du,u,p,args..)` for the in-place version.

```julia
build_function(ex, args...;
               conv = simplified_expr, expression = Val{true},
               checkbounds = false, convert_oop = true,
			   force_SA = false,
               linenumbers = false, target = JuliaTarget())
```

Arguments:

- `ex`: The `Expression` to compile
- `vs`: The variables of the expression
- `ps`: The parameters of the expression
- `args`: Extra arguments to the function
- `conv`: The conversion function of the Operation to Expr. By default this uses
  the `simplified_expr` function utilized in `convert(Expr,x)`.
- `expression`: Whether to generate code or whether to generate the compiled form.
  By default, `expression = Val{true}`, which means that the code for the
  function is returned. If `Val{false}`, then the returned value is a compiled
  Julia function, which utilizes GeneralizedGenerated.jl in order to world-age
  free.

Keyword Arguments:

- `checkbounds`: For whether to enable bounds checking inside of the generated
  function. Defaults to false, meaning that `@inbounds` is applied.
- `linenumbers`: Determines whether the generated function expression retains
  the line numbers. Defaults to true.
- `convert_oop`: Determines whether the OOP version should try to convert
  the output to match the type of the first input. This is useful for
  cases like LabelledArrays or other array types that carry extra
  information. Defaults to true.
- `force_SA`: Forces the output of the OOP version to be a StaticArray.
  Defaults to `false`, and outputs a static array when the first argument
  is a static array.
- `target`: The output target of the compilation process. Possible options are:
    - `JuliaTarget`: Generates a Julia function
    - `CTarget`: Generates a C function
    - `StanTarget`: Generates a function for compiling with the Stan probabilistic
      programming language
    - `MATLABTarget`: Generates an anonymous function for use in MATLAB and Octave
      environments
"""
function build_function(args...;target = JuliaTarget(),kwargs...)
  _build_function(target,args...;kwargs...)
end

function addheader(ex, fargs, iip; X=gensym(:MTIIPVar))
  if iip
    wrappedex = :(
        ($X,$(fargs.args...)) -> begin
        $ex
        nothing
      end
    )
  else
    wrappedex = :(
      ($(fargs.args...),) -> begin
        $ex
      end
    )
  end
  wrappedex
end

function add_integrator_header(ex, fargs, iip; X=gensym(:MTIIPVar))
  integrator = gensym(:MTKIntegrator)
  if iip
    wrappedex = :(
        $integrator -> begin
        ($X,$(fargs.args...)) = (($integrator).u,($integrator).u,($integrator).p,($integrator).t)
        $ex
        nothing
      end
    )
  else
    wrappedex = :(
      $integrator -> begin
      ($(fargs.args...),) = (($integrator).u,($integrator).p,($integrator).t)
        $ex
      end
    )
  end
  wrappedex
end

# Scalar output
function _build_function(target::JuliaTarget, op::Operation, args...;
                         conv = simplified_expr, expression = Val{true},
                         checkbounds = false,
                         linenumbers = true, headerfun=addheader)

    argnames = [gensym(:MTKArg) for i in 1:length(args)]
    arg_pairs = map(vars_to_pairs,zip(argnames,args))
    ls = reduce(vcat,first.(arg_pairs))
    rs = reduce(vcat,last.(arg_pairs))
    var_eqs = Expr(:(=), ModelingToolkit.build_expr(:tuple, ls), ModelingToolkit.build_expr(:tuple, rs))

    fname = gensym(:ModelingToolkitFunction)
    out_expr = conv(op)
    let_expr = Expr(:let, var_eqs, Expr(:block, out_expr))
    bounds_block = checkbounds ? let_expr : :(@inbounds begin $let_expr end)

    fargs = Expr(:tuple,argnames...)
    oop_ex = headerfun(bounds_block, fargs, false)

    if !linenumbers
        oop_ex = striplines(oop_ex)
    end

    if expression == Val{true}
        return ModelingToolkit.inject_registered_module_functions(oop_ex)
    else
        _build_and_inject_function(@__MODULE__, oop_ex)
    end
end

function _build_and_inject_function(mod::Module, ex)
    # Generate the function, which will process the expression
    runtimefn = GeneralizedGenerated.mk_function(mod, ex)

    # Extract the processed expression of the function body
    params = typeof(runtimefn).parameters
    fn_expr = GeneralizedGenerated.NGG.from_type(params[3])

    # Inject our externally registered module functions
    new_expr = ModelingToolkit.inject_registered_module_functions(fn_expr)

    # Reconstruct the RuntimeFn's Body
    new_body = GeneralizedGenerated.NGG.to_type(new_expr)
    return GeneralizedGenerated.RuntimeFn{params[1:2]..., new_body, params[4]}()
end

# Detect heterogeneous element types of "arrays of matrices/sparce matrices"
function is_array_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray), F)
end
function is_array_sparse_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractSparseMatrix), F)
end
# Detect heterogeneous element types of "arrays of arrays of matrices/sparce matrices"
function is_array_array_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray{<:AbstractMatrix}), F)
end
function is_array_array_sparse_matrix(F)
    return isa(F, AbstractVector) && all(x->isa(x, AbstractArray{<:AbstractSparseMatrix}), F)
end

function fill_array_with_zero!(x::AbstractArray)
    if eltype(x) <: AbstractArray
        foreach(fill_array_with_zero!, x)
    else
        fill!(x, false)
    end
    return x
end

function _build_function(target::JuliaTarget, rhss, args...;
                         conv = simplified_expr, expression = Val{true},
                         checkbounds = false,
                         linenumbers = false, multithread=nothing,
                         headerfun=addheader, outputidxs=nothing,
						 convert_oop = true, force_SA = false,
                         skipzeros = outputidxs===nothing, fillzeros = skipzeros, parallel=SerialForm(), kwargs...)
	if multithread isa Bool
		@warn("multithraded is deprecated for the parallel argument. See the documentation.")
		parallel = multithread ? MultithreadedForm() : SerialForm()
	end

    argnames = [gensym(:MTKArg) for i in 1:length(args)]
    arg_pairs = map(vars_to_pairs,zip(argnames,args))
    ls = reduce(vcat,first.(arg_pairs))
    rs = reduce(vcat,last.(arg_pairs))
    var_eqs = Expr(:(=), ModelingToolkit.build_expr(:tuple, ls), ModelingToolkit.build_expr(:tuple, rs))

    fname = gensym(:ModelingToolkitFunction)
    fargs = Expr(:tuple,argnames...)

    oidx = isnothing(outputidxs) ? (i -> i) : (i -> outputidxs[i])
    X = gensym(:MTIIPVar)

	rhs_length = rhss isa SparseMatrixCSC ? length(rhss.nzval) : length(rhss)

	if parallel isa DistributedForm
		numworks = Distributed.nworkers()
		reducevars = [Variable(gensym(:MTReduceVar))() for i in 1:numworks]
		lens = Int(ceil(rhs_length/numworks))
		finalsize = rhs_length - (numworks-1)*lens
		_rhss = vcat(reduce(vcat,[[getindex(reducevars[i],j) for j in 1:lens] for i in 1:numworks-1],init=Expr[]),
						 [getindex(reducevars[end],j) for j in 1:finalsize])
    elseif parallel isa DaggerForm
		computevars = [Variable(gensym(:MTComputeVar))() for i in axes(rhss,1)]
        reducevar = Variable(gensym(:MTReduceVar))()
        _rhss = [getindex(reducevar,i) for i in axes(rhss,1)]
	elseif rhss isa SparseMatrixCSC
		_rhss = rhss.nzval
	else
		_rhss = rhss
	end

    ip_sys_exprs = Expr[]
    # we cannot reliably fill the array with the presence of index translation
    if fillzeros && outputidxs === nothing
        push!(ip_sys_exprs, :($fill_array_with_zero!($X)))
    end
    if is_array_array_sparse_matrix(rhss) # Array of arrays of sparse matrices
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhsel2) ∈ enumerate(rhsel)
                for (k, rhs) ∈ enumerate(rhsel2.nzval)
                    rhs′ = conv(rhs)
                    (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                    push!(ip_sys_exprs, :($X[$i][$j].nzval[$k] = $rhs′))
                end
            end
        end
    elseif is_array_array_matrix(rhss) # Array of arrays of arrays
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhsel2) ∈ enumerate(rhsel)
                for (k, rhs) ∈ enumerate(rhsel2)
                    rhs′ = conv(rhs)
                    (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                    push!(ip_sys_exprs, :($X[$i][$j][$k] = $rhs′))
                end
            end
        end
    elseif is_array_sparse_matrix(rhss) # Array of sparse matrices
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhs) ∈ enumerate(rhsel.nzval)
                rhs′ = conv(rhs)
                (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                push!(ip_sys_exprs, :($X[$i].nzval[$j] = $rhs′))
            end
        end
    elseif is_array_matrix(rhss) # Array of arrays
        for (i, rhsel) ∈ enumerate(_rhss)
            for (j, rhs) ∈ enumerate(rhsel)
                rhs′ = conv(rhs)
                (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
                push!(ip_sys_exprs, :($X[$i][$j] = $rhs′))
            end
        end
    elseif rhss isa SparseMatrixCSC
        for (i, rhs) ∈ enumerate(_rhss)
            rhs′ = conv(rhs)
            (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
            push!(ip_sys_exprs, :($X.nzval[$i] = $rhs′))
        end
    else
        for (i, rhs) ∈ enumerate(_rhss)
            rhs′ = conv(rhs)
            (skipzeros && rhs′ isa Number && iszero(rhs′)) && continue
            push!(ip_sys_exprs, :($X[$(oidx(i))] = $rhs′))
        end
    end

    ip_let_expr = Expr(:let, var_eqs, build_expr(:block, ip_sys_exprs))

    if parallel isa MultithreadedForm
        lens = Int(ceil(length(ip_let_expr.args[2].args)/Threads.nthreads()))
        threaded_exprs = vcat([quote
           Threads.@spawn begin
              $(ip_let_expr.args[2].args[((i-1)*lens+1):i*lens]...)
           end
        end for i in 1:Threads.nthreads()-1],
           quote
              Threads.@spawn begin
                 $(ip_let_expr.args[2].args[((Threads.nthreads()-1)*lens+1):end]...)
              end
           end)
        ip_let_expr.args[2] =  ModelingToolkit.build_expr(:block, threaded_exprs)
		ip_let_expr = :(@sync begin $ip_let_expr end)
    elseif parallel isa DistributedForm
		numworks = Distributed.nworkers()
		lens = Int(ceil(length(ip_let_expr.args[2].args)/numworks))
		spawnvars = [gensym(:MTSpawnVar) for i in 1:numworks]
		rhss_flat = rhss isa SparseMatrixCSC ? rhss.nzval : rhss
		spawnvectors = vcat(
					   [build_expr(:vect, [conv(rhs) for rhs ∈ rhss_flat[((i-1)*lens+1):i*lens]]) for i in 1:numworks-1],
					   build_expr(:vect, [conv(rhs) for rhs ∈ rhss_flat[((numworks-1)*lens+1):end]]))

        spawn_exprs = [quote
           $(spawnvars[i]) = ModelingToolkit.Distributed.remotecall($(i+1)) do
              $(spawnvectors[i])
           end
        end for i in 1:numworks]
        spawn_exprs = ModelingToolkit.build_expr(:block, spawn_exprs)
		resunpack_exprs = [:($(Symbol(reducevars[iter])) = fetch($(spawnvars[iter]))) for iter in 1:numworks]

		ip_let_expr.args[2] = quote
			@sync begin
				$spawn_exprs
				$(resunpack_exprs...)
				$(ip_let_expr.args[2])
			end
		end
    elseif parallel isa DaggerForm
        @assert HAS_DAGGER[] "Dagger.jl is not loaded; please do `using Dagger`"
        dagwrap(x) = x
        dagwrap(ex::Expr) = dagwrap(ex, Val(ex.head))
        dagwrap(ex::Expr, ::Val) = ex
        dagwrap(ex::Expr, ::Val{:call}) = :(Dagger.delayed($(ex.args[1]))($(dagwrap.(ex.args[2:end])...)))
        new_rhss = dagwrap.(conv.(rhss))
        delayed_exprs = build_expr(:block, [:($(Symbol(computevars[i])) = Dagger.delayed(identity)($(new_rhss[i]))) for i in axes(computevars,1)])
        # TODO: treereduce?
        reduce_expr = quote
            $(Symbol(reducevar)) = collect(Dagger.delayed(vcat)($(computevars...)))
        end
        ip_let_expr.args[2] = quote
			@sync begin
	            $delayed_exprs
	            $reduce_expr
	            $(ip_let_expr.args[2])
			end
        end
	end

    tuple_sys_expr = build_expr(:tuple, [conv(rhs) for rhs ∈ rhss])

    if rhss isa Matrix
        arr_sys_expr = build_expr(:vcat, [build_expr(:row,[conv(rhs) for rhs ∈ rhss[i,:]]) for i in 1:size(rhss,1)])
    elseif typeof(rhss) <: Array && !(typeof(rhss) <: Vector)
        vector_form = build_expr(:vect, [conv(rhs) for rhs ∈ rhss])
        arr_sys_expr = :(reshape($vector_form,$(size(rhss)...)))
    elseif rhss isa SparseMatrixCSC
        vector_form = build_expr(:vect, [conv(rhs) for rhs ∈ nonzeros(rhss)])
        arr_sys_expr = :(SparseMatrixCSC{eltype($(first(argnames))),Int}($(size(rhss)...), $(rhss.colptr), $(rhss.rowval), $vector_form))
    else # Vector
        arr_sys_expr = build_expr(:vect, [conv(rhs) for rhs ∈ rhss])
    end

	xname = gensym(:MTK)

	arr_sys_expr = (typeof(rhss) <: Vector || typeof(rhss) <: Matrix) && !(eltype(rhss) <: AbstractArray) ? quote
		if $force_SA || typeof($(fargs.args[1])) <: Union{ModelingToolkit.StaticArrays.SArray,ModelingToolkit.LabelledArrays.SLArray}
			$xname = ModelingToolkit.StaticArrays.@SArray $arr_sys_expr
			if $convert_oop && !(typeof($(fargs.args[1])) <: Number) && $(typeof(rhss) <: Vector) # Only try converting if it should match `u`
				return similar_type($(fargs.args[1]),eltype($xname))($xname)
			else
				return $xname
			end
		else
			$xname = $arr_sys_expr
			if $convert_oop && $(typeof(rhss) <: Vector)
				if !(typeof($(fargs.args[1])) <: Array) && !(typeof($(fargs.args[1])) <: Number) && eltype($(fargs.args[1])) <: eltype($xname)
					# Last condition: avoid known error because this doesn't change eltypes!
					return convert(typeof($(fargs.args[1])),$xname)
				elseif typeof($(fargs.args[1])) <: ModelingToolkit.LabelledArrays.LArray
					# LArray just needs to add the names back!
					return ModelingToolkit.LabelledArrays.LArray{ModelingToolkit.LabelledArrays.symnames(typeof($(fargs.args[1])))}($xname)
				else
					return $xname
				end
			else
				return $xname
			end
		end
	end : arr_sys_expr

    let_expr = Expr(:let, var_eqs, tuple_sys_expr)
    arr_let_expr = Expr(:let, var_eqs, arr_sys_expr)
    bounds_block = checkbounds ? let_expr : :(@inbounds begin $let_expr end)
    oop_bounds_block = checkbounds ? arr_let_expr : :(@inbounds begin $arr_let_expr end)
    ip_bounds_block = checkbounds ? ip_let_expr : :(@inbounds begin $ip_let_expr end)

    oop_ex = headerfun(oop_bounds_block, fargs, false)
    iip_ex = headerfun(ip_bounds_block, fargs, true; X=X)

    if !linenumbers
        oop_ex = striplines(oop_ex)
        iip_ex = striplines(iip_ex)
    end

    if expression == Val{true}
        return ModelingToolkit.inject_registered_module_functions(oop_ex), ModelingToolkit.inject_registered_module_functions(iip_ex)
    else
        return _build_and_inject_function(@__MODULE__, oop_ex), _build_and_inject_function(@__MODULE__, iip_ex)
    end
end

vars_to_pairs(args) = vars_to_pairs(args[1],args[2])
function vars_to_pairs(name,vs::AbstractArray)
	_vs = convert.(Variable,vs)
	names = [Symbol(u) for u ∈ _vs]
	exs = [:($name[$i]) for (i, u) ∈ enumerate(_vs)]
	names,exs
end

function vars_to_pairs(name,vs)
	_vs = convert(Variable,vs)
	names = [Symbol(_vs)]
	exs = [name]
	names,exs
end

get_varnumber(varop::Operation,vars::Vector{Operation}) =  findfirst(x->isequal(x,varop),vars)
get_varnumber(varop::Operation,vars::Vector{<:Variable})  =  findfirst(x->isequal(x,varop.op),vars)

function numbered_expr(O::Equation,args...;kwargs...)
  :($(numbered_expr(O.lhs,args...;kwargs...)) = $(numbered_expr(O.rhs,args...;kwargs...)))
end

function numbered_expr(O::Operation,vars,parameters;
                       derivname=:du,
                       varname=:u,paramname=:p)
  if isa(O.op, ModelingToolkit.Differential)
    varop = O.args[1]
    i = get_varnumber(varop,vars)
    return :($derivname[$i])
  elseif isa(O.op, ModelingToolkit.Variable)
    i = get_varnumber(O,vars)
    if i == nothing
      i = get_varnumber(O,parameters)
      return :($paramname[$i])
    else
      return :($varname[$i])
    end
  end
  return Expr(:call, Symbol(O.op),
         [numbered_expr(x,vars,parameters;derivname=derivname,
                        varname=varname,paramname=paramname) for x in O.args]...)
end

function numbered_expr(de::ModelingToolkit.Equation,vars::Vector{<:Variable},parameters;
                       derivname=:du,varname=:u,paramname=:p)
    i = findfirst(x->isequal(x.name,var_from_nested_derivative(de.lhs)[1].name),vars)
    :($derivname[$i] = $(numbered_expr(de.rhs,vars,parameters;
                                     derivname=derivname,
                                     varname=varname,paramname=paramname)))
end
function numbered_expr(de::ModelingToolkit.Equation,vars::Vector{Operation},parameters;
                       derivname=:du,varname=:u,paramname=:p)
    i = findfirst(x->isequal(x.op.name,var_from_nested_derivative(de.lhs)[1].name),vars)
    :($derivname[$i] = $(numbered_expr(de.rhs,vars,parameters;
                                     derivname=derivname,
                                     varname=varname,paramname=paramname)))
end
numbered_expr(c::ModelingToolkit.Constant,args...;kwargs...) = c.value

function _build_function(target::StanTarget, eqs, vs, ps, iv,
                         conv = simplified_expr, expression = Val{true};
                         fname = :diffeqf, derivname=:internal_var___du,
                         varname=:internal_var___u,paramname=:internal_var___p)
    differential_equation = string(join([numbered_expr(eq,vs,ps,derivname=derivname,
                                   varname=varname,paramname=paramname) for
                                   (i, eq) ∈ enumerate(eqs)],";\n  "),";")
    """
    real[] $fname(real $iv,real[] $varname,real[] $paramname,real[] x_r,int[] x_i) {
      real $derivname[$(length(eqs))];
      $differential_equation
      return $derivname;
    }
    """
end

function _build_function(target::CTarget, eqs, vs, ps, iv;
                         conv = simplified_expr, expression = Val{true},
                         fname = :diffeqf, derivname=:internal_var___du,
                         varname=:internal_var___u,paramname=:internal_var___p)
    differential_equation = string(join([numbered_expr(eq,vs,ps,derivname=derivname,
                                  varname=varname,paramname=paramname) for
                                  (i, eq) ∈ enumerate(eqs)],";\n  "),";")
    """
    void $fname(double* $derivname, double* $varname, double* $paramname, $iv) {
      $differential_equation
    }
    """
end

function _build_function(target::MATLABTarget, eqs, vs, ps, iv;
                         conv = simplified_expr, expression = Val{true},
                         fname = :diffeqf, derivname=:internal_var___du,
                         varname=:internal_var___u,paramname=:internal_var___p)
    matstr = join([numbered_expr(eq.rhs,vs,ps,derivname=derivname,
                                  varname=varname,paramname=paramname) for
                                  (i, eq) ∈ enumerate(eqs)],"; ")

    matstr = replace(matstr,"["=>"(")
    matstr = replace(matstr,"]"=>")")
    matstr = "$fname = @(t,$varname) ["*matstr*"];"
    matstr
end


text:
export cone_positiveorthant,
    cone_vectorspace,
    cone_positivenullspace,
    Newtonpolytope,
    verticesof,
    raysof,
    linearcombination,
    rays_outernormalcone,
    findnegativepoint,
    printfound,
    pRoots_qPossitive

## Nemo: matrix, FlintZZ, nullspace
## Polymake: polytope:cone,intersection
## LinearAlgebra: I

"""

    cone_positiveorthant(n)

Return the cone (Polymake big object) corresponding to the nonnegative orthant of R^n.

"""
function cone_positiveorthant(n)
    rays = cat(ones(Int, n)..., dims=(1, 2))
    return Polymake.polytope.Cone(INPUT_RAYS=rays)
end

"""

    cone_vectorspace(M::AbstractMatrix{T}) where {T <: Integer}

Return the cone (Polymake big object) corresponding to vector space generated by the columns of M.

"""
function cone_vectorspace(M::AbstractMatrix{T}) where {T <: Integer}
    rays = transpose(hcat(M, -M))
    return Polymake.polytope.Cone(INPUT_RAYS=rays)
end

"""

    cone_positivenullspace(N::AbstractMatrix{T}) where {T<:Integer}

Return the cone (Polymake big object) intersection of the nonnegative orthant and the nullspace of `N`.

"""
function cone_positivenullspace(N::AbstractMatrix{T}) where {T <: Integer}
    nullspace = Nemo.nullspace_right_rational(N)
    return Polymake.polytope.intersection(cone_positiveorthant(size(N, 2)), cone_vectorspace(nullspace))
end

function Newtonpolytope(p)
    points = transpose(collect_homoexponent_vectors(p))
    return Polymake.polytope.Polytope(POINTS=points)
end

function verticesof(polyt, polytname="")
    length(polytname) > 0 && print("Computing vertices of $(polytname), it may take some time.\n")
    return convert_to_array(polyt.VERTICES[:,2:end])
end

function raysof(cone, conename="")
    length(conename) > 0 && print("Computing rays of $(conename), it may take some time.\n")
    return Int.(integermultiple(Rational.(transpose(Array(cone.RAYS)))))
end

## Input:  polytope P, vertex v::Int.
## Output: point/vector in the outer normal cone
function linearcombination(A, coeffs=ones(Int, size(A, 2)))
    return A * coeffs
end

function rays_outernormalcone(polyt, vertex)
    cone = Polymake.polytope.normal_cone(polyt, vertex - 1, outer=1)
    return raysof(cone)
end

rays_outernormalcone(polyt) =  vertex -> rays_outernormalcone(polyt, vertex)

function findnegativepoint(p)
    Newtonp = Newtonpolytope(p)
    ver = verticesof(Newtonp, "Newton polytope")
    negver = Findallcols(isnegative, p, ver)
    points = []
    for v in negver
        rays = rays_outernormalcone(Newtonp, v)
        exponents = integermultiple(rays * ones(Int, size(rays, 2)))
        push!(points, collect_realpositiveroots(p, exponents))
    end
    return points
end

## Find roots of p for which q is possitive
function pRoots_qPossitive(p, q, nattemps::Integer=10, randbound::Integer=50)
    Newp = Newtonpolytope(p)
    Newq = Newtonpolytope(q)
    Vp = verticesof(Newp)
    Vq = verticesof(Newq, "q")
    println("Vertices of the Newton polytopes computed!!!!\n")
    #     P = NPq
    #     idx = posVq
    # else
    #     Cones = [Polymake.polytope.normal_cone(NPq, v-1, outer=1) for v in posVq]
    #     P = NPp
    #     idx = negVp
    #     ## Flip matrices of vertices to homogenize printing (the indices are flip)
    #     aux = VNPp
    #     VNPp = VNPq
    #     VNPq = aux
    #     end
    for i in posVq
        println("Computing cone to the $(i)th positive vertex of q\n")
        c1 = Polymake.polytope.normal_cone(NPq, i - 1, outer=1)
        for (nc2, c2) in enumerate(Conesp)
            rays = raysof(Polymake.polytope.intersection(c1, c2))
            r = size(rays, 1)
            if r > 0
                println("Intersecting cones found :)")
                found = false
                point = IntPointIncone(rays, ones(Int, 1, r))
                println("Computing real positive roots\n")
                proots = CollectallPossitiveRoots(p, point)
                qvals = [Float64(evaluate(q, (Nemo.fmpq).(pt))) for pt in proots]
                if findfirst(>(0), qvals) == nothing
                    j = 1;
                    while !found && j < nattemps
                        j += 1
                        println("Computing point in cone\n")
                        point = IntPointIncone(rays, rand(1:randbound, 1, r))
                        println("Computing real positive roots\n")
                        proots = CollectallPossitiveRoots(p, point)
                        qvals = [Float64(evaluate(q, (Nemo.fmpq).(pt))) for pt in proots]
                        found = findfirst(>(0), qvals) != nothing
                    end
                end
                if found
                    printfound(VNPp[negVp[nc2],:], VNPq[i,:], point, proots, qvals)
                end
            end
        end
    end
end

function printfound(Vp, Vq, point, proots, qvals)
    print("==============================================\n")
    print("==============================================\n")
    print("==============  Points found  ================\n")
    print("Negative Vertex p\n")
    print("$(Vp)\n")
    print("Possitive Vertex q\n")
    print("$(Vq)\n")
    print("Point in both outer normal cones\n")
    print("$(point)\n")
    for j in findall(>(0), qvals)
        print("$(proots[j])\n")
    end
end

function pRoots_qPossitive(p, NPp, VNPp, q, NPq, VNPq, nattemps::Integer=10, randbound::Integer=50)
    println("Finding negative/possitive vertices\n")
    isneg(x) = x < 0
    ispos(x) = 0 < x
    negVp = Findallrows(isneg, p, VNPp)
    posVq = Findallrows(ispos, q, VNPq)
    println("Computing cones to negative vertices of p\n")
    Conesp = [Polymake.polytope.normal_cone(NPp, v - 1, outer=1) for v in negVp]
    Listpt = []
    for i in posVq
        print("==============================================\n")
        print("==============================================\n")
        println("Computing the cone to the $(i)th positive vertex of q\n")
        c1 = Polymake.polytope.normal_cone(P, i - 1, outer=1)
        for (nc2, c2) in enumerate(Conesp)
            # print("==============================================\n")
            # print("==============================================\n")
            # println("Intersecting the $(nc2) first cone and the $i second\n")
            rays = raysof(Polymake.polytope.intersection(c1, c2))
            r = size(rays, 1)
            if r > 0
                # println("Intersecting cones found :)")
                found = false
                # println("Computing point in cone (coeffs 1)\n")
                point = IntPointIncone(rays, ones(Int, 1, r))
                # print("Point in both outer normal cones\n")
                # print("$(point)\n")
                # println("Computing real positive roots\n")
                proots = CollectallPossitiveRoots(p, point)
                qvals = [Float64(Nemo.evaluate(q, (Nemo.fmpq).(vec(pt)))) for pt in proots]
                pvals = [Float64(Nemo.evaluate(p, (Nemo.fmpq).(vec(pt)))) for pt in proots]
                found = findfirst(>(0), qvals) != nothing
                if !found
                    j = 1;
                    while !found && j < nattemps
                        j += 1
                        # println("Computing point in cone (rand coeffs, attempt $j)\n")
                        point = IntPointIncone(rays, rand(1:randbound, 1, r))
                        # print("Point in both outer normal cones\n")
                        # print("$(point)\n")
                        # println("Computing real positive roots\n")
                        proots = CollectallPossitiveRoots(p, point)
                        qvals = [Float64(Nemo.evaluate(q, (Nemo.fmpq).(vec(pt)))) for pt in proots]
                        pvals = [Float64(Nemo.evaluate(p, (Nemo.fmpq).(vec(pt)))) for pt in proots]
                        # print("Values of p for the roots\n")
                        # print("$(pvals)\n")
                        # print("Values of q for the roots\n")
                        # print("$(qvals)\n")
                        found = findfirst(>(0), qvals) != nothing
                    end
                end
                if found
                    printfound(VNPp[negVp[nc2],:], VNPq[i,:], point, proots, qvals)
                    Listpt = [Listpt; [[proots, qvals, pvals, point, VNPp[negVp[nc2],:], VNPq[i,:]]]]
                end
            end
        end
    end
    return Listpt
end


text:
using Statistics
using StatsBase:sample
using Revise
using EvoTrees
using BenchmarkTools
using CUDA

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

###########################
# Tree CPU
###########################
params_c = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32);

model_c, cache_c = EvoTrees.init_evotree(params_c, X_train, Y_train);

# initialize from cache
params_c = model_c.params
train_nodes = cache_c.train_nodes
splits = cache_c.splits
X_size = size(cache_c.X_bin)

# select random rows and cols
𝑖 = cache_c.𝑖_[sample(params_c.rng, cache_c.𝑖_, ceil(Int, params_c.rowsample * X_size[1]), replace=false, ordered=true)]
𝑗 = cache_c.𝑗_[sample(params_c.rng, cache_c.𝑗_, ceil(Int, params_c.colsample * X_size[2]), replace=false, ordered=true)]
# reset gain to -Inf
for feat in cache_c.𝑗_
    splits[feat].gain = -Inf
end

# build a new tree
# 897.800 μs (6 allocations: 736 bytes)
@time EvoTrees.update_grads!(params_c.loss, params_c.α, cache_c.pred_cpu, cache_c.Y_cpu, cache_c.δ, cache_c.δ², cache_c.𝑤)
∑δ, ∑δ², ∑𝑤 = sum(cache_c.δ[𝑖]), sum(cache_c.δ²[𝑖]), sum(cache_c.𝑤[𝑖])
gain = EvoTrees.get_gain(params_c.loss, ∑δ, ∑δ², ∑𝑤, params_c.λ)
# assign a root and grow tree
train_nodes[1] = EvoTrees.TrainNode(0, 1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)
# 69.247 ms (1852 allocations: 38.41 MiB)
@time tree = grow_tree(cache_c.δ, cache_c.δ², cache_c.𝑤, cache_c.hist_δ, cache_c.hist_δ², cache_c.hist_𝑤, params_c, train_nodes, splits, cache_c.edges, cache_c.X_bin);
@btime tree = grow_tree($cache_c.δ, $cache_c.δ², $cache_c.𝑤, $cache_c.hist_δ, $cache_c.hist_δ², $cache_c.hist_𝑤, $params_c, $train_nodes, $splits, $cache_c.edges, $cache_c.X_bin);
push!(model_c.trees, tree)
@btime EvoTrees.predict!(cache_c.pred_cpu, tree, cache_c.X)

δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, edges, X_bin = cache_c.δ, cache_c.δ², cache_c.𝑤, cache_c.hist_δ, cache_c.hist_δ², cache_c.hist_𝑤, cache_c.edges, cache_c.X_bin;

T = Float32
L = 1
active_id = ones(Int, 1)
leaf_count = one(Int)
tree_depth = one(Int)
tree = EvoTrees.Tree(Vector{EvoTrees.TreeNode{L,T,Int,Bool}}())

id = 1
node = train_nodes[id]
# 9.613 ms (81 allocations: 13.55 KiB)
@time EvoTrees.update_hist!(hist_δ[id], hist_δ²[id], hist_𝑤[id], δ, δ², 𝑤, X_bin, node)
@btime EvoTrees.update_hist!($hist_δ[id], $hist_δ²[id], $hist_𝑤[id], $δ, $δ², $𝑤, $X_bin, $node)

j = 1
# 601.685 ns (6 allocations: 192 bytes) 8 100 feat ~ 60us
@btime EvoTrees.find_split!(view(hist_δ[id], :, j), view(hist_δ²[id], :, j), view(hist_𝑤[id], :, j), params_c, node, splits[j], edges[j])

for j in node.𝑗
    splits[j].gain = node.gain
    EvoTrees.find_split!(view(hist_δ[id],:,j), view(hist_δ²[id],:,j), view(hist_𝑤[id],:,j), params_c, node, splits[j], edges[j])
end
best_cpu = EvoTrees.get_max_gain(splits)

set = node.𝑖
best = X_bin[3]
@btime EvoTrees.update_set(set, best, X_bin[:,1])
@btime EvoTrees.update_set(node.𝑖, best, view(X_bin, :, 1))

###################################################
# GPU
###################################################
params_g = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=32);

model_g, cache_g = EvoTrees.init_evotree_gpu(params_g, X_train, Y_train);

params_g = model_g.params;
train_nodes = cache_g.train_nodes;
# splits = cache_g.splits;
X_size = size(cache_g.X_bin);

# select random rows and cols
𝑖 = CuVector(cache_g.𝑖_[sample(params_g.rng, cache_g.𝑖_, ceil(Int, params_g.rowsample * X_size[1]), replace=false, ordered=true)])
𝑗 = CuVector(cache_g.𝑗_[sample(params_g.rng, cache_g.𝑗_, ceil(Int, params_g.colsample * X_size[2]), replace=false, ordered=true)])
# reset gain to -Inf
# splits.gains .= -Inf

# build a new tree
# 144.600 μs (23 allocations: 896 bytes) - 5-6 X time faster on GPU
@time CUDA.@sync EvoTrees.update_grads_gpu!(params_g.loss, cache_g.δ, cache_g.pred_gpu, cache_g.Y)
# sum Gradients of each of the K parameters and bring to CPU
∑δ = Array(vec(sum(cache_g.δ[𝑖,:], dims=1)))
gain = EvoTrees.get_gain_gpu(params_g.loss, ∑δ, params_g.λ)
# assign a root and grow tree
train_nodes[1] = EvoTrees.TrainNodeGPU(UInt32(0), UInt32(1), ∑δ, gain, 𝑖, 𝑗)
# 60.736 ms (108295 allocations: 47.95 MiB) - only 15% faster than CPU

EvoTrees.grow_tree_gpu(cache_g.δ, cache_g.hist, params_g, cache_g.K, train_nodes, cache_g.edges, cache_g.X_bin);
@btime CUDA.@sync tree = EvoTrees.grow_tree_gpu(cache_g.δ, cache_g.δ², cache_g.𝑤, cache_g.hist_δ, cache_g.hist_δ², cache_g.hist_𝑤, params_g, cache_g.K, train_nodes, splits, cache_g.edges, cache_g.X_bin, cache_g.X_bin_cpu);
push!(model_g.trees, tree);
# 2.736 ms (93 allocations: 13.98 KiB)
@btime CUDA.@sync EvoTrees.predict_gpu!(cache_g.pred_cpu, tree, cache_g.X)
# 1.013 ms (37 allocations: 1.19 KiB)
@btime CUDA.@sync cache_g.pred .= CuArray(cache_g.pred_cpu);

###########################
# Tree GPU
###########################
δ, hist, K, edges, X_bin = cache_g.δ, cache_g.hist, cache_g.K, cache_g.edges, cache_g.X_bin;
T = Float32
S = UInt32
active_id = ones(S, 1)
leaf_count = one(S)
tree_depth = one(S)
tree = EvoTrees.TreeGPU(Vector{EvoTrees.TreeNodeGPU{T,S,Bool}}())

id = S(1)
node = train_nodes[id];
# 2.930 ms (24 allocations: 656 bytes)
@time CUDA.@sync EvoTrees.update_hist_gpu!(hist[1], δ, X_bin, node.𝑖, node.𝑗, K);
@btime CUDA.@sync EvoTrees.update_hist_gpu!($hist[1], $δ, $X_bin, $node.𝑖, $node.𝑗, $K, MAX_THREADS=128);

j = 1
# 2.925 μs (78 allocations: 6.72 KiB) * 100 features ~ 300us
EvoTrees.find_split_gpu_v1!(hist[j], params_g, node, splits, edges, node.𝑗, K)
@btime CUDA.@sync EvoTrees.find_split_gpu_v1!(hist[j], edges, params_g)

j = 1
# 347.199 μs (403 allocations: 13.31 KiB)
best_g = EvoTrees.find_split_gpu_v2!(hist[j], edges, params_g);
best_g[1]
best_g

@btime CUDA.@sync EvoTrees.find_split_gpu_v2!(hist[j], edges, params_g);

# 673.900 μs (600 allocations: 29.39 KiB)
left, right = EvoTrees.update_set_gpu(node.𝑖, 16, X_bin[:,1], MAX_THREADS=1024);
@btime CUDA.@sync EvoTrees.update_set_gpu(node.𝑖, 16, X_bin[:,1], MAX_THREADS=1024);


text:
using SplineMakie
using GeometryBasics
using Test
using Makie
using GLMakie

@testset "SplineMakie.jl" begin
    @testest "scatterlines of markers" begin
        path = SplineMakie.BezierCircle
        scatterlines(discretize(path))

        path = SplineMakie.BezierUTriangle
        scatterlines(discretize(path))

        path = SplineMakie.BezierDTriangle
        scatterlines(discretize(path))

        path = SplineMakie.BezierSquare
        scatterlines(discretize(path))

        path = SplineMakie.BezierCross
        scatterlines(discretize(path))
    end

    @testest "convert arguments for lines" begin
        path = SplineMakie.BezierCircle
        lines(path)
    end

    @testest "cubic spline" begin
        # 2d
        using SplineMakie: splinepath
        p = Point2f0[]
        for x in 0:10
            push!(p, Point2f0(x, rand()))
        end

        scatter(p)
        lines!(splinepath(p))

        # 3d
        p = rand(Point3f0, 5)
        scatter(p)
        lines!(splinepath(p))

    end

    @testset "recipe" begin
        p = rand(Point2f0, 5)
        spline(p)
        scatter!(p)

        p = rand(Point3f0, 5)
        spline(p)
        scatter!(p, markersize=20)

        x = [1.0,2,3]
        y = [1.0,-1,5]
        z = [1.0,-4.5,0.2]
        spline(x, y)
        spline(x, y, z)

        p = [(0,0), (-.2,.5), (0,1)]
        spline(p)

        p = [(0,0), (1,1)]
        f, ax, p = spline(p)
        p.plots[1][1][]
    end

end


text:
export selector, clearconsole, input, structure, @sh, @profiler

"""
    selector([xs...]) -> x

Allow the user to select one of the `xs`.

`xs` should be an iterator of strings. Currently there is no fallback in other
environments.
"""
selector(xs) = Main.Atom.selector(xs)

"""
    clearconsole()

Clear the console if Juno is used; does nothing otherwise.
"""
clearconsole() = isactive() && Main.Atom.clearconsole()

"""
    input(prompt = "") -> "..."

Prompt the user to input some text, and return it. Optionally display a prompt.
"""
input(prompt = "") = (print(prompt); isactive() ? Main.Atom.input() : readline())

"""
    info(msg)

Show the given message in Juno's console using blue styling, or fall back to
`Base.info`.

In a package, you can use `import Juno: info` to replace the default version
with this one.
"""
info(msg) = (isactive() ? Main.Atom : Base).info(msg)

"""
    notify(msg)

Display `msg` as an OS specific notification.

Useful for signaling the end of a long running computation or similar. This
disregards the `Notifications` setting in `julia-client`. Falls back to
`info(msg)` in other environments.
"""
notify(msg::AbstractString) = isactive() ? Main.Atom.sendnotify(msg) : @info(msg)

"""
    plotsize()

Get the size of Juno's plot pane in `px`. Returns `[100, 100]` as a fallback value.
"""
plotsize() = isactive() ? Main.Atom.plotsize() : [100, 100]

"""
    syntaxcolors(selectors = Atom.SELECTORS)::Dict{String, UInt32}

Get the colors used by the current Atom theme.
`selectors` should be a `Dict{String, Vector{String}}` which assigns a css
selector (e.g. `syntax--julia`) to a name (e.g. `variable`).
"""
syntaxcolors(selectors) = isactive() ? Main.Atom.syntaxcolors(selectors) : Dict{String, UInt32}()
syntaxcolors() = isactive() ? Main.Atom.syntaxcolors() : Dict{String, UInt32}()

"""
    profiler(data=Profile.fetch(); lidict=nothing, C=false, combine=true, recur=:off, pruned=[])

Show profile information as an in-editor flamechart.
Any keyword argument that [`FlameGraphs.flamegraph`](@ref) can take could be given.
"""
profiler(args...; kwargs...) = isactive() && Main.Atom.Profiler.profiler(args...; kwargs...)

"""
    @profiler exp [kwargs...]

Clear currently collected profile traces, profile the provided expression and show
  it via `Juno.profiler()`.

Any keyword argument that [`FlameGraphs.flamegraph`](@ref) can take could be given
  as optional arugments `kwargs...`

```julia
# profile a function call
@profiler fname(fargs)

# include ccalls and compress recursive calls
@profiler fname(fargs) C = true recur = :flat
```
"""
macro profiler(exp, kwargs...)
  quote
    let
      $(Profile).clear()
      res = $(Profile).@profile $(esc(exp))
      profiler(; $(map(esc, kwargs)...))
      res
    end
  end
end

"""
    profiletree()

Show currently collected profile information in tree-form. Falls back to `Profile.print()`.
"""
profiletree() = isactive() ? Main.Atom.Profiler.tree() : Profile.print()

"""
    structure(x)

Display `x`'s underlying representation, rather than using its normal display
method.

For example, `structure(:(2x+1))` displays the `Expr` object with its
`head` and `args` fields instead of printing the expression.
"""
structure(args...) = Main.Atom.structure(args...)


text:
using KernelAbstractions, Test

@kernel function kernel_print()
  I = @index(Global)
  @print("Hello from thread ", I, "!\n")
end

function printing_testsuite(backend)
    @testset "print test" begin
        kernel = kernel_print(backend(), 4)
        wait(kernel(ndrange=(4,)))
        @test true

        @print("Why this should work\n")
        @test true
    end
end


text:
function predict_f(gp::GPBase, x::AbstractMatrix; full_cov::Bool=false)
    size(x,1) == gp.dim || throw(ArgumentError("Gaussian Process object and input observations do not have consistent dimensions"))
    if full_cov
        return predict_full(gp, x)
    else
        ## Calculate prediction for each point independently
        μ = Array{eltype(x)}(undef, size(x,2))
        σ2 = similar(μ)
        for k in 1:size(x,2)
            m, sig = predict_full(gp, x[:,k:k])
            μ[k] = m[1]
            σ2[k] = max(diag(sig)[1], 0.0)
        end
        return μ, σ2
    end
end

predict_full(gp::GPA, xpred::AbstractMatrix, Q::Approx) = predictMVN(gp,xpred, gp.x, gp.y, gp.kernel, gp.mean, gp.v, gp.covstrat, Q)

function predictMVN(xpred::AbstractMatrix, xtrain::AbstractMatrix, ytrain::AbstractVector,
    kernel::Kernel, meanf::Mean, alpha::AbstractVector,
    covstrat::CovarianceStrategy, Ktrain::AbstractPDMat)
crossdata = KernelData(kernel, xtrain, xpred)
priordata = KernelData(kernel, xpred, xpred)
Kcross = cov(kernel, xtrain, xpred, crossdata)
Kpred = cov(kernel, xpred, xpred, priordata)
mx = mean(meanf, xpred)
mu, Sigma_raw = predictMVN!(Kpred, Ktrain, Kcross, mx, alpha)
return mu, Sigma_raw
end

function predictMVN!(Kxx, Kff, Kfx, mx, αf)
    mu = mx + Kfx' * αf
    Lck = whiten!(Kff, Kfx)
    subtract_Lck!(Kxx, Lck)
    return mu, Kxx
end

@inline function subtract_Lck!(Sigma_raw::AbstractArray{<:AbstractFloat}, Lck::AbstractArray{<:AbstractFloat})
    LinearAlgebra.BLAS.syrk!('U', 'T', -1.0, Lck, 1.0, Sigma_raw)
    LinearAlgebra.copytri!(Sigma_raw, 'U')
end

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
pyplot()
cd(pwd)
dir = "results/"
font_size = 18

# Plot mean run time as a function of number of data points (Nd) for each sampler
meantimePlot = plotsummary(results, :Nd, :time, (:sampler,); save=true, dir=dir, yscale=:log10,
  xaxis=font(font_size), yaxis=font(font_size), legendfont = 12, size=(500,500), width=3
  )

# Plot mean allocations as a function of number of data points (Nd) for each sampler
meanallocPlot = plotsummary(results, :Nd, :allocations, (:sampler,); save=true, dir=dir, yscale=:log10,
  ylabel="Allocations (log scale)", xaxis=font(font_size), yaxis=font(font_size), legendfont = 12, size=(500,500), width=3
)

# Plot density of effective sample size as function of number of data points (Nd) for each sampler
essPlots = plotdensity(results, :ess, (:sampler, :Nd); save=true, dir=dir,
xaxis=font(font_size), yaxis=font(font_size), legendfont = 12, size=(700,700), width=3
)