text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks

using ChainRulesCore
import ForwardDiff, Zygote

# helper to collect variable IdSet
function indiciesForValueReferences(fmu::FMU2,
                                    refs::Array{fmi2ValueReference})
    array = fmi2ValueReference[]
    for i in 1:length(fmu.modelDescription.valueReferences)
        if fmu.modelDescription.valueReferences[i] in refs
            push!(array, i)
        end
    end
    array
end

"""
Performs something similar to `fmiDoStep` for ME-FMUs (note, that fmiDoStep is for CS-FMUs only).
Event handling (state- and time-events) is supported. If you don't want events to be handled, you can disable event-handling for the NeuralFMU `nfmu` with the attribute `eventHandling = false`.

Optional, additional FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, additional FMU-values can be retrieved by keyword argument `getValueReferences`.

Function takes the current system state array ("x") and returns an array with state derivatives ("x dot") and optionally the FMU-values for `getValueReferences`.
Setting the FMU time via argument `t` is optional, if not set, the current time of the ODE solver around the NeuralFMU is used.
"""
function fmi2EvaluateME(comp::FMU2Component,
        x::Array{<:Real},
        t = -1.0,#::Real,
        setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
        setValues::Array{<:Real} = zeros(Real, 0),
        getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    _fmi2EvaluateME(comp, x, t, setValueReferences, setValues, getValueReferences)
end
function fmi2EvaluateME(comp::FMU2Component,
    x::Array{<:ForwardDiff.Dual{Tx, Vx, Nx}},
    t = -1.0,#::Real,
    setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
    setValues::Array{<:Real} = zeros(Real, 0),
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tx, Vx, Nx}

    _fmi2EvaluateME_fd((Tx, Vx, Nx), (Tx, Vx, Nx), comp, x, t, setValueReferences, setValues, getValueReferences)
end
function fmi2EvaluateME(comp::FMU2Component,
    x::Array{<:ForwardDiff.Dual{Tx, Vx, Nx}},
    t,#::Real,
    setValueReferences::Array{fmi2ValueReference},
    setValues::Array{<:ForwardDiff.Dual{Tu, Vu, Nu}},
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tx, Vx, Nx, Tu, Vu, Nu}

    _fmi2EvaluateME_fd((Tx, Vx, Nx), (Tu, Vu, Nu), comp, x, t, setValueReferences, setValues, getValueReferences)
end

# ForwardDiff backend using the existing ChainRulesCore.frule
# adapted from: https://discourse.julialang.org/t/chainrulescore-and-forwarddiff/61705/8
function _fmi2EvaluateME_fd(TVNx, TVNu, comp, x, t, setValueReferences, setValues, getValueReferences) #where {T, V, N} #where {T} # <:ForwardDiff.Dual}

    Tx, Vx, Nx = TVNx
    Tu, Vu, Nu = TVNu

    ȧrgs = [NoTangent(), NoTangent(), collect(ForwardDiff.partials(e) for e in x), ForwardDiff.partials(t), NoTangent(), collect(ForwardDiff.partials(e) for e in setValues), NoTangent()]
    args = [fmi2EvaluateME, comp, collect(ForwardDiff.value(e) for e in x), ForwardDiff.value(t), setValueReferences, collect(ForwardDiff.value(e) for e in setValues), getValueReferences]

    # ToDo: Find a good fix!
    #Zygote.@ignore @debug "From $(typeof(args[6]))"
    if typeof(args[6]) == Vector{Any}
        args[6] = convert(Vector{Float64}, args[6])
        #Zygote.@ignore @debug "To $(typeof(args[6]))"
    end

    ȧrgs = (ȧrgs...,)
    args = (args...,)

    @assert (typeof(args[3]) == Vector{Float64}) || (typeof(args[3]) == Vector{Float32}) "After conversion, `x` is still an invalid type `$(typeof(args[3]))`."

    y, _, dx, _, _, du, _ = ChainRulesCore.frule(ȧrgs, args...)

    if Vx != Float64
        Vx = Float64
    end

    if Vu != Float64
        Vu = Float64
    end

    # original function returns [dx_1, ..., dx_n, y_1, ..., y_m]
    # ToDo: Add sensitivities (partials) from du -> (dx, y)
    [collect( ForwardDiff.Dual{Tx, Vx, Nx}(y[i], dx[i]) for i in 1:length(dx) )...]
end

function _fmi2EvaluateME(comp::FMU2Component,
                      x::Array{<:Real},
                      t::Real,
                      setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                      setValues::Array{<:Real} = zeros(Real, 0),
                      getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    @assert fmi2IsModelExchange(comp.fmu) ["fmi2EvaluateME(...): As in the name, this function only supports ME-FMUs."]

    setter = (length(setValueReferences) > 0)
    getter = (length(getValueReferences) > 0)

    if setter
        @assert length(setValueReferences) == length(setValues) ["fmi2EvaluateME(...): `setValueReferences` and `setValues` need to be the same length!"]
    end

    if t >= 0.0
        fmi2SetTime(comp, t)
    end

    if setter
        fmi2SetReal(comp, setValueReferences, setValues)
    end

    fmi2SetContinuousStates(comp, x)

    y = []
    if getter
        y = fmi2GetReal(comp, getValueReferences)
    end

    dx = fmi2GetDerivatives(comp)

    [dx..., y...]
end

function ChainRulesCore.rrule(::typeof(fmi2EvaluateME),
                              comp::FMU2Component,
                              x::Array{<:Real},
                              t::Real = -1.0,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2EvaluateME(comp, x, t, setValueReferences, setValues, getValueReferences)
    function fmi2EvaluateME_pullback(ȳ)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.rrule(fmi2EvaluateME, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        if t >= 0.0
            fmi2SetTime(comp, t)
        end

        fmi2SetContinuousStates(comp, x)

        rdx = vcat(comp.fmu.modelDescription.derivativeValueReferences, getValueReferences)
        rx = comp.fmu.modelDescription.stateValueReferences
        ru = setValueReferences

        n_dx_x = NoTangent()
        n_dx_u = NoTangent()

        if comp.senseFunc == :full || comp.senseFunc == :directionalDerivatives || comp.senseFunc == :auto
            # OPTIMIZATION: compute new jacobians only if system state or time changed, otherwise return the cached one
            if comp.jac_x != x || comp.jac_t != t

                if size(comp.jac_dxy_x) != (length(rdx), length(rx))
                    comp.jac_dxy_x = zeros(length(rdx), length(rx))
                end
                comp.jacobianFct(comp.jac_dxy_x, comp, rdx, rx)

                if size(comp.jac_dxy_u) != (length(rdx), length(ru))
                    comp.jac_dxy_u = zeros(length(rdx), length(ru))
                end
                comp.jacobianFct(comp.jac_dxy_u, comp, rdx, ru)

                comp.jac_x = x
                comp.jac_t = t
            end

            n_dx_x = @thunk(comp.jac_dxy_x' * ȳ)

            if setter
                n_dx_u = @thunk(comp.jac_dxy_u' * ȳ)
            end
        elseif comp.senseFunc == :adjointDerivatives
            @assert false "Adjoint Derivatives not supported by FMI2."
        else
            @assert false "`senseFunc=$(comp.senseFunc)` unknown value for `senseFunc`."
        end

        f̄ = NoTangent()
        c̄omp = ZeroTangent()
        x̄ = n_dx_x
        t̄ = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return f̄, c̄omp, x̄, t̄, s̄etValueReferences, s̄etValues, ḡetValueReferences
    end
    return y, fmi2EvaluateME_pullback
end

function ChainRulesCore.frule((Δself, Δcomp, Δx, Δt, ΔsetValueReferences, ΔsetValues, ΔgetValueReferences),
                              ::typeof(fmi2EvaluateME),
                              comp, #::FMU2Component,
                              x,#::Array{<:Real},
                              t,#::Real = -1.0,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2EvaluateME(comp, x, t, setValueReferences, setValues, getValueReferences)
    function fmi2EvaluateME_pullforward(Δx, ΔsetValues)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.frule(fmi2EvaluateME, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        if t >= 0.0
            fmi2SetTime(comp, t)
        end

        if all(isa.(x, ForwardDiff.Dual))
            xf = collect(ForwardDiff.value(e) for e in x)
            fmi2SetContinuousStates(comp, xf)
        else
            fmi2SetContinuousStates(comp, x)
        end

        rdx = vcat(comp.fmu.modelDescription.derivativeValueReferences, getValueReferences)
        rx = comp.fmu.modelDescription.stateValueReferences
        ru = setValueReferences

        n_dx_x = NoTangent()
        n_dx_u = NoTangent()

        if comp.senseFunc == :auto || comp.senseFunc == :full
            # OPTIMIZATION: compute new jacobians only if system state or time changed, otherwise return the cached one
            # ToDo: Optimize for getDirectionalDerivatives with seed vector Δx
            if comp.jac_x != x || comp.jac_t != t

                if size(comp.jac_dxy_x) != (length(rdx), length(rx))
                    comp.jac_dxy_x = zeros(length(rdx), length(rx))
                end
                comp.jacobianFct(comp.jac_dxy_x, comp, rdx, rx)

                if size(comp.jac_dxy_u) != (length(rdx), length(ru))
                    comp.jac_dxy_u = zeros(length(rdx), length(ru))
                end
                comp.jacobianFct(comp.jac_dxy_u, comp, rdx, ru)

                comp.jac_x = x
                comp.jac_t = t
            end

            n_dx_x = comp.jac_dxy_x * Δx

            if setter
                n_dx_u = comp.jac_dxy_u * ΔsetValues
            end
        elseif comp.senseFunc == :directionalDerivatives
            n_dx_x = fmi2GetDirectionalDerivative(comp, rdx, rx, Δx)
            if setter
                n_dx_u = fmi2GetDirectionalDerivative(comp, rdx, ru, ΔsetValues)
            end
        elseif comp.senseFunc == :adjointDerivatives
            @assert false "Adjoint Derivatives not supported by FMI2."
        else
            @assert false "`senseFunc=$(comp.senseFunc)` unknown value for `senseFunc`."
        end

        c̄omp = ZeroTangent()
        x̄ = n_dx_x
        t̄ = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (c̄omp, x̄, t̄, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return (y, fmi2EvaluateME_pullforward(Δx, ΔsetValues)...)
end

"""
Performs a fmiDoStep for CS-FMUs (note, that fmiDoStep is for CS-FMUs only).

Optional, FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, FMU-values can be retrieved by keyword argument `getValueReferences`.

Function returns the FMU-values for the optional keyword argument `getValueReferences`.
The CS-FMU performs one macro step with step size `dt`. Dependent on the integrated numerical solver, the FMU may perform multiple (internal) micro steps if needed to meet solver requirements (stability/accuracy). These micro steps are hidden by FMI2.
"""
function fmi2DoStepCS(comp::FMU2Component,
                      dt::Real,
                      setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                      setValues::Array{<:Real} = zeros(Real, 0),
                      getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))
    _fmi2DoStepCS(comp, dt, setValueReferences, setValues, getValueReferences)
end
function fmi2DoStepCS(comp::FMU2Component,
    dt,#::Real,
    setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
    setValues::Array{<:ForwardDiff.Dual{Tu, Vu, Nu}} = Array{ForwardDiff.Dual{Tu, Vu, Nu}}(),
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tu, Vu, Nu}

    _fmi2DoStepCS_fd((Td, Vd, Nd), comp, dt, setValueReferences, setValues, getValueReferences)
end

# Helper because keyword arguments are (currently) not supported by Zygote.
function _fmi2DoStepCS(comp::FMU2Component,
                       dt::Real,
                       setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                       setValues::Array{<:Real} = zeros(Real, 0),
                       getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    @assert fmi2IsCoSimulation(comp.fmu) ["fmi2DoStepCS(...): As in the name, this function only supports CS-FMUs."]
    @assert length(setValueReferences) == length(setValues) ["fmi2DoStepCS(...): `setValueReferences` ($(length(setValueReferences))) and `setValues` ($(length(setValues))) need to be the same length!"]

    if length(setValueReferences) > 0
        fmi2SetReal(comp, setValueReferences, setValues)
    end

    fmi2DoStep(comp, dt)

    y = zeros(Float64, 0)

    if length(getValueReferences) > 0
        y = fmi2GetReal(comp, getValueReferences)
    end

    y
end

# ForwardDiff backend using the existing ChainRulesCore.frule
# adapted from: https://discourse.julialang.org/t/chainrulescore-and-forwarddiff/61705/8
function _fmi2DoStepCS_fd(TVNu,
    comp,
                          dt,
                          setValueReferences,
                          setValues,
                          getValueReferences)

    Tu, Vu, Nu = TVNu

    ȧrgs = [NoTangent(), NoTangent(), ForwardDiff.partials(dt), NoTangent(), collect(ForwardDiff.partials(e) for e in setValues), NoTangent()]
    args = [fmi2DoStepCS, comp, ForwardDiff.value(dt), setValueReferences, collect(ForwardDiff.value(e) for e in setValues), getValueReferences]

    # ToDo: Find a good fix!
    if typeof(args[5]) == Vector{Any}
        args[5] = convert(Vector{Float64}, args[5])
    end

    ȧrgs = (ȧrgs...,)
    args = (args...,)

    y, _, _, _, du, _ = ChainRulesCore.frule(ȧrgs, args...)

    if Vu != Float64
        Vu = Float64
    end

    # original function returns [y_1, ..., y_m]
    [collect( ForwardDiff.Dual{Tu, Vu, Nu}(y[i], du[i]) for i in 1:length(du) )...]
end

function ChainRulesCore.rrule(::typeof(fmi2DoStepCS),
                              comp::FMU2Component,
                              dt::Real = -1.0,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2DoStepCS(comp, dt, setValueReferences, setValues, getValueReferences)
    function fmi2DoStepCS_pullback(ȳ)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.rrule(fmi2DoStepCS, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        rdx = getValueReferences
        rx = setValueReferences

        n_dx_u = ZeroTangent()

        if getter
            mat = zeros(length(rdx), length(rx))
            comp.jacobianFct(mat, comp, rdx, rx)
            n_dx_u = @thunk(mat' * ȳ)
        end

        f̄ = NoTangent()
        f̄mu = ZeroTangent()
        d̄t = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (f̄, f̄mu, d̄t, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return y, fmi2DoStepCS_pullback
end

function ChainRulesCore.frule((Δself, Δcomp, Δdt, ΔsetValueReferences, ΔsetValues, ΔgetValueReferences),
                              ::typeof(fmi2DoStepCS),
                              comp, #::FMU2,
                              dt,#::Real = -1.0,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2DoStepCS(comp, dt, setValueReferences, setValues, getValueReferences)
    function fmi2DoStepCS_pullforward(ΔsetValues)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.frule(fmi2DoStepCS, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        rdx = getValueReferences
        rx = setValueReferences

        n_dx_u = ZeroTangent()

        if getter
            mat = zeros(length(rdx), length(rx))
            comp.jacobianFct(mat, comp, rdx, rx)
            n_dx_u = mat * ΔsetValues
        end

        f̄mu = ZeroTangent()
        d̄t = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (f̄mu, d̄t, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return (y, fmi2DoStepCS_pullforward(ΔsetValues)...)
end

"""
    fmi2InputDoStepCSOutput(comp::FMU2Component,
                            dt::Real,
                            u::Array{<:Real})

Sets all FMU inputs to `u`, performs a ´´´fmi2DoStep´´´ and returns all FMU outputs.
"""
function fmi2InputDoStepCSOutput(comp::FMU2Component,
                                 dt::Real,
                                 u::Array{<:Real})

    @assert fmi2IsCoSimulation(comp.fmu) ["fmi2InputDoStepCSOutput(...): As in the name, this function only supports CS-FMUs."]

    fmi2DoStepCS(comp, dt,
                 comp.fmu.modelDescription.inputValueReferences,
                 u,
                 comp.fmu.modelDescription.outputValueReferences)
end

# FMU wrappers

function fmi2EvaluateME(fmu::FMU2, args...; kwargs...)
    fmi2EvaluateME(fmu.components[end], args...; kwargs...)
end

function fmi2DoStepCS(fmu::FMU2, args...; kwargs...)
    fmi2DoStepCS(fmu.components[end], args...; kwargs...)
end

function fmi2InputDoStepCSOutput(fmu::FMU2, args...; kwargs...)
    fmi2InputDoStepCSOutput(fmu.components[end], args...; kwargs...)
end


text:
const RSTRINGMACRO = r"@([^\s\(\)]*)_str\((\"[^\"]*\")\)"  # care about triple quotes?
function sugar_stringmacro(s)
    m = match(RSTRINGMACRO, s)
    string(m[1], m[2])
end
function sugar_stringmacros(s)
    replace(s, RSTRINGMACRO, sugar_stringmacro)
end

function sugar(s)
    sugar_stringmacros(s)
end

ext: jl


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

JLLWrappers.@generate_main_file_header("MMseqs2")
JLLWrappers.@generate_main_file("MMseqs2", UUID("2f0138a6-2765-53ae-bca2-2fe9d5077335"))
end  # module MMseqs2_jll


text:
@testset "Statistics" begin
    @testset "Statistics: Ordinary Least Squares" begin
        a = [10, 14,17,21,20, 18, 42, 51, 77, 11, 91]
        b = [0.11,0.7,0.2,0.19,0.09,0.8, 0.71,0.1, 0.6,0.3,0.81]
        OLSbeta(a,b)

        df = DataFrame(a=a, b=b)
        ols = GLM.lm(@formula(a~b), df)
        coef(ols) # print coefficients

        # Predict
        predict(ols) # predict the model
    end

end

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
#Creates a bifurcation diagram of a specific reaction network and for the given parameter and range.
function bifurcation_diagram(model::DiffEqBase.AbstractReactionNetwork, pOrg::AbstractVector, parameter::Symbol, range::AbstractVector, replacements::Expr...;mult=false::Bool,progress_update=2147483640.::Float64)
    p_idx = findfirst(model.params.==parameter)
    mult && (range*=pOrg[p_idx])
    params_const = Dict(zip(deleteat!(deepcopy(model.params),p_idx), deleteat!(deepcopy(pOrg),p_idx)))
    expr = meta_program(deepcopy(model.f_func), model.syms, params_const,replacements...,allowed_exp_params=Set{Symbol}([parameter]))
    jacobian = map(val->recursive_replace!(val,params_const),deepcopy(model.symjac))
    fps = Vector{Vector{Vector{Float64}}}(undef,length(range))
    stabilities = Vector{Vector{Tuple{Bool,Bool}}}(undef,length(range))
    colors = Vector{Vector{Symbol}}(undef,length(range))
    @showprogress progress_update for i = 1:length(range)
        fps[i] = eval(recursive_replace!(deepcopy(expr),Dict(parameter=>range[i])))
        stabilities[i] = map(fp->stability_internal(jacobian,Dict(zip(model.syms,fp)),Dict(parameter=>range[i])),fps[i])
        colors[i] = bifurcation_color.(stabilities[i])
    end
    return bif_diagram(parameter,range,fps,stabilities,colors,deepcopy(model.syms))
end

#Bifurcation diagram structure, contains all the information in a bifurcation diagram.
struct bif_diagram
    param::Symbol
    param_values::AbstractVector
    fixed_points::Vector{Vector{Vector{Float64}}}
    stabilities::Vector{Vector{Tuple{Bool,Bool}}}
    colors::Vector{Vector{Symbol}}
    reactant_syms::Vector{Symbol}
end

#Creates a two dimensional bifurcation diagram of a specific reaction network and for the given parameters and ranges.
function bifurcation_diagram_2d(model::DiffEqBase.AbstractReactionNetwork, pOrg::AbstractVector, parameter1::Symbol, parameter2::Symbol, range1::AbstractVector, range2::AbstractVector, replacements::Expr...;mult=false::Bool,progress_update=2147483640.::Float64)
    p_idx1 = findfirst(model.params.==parameter1); p_idx2 = findfirst(model.params.==parameter2);
    mult && (range1*=pOrg[p_idx1]; range2*=pOrg[p_idx2];)
    params_const = Dict(zip(deleteat!(deepcopy(model.params),sort([p_idx1,p_idx2])), deleteat!(deepcopy(pOrg),sort([p_idx1,p_idx2]))))
    expr = meta_program(deepcopy(model.f_func), model.syms, params_const,replacements...,allowed_exp_params=Set{Symbol}([parameter1,parameter2]))
    jacobian = map(val->recursive_replace!(val,params_const),deepcopy(model.symjac))
    fps = Matrix{Vector{Vector{Float64}}}(undef,length(range1),length(range2))
    stabilities = Matrix{Vector{Tuple{Bool,Bool}}}(undef,length(range1),length(range2))
    colors = Matrix{Vector{Symbol}}(undef,length(range1),length(range2))
    @showprogress progress_update for i = 1:length(range1), j = 1:length(range2)
        fps[i,j] = eval(recursive_replace!(deepcopy(expr),Dict(parameter1=>range1[i],parameter2=>range2[j])))
        stabilities[i,j] = map(fp->stability_internal(jacobian,Dict(zip(model.syms,fp)),Dict(parameter1=>range1[i],parameter2=>range2[j])),fps[i,j])
        colors[i,j] = bifurcation_color.(stabilities[i,j])
    end
    return bif_diagram_2d(parameter1,parameter2,range1,range2,fps,stabilities,colors,deepcopy(model.syms))
end

#Two dimensional bifurcation diagram structure, contains all the information in a two dimensional bifurcation diagram.
struct bif_diagram_2d
    param1::Symbol
    param2::Symbol
    param_values1::AbstractVector
    param_values2::AbstractVector
    fixed_points::Matrix{Vector{Vector{Float64}}}
    stabilities::Matrix{Vector{Tuple{Bool,Bool}}}
    colors::Matrix{Vector{Symbol}}
    reactant_syms::Array{Symbol}
end

#Function for determening stability of a reaction network. Only used internally by bifurcation diagram functions. Not to be used by the user.
function stability_internal(jacobian::Matrix{Expr},var_vals::Dict,par_vals::Dict)
    jac = map(val->recursive_replace!(recursive_replace!(val,var_vals),par_vals),deepcopy(jacobian))
    eigenvalues = eigvals(eval.(jac))
    return (maximum(real(eigenvalues))<1e-6,any(imag(eigenvalues).>1e-6))
end

#Determines the color of a point in the bifurcation diagram for the given stability information.
function bifurcation_color(stability::Tuple{Bool,Bool})
    stability[1]&&stability[2]&&(return :green)
    stability[1]&&(return :blue)
    stability[2]&&(return :yellow)
    return :red
end

max_stars_repo_name: TorkelE/BRNEquilibrate

max_issues_repo_name: TorkelE/BRNEquilibrate

max_forks_repo_name: TorkelE/BRNEquilibrate


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module TrajSaveMCTS

export MCTSStudy

using AdaptiveStressTesting
using SISLES: GenerativeModel

using CPUTime
using RLESUtils, Obj2Dict, RunCases, Observers, Loggers, Obj2DataFrames
using DataFrames

using ..Config_ACASX_GM
using ..ConfigAST
using ..ConfigMCTS
using ..DefineSave
using ..TrajSaveCommon
using ..DefineLog
using ..SaveTypes
using ..SaveHelpers
using ..PostProcess

import ..DefineSave.trajSave
import ..TrajSaveCommon.get_study_params

type MCTSStudy
    fileroot::String
end

function MCTSStudy(;
    fileroot::AbstractString="trajSaveMCTS")
    MCTSStudy(fileroot)
end

function trajSave(study_params::MCTSStudy,
                  cases::Cases=Cases(Case());
                  outdir::AbstractString="./",
                  postproc::PostProcessing=StandardPostProc())

    println("Starting MCTS Study...")
    pmap(case -> begin
        starttime_us = CPUtime_us()
        startnow = string(now())

        sim_params = defineSimParams(; extract_params(case, "sim_params")...)
        ast_params = defineASTParams(; extract_params(case, "ast_params")...)
        mcts_params = defineMCTSParams(; extract_params(case, "mcts_params")...)
        study_params = extract_params!(study_params, case, "study_params")

        sim = defineSim(sim_params)
        ast = defineAST(sim, ast_params)

        results = stress_test2(ast, mcts_params)

        compute_info = ComputeInfo(startnow, string(now()), gethostname(),
            (CPUtime_us() - starttime_us) / 1e6)

        save_k_logs(sim, mcts_params, ast, results, compute_info,
            sim_params, ast_params, study_params, postproc, outdir)

        results.rewards[1]
    end, cases)

end

function get_study_params(d::TrajLog, ::Type{Val{:MCTS}})
    study = MCTSStudy()
    Obj2DataFrames.set!(study, ObjDataFrame(d["study_params"]))
end

function save_k_logs(sim::Any, mcts_params::DPWParams, ast::AdaptiveStressTest,
    results::StressTestResults, compute_info::ComputeInfo, sim_params,
    ast_params::ASTParams, study_params::MCTSStudy, postproc::PostProcessing,
    outdir::AbstractString)
    for k = 1:mcts_params.top_k
        fileroot_ = "$(study_params.fileroot)_$(sim.string_id)_k$k"
        outfileroot = joinpath(outdir, fileroot_)
        save_mcts_log(mcts_params, ast, results, k, compute_info, sim_params,
            ast_params, study_params, postproc, outfileroot)
    end
end

function save_k_logs(sim::DualSim, mcts_params::DPWParams, ast::AdaptiveStressTest,
    results::StressTestResults, compute_info::ComputeInfo, sim_params,
    ast_params::ASTParams, study_params::MCTSStudy, postproc::PostProcessing,
    outdir::AbstractString)
    #save
    ds_transition_model = ast.transition_model
    ds_sim = ast.sim

    #sim1
    ast.sim = ds_sim.sim1
    ast.transition_model = transition_model(ast)
    for k = 1:mcts_params.top_k
        fileroot_ = "$(study_params.fileroot)_$(ast.sim.string_id)_k$(k)_sim1"
        outfileroot = joinpath(outdir, fileroot_)
        save_mcts_log(mcts_params, ast, results, k, compute_info, sim_params[1],
            ast_params, study_params, postproc, outfileroot; suppress_warn=true)
    end

    #sim2
    ast.sim = ds_sim.sim2
    ast.transition_model = transition_model(ast)
    for k = 1:mcts_params.top_k
        fileroot_ = "$(study_params.fileroot)_$(ast.sim.string_id)_k$(k)_sim2"
        outfileroot = joinpath(outdir, fileroot_)
        save_mcts_log(mcts_params, ast, results, k, compute_info, sim_params[2],
            ast_params, study_params, postproc, outfileroot; suppress_warn=true)
    end

    #restore
    ast.sim = ds_sim
    ast.transition_model = ds_transition_model

    ast
end

function save_mcts_log(mcts_params::DPWParams, ast::AdaptiveStressTest,
    results::StressTestResults, k::Int64, compute_info::ComputeInfo,
    sim_params, ast_params::ASTParams, study_params::MCTSStudy, postproc,
    outfileroot::AbstractString; suppress_warn::Bool=false)

    log = trajLoggedPlay(ast, results.rewards[k], results.action_seqs[k],
        compute_info, sim_params, ast_params; suppress_warn=suppress_warn)

    set_run_type!(log, "MCTS")
    set_mcts_params!(log, mcts_params)
    set_q_values!(log, results.q_values[k])
    set_study_params!(log, study_params)
    set_action_seq!(log, results.action_seqs[k])
    set_q_values!(log, results.q_values[k])

    outfile = trajSave(outfileroot, log)

    #callback for postprocessing
    postprocess(outfile, postproc)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using KiteUtils
using Documenter

DocMeta.setdocmeta!(KiteUtils, :DocTestSetup, :(using KiteUtils); recursive=true)

makedocs(;
    modules=[KiteUtils],
    authors="Uwe Fechner <uwe.fechner.msc@gmail.com> and contributors",
    repo="https://github.com/ufechner7/KiteUtils.jl/blob/{commit}{path}#{line}",
    sitename="KiteUtils.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ufechner7.github.io/KiteUtils.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Exported Functions" => "functions.md",
        "Exported Types" => "types.md",
        "Examples" => "examples.md",
    ],
)

deploydocs(;
    repo="github.com/ufechner7/KiteUtils.jl",
    devbranch="main",
)


text:
"""
Hold state for use with `exec` to determine the objects referenced by symbols.
"""
type State
    mod    :: Module
    scopes :: Vector{Dict}
    refs   :: Vector

    State(mod, scopes = Dict[]) = new(mod, scopes, Any[])
end

# Scope manipulations.

pushscope!(state, scope) = push!(state.scopes, scope)
popscope!(state)         = pop!(state.scopes)

# Ref manipulations.

pushref!(state, ref) = push!(state.refs, ref)
popref!(state)       = pop!(state.refs)

"""
Push reference onto `state`, run function block, and pop reference afterwards.
"""
function withref(fn, state, ref)
    pushref!(state, ref)
    result = fn()
    popref!(state)
    result
end

function getvar(state, var)
    for scope in reverse(state.scopes)
        haskey(scope, var) && return scope[var]
    end
    var
end

"""
Add new variable and it's value to topmost scope.
"""
function addtoscope!(state, var, value)
    isempty(state.scopes) && pushscope!(state, Dict())
    setindex!(state.scopes[end], value, var)
end


text:
module CUDAAtomicUtils
using CUDA
using CUDA: @atomic

export atomicallyAddOne,atomicallyAddToSpot, atomicAdd, atomicMinSet, atomicMaxSet
"""
atomically add to given 1 length array 1
"""
function atomicallyAddOne(arr)
   return  CUDA.@atomic arr[1]+=1
end

"""
atomically add given value to the coordinate (linear) of supplied array
"""
function atomicallyAddToSpot( arr,coord,value)
#    return CUDA.atomic_add!(pointer(arr, coord),typ(value))
return  CUDA.@atomic arr[coord]+=value
end

"""
adds atomically number to target and return old value
"""
function atomicAdd(target, number)
    return  CUDA.@atomic target[1]+=number

end

"""
given array (or some pointer ...)  and value to set it will set target to miniumum of targets value and supplied value and return old value
target - what we want to updated
value - new value if is the minimum
coord - coordinate of a pointer in array
"""
function atomicMinSet( arr,value,coord=1)

    #return CUDA.atomic_add!(pointer(arr, coord),typ(value))
  return  @inbounds CUDA.@atomic arr[coord]=min(arr[coord],value)
end#atomicMinSet

"""
given array (or some pointer ...)  and value to set it will set target to maximum of targets value and supplied value  and return old value
"""
function atomicMaxSet( arr,value,coord=1)
    return @inbounds CUDA.@atomic arr[coord]=max(arr[coord],value)

end#atomicMaxSet

 end#module CUDAAtomicUtils


text:
function checkhash(filename)
    filename |> read |> sha256 |> bytes2hex
end
function verifyfiles(allfilesjson; skip = false)
    allfiles_dict = JSON3.read(read(allfilesjson, String), SortedDict{String, String})
    for (f, v) in allfiles_dict
        if !isfile(f)
            if !skip
                error("$f not found and is required.")
            else
                @warn "$f not found, skipping"
                continue
            end
        end
        fh = checkhash(f)
        if fh == v
            println("Found file $f, hash OK")
        else
            @warn "Hash for file $f is $fh, expected $v"
        end
    end
end


text:
using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

df = DataFrame(CSV.read(joinpath("..", "..", "data", "chimpanzees.csv"), delim=';'))
df[!, :pulled_left] = convert(Array{Int64}, df[:, :pulled_left])
df[!, :prosoc_left] = convert(Array{Int64}, df[:, :prosoc_left])
first(df, 5)

Base.@kwdef mutable struct Chimpanzees{Ty <: AbstractVector,
  Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
    "Number of observations"
    N::Int
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::Chimpanzees)
  as( (β = as(Array, size(model.x, 2)), ) )
end

# Instantiate the model with data and inits.

N = size(df, 1)
x = hcat(ones(Int64, N), df[:, :prosoc_left]);
y = df[:, :pulled_left]
model = Chimpanzees(;y=y, x=x, N=N);

# Make the model callable with a single argument.

function (model::Chimpanzees)(θ)
    @unpack y, x, N = model   # extract the data
    @unpack β = θ  # works on the named tuple too
    ll = 0.0
    ll += sum(logpdf.(Normal(0, 10), β)) # a & bp
    ll += sum([loglikelihood(Binomial(1, logistic(dot(x[i, :], β))), [y[i]]) for i in 1:N])
    ll
end

println()
θ = (β = [1.0, 2.0],)
model(θ)
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

DynamicHMC.Diagnostics.EBFMI(results.tree_statistics)

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)

a3d = Array{Float64, 3}(undef, 1000, 2, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[i].β)
  end
end

# Create MCMCChains object

cnames = ["bp", "bpC"]
sections =   Dict(
  :parameters =>  ["bp", "bpC"]
)
chns = create_mcmcchains(a3d, cnames, sections, start=1)

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
      Mean        SD       Naive SE       MCSE      ESS
 a 0.05103234 0.12579086 0.0019889282 0.0035186307 1000
bp 0.55711212 0.18074275 0.0028577937 0.0040160451 1000

Quantiles:
       2.5%        25.0%       50.0%      75.0%      97.5%
 a -0.19755400 -0.029431425 0.05024655 0.12978825 0.30087758
bp  0.20803447  0.433720250 0.55340400 0.67960975 0.91466915
";

# Summarize results

describe(chns)

# End of `10/m10.2d.jl`


text:
using CombinedParsers
using CombinedParsers.Regexp
using CombinedParserTools
using CombinedParserTools.Tokens
using Test

@testset "CombinedParserTools.jl" begin
    import CombinedParserTools.Tokens: attributes
    @testset "html attributes" begin
        @test parse(attributes,"a = 1 b=6% font=\"+1asd\"") == [T"1"a, T"6%"b, T"+1asd"font]
        @test parse(attributes, "size=10% class=1") == [ T"10%"size, T"1"class ]
    end

    @testset "html" begin
        inner = Either{Any}(Any[!Repeat(CharNotIn("<>"))]);
        pushfirst!(inner,html(!re"[[:alpha:]]+",inner,attributes));
        parse(inner,"ba")
        @test parse(inner,"iboldj") ==
            Node("a",[T"+1"font], ["i",Node("b", [], ["bold"]),"j"])

        @test parse(inner,"") ==
            Node("a", [T"+1"font], [])
    end
    # Write your own tests here.
end


text:
#=
A lot of visualization forms in GLVisualize are realised in the form of instanced
particles. This is because they can be handled very efficiently by OpenGL.
There are quite a few different ways to feed instances with different attributes.
The main constructor for particles is a tuple of (Primitive, Position), whereas
position can come in all forms and shapes. You can leave away the primitive.
In that case, GLVisualize will fill in some default that is anticipated to make
the most sense for the datatype.
=#

#3D primitives
const Primitives3D = Union{AbstractGeometry{3}, AbstractMesh}
#2D primitives AKA sprites, since they are shapes mapped onto a 2D rectangle
const Sprites = Union{AbstractGeometry{2}, Shape, Char, Type}
const AllPrimitives = Union{AbstractGeometry, Shape, Char, AbstractMesh}

"""
We plot simple Geometric primitives as particles with length one.
At some point, this should all be appended to the same particle system to increase
performance.
"""
function _default(
        geometry::TOrSignal{G}, s::Style, data::Dict
    ) where G <: GeometryPrimitive{2}
    data[:offset] = Vec2f0(0)
    _default((geometry, const_lift(x-> Point2f0[minimum(x)], geometry)), s, data)
end

"""
Vectors of floats are treated as barplots, so they get a HyperRectangle as
default primitive.
"""
function _default(main::VectorTypes{T}, s::Style, data::Dict) where T <: AbstractFloat
    _default((centered(HyperRectangle{2, Float32}), main), s, data)
end
"""
Matrices of floats are represented as 3D barplots with cubes as primitive
"""
function _default(main::MatTypes{T}, s::Style, data::Dict) where T <: AbstractFloat
    _default((AABB(Vec3f0(-0.5,-0.5,0), Vec3f0(1.0)), main), s, data)
end
"""
Vectors of n-dimensional points get ndimensional rectangles as default
primitives. (Particles)
"""
function _default(main::VectorTypes{P}, s::Style, data::Dict) where P <: Point
    N = length(P)
    @gen_defaults! data begin
        scale = N == 2 ? Vec2f0(30) : Vec3f0(0.03) # for 2D points we assume they're in pixels
    end
    _default((centered(HyperRectangle{N, Float32}), main), s, data)
end

"""
3D matrices of vectors are 3D vector field with a pyramid (arrow) like default
primitive.
"""
function _default(main::ArrayTypes{T, 3}, s::Style, data::Dict) where T<:Vec
    _default((Pyramid(Point3f0(0, 0, -0.5), 1f0, 1f0), main), s, data)
end
"""
2D matrices of vectors are 2D vector field with a an unicode arrow as the default
primitive.
"""
function _default(main::ArrayTypes{T, 2}, s::Style, data::Dict) where T<:Vec
    _default(('⬆', main), s, data)
end

"""
Vectors with `Vec` as element type are treated as vectors of rotations.
The position is assumed to be implicitely on the grid the vector defines (1D,2D,3D grid)
"""
function _default(
        main::Tuple{P, ArrayTypes{T, N}}, s::Style, data::Dict
    ) where {P <: AllPrimitives, T <: Vec, N}
    primitive, rotation_s = main
    rotation_v = to_value(rotation_s)
    @gen_defaults! data begin
        color_norm = const_lift(extrema2f0, rotation_s)
        ranges = ntuple(i->LinRange(0f0, 1f0, size(rotation_v, i)), N)
    end
    grid = Grid(rotation_v, ranges)

    if N == 1
        scalevec = Vec2f0(step(grid.dims[1]), 1)
    elseif N == 2
        scalevec = Vec2f0(step(grid.dims[1]), step(grid.dims[2]))
    else
        scalevec = Vec3f0(ntuple(i->step(grid.dims[i]), 3)).*Vec3f0(0.4,0.4, 1/to_value(color_norm)[2]*4)
    end
    if P <: Char # we need to preserve proportion of the glyph
        scalevec = Vec2f0(glyph_scale!(primitive, scalevec[1]))
        @gen_defaults! data begin # for chars we need to make sure they're centered
            offset = -scalevec/2f0
        end
    end
    @gen_defaults! data begin
        rotation   = const_lift(vec, rotation_s)
        color_map  = default(Vector{RGBA})
        scale      = scalevec
        color      = nothing
    end
    _default((primitive, grid), s, data)
end

"""
arrays of floats with any geometry primitive, will be spaced out on a grid defined
by `ranges` and will use the floating points as the
height for the primitives (`scale_z`)
"""
function _default(
        main::Tuple{P, ArrayTypes{T,N}}, s::Style, data::Dict
    ) where {P<:AbstractGeometry, T<:AbstractFloat, N}
    primitive, heightfield_s = main
    heightfield = to_value(heightfield_s)
    @gen_defaults! data begin
        ranges = ntuple(i->LinRange(0f0, 1f0, size(heightfield, i)), N)
    end
    grid = Grid(heightfield, ranges)
    @gen_defaults! data begin
        scale = nothing
        scale_x::Float32 = step(grid.dims[1])
        scale_y::Float32 = N == 1 ? 1f0 : step(grid.dims[2])
        scale_z = const_lift(vec, heightfield_s)
        color = nothing
        color_map  = color == nothing ? default(Vector{RGBA}) : nothing
        color_norm = color == nothing ? const_lift(extrema2f0, heightfield_s) : nothing
    end
    _default((primitive, grid), s, data)
end

"""
arrays of floats with the sprite primitive type (2D geometry + picture like),
will be spaced out on a grid defined
by `ranges` and will use the floating points as the
z position for the primitives.
"""
function _default(
        main::Tuple{P, ArrayTypes{T,N}}, s::Style, data::Dict
    ) where {P <: Sprites, T <: AbstractFloat, N}
    primitive, heightfield_s = main
    heightfield = to_value(heightfield_s)
    @gen_defaults! data begin
        ranges = ntuple(i->LinRange(0f0, 1f0, size(heightfield, i)), N)
    end
    grid = Grid(heightfield, ranges)
    @gen_defaults! data begin
        position_z = const_lift(vec, heightfield_s)
        scale      = Vec2f0(step(grid.dims[1]), N>=2 ? step(grid.dims[2]) : 1f0)
        color_map  = default(Vector{RGBA})
        color      = nothing
        color_norm = const_lift(extrema2f0, heightfield_s)
    end
    _default((primitive, grid), s, data)
end

"""
Sprites primitives with a vector of floats are treated as something barplot like
"""
function _default(
        main::Tuple{P, VectorTypes{T}}, s::Style, data::Dict
    ) where {P <: AllPrimitives, T <: AbstractFloat}
    primitive, heightfield_s = main
    heightfield = to_value(heightfield_s)
    @gen_defaults! data begin
        ranges = range(0f0, stop = 1f0, length = length(heightfield))
    end
    grid = Grid(heightfield, ranges)
    delete!(data, :ranges)
    @gen_defaults! data begin
        scale = nothing
        scale_x::Float32 = step(grid.dims[1])
        scale_y = heightfield_s
        scale_z::Float32 = 1f0
    end
    _default((primitive, grid), s, data)
end

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to meshparticle
function _default(
        p::Tuple{TOrSignal{Pr}, VectorTypes{P}}, s::Style, data::Dict
    ) where {Pr <: Primitives3D, P <: Point}
    meshparticle(p, s, data)
end

function _default(
        p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
    ) where {Pr <: Primitives3D, G <: Tuple}
    @gen_defaults! data begin
        primitive = p[1]
        position         = nothing => TextureBuffer
        position_x       = p[2][1] => TextureBuffer
        position_y       = p[2][2] => TextureBuffer
        position_z       = length(p[2]) > 2 ? p[2][3] : 0f0 => TextureBuffer
        instances        = const_lift(length, position_x)
    end
    meshparticle(p, s, data)
end

function _default(
        p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
    ) where {Pr <: Primitives3D, G <: Grid}
    meshparticle(p, s, data)
end

# make conversion of mesh signals work. TODO move to GeometryTypes?
function Base.convert(::Type{T}, mesh::Node) where T<:GeometryTypes.HomogenousMesh
    lift(T, mesh)
end

function to_meshcolor(color::TOrSignal{Vector{T}}) where T <: Colorant
    TextureBuffer(color)
end

function to_meshcolor(color::TOrSignal{Matrix{T}}) where T <: Colorant
    Texture(color)
end
function to_meshcolor(color)
    color
end

function to_mesh(mesh::TOrSignal{<: GeometryPrimitive})
    gl_convert(const_lift(GLNormalMesh, mesh))
end

function to_mesh(mesh::TOrSignal{<: HomogenousMesh})
    gl_convert(to_value(mesh))
end

function orthogonal(v::T) where T <: StaticVector{3}
    x, y, z = abs.(v)
    other = x < y ? (x < z ? unit(T, 1) : unit(T, 3)) : (y < z ? unit(T, 2) : unit(T, 3))
    return cross(v, other)
end

using AbstractPlotting
using AbstractPlotting: get_texture_atlas

vec2quaternion(rotation::StaticVector{4}) = rotation

function vec2quaternion(r::StaticVector{2})
    vec2quaternion(Vec3f0(r[1], r[2], 0))
end
function vec2quaternion(rotation::StaticVector{3})
    AbstractPlotting.rotation_between(Vec3f0(0, 0, 1), Vec3f0(rotation))
end

vec2quaternion(rotation::Vec4f0) = rotation
vec2quaternion(rotation::VectorTypes) = const_lift(x-> vec2quaternion.(x), rotation)
vec2quaternion(rotation::Node) = lift(vec2quaternion, rotation)
vec2quaternion(rotation::AbstractPlotting.Quaternion)= Vec4f0(rotation.data)
"""
This is the main function to assemble particles with a GLNormalMesh as a primitive
"""
function meshparticle(p, s, data)
    rot = get!(data, :rotation, Vec4f0(0, 0, 0, 1))
    rot = vec2quaternion(rot)
    delete!(data, :rotation)
    @gen_defaults! data begin
        primitive = p[1] => to_mesh
        position = p[2] => TextureBuffer
        position_x = nothing => TextureBuffer
        position_y = nothing => TextureBuffer
        position_z = nothing => TextureBuffer

        scale = Vec3f0(1) => TextureBuffer
        scale_x = nothing => TextureBuffer
        scale_y = nothing => TextureBuffer
        scale_z = nothing => TextureBuffer

        rotation = rot => TextureBuffer
        texturecoordinates = nothing
    end

    @gen_defaults! data begin
        color_map  = nothing => Texture
        color_norm = nothing
        intensity  = nothing
        color      = if color_map == nothing
            default(RGBA{Float32}, s)
        else
            nothing
        end => to_meshcolor

        instances = const_lift(length, position)
        shading = true
        shader = GLVisualizeShader(
            "util.vert", "particles.vert", "fragment_output.frag", "standard.frag",
            view = Dict(
                "position_calc" => position_calc(position, position_x, position_y, position_z, TextureBuffer),
                "light_calc" => light_calc(shading)
            )
        )
    end
    if AbstractPlotting.to_value(intensity) != nothing
        if AbstractPlotting.to_value(position) != nothing
            data[:intensity] = intensity_convert_tex(intensity, position)
            data[:len] = const_lift(length, position)
        else
            data[:intensity] = intensity_convert_tex(intensity, position_x)
            data[:len] = const_lift(length, position_x)
        end
    end
    data
end

"""
This is the most primitive particle system, which uses simple points as primitives.
This is supposed to be the fastest way of displaying particles!
"""
_default(position::VectorTypes{T}, s::style"speed", data::Dict) where {T <: Point} = @gen_defaults! data begin
    vertex       = position => GLBuffer
    color_map    = nothing  => Vec2f0
    color        = (color_map == nothing ? default(RGBA{Float32}, s) : nothing) => GLBuffer

    color_norm   = nothing  => Vec2f0
    intensity = nothing  => GLBuffer
    point_size   = 2f0
    prerender    = ()->glPointSize(point_size)
    shader       = GLVisualizeShader("fragment_output.frag", "dots.vert", "dots.frag")
    gl_primitive = GL_POINTS
end

"""
returns the Shape for the distancefield algorithm
"""
primitive_shape(::Char) = DISTANCEFIELD
primitive_shape(x::X) where {X} = primitive_shape(X)
primitive_shape(::Type{T}) where {T <: Circle} = CIRCLE
primitive_shape(::Type{T}) where {T <: SimpleRectangle} = RECTANGLE
primitive_shape(::Type{T}) where {T <: HyperRectangle{2}} = RECTANGLE
primitive_shape(x::Shape) = x

"""
Extracts the scale from a primitive.
"""
primitive_scale(prim::GeometryPrimitive) = Vec2f0(widths(prim))
primitive_scale(::Union{Shape, Char}) = Vec2f0(40)
primitive_scale(c) = Vec2f0(0.1)

"""
Extracts the offset from a primitive.
"""
primitive_offset(x, scale::Nothing) = Vec2f0(0) # default offset
primitive_offset(x, scale) = const_lift(/, scale, -2f0)  # default offset

"""
Extracts the uv offset and width from a primitive.
"""
primitive_uv_offset_width(c::Char) = glyph_uv_width!(c)
primitive_uv_offset_width(x) = Vec4f0(0,0,1,1)

"""
Gets the texture atlas if primitive is a char.
"""
primitive_distancefield(x) = nothing
primitive_distancefield(::Char) = get_texture!(get_texture_atlas())
primitive_distancefield(::Node{Char}) = get_texture!(get_texture_atlas())

function _default(
        p::Tuple{TOrSignal{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: Colorant, P <: Point}
    data[:image] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = lift(x-> Vec2f0(size(x)), p[1])
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end
function _default(
        p::Tuple{TOrSignal{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: AbstractFloat, P <: Point}
    data[:distancefield] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = lift(x-> Vec2f0(size(x)), p[1])
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

function _default(
        p::Tuple{VectorTypes{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: Colorant, P <: Point}
    images = to_value(p[1])
    isempty(images) && error("Can not display empty vector of images as primitive")
    images = sort(images, by=size)
    sizes = map(size, images)
    scale = Vec2f0(first(sizes))
    if !all(x-> x == sizes[1], sizes) # if differently sized
        # create texture atlas
        maxdims = sum(map(Vec{2, Int}, sizes))
        rectangles = map(x->SimpleRectangle(0,0,x...), sizes)
        rpack = RectanglePacker(SimpleRectangle(0, 0, maxdims...))
        uv_coordinates = [push!(rpack, rect).area for rect in rectangles]
        max_xy = mapreduce(maximum, max, uv_coordinates)
        texture_atlas = Texture(C, tuple(max_xy...))
        for (area, img) in zip(uv_coordinates, images)
            texture_atlas[area] = img #transfer to texture atlas
        end
        scale = Vec2f0.(widths.(uv_coordinates))
        data[:uv_offset_width] = map(uv_coordinates) do uv
            mini = minimum(uv) ./ max_xy
            maxi = maximum(uv) ./ max_xy
            Vec4f0(mini..., maxi...)
        end
        images = texture_atlas
    end
    data[:image] = images # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = scale
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to sprites
_default(p::Tuple{TOrSignal{Pr}, VectorTypes{P}}, s::Style, data::Dict) where {Pr <: Sprites, P<:Point} =
    sprites(p,s,data)

_default(p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict) where {Pr <: Sprites, G<:Grid} =
    sprites(p,s,data)

function _default(
            p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
        ) where {Pr <: Sprites, G <: Tuple}
    @gen_defaults! data begin
        shape      = const_lift(primitive_shape, p[1])
        position   = nothing => GLBuffer
        position_x = p[2][1] => GLBuffer
        position_y = p[2][2] => GLBuffer
        position_z = length(p[2]) > 2 ? p[2][3] : 0f0 => GLBuffer
    end
    sprites(p, s, data)
end

"""
Main assemble functions for sprite particles.
Sprites are anything like distance fields, images and simple geometries
"""
function sprites(p, s, data)
    rot = get!(data, :rotation, Vec4f0(0, 0, 0, 1))
    rot = vec2quaternion(rot)
    delete!(data, :rotation)

    @gen_defaults! data begin
        shape       = const_lift(x-> Int32(primitive_shape(x)), p[1])
        position    = p[2]    => GLBuffer
        position_x  = nothing => GLBuffer
        position_y  = nothing => GLBuffer
        position_z  = nothing => GLBuffer

        scale       = const_lift(primitive_scale, p[1]) => GLBuffer
        scale_x     = nothing                => GLBuffer
        scale_y     = nothing                => GLBuffer
        scale_z     = nothing                => GLBuffer

        rotation    = rot => GLBuffer
        image       = nothing => Texture
    end
    # TODO don't make this dependant on some shady type dispatch
    if isa(to_value(p[1]), Char) && !isa(to_value(scale), Vec) # correct dimensions
        scale = const_lift(s-> Vec2f0(glyph_scale!(p[1], s)), scale)
        data[:scale] = scale
    end

    @gen_defaults! data begin
        offset          = primitive_offset(p[1], scale) => GLBuffer
        intensity       = nothing => GLBuffer
        color_map       = nothing => Texture
        color_norm      = nothing
        color           = (color_map == nothing ? default(RGBA, s) : nothing) => GLBuffer

        glow_color      = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_color    = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_width    = 0f0
        glow_width      = 0f0
        uv_offset_width = const_lift(primitive_uv_offset_width, p[1]) => GLBuffer

        distancefield   = primitive_distancefield(p[1]) => Texture
        indices         = const_lift(length, p[2]) => to_index_buffer
        # rotation and billboard don't go along
        billboard        = rotation == Vec4f0(0,0,0,1) => "if `billboard` == true, particles will always face camera"
        preferred_camera = :orthographic_pixel
        fxaa             = false
        shader           = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "sprites.geom",
            "sprites.vert", "distance_shape.frag",
            view = Dict("position_calc"=>position_calc(position, position_x, position_y, position_z, GLBuffer))
        )
        gl_primitive = GL_POINTS
    end
    # Exception for intensity, to make it possible to handle intensity with a
    # different length compared to position. Intensities will be interpolated in that case
    if position != nothing
        data[:intensity] = intensity_convert(intensity, position)
        data[:len] = const_lift(length, position)
    else
        data[:intensity] = intensity_convert(intensity, position_x)
        data[:len] = const_lift(length, position_x)
    end
    data
end

"""
Transforms text into a particle system of sprites, by inferring the
texture coordinates in the texture atlas, widths and positions of the characters.
"""
function _default(main::Tuple{TOrSignal{S}, P}, s::Style, data::Dict) where {S <: AbstractString, P}
    data[:position] = main[2]
    _default(main[1], s, data)
end

import ..AbstractPlotting: to_font, glyph_uv_width!, glyph_scale!
import ..Makie: get_texture!

function _default(main::TOrSignal{S}, s::Style, data::Dict) where S <: AbstractString
    @gen_defaults! data begin
        relative_scale  = 4 #
        start_position  = Point2f0(0)
        atlas           = get_texture_atlas()
        distancefield   = get_texture!(atlas)
        stroke_width    = 0f0
        glow_width      = 0f0
        font            = to_font("default")
        scale_primitive = true
        position        = const_lift(calc_position, main, start_position, relative_scale, font, atlas)
        offset          = const_lift(calc_offset, main, relative_scale, font, atlas)
        prerender       = () -> begin
            glDisable(GL_DEPTH_TEST)
            glDepthMask(GL_TRUE)
            glDisable(GL_CULL_FACE)
            enabletransparency()
        end
        uv_offset_width = const_lift(main) do str
            Vec4f0[glyph_uv_width!(atlas, c, font) for c = str]
        end
        scale = const_lift(main, relative_scale) do str, s
            Vec2f0[glyph_scale!(atlas, c, font, s) for c = str]
        end
    end
    delete!(data, :font)
    _default((DISTANCEFIELD, position), s, data)
end


text:
# "smart" summary and REPL printing

function Base.summary(sys::System)
    return "System (base power $(get_base_power(sys)))"
end

function Base.show(io::IO, sys::System)
    println(io, "base_power=$(get_base_power(sys))")
    show(io, sys.data)
end

function Base.show(io::IO, ::MIME"text/plain", sys::System)
    println(io, "System")
    println(io, "======")
    println(io, "System Units Base: $(get_units_base(sys))")
    println(io, "Base Power: $(get_base_power(sys))")
    println(io, "Base Frequency: $(get_frequency(sys))\n")
    show(io, MIME"text/plain"(), sys.data)
end

function Base.show(io::IO, ::MIME"text/html", sys::System)
    println(io, "
System
")
    println(io, "
Base Power: $(get_base_power(sys))
")
    show(io, MIME"text/html"(), sys.data)
end

function Base.summary(tech::DeviceParameter)
    return "$(typeof(tech))"
end

function Base.show(io::IO, ::MIME"text/plain", data::PowerSystemTableData)
    println(io, "$(typeof(data)):")
    println(io, "  directory:  $(data.directory)")
    if !isnothing(data.timeseries_metadata_file)
        println(io, "  timeseries_metadata_file:  $(data.timeseries_metadata_file)")
    end
    println(io, "  base_power:  $(data.base_power)")
    for (field, df) in data.category_to_df
        print(io, "  $field:  ")
        println(io, "$(summary(df))")
    end
end

function Base.show(io::IO, ist::Component)
    print(io, string(nameof(typeof(ist))), "(")
    is_first = true
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        getter_name = Symbol("get_$name")
        if field_type <: IS.TimeSeriesContainer ||
           field_type <: InfrastructureSystemsInternal
            continue
        elseif hasproperty(PowerSystems, getter_name)
            getter_func = getproperty(PowerSystems, getter_name)
            val = getter_func(ist)
        else
            val = getproperty(ist, name)
        end
        if is_first
            is_first = false
        else
            print(io, ", ")
        end
        print(io, val)
    end
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", ist::Component)
    default_units = false
    if !has_units_setting(ist)
        print(io, "\n")
        @warn(
            "SystemUnitSetting not defined, using NATURAL_UNITS for displaying device specification."
        )
        set_units_setting!(
            ist,
            SystemUnitsSettings(100.0, UNIT_SYSTEM_MAPPING["NATURAL_UNITS"]),
        )
        default_units = true
    end
    try
        print(io, summary(ist), ":")
        for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
            obj = getproperty(ist, name)
            getter_name = Symbol("get_$name")
            if (obj isa InfrastructureSystemsInternal) && !default_units
                print(io, "\n   ")
                show(io, MIME"text/plain"(), obj.units_info)
                continue
            elseif obj isa IS.TimeSeriesContainer ||
                   obj isa InfrastructureSystemsType ||
                   obj isa Vector{<:InfrastructureSystemsComponent}
                val = summary(getproperty(ist, name))
            elseif hasproperty(PowerSystems, getter_name)
                getter_func = getproperty(PowerSystems, getter_name)
                val = getter_func(ist)
            else
                val = getproperty(ist, name)
            end
            # Not allowed to print `nothing`
            if isnothing(val)
                val = "nothing"
            end
            print(io, "\n   ", name, ": ", val)
        end
    finally
        if default_units
            set_units_setting!(ist, nothing)
        end
    end
end


text:
"""
- x0: [gene, P2_gene, RNA, P, P2]
- k1: repression binding rate
- k1r: reverse repression binding rate
- k2: transcription rate
- k3: translation rate
- k4: dimerization rate
- k4r: dissociation rate
- k5: RNA degradation
- k6: protein degradation rate
"""
function neg_autoreg(;x0 = [10, 0, 0, 0, 0], k1=1.0, k1r=10.0, k2=0.01, k3=10.0, k4=1.0, k4r=1.0, k5=0.1, k6=0.01)
  m = Network("auto-regulation")

  m <= Species("gene",    x0[1])
  m <= Species("P2_gene", x0[2])
  m <= Species("RNA",     x0[3])
  m <= Species("P",       x0[4])
  m <= Species("P2",      x0[5])

  m <= Reaction("repression binding",         k1,  "gene + P2 --> P2_gene")
  m <= Reaction("reverse repression binding", k1r, "P2_gene --> gene + P2")
  m <= Reaction("transcription",              k2,  "gene --> gene + RNA")
  m <= Reaction("translation",                k3,  "RNA --> RNA + P")
  m <= Reaction("dimerization",               k4,  "P + P --> P2")
  m <= Reaction("dissociation",               k4r, "P2 --> P + P")
  m <= Reaction("RNA degradation",            k5,  "RNA --> 0")
  m <= Reaction("protein degradation",        k6,  "P --> 0")

  return m
end


text:
### sukkiri_db_sqlite.jl -- SQLite3 implementation for Sukkiri database layer
###   Copyright © 2014 by Matthew C. Gushee <matt@gushee.net>
###   This program is open-source software, released under the BSD license.
###   See the accompanying LICENSE file for details.

module SukkiriDBSqlite
using SQLite

export SukkiriSqliteDB, connect, disconnect, setup_db

include("sukkiri_common_sql.jl")

type SukkiriSqliteDB
    file::String
    connection::Union(Nothing,SQLiteDB)
    SukkiriSqliteDB(filename) = new(filename, nothing)
end

function connect(db::SukkiriSqliteDB)
    if db.connection == nothing
        db.connection = SQLiteDB(db.file)
    end
end

function disconnect(db::SukkiriSqliteDB)
    if db.connection != nothing
        close(db.connection)
    end
    db.connection = nothing
end

### IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
### ----  DATABASE SETUP  --------------------------------------------------

function setup_db(db::SukkiriSqliteDB)
    connect(db)
    execute(db.connection, "BEGIN")
    for q = [create_primitive_table_query, create_string_type_table_query,
             create_number_type_table_query, create_vocab_table_query,
             create_cardinality_table_query, create_struct_type_table_query,
             create_type_class_table_query, create_types_table_query,
             create_union_type_table_query, create_struct_members_table_query,
             create_statement_table_query]
        execute(db.connection, q)
    end
    for qlist = (populate_primitive_table_queries, populate_cardinality_table_queries,
                 populate_type_class_table_queries, populate_types_table_queries,
                 populate_union_type_table_queries)
        for q = qlist
            execute(db.connection, q)
        end
    end
    execute(db.connection, "COMMIT")
    disconnect(db)
end

### OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

end

max_stars_repo_name: mgushee/sukkiri

max_issues_repo_name: mgushee/sukkiri

max_forks_repo_name: mgushee/sukkiri


text:
"""
    smooth(z, sdim)

Smoothen field `z` with a circular kernel of radius `sdim` using convolution
"""
function smooth(zi, sdim)

    println("Smoothin ....")
    iinan = findall(isnan.(zi))
    iinotnan = findall(.~isnan.(zi))
    zi[iinan] .= 0.

    (nsx,nsy) = size(zi)
    #sdim=div(div(ns,nas),2)

    nsx2 = div(nsx, 2)
    nsy2 = div(nsy, 2)
    mask = zeros(nsx, nsy)

    for i = 1:nsx
        for j = 1:nsy
            kx = i-1
            ky = j-1
            if (i>nsx2+1)
                kx = i-nsx-1
            end
            if (j>nsy2+1)
                ky = j-nsy-1
            end
            r2 = kx*kx+ky*ky
            if (r2<=sdim*sdim)
                mask[i,j] = 1.0
            end
        end
    end
    fm = fft(mask)
    zf = real(ifft(fm.*fft(zi)))/sum(mask)
    if length(iinan)>0
        zi1 = deepcopy(zi)
        zi1[iinotnan] .= 1.0
        zf = zf./(real(ifft(fm.*fft(zi1)))/sum(mask))
    end
    zf[iinan] .= NaN

    return zf
end


text:
module TestConsoleProgressMonitor

using ConsoleProgressMonitor
using Logging
using Test

include("demos.jl")

@testset "with_progresslogger" begin
    @test ConsoleProgressMonitor.with_progresslogger(demo1, dt=0) isa Any
    @test ConsoleProgressMonitor.with_progresslogger(demo2, dt=0) isa Any
end

@testset "(un)install_logger" begin
    if get(ENV, "CI", "false") == "true"
        orig = global_logger()
        @test ConsoleProgressMonitor.install_logger() === orig
        @test global_logger() !== orig
        demo1()
        demo2()
        @test ConsoleProgressMonitor.uninstall_logger() isa AbstractLogger
        @test global_logger() === orig
    end
end

end  # module


text:
"""
Make a grid with the given size, for which the given domain is the covering domain.
"""
function discretize_togrid end

discretize_togrid(d::AbstractInterval, n::Int) = EquispacedGrid(n, d)
discretize_togrid(d::Interval{:closed,:open}, n::Int) = PeriodicEquispacedGrid(n, d)
discretize_togrid(d::Interval{:open,:open}, n::Int) = MidpointEquispacedGrid(n, d)
discretize_togrid(d::UnitInterval{T}, n::Int) where {T} = UnitEquispacedGrid{T}(n)
discretize_togrid(d::ChebyshevInterval{T}, n::Int) where {T} = ChebyshevNodes{T}(n)

discretize_togrid(d::ProductDomain, dims) = productgrid(map(discretize_togrid, components(d), dims)...)

discretize_togrid(d::MappedDomain, dims) =
    MappedGrid(forward_map(d), discretize_togrid(superdomain(d), dims))

discretize_togrid(d::UnitCircle{T}, n::Int) where {T} =
    MappedGrid(mapfrom_parameterdomain(d), UnitPeriodicEquispacedGrid{T}(n))


text:
# SQL operators.

mutable struct OperatorClause <: AbstractSQLClause
    name::Symbol
    args::Vector{SQLClause}

    OperatorClause(;
                   name::Union{Symbol, AbstractString},
                   args = SQLClause[]) =
        new(Symbol(name), args)
end

OperatorClause(name; args = SQLClause[]) =
    OperatorClause(name = name, args = args)

OperatorClause(name, args...) =
    OperatorClause(name, args = SQLClause[args...])

"""
    OP(; name, args = [])
    OP(name; args = [])
    OP(name, args...)

An application of a SQL operator.

# Examples

```jldoctest
julia> c = OP("NOT", OP("=", :zip, "60614"));

julia> print(render(c))
(NOT ("zip" = '60614'))
```
"""
OP(args...; kws...) =
    OperatorClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(OP), pats::Vector{Any}) =
    dissect(scr, OperatorClause, pats)

Base.convert(::Type{AbstractSQLClause}, ::typeof(*)) =
    OperatorClause(:*)

PrettyPrinting.quoteof(c::OperatorClause, ctx::QuoteContext) =
    Expr(:call, nameof(OP), string(c.name), quoteof(c.args, ctx)...)


text:
n = 6;
ord = 5;

init_femshop("reordering");

@domain(3)
@functionSpace(LEGENDRE, ord)

@mesh(HEXMESH, n)

@variable(u)
@variable(q1)
@variable(q2)
@variable(q3)
@variable(q4)
@variable(q5)
@variable(q6)

@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)
@boundary(q1, 1, DIRICHLET, 0)
@boundary(q2, 1, DIRICHLET, 0)
@boundary(q3, 1, DIRICHLET, 0)
@boundary(q4, 1, DIRICHLET, 0)
@boundary(q5, 1, DIRICHLET, 0)
@boundary(q6, 1, DIRICHLET, 0)

@coefficient(f, "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

@coefficient(g, "(-3*pi*pi + 5)*sin(pi*x)*sin(pi*y)*sin(pi*z)")
@weakForm([q1,q2,q3,q4,q5,q6], ["-dot(grad(q1), grad(v)) + q2*v + q3*v + q4*v + q5*v + q6*v - g*v",
                                "-dot(grad(q2), grad(v)) + q3*v + q4*v + q5*v + q6*v + q1*v - g*v",
                                "-dot(grad(q3), grad(v)) + q4*v + q5*v + q6*v + q1*v + q2*v - g*v",
                                "-dot(grad(q4), grad(v)) + q5*v + q6*v + q1*v + q2*v + q3*v - g*v",
                                "-dot(grad(q5), grad(v)) + q6*v + q1*v + q2*v + q3*v + q4*v - g*v",
                                "-dot(grad(q6), grad(v)) + q1*v + q2*v + q3*v + q4*v + q5*v - g*v"])

times = 2;
timings = zeros(3);


text:
Random.seed!(223)

@testset "data_collector" begin
  forest = model_initiation(f=0.05, d=0.8, p=0.01, griddims=(20, 20), seed=2);
  agent_properties = [:status, :pos]
  aggregators = [length, count]
  steps_to_collect_data = collect(1:10);
  data = step!(dummy_agent_step, forest_step!, forest, 10, agent_properties, steps_to_collect_data);
  @test size(data)[2] == 21
  colnames = names(data)
  @test colnames[1] == :id
  @test colnames[end] == :pos_10

  agent_properties = [:status]
  forest = model_initiation(f=0.05, d=0.8, p=0.01, griddims=(20, 20), seed=2);
  data = step!(dummy_agent_step, forest_step!, forest, 10, agent_properties, aggregators, steps_to_collect_data);
  @test size(data) == (10,3)
end


text:
using Test

include("fridgeTest.jl")
# add here the file with your unit tests


text:
mutable struct Lorentz
    σ::Float64
    ρ::Float64
    β::Float64
    x::Float64
    y::Float64
    z::Float64
end

function step!(l::Lorentz, dt::Float64)
    dx = l.σ*(l.y - l.x)
    dy = l.x*(l.ρ - l.z) - l.y
    dz = l.x*l.y - l.β*l.z
    l.x += dt*dx
    l.y += dt*dy
    l.z += dt*dz
end

attractor = Lorentz(10.0,28.0,8/3,1,1,1)

using Plots

default(legendfontsize=16, dpi=150)
fig = plot3d(1,xlim=(-20.0,20.0),ylim=(-20.0,20.0),zlim=(0,50.0),
            title="Now this is the nice Lorentz",color="purple",
            marker=2,linestyle=:dot,linewidth=0.5,markerstrokewidth=0.0,
            legend=false)

@gif for i = 1:1500
    step!(attractor,0.01)
    push!(fig,attractor.x,attractor.y,attractor.z)
end every 10

max_stars_repo_name: oist/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

# """
#     AbstractSolverData

# An abstract type for the internal data structure of a dynamic solver.

# ### Implementation (for developers)

# To implement a derived class, one must specialize the following methods. Note
# that in the signatures of most of these methods, the `sd` argument of the new
# type is the last argument..

#   * Constructor - normally it would take a model and a plan together with a
#     final condition type.
#   * `global_RJ(x, exog, sd)` - compute the residual and the Jacobian (return a
#     2-tuple) of the system at the given point `x` and exogenous data `exog`.
#   * `global_R!(res, x, exog, sd)` (optional) - compute the residual in place.
#   * `assign_exog_data!(x, exog, sd)` - this function applies the exogenous
#     constraints as per the plan and the given exogenous data.
#   * `assign_final_condition!(x, exog, sd, ::Val{T})` where T is one of the
#     FCType constants. This method applies the given type of final condition. If
#     T is fcgiven, the values of the final conditions are taken from exog.
#   * `assign_update_step(x, lambda, dx, sd)` - this method must implement x = x +
#     lambda * dx
#   * `update_plan!(sd, model, plan)` - update the solver data to reflect the given
#     simulation plan.

# """
# abstract type AbstractSolverData end

# """
#     global_RJ(x::Array, exog::Array, sd::AbstractSolverData)

# Compute the residual and the Jacobian of the dynamic system of equations at the
# given point `x` using the given exogenous data `exog`. Returns a tuple of two
# elements. The first is a vector of residuals and the second is the Jacobian
# matrix.

# """
# global_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# # export global_RJ

# """
#     global_R!(res::Vector, x::Array, exog::Array, sd::AbstractSolverData)

# Compute the residual of the dynamic system of equations at the given point `x` using the given exogenous data `exog`.
# The residual is updated in place.
# """
# global_R!(res::AbstractArray{Float64}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# # export global_R!

# """
#     assign_exog_data(x::Array, exog::Array, sd::AbstractSolverData)

# Apply the exogenous constraints according to the simulation plan (passed to the
# constructor of sd) with exogenous values taken from `exog`.

# """
# assign_exog_data!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# # export assign_exog_data!

# """
#     assign_final_condition!(x::Array, exog::Array, sd::AbstractSolverData, ::Val{T})

# Apply final conditions of the given type. The type parameter `T` must be one of
# the constants of type `FCType`. If `T` is fcgiven, the final condition values
# are taken from `exog`. If `T` is one of `fclevel` or `fcslope`, then the final
# conditions are taken from the steady state solution in the model instance passed
# to the constructor of `sd`.

# """
# assign_final_condition!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, g::AbstractSolverData, ::Val{T}) where T = error("Not implemented for $(typeof(g)) with final condition $(T).")
# # export assign_final_condition!

# """
#     assign_update_step!(x::Array, lambda::Real, dx::Array, sd::AbstractSolverData)

# Assign the iterative update. Equivalent to `x = x + lambda * dx`.
# """
# function assign_update_step! end
# @inline assign_update_step!(x::AbstractArray{Float64}, lambda::Float64, Δx::AbstractArray{Float64}, ::AbstractSolverData) = BLAS.axpy!(lambda, Δx, x)  # axpy!(a,x,y) <=> y .= a .* x .+ y
# # export assign_update_step!

# """
#     update_plan!(sd, model, plan)

# Update the solver data structure for the given simulation plan. The `model` must
# be the same instance that was used to create `sd`. Typically, the plan should
# span the same range as the one used when `sd` was created, although it might
# have different exogenous/endogenous assignments.

# """
# update_plan!(sd::AbstractSolverData, m::Model, p::Plan) = error("Not implemented for $(typeof(sd)).")
# # export update_plan!


text:
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../"))  # CytofRepFam
include(joinpath(@__DIR__, "ParTemp.jl"))
using Distributions
using LinearAlgebra
import Random

function pretty(x; digits=nothing)
  if digits == nothing
    show(stdout, "text/plain", x)
  else
    show(stdout, "text/plain", round.(x, digits=digits))
  end
  println()
end

Random.seed!(0)

eye(n) = Matrix{Float64}(LinearAlgebra.I, n, n)

function gen_distribution(P)
  S = randn(P, P)
  SST = S * S'
  m = randn(P)
  (m, SST)
end

# Dimensions
P = 10
m, S = gen_distribution(P)
mvn = MvNormal(S)

function gensamples(init::T, propCov, iters) where T
  logprob(x) = logpdf(mvn, x)
  out = T[]
  x = copy(init)
  for i in 1:iters
    x, _ = ParTemp.metropolisMV(x, logprob, propCov)
    append!(out, [copy(x)])
  end
  return out
end

delta(n) = min(n^(-0.5), .01)
compute_acceptance_rate(x) = length(unique(x)) / length(x)

function adapt(init, epochs, batchsize; target_acceptance_rate=.23)
  x = copy(init)
  P = length(init)
  propCov = 0.01 * eye(P) / P
  for i in 1:epochs
    print("\r$(i)/$(epochs)")
    samps = gensamples(x, propCov, batchsize)
    acceptance_rate = compute_acceptance_rate(samps)
    # Adapt s, propCov
    propCov = let
      # See: http://probability.ca/jeff/ftpdir/adaptex.pdf
      Cest = cov(samps) # * (1 / i) + propCov * (i - 1) / i
      C1 = Cest * 5.66 / P + 1e-6 * eye(P)
      C0 = 0.01 * eye(P) / P
      if i <= 2 * P
        C = C0
      else
        C = 0.95 > rand() ? C1 : C0
      end
      C
    end
    x .= copy(samps[end])

    if i == epochs
      println()
      return (samps, propCov)
    end
  end
end

@time samps, propCov = adapt(randn(P), 1000, 500);
println("Dims: $P")
println("Accept rate: $(compute_acceptance_rate(samps))")
K = 5
println("Est. Cov.:")
pretty(cov(samps)[1:K, 1:K], digits=3)
println("True Cov.:")
pretty(cov(mvn)[1:K, 1:K], digits=3)


text:
using CircGeometry, Test

@testset "circle example test" begin
    vf = 0.05
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    radius = 1.5
    center = Point(-0.5,1.0)
    outline = OutlineCircle(radius,center)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(outline,material,between_buffer;log=true)

    write_circ("circle.circ",ps)
    save_image("circle.svg",ps,outline)
end

@testset "rectangle example test" begin
    vf = 0.15
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    p1 = Point(-1.0,0.0)
    p2 = Point(1.0,2.0)
    outline = OutlineRectangle(p1,p2)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(outline,material,between_buffer;log=false)

    translate!(outline,1.0,-1.2)

    write_circ("rectangle.circ",ps)
    save_image("rectangle.svg",ps,outline)
end

@testset "polygon car example test" begin
    vf = 0.25
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    polygon = csv_to_polygon("test-data/car.txt")
    outline = OutlinePolygon(polygon)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(outline,material,between_buffer;log=false)

    translate!(outline,1.0,-1.2)

    write_circ("car.circ",ps)
    save_image("car.svg",ps,outline)
end

@testset "polygon car example test" begin
    vf = 0.35
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    polygon = csv_to_polygon("test-data/car.txt")
    outline = OutlinePolygon(polygon)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(outline,material,between_buffer;log=false)

    translate!(outline,1.0,-1.2)

    write_circ("car.circ",ps)
    save_image("car.svg",ps,outline)
end

# other compute_between_buffer tests
vf = 0.55
n_bodies = 400
material = MaterialParameters(vf,n_bodies)
polygon = csv_to_polygon("test-data/car.txt")
outline = OutlinePolygon(polygon)
between_buffer = compute_between_buffer(outline,material)

vf = 0.65
material = MaterialParameters(vf,n_bodies)
between_buffer = compute_between_buffer(outline,material)

rm("circle.circ")
rm("circle.svg")
rm("rectangle.circ")
rm("rectangle.svg")
rm("car.circ")
rm("car.svg")

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Test conversion from osculating to mean elements using SGP4.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/orbit/mean_elements/sgp4.jl
# =======================================

# Function rv_to_mean_elements_sgp4
# ---------------------------------

@testset "Function rv_to_mean_elements_sgp4" begin
    # We just need to run the SGP4, obtain the osculating elements, convert them
    # to mean elements, and compare with the original TLE.

    tle = tle"""
    AMAZONIA 1
    1 47699U 21015A   21270.48626105 -.00000044  00000-0  19860-2 0  9993
    2 47699  98.4889 344.6059 0001597  74.4244 285.7135 14.40801240 30436"""[1]

    # Generate the osculating elements (TEME).
    orbp   = init_orbit_propagator(Val(:sgp4), tle)
    vr, vv = propagate!(orbp, 0:10:12000)
    vjd    = get_epoch(orbp) .+ (0:10:12000) ./ 86400

    # Obtain the mean elements.
    epoch, xo, ~ = rv_to_mean_elements_sgp4(
        vjd,
        vr,
        vv;
        mean_elements_epoch = :begin,
        print_debug         = false,
        max_iterations      = 1000,
        atol                = 1e-10,
        rtol                = 1e-10
    )

    # Convert some elements to compare with the TLE.
    a₀ = xo[1] / (1000sgp4_gc_wgs84.R0)
    e₀ = xo[2]
    i₀ = rad2deg(xo[3])
    Ω₀ = rad2deg(xo[4])
    ω₀ = rad2deg(xo[5])
    M₀ = rad2deg(f_to_M(e₀, xo[6]))

    # Obtain the mean motion [rev / day].
    n₀ = 720 * sgp4_gc_wgs84.XKE / (sqrt(a₀ * a₀ * a₀) * π)

    # Compare.
    @test epoch == tle.epoch
    @test n₀    ≈  tle.n     atol = 1e-7
    @test e₀    ≈  tle.e     atol = 1e-7
    @test i₀    ≈  tle.i     atol = 1e-5
    @test Ω₀    ≈  tle.Ω     atol = 1e-5
    @test ω₀    ≈  tle.ω     atol = 1e-5
    @test M₀    ≈  tle.M     atol = 1e-5
    @test xo[7] ≈  tle.bstar atol = 1e-6
end

# Function rv_to_tle
# ------------------

@testset "Function rv_to_tle" begin
    # We just need to run the SGP4, obtain the osculating elements, convert them
    # to TLE, and compare with the original TLE.

    tle_input = tle"""
    AMAZONIA 1
    1 47699U 21015A   21270.48626105 -.00000044  00000-0  19860-2 0  9993
    2 47699  98.4889 344.6059 0001597  74.4244 285.7135 14.40801240 30436"""[1]

    # Generate the osculating elements (TEME).
    orbp   = init_orbit_propagator(Val(:sgp4), tle_input)
    vr, vv = propagate!(orbp, 0:10:12000)
    vjd    = get_epoch(orbp) .+ (0:10:12000) ./ 86400

    # Obtain the mean elements.
    tle, ~ = rv_to_tle(
        vjd,
        vr,
        vv;
        elem_set_number     = 999,
        int_designator      = "21015A  ",
        mean_elements_epoch = :begin,
        name                = "AMAZONIA 1",
        rev_num             = 3043,
        sat_num             = 47699,
        print_debug         = false,
        max_iterations      = 1000,
        atol                = 1e-10,
        rtol                = 1e-10
    )

    # Compare.
    @test tle.classification  == tle_input.classification
    @test tle.elem_set_number == tle_input.elem_set_number
    @test tle.epoch           == tle_input.epoch
    @test tle.epoch_year      == tle_input.epoch_year
    @test tle.int_designator  == tle_input.int_designator
    @test tle.name            == tle_input.name
    @test tle.rev_num         == tle_input.rev_num
    @test tle.sat_num         == tle_input.sat_num

    @test tle.bstar     ≈  tle_input.bstar     atol = 1e-6
    @test tle.e         ≈  tle_input.e         atol = 1e-7
    @test tle.epoch_day ≈  tle_input.epoch_day atol = 1e-8
    @test tle.i         ≈  tle_input.i         atol = 1e-4
    @test tle.M         ≈  tle_input.M         atol = 1e-4
    @test tle.n         ≈  tle_input.n         atol = 1e-7
    @test tle.Ω         ≈  tle_input.Ω         atol = 1e-4
    @test tle.ω         ≈  tle_input.ω         atol = 1e-4
end


text:
module CameraManager

export Camera

using CoreExtended
using MatrixMath

using ..EntityManager

const Entity = EntityManager.Typ

type Camera
	id::Symbol
	entity::Entity

	fov::Float32
	ratio::Float32
	depth::Vec2f

  viewport::Vec4f

  position::Vec3f
  rotation::Vec3f #quat
	scaling::Vec3f

  viewMat::Mat4x4f
  projectionMat::Mat4x4f
	modelMat::Mat4x4f

	translateMat::Mat4x4f
	rotationMat::Mat4x4f
	scalingMat::Mat4x4f

	mvpMat::Mat4x4f

	Camera(id=:_) = new(id,nothing, 0,0,Vec2f(0,1),zeros(Vec4f),zeros(Vec3f),zeros(Vec3f),ones(Vec3f),IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f)
end

# preset for manager: includes code here
presetManager(CameraManager,Camera)

function zoom(this::Camera, zf::Float32)
	this.scaling = ones(Vec3f) * (1+zoom*0.01f0)
end

function updateTranslation(this::Camera)
	this.translateMat = MatrixMath.translationmatrix(this.position)
end

function updateRotation(this::Camera)
	this.rotationMat = MatrixMath.computeRotation(this.rotation)
end

function updateScaling(this::Camera)
	this.scalingMat = MatrixMath.scalingmatrix(this.scaling)
end

function updateView(this::Camera)
	updateTranslation(this)
	updateRotation(this)
	updateScaling(this)
	this.viewMat = this.rotationMat * this.translateMat * this.scalingMat
end

function setPerspective(this::Camera, fov::Float32,ratio::Float32,near::Float32,far::Float32) #(60.0f0, Float32(windowsize[1]/windowsize[2]), 0.1f0, 10.0f0)
	this.fov=fov
	this.ratio=ratio
	this.depth=Vec2f(near,far)
end

function updateProjection(this::Camera)
	this.projectionMat = MatrixMath.perspectiveprojection(this.fov, this.ratio, this.depth[1], this.depth[2])
end

function updateMVP(this::Camera)
	this.mvpMat = this.projectionMat * this.viewMat * this.modelMat
end

function updateModel(this::Camera)
	if this.entity == nothing return end
	this.modelMat = MatrixMath.translationmatrix(this.entity.position) * MatrixMath.computeRotation(this.entity.rotation) * MatrixMath.scalingmatrix(this.entity.scaling)
end

function update(this::Camera)
	updateProjection(this)
	updateView(this)
	updateModel(this)
	updateMVP(this)
end

end #CameraManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
couch() = 2
apple() = :orange

begin
    plant() = "happy"
    """ good boy """
    dog() = 3
    @inline inlined() = 5
end

################################################################################

using SimpleTraits
isnice(::Type{Int}) = true
isnice(::Type{Float64}) = false
@traitdef IsNice{X}
@traitdef BelongTogether{X,Y} # traits can have several parameters

@traitimpl IsNice{Int}
@traitimpl BelongTogether{Int,String}

@traitimpl IsNice{X} <- isnice(X)

@traitfn f{X; IsNice{X}}(x::X) = "Very nice!"
@traitfn f{X; !IsNice{X}}(x::X) = "Not so nice!"

################################################################################

@generated generated_mouse(x::String) = :(string(x, " ", $(typeof(x))))
@generated generated_missing(::Any) = 42  # JuliaLang#23809

################################################################################

module Moo
# Unsupported. FIXME? - September'17
export cow
cow() = 1
end


text:
function kpm( A, Vinit, M::Integer )

    nvec = size(Vinit,2)
    n = size(A,1);
    zeta = zeros(M,1);
    mu = zeros(M,1);
    V = zeros(n,M+1) + zeros(n,M+1)*im;

    for l = 1:nvec
        V[:,1] = Vinit[:,l];
        for k = 1:M
            zeta[k] = zeta[k] + real(dot(V[:,1],V[:,k]));
            if ( k == 1 )
                V[:,k+1] = A*V[:,k];
            else
                V[:,k+1] = 2*(A*V[:,k]) - V[:,k-1];
            end
        end
    end

    zeta = zeta ./ nvec;
    for k = 1:M
        kd = 2;
        if( k == 1 ) kd = 1; end
        mu[k] = kd / n / pi * zeta[k];
    end

    p = Polynomials.Polynomial(0.0)
    for k=1:M
        p += chebyshevpoly(k-1) * mu[k]
    end

    return p

end


text:
using Dates
using DataStructures

function headline(title, date, tags)
    tag_page = globvar(:tag_page_path)
    io = IOBuffer()
    write(io, """

$title

    """)
    date == Date(1) || write(io, "
$date
")
    write(io, "")
    for tag in tags
        write(io, """
        #$tag
        """)
    end
    write(io, "
")
    return String(take!(io))
end

function bydate!(pagelist)
    sorter(p) = begin
        pvd = pagevar(p, :date)
        if isnothing(pvd)
            return Date(Dates.unix2datetime(stat(p * ".md").ctime))
        end
        return pvd
    end
    sort!(pagelist, by = sorter, rev = true)
end

function postlist(rpaths)
    io = IOBuffer()
    for post in rpaths
        write(io, "
")
        url = get_url(post)
        title = pagevar(post, :title)
        date = pagevar(post, :date)
        tags = pagevar(post, :tags)
        rss = pagevar(post, :rss_description)
        linktitle = "$title"
        write(io, headline(linktitle, date, tags))
        write(io, """

$rss

        Read more →

        """)
    end
    return String(take!(io))
end

function getpostpaths(path = "posts")
    posts = readdir(path)
    return path .* rstrip.(get_url.(posts), '/')
end

hfun_year() = year(now())

function hfun_headline()
    title = locvar(:title)
    date = locvar(:date)
    tags = locvar(:tags)

    io = IOBuffer()
    write(io, headline(title, date, tags))
    toc = Franklin.hfun_toc(
        string.(locvar.([:mintoclevel, :maxtoclevel]))
    )
    write(io, toc)

    return String(take!(io))
end

function hfun_allposts()
    rpaths = getpostpaths()
    bydate!(rpaths)
    return postlist(rpaths)
end

function hfun_taglist()
    tag = locvar(:fd_tag)
    rpaths = globvar(:fd_tag_pages)[tag]
    bydate!(rpaths)
    return postlist(rpaths)
end

function hfun_tagpage()
    rpaths = getpostpaths()
    tags = Iterators.flatten(pagevar.(rpaths, :tags))
    namesortedtags = collect(SortedDict(counter(tags)))
    countsortedtags = sort(namesortedtags; by = x -> x[2], rev = true)
    count = countsortedtags[1][2]
    io = IOBuffer()
    write(io, """

count	name

$count

    """)
    for (tag, c) in countsortedtags
        if c < count
            write(io, """

$c

            """)
            count = c
        end
        write(io, """
            #$tag
        """)
    end
    write(io, "

")
    return String(take!(io))
end