text:
module scenario
using  ..Types, ..materials
using Distributions, LinearAlgebra
#import generic functions
#import .Ballistics: InternalBallistics!, ExternalBallistics!, Vulnerability!

export createProjectile, createFragment

function createProjectile(calibre::Float64,mass::Float64, velocity::Float64,  position::Float64 ;ev = nothing, spin = nothing, bf = nothing, type=nothing, mat = nothing, yaw=nothing )
            Projectile1D(calibre, type, mat, mass, velocity, ev, position, spin, yaw, bf)
end

function createProjectile(calibre::Float64, mass::Float64, velocity::Array{Float64,1},  position::Array{Float64,1} ;ev = nothing, spin = nothing, bf = nothing, type=nothing, mat = nothing, yaw=nothing, inertia=nothing, tof=nothing )
            Projectile2D(calibre, type, mat, mass, velocity, ev, position, spin, yaw, bf, inertia, tof)
end

function createTarget(size::Float64, position::Float64)
    Target1D(size,position)
end

function createTarget(size::Float64, position::Array{Float64,1}, velocity::Array{Float64,1})
    Target2D(size,position,velocity)
end

function createCanon(tc::Float64, lw::Float64, X2w::Float64)
    Canon(tc, lw, X2w)
end

function createFragment(z::Zone, s::Array{FragShapes,1}, cf::ShapeFactor,proj::Projectile2D,αₑ_bar::Array{Float64,1})
    number = Int64(rand(Normal(z.μ_number,z.σ_number)))
    rads = rand(Uniform(0,360), number)
    axs = rand(Uniform(z.angleBegin, z.angleEnd), number)
    f = Array{Fragment}(undef,number)
    c_bar = proj.velocity/norm(proj.velocity)*cos(norm(αₑ_bar))+αₑ_bar
    y_bar = [0.0, 1.0, 0.0]
    e_bar = cross(c_bar,y_bar)/norm(cross(c_bar,y_bar))
    d_bar = cross(e_bar,c_bar)
    for i=1:number
        velocity = rand(Normal(z.μ_velocity,z.σ_velocity))
        mass = abs(rand(Normal(z.μ_mass,z.σ_mass)))
        shape = rand(s)
        cd = rand(Normal(shape.μ_cd,shape.σ_cd))
        γ = rand(LogNormal(cf.μ, cf.σ))
        Aₚ = γ*(mass/z.density)^(2/3)
        #vx = velocity*cosd(axs[i])*cosd(rads[i])
        #vy = velocity*sind(axs[i])
        #vz = velocity*cosd(axs[i])*sind(rads[i])
        f_velocity = proj.velocity + velocity*c_bar + proj.spin * proj.calibre/2*(-d_bar*sind(rads[i])+e_bar*cosd(rads[i]))
        f_position = proj.position*1
        f[i] = Fragment(rads[i], axs[i], mass, z.density,cd, Aₚ, f_position,f_velocity,0.0)
    end
    return f
end

end


text:
# """
# fetchSHparameters(molecule_name::String, record_name::String, base_path::String)
#
# Get the chemical shift and J-coupling parameters from the locally-stored GISSMO record.
# """
# function fetchSHparameters(molecule_name::String, record_name::String,
#     base_path::String)
#
#     file_name = "$(molecule_name)_$(record_name)"
#
#     css_sys, J_IDs_sys, J_vals_sys, cs_singlets,
#     J_lb_sys, J_ub_sys, cs_lb_sys, cs_ub_sys,
#     cs_LUT, reference_concentration, solute_concentration,
#     _ = GISSMOReader.loadGISSMOmolecule(base_path, file_name)
#
#     #N_groups = length(cs_group)
#     p_cs_sys = cstopcs(css_sys, cs_LUT)
#
#     N_spins_singlet = Vector{Int}(undef, length(cs_singlets))
#     cs_singlets_compact = Vector{Float64}(undef, length(cs_singlets))
#     for m = 1:length(cs_singlets)
#         cs_singlets_compact[m] = cs_singlets[m][1]
#         N_spins_singlet[m] = length(cs_singlets[m])
#     end
#
#     N_sys = length(css_sys)
#     @assert length(p_cs_sys) == N_sys
#
#     N_spins_sys = collect( length(css_sys[m]) for m = 1:length(css_sys) )
#     intermediates_sys = prepcouplingalgorithm(N_spins_sys)
#
#     # prepare initial guess and bounds.
#     #tmp = collect( [cs_group[i]; (J_lb_group[i] + J_ub_group[i]) ./ 2] for i = 1:N_groups )
#     css = combinevectors(p_cs_sys)
#     push!(css, cs_singlets_compact...)
#
#     N_p_cs_sys = collect( length(p_cs_sys[i]) for i = 1:N_sys )
#     p0 = copy(css)
#
#     # could be used for dim-reduced p.
#     p_compound = Vector{Vector{Float64}}(undef, N_sys + length(cs_singlets_compact))
#     for i = 1:N_sys
#         #p_compound[i] = p_cs_sys[i]
#         p_compound[i] = copy(p_cs_sys[i])
#     end
#     for i = 1:length(cs_singlets_compact)
#         p_compound[N_sys+i] = [ cs_singlets_compact[i] ]
#     end
#
#     # the full cs version of p_compound.
#     cs_compound = Vector{Vector{Float64}}(undef, N_sys + length(cs_singlets))
#     for i = 1:N_sys
#         cs_compound[i] = copy(p_cs_sys[i])
#     end
#     for i = 1:length(cs_singlets)
#         cs_compound[N_sys+i] = copy(cs_singlets[i])
#     end
#
#     return p_cs_sys, css_sys, cs_singlets, J_vals_sys, J_IDs_sys, intermediates_sys,
#     cs_LUT, cs_singlets_compact, N_spins_singlet, css, p_compound, cs_compound
# end

ext: jl


text:
using OrdinaryDiffEq

# https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/390
solve( ODEProblem( (x,p,t)->-x, 1.0, (0.0, 50.0) ), Rosenbrock23(autodiff=false) )


text:
module ShowLint

using Dates
using DocStringExtensions
using Formatting
using Serialization

export
    apply,
    Pattern,
    Repo,
    project_root,
    ansi2html,
    clone_repositories,
    target_dir,
    patterns,
    repositories

include("patterns.jl")
include("repositories.jl")

const project_root = dirname(dirname(pathof(ShowLint)))
const configs_dir = joinpath(project_root, "configs")

"""
Remove ANSI color codes from `text`.

$(TYPEDSIGNATURES)
"""
function ansi2html(text::String)::String
  escape_codes = Dict(
    "0;30" => "inherit",
    "0;31" => "red",
    "0;32" => "green",
    "0;34" => "blue",
    "0;1" => "inherit",
    "0;100;30" => "inherit",
    "0;41;30" => "red",
    "0;42;30" => "green",
    "0;43;30" => "orange",
    "0;2" => "gray"
  )
  for code in keys(escape_codes)
    before = "[" * code * "m"
    color = escape_codes[code]
    after = ""
    text = replace(text, before => after)
  end
  no_color = "[0m"
  text = replace(text, no_color => "")
  text = replace(text, r"WARNING: [^\n]*$" => "")
  text = replace(text, "/home/runner/clones/github/" => "")
  text = replace(text, "" => "")
  text = replace(text, "<!" => "<!")
  text = replace(text, "[0;7;2m" => "")
  text = text[1:end-6]
  text = strip(text)
  return length(text) < 6 ?
        "" :
        """
        ~~~

$text

        ~~~
        """
end

const clones_dir = joinpath(homedir(), "clones")

"""
Directory of the host of `repo`.

$(TYPEDSIGNATURES)

### Example
```
julia> host_dir("https://github.com")
"github"
```
"""
function host_dir(repo::Repo)::String
    host = repo.host
    host_dir = startswith(host, "https://") ? host[9:end] : host
    host_dir = lowercase(host_dir) == "github.com" ? "github" : host_dir
end

function target_dir(repo::Repo)
    name = repo.name
    joinpath(clones_dir, host_dir(repo), name)
end

function page_path(repo::Repo)
    name = repo.name
    lowercase(joinpath(host_dir(repo), "$name.md"))
end

function clone_repositories(; production=is_production())
    if !isdir(clones_dir)
        mkdir(clones_dir)
    end
    for repo in repositories(; production)
        println("Cloning $(repo.name)")
        name = repo.name
        host = repo.host
        dest = target_dir(repo)
        if isdir(dest)
            cd(dest)
            try
                run(`git pull --ff-only`)
            catch
                @warn "git pull failed for $(repo.name)"
            end
        else
            run(`git clone
                --depth=1
                --no-tags
                --single-branch
                $host/$name $dest`)
        end
        println()
    end
end

function config_path(pat::Pattern)::String
    id = pat.id
    name = "p$id"
    joinpath(configs_dir, "$name.toml")
end

"""
Create a Comby configuration file for `pat` at `dir`.

$(TYPEDSIGNATURES)
"""
function write_config(pat::Pattern, dir::String)
    open(config_path(pat), "w") do io
        write(io, pat.toml)
    end
end

"""
Write Comby configurations from `src/patterns.jl` to `project_root/configs`.

$(TYPEDSIGNATURES)
"""
function write_configs()
    rm(configs_dir; force=true, recursive=true)
    mkpath(configs_dir)
    for pat in patterns
        write_config(pat, configs_dir)
    end
end

"""
Obtains the number of matches from Comby's `-stats` output
on stderr.

$(TYPEDSIGNATURES)
"""
function number_of_matches(stderr::String)::Int
    regex = r"\"number_of_matches\": ([0-9]*)"
    m = match(regex, stderr)
    n_str = m[1]
    parse(Int, n_str)
end

"""
Apply `pattern` to `repo`.

$(TYPEDSIGNATURES)
"""
function apply(pat::Pattern, repo::Repo;
    file_extensions="jl", in_place=false)

    start = Dates.now()
    repo_path = target_dir(repo)

    exclude_prefixes = join(repo.exclude_prefixes, ',')
    stats_path = joinpath(tempdir(), "comby_stats.log")

    cmd = !in_place ?
        `comby
            -exclude-dir $exclude_prefixes
            -config $(config_path(pat))
            -directory $repo_path
            -file-extensions $file_extensions
            -stats
        ` :
        `comby
            -exclude-dir $exclude_prefixes
            -config $(config_path(pat))
            -directory $repo_path
            -file-extensions $file_extensions
            -in-place
        `

    stdout = IOBuffer()
    stderr = IOBuffer()
    run(pipeline(cmd; stdout, stderr))
    out = String(take!(stdout))
    err = String(take!(stderr))
    function avoid_franklin_parse_errors(out)
        out = replace(out, '`' => "`")
        out = replace(out, '"' => """)
        out
    end
    if !in_place
        out = avoid_franklin_parse_errors(out)
        if !isnothing(out)
            out = ansi2html(out)
        end
    end
    time = Dates.now() - start
    println("$(pat.id): $time")
    return (out, err)
end

"""
Apply `pattern` to `code`.

$(TYPEDSIGNATURES)
"""
function apply(pat::Pattern, code::String; file_extension="jl")
    test_dir = joinpath(clones_dir, "tmp", "tmp")
    rm(test_dir; recursive=true, force=true)
    mkpath(test_dir)

    code_file = joinpath(test_dir, "code.$file_extension")
    open(code_file, "w") do io
        write(io, code)
    end

    repo = Repo("tmp", "tmp"; exclude=["nah"])

    apply(pat, repo; in_place=true)

    read(code_file, String)
end

function filter_patterns(repo::Repo, patterns::Vector{Pattern})
    predicates = repo.tags_predicates
    all_predicates_hold(tags) = all([P(tags) for P in predicates])
    filtered_patterns = filter(pat -> all_predicates_hold(pat.tags), patterns)
end

function repo_page(repo::Repo)
    predicates_str = join(repo.tags_predicates, " && ")
    headers = []
    repo_url = joinpath(repo.host, repo.name)
    head = """
        ~~~

$(repo.name)

        ~~~

        Showing patterns for which the `tags` satisfy:
        *$predicates_str*, and the pattern resulted in at least one match.

        \\toc

        """

    function pattern_section(pat)
        diff, err = ShowLint.apply(pat, repo)
        title = pat.title
        id = pat.id
        n_matches = number_of_matches(err)
        n_text = n_matches == 1 ? "hit" : "hits"
        unicode_em_space = " "
        header = "$title$(unicode_em_space)➤$(unicode_em_space)$n_matches $n_text"

        if diff == ""
            return ""
        else
            push!(headers, header)
            return """
            ### $header
            [Pattern #$id](/patterns/#$id)

            $diff

            """
        end
    end
    filtered_patterns = filter_patterns(repo, patterns)
    sections = pattern_section.(filtered_patterns)

    return (
        headers = headers,
        text = join([head, sections...], '\n')
    )
end

"""
Create one webpage per repository.
We could process all the diffs when this function is called or when `serve` runs.
It seems more flexible to do it as early as possible.

$(TYPEDSIGNATURES)
"""
function create_repo_pages(; production=is_production())
    write_configs()
    pages_headers = []

    for repo in repositories(; production)
        if !isdir(host_dir(repo))
            mkdir(host_dir(repo))
        end
        franklin_file = page_path(repo)
        println("Creating repository page at $franklin_file")
        if !isdir(dirname(franklin_file))
            mkdir(dirname(franklin_file))
        end

        page_headers, page_text = repo_page(repo)
        push!(pages_headers, page_headers)

        open(franklin_file, "w") do io
            write(io, """
                +++
                title = "$(repo.name)"
                +++

                [//]: # (Generated file. Do not modify.)

                $page_text
                """
            )
        end
    end
    headers_path = joinpath(project_root, "__site", "pages-headers.txt")
    mkpath(dirname(headers_path))
    Serialization.serialize(headers_path, pages_headers)
    pages_headers
end

function count_without_comments(path)::Int
    content = read(path, String)
    lines = split(content, '\n')
    is_comment_line(line) = startswith(lstrip(line), '#')
    non_comment_lines = filter(!is_comment_line, lines)
    length(non_comment_lines)
end

"""
Counts the number of lines of codes in files with `extension` in
`dir`.

$(TYPEDSIGNATURES)
"""
function cloned_loc(; extension=".jl")::Int
    counts = []

    for (root, dirs, files) in walkdir(clones_dir)
        for file in files
            ext = last(splitext(file))
            if ext == extension
                path = joinpath(root, file)
                count = count_without_comments(path)
                push!(counts, count)
            end
        end
    end
    sum(counts)
end

prettify_loc(n::Number)::String = format(n, commas=true)

"""
Helper function to build the site locally.
`Franklin.serve` can be called after this function, but it is also possible to have Franklin already running.

$(TYPEDSIGNATURES)
"""
function build_site(; production=true)
    clone_repositories()
    create_repo_pages()
end

end # module

max_stars_repo_name: rikhuijzer/display-lint

max_issues_repo_name: rikhuijzer/display-lint

max_forks_repo_name: rikhuijzer/ShowLint


text:
using Logging
using Random
using LinearAlgebra
using StatsBase
using Distributions
using CubicSplines
using LightGraphs
using Plots
using DelimitedFiles
using GraphPlot
using Compose
using Cairo
using HypothesisTests
using Clustering

include("./logger.jl")
ConsoleLogger() |> logger.ElapsedTimeLoggerDecorator |> global_logger
findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2] # find the nearest absolute value `t` in `A`

# general parameters

stim_start = 1 # ms
stim_end_time = 80 # ms
inhib_stim_end_time = 80 # ms
tau_exc = 10 # ms (E current decay time constant)
tau_inh_global = 2 # ms (global I current decay time constant)
tau_inh_local = 2 # ms (local I current decay time constant)
E_to_GI_conn_prob = 0.1 # connection probability of E neurons to randomly synapse to global I neurons
GI_to_E_conn_prob = 0.5 # connection probability of global I neurons to randomly synapse to E neurons

phi_x_train = [-0.015,      0, 0.025, 0.05, 0.075,  0.1, 0.15]
phi_y_train = [     0,  0.005, 0.033, 0.05, 0.06, 0.068, 0.08]
phi_spline = CubicSpline(phi_x_train, phi_y_train)
phi_curve = phi_spline[range(phi_x_train[1], stop=phi_x_train[end], length=1_000)]
phi_curve[phi_curve.<0.0] .= 0.0
xs = range(phi_x_train[1], stop=phi_x_train[end], length=1_000)

# graph generator functions

function generate_triK5_graph()
	triK5_adjacency_matrix = readdlm("triK5_adjacency_matrix.csv", ',', Int)
	g = SimpleGraph(triK5_adjacency_matrix)
	return g
end

function generate_karate_club_graph()
	karate_edges = readdlm("karate_edges.csv", ',', Int)
	g = SimpleGraph(maximum(karate_edges))
	for i in 1:size(karate_edges)[1]
		add_edge!(g, karate_edges[i,1], karate_edges[i,2])
	end
	return g
end

function generate_tutte_graph()
	return smallgraph(:tutte)
end

function generate_gridworld_graph()
	g_grid = LightGraphs.grid([5,5])
	h_grid = LightGraphs.grid([5,5])
	gh_grid = blockdiag(g_grid, h_grid)
	i_grid = LightGraphs.grid([5,5])
	j_grid = LightGraphs.grid([5,5])
	ij_grid = blockdiag(i_grid, j_grid)
	ghij_grid = blockdiag(gh_grid, ij_grid)
	add_edge!(ghij_grid, 23, 28)
	add_edge!(ghij_grid, 15, 61)
	add_edge!(ghij_grid, 73, 78)
	add_edge!(ghij_grid, 86, 40)
	return ghij_grid
end

# weight connection functions

function create_XX_graphical_grinasty_connection_matrix(NE, E_patterns, g)
	E_assembly_length = length(findall(isequal(1), E_patterns))
	n_assemblies = Int(length(findall(>(0), E_patterns))/E_assembly_length)
	@assert length(vertices(g)) == n_assemblies
	J = zeros(NE,NE)
	# auto-association
	for a in 1:n_assemblies
		J[findall(isequal(a), E_patterns),findall(isequal(a), E_patterns)] .+= 1
	end
	# hetero-association
	for v in vertices(g)
		for n in neighbors(g, v)
			J[findall(isequal(v), E_patterns),findall(isequal(n), E_patterns)] .+= 1
			J[findall(isequal(n), E_patterns),findall(isequal(v), E_patterns)] .+= 1
		end
	end
	J .-= Diagonal(J) # no self-edges/autopases
	return J
end

function create_XX_graphical_amit_connection_matrix(NE, E_patterns, g, a)
	E_assembly_length = length(findall(isequal(1), E_patterns))
	n_assemblies = Int(length(findall(>(0), E_patterns))/E_assembly_length)
	@assert length(vertices(g)) == n_assemblies
	J = zeros(NE,NE)
	# auto-association
	for a in 1:n_assemblies
		J[findall(isequal(a), E_patterns),findall(isequal(a), E_patterns)] .= 1
	end
	# hetero-association
	for v in vertices(g)
		for n in neighbors(g, v)
			J[findall(isequal(v), E_patterns),findall(isequal(n), E_patterns)] .= a
			J[findall(isequal(n), E_patterns),findall(isequal(v), E_patterns)] .= a
		end
	end
	J .-= Diagonal(J) # no self-edges/autopases
	return J
end

function create_XX_amit_connection_matrix(patterns, a)
	patterns *= 1.0
	pattern_step_idx = circshift(patterns, 1)
	J_same = transpose(patterns)*patterns
	J_same[J_same.>0] .= 1.0
	J_ahead = transpose(pattern_step_idx)*patterns
	J_ahead[J_ahead.>0] .= a
	J_behind = transpose(patterns)*pattern_step_idx
	J_behind[J_behind.>0] .= a
	J = max.(J_same, max.(J_ahead, J_behind))
    J .-= Diagonal(J) # no self-edges/autopases
	return J
end

function create_XX_grinasty_connection_matrix(patterns)
	patterns *= 1.0
	pattern_step_idx = circshift(patterns, 1)
	J_same = transpose(patterns)*patterns
	J_ahead = transpose(pattern_step_idx)*patterns
	J_behind = transpose(patterns)*pattern_step_idx
	J = J_same + J_ahead + J_behind
    J .-= Diagonal(J) # no self-edges/autopases
	return J
end

function create_graphical_local_EI_IE_matrices(E_patterns, I_assemblies, NE, NLI)
	E_pattern_length = length(findall(isequal(1), E_patterns))
	I_assembly_length = length(findall(isequal(1), I_assemblies))
	n_patterns = Int(length(E_patterns)/E_pattern_length)

	E_patterns_idxs = zeros(Int64, E_pattern_length, length(E_patterns))
	I_assemblies_idxs = zeros(Int64, I_assembly_length, length(E_patterns))

	for pattern = 1:n_patterns
		E_patterns_idxs[:,pattern] = findall(isequal(pattern), E_patterns)
		I_assemblies_idxs[:,pattern] = findall(isequal(pattern), I_assemblies)
	end

	J_EI = zeros(NE, NLI)

	for pattern = 1:n_patterns
		J_EI[E_patterns_idxs[:,pattern],I_assemblies_idxs[:,pattern]] = ones(E_pattern_length, I_assembly_length)
	end

	J_IE = deepcopy(J_EI)'

	return J_EI, J_IE
end

function create_local_EI_IE_matrices(E_patterns, I_assemblies)
	E_pattern_length = length(findall(isone.(E_patterns[1,:])))
	I_assembly_length = length(findall(isone.(I_assemblies[1,:])))
	n_patterns = size(E_patterns)[1]

	E_patterns_idxs = zeros(Int64, E_pattern_length,size(E_patterns)[1])
	I_assemblies_idxs = zeros(Int64, I_assembly_length,size(E_patterns)[1])

	for pattern = 1:n_patterns
		E_patterns_idxs[:,pattern] = findall(isone.(E_patterns[pattern,:]))
		I_assemblies_idxs[:,pattern] = findall(isone.(I_assemblies[pattern,:]))
	end

	J_EI = zeros(size(E_patterns)[2], size(I_assemblies)[2])

	for pattern = 2:n_patterns-1
		J_EI[E_patterns_idxs[:,pattern],I_assemblies_idxs[:,pattern]] = ones(E_pattern_length, I_assembly_length)
	end

	pattern = 1
	J_EI[E_patterns_idxs[:,pattern],I_assemblies_idxs[:,pattern]] = ones(E_pattern_length, I_assembly_length)
	pattern = n_patterns
	J_EI[E_patterns_idxs[:,pattern],I_assemblies_idxs[:,pattern]] = ones(E_pattern_length, I_assembly_length)

	J_IE = deepcopy(J_EI)'

	return J_EI, J_IE
end

# neural response functions

function phi(x)
	if x < -0.05
		return 0
	elseif x > 0.15
		return 0.08
	else
		return phi_curve[findnearest(xs, x)]
	end
end

function psi(x)
	return max(0, (x-0.05)*0.1)
end

# main simulation function

function simulate(graphical, graph_name, grinasty, gamma, c, alpha_var, J_EE, J_ELI, J_LIE, J_GIE, J_EGI, E_patterns, I_assemblies, NE, NLI, NGI, sim_length, stim_target, f, target_dishib=false)
	E_current = zeros(NE)
	E_rate = zeros(NE)
	LI_current = zeros(NLI)
	LI_rate = zeros(NLI)
	GI_current = zeros(NGI)
	GI_rate = zeros(NGI)

	E_rates = zeros(Int(sim_length/0.1), NE)
	LI_rates = zeros(Int(sim_length/0.1), NLI)
	GI_rates = zeros(Int(sim_length/0.1), NGI)

	# These factors are to give an 'average' response population-to-population
	# In Amit et al. 1994:
	# - E to E factor given is in eq. 1: 1/(f*N), where f*N=40
	# - E to GI factor given is in eq. 5: 1/N_inh, where N_inh=f*N=40
	# Notice: in Amit et al. 1994: N_exc/p=40
	# These factors were hand-tuned and the tuning procedure is not described.
	# However, the given factors seem to rely on a constant size of E assemblies (40 neurons).
	# currently, these factors rely on p=100
	# In the grinasty case, we multiply f by 1.5 due to Amit rule discounting (1*40)+(0.5*40)+(0.5*40)=80 c.f. Grinasty rule (1*40)+(1*40)+(1*40)=120
	# this multiplication seems to make Amit and Grinasty rules more comparable with the same values of c
	# the fact this multiplication works suggests the original factors are based on the maximum possible sum of indegree weights
	if graphical
		if graph_name == "triK5"
			if grinasty
				E_to_E_factor = 1/(NE*f*2.5) # triK5 has degree 4, so max input is 40*4 (Grinasty heteroassociation) + 40 (autoassociation) = 200, and the Amit rule max is 120, so 200/80=2.5
			else
				E_to_E_factor = 1/(NE*f*1.5) # triK5 has degree 4, so max input is 20*4 (Amit heteroassociation) + 40 (autoassociation) = 120
			end
		elseif graph_name == "karate"
			if grinasty
				E_to_E_factor = 1/(NE*f*2.794) # karate club has avg degree 4.58824, so max input is 40*4.58824 (Grinasty heteroassociation) + 40 (autoassociation) = 223.5296, and original Amit model max is 80, so 223.5296/80=2.794
			else
				E_to_E_factor = 1/(NE*f*1.647) # karate club has avg degree 4.58824, so max input is 20*4.58824 (Amit heteroassociation) + 40 (autoassociation) = 131.764798, 131.764798/80=1.647
			end
		elseif graph_name == "tutte"
			if grinasty
				E_to_E_factor = 1/(NE*f*2) # the tutte graph has degree 3, so max input is 40*3 (Grinasty heteroassociation) + 40 (autoassociation) = 160, and original Amit model max is 80, so 160/80=2
			else
				E_to_E_factor = 1/(NE*f*1.25) # the tutte graph has degree, so max input is 20*3 (Amit heteroassociation) + 40 (autoassociation) = 100, 100/80=1.25
			end
		elseif graph_name == "gridworld"
			if grinasty
				E_to_E_factor = 1/(NE*f*2.14) # the gridworld graph has avg degree 3.28, so max input is 40*3.28 (Grinasty heteroassociation) + 40 (autoassociation) = 171.2, and original Amit model max is 80, so 171.2/80=2.14
			else
				E_to_E_factor = 1/(NE*f*1.32) # the gridworld graph has degree, so max input is 20*3.28 (Amit heteroassociation) + 40 (autoassociation) = 105.6, 105.6/80=1.32
			end
		end
	else
		if grinasty
			E_to_E_factor = 1/(NE*f*1.5)
		else
			E_to_E_factor = 1/(NE*f)
		end
	end
	E_to_GI_factor = 1/(NE*f*E_to_GI_conn_prob)
	E_to_LI_factor = 1/(NLI*f)
	GI_to_E_factor = 1/(NGI*GI_to_E_conn_prob)
	LI_to_E_factor = 1/(NLI*f)

	# we use `gamma` to scale the overall level of inhibition (for testing only -- in the paper, gamma=1)
	# we use `c` to balance between global (c=0) and local (c=1) inhibition
	# we use `alpha_var` to further rescale global inhibition independently of local inhibition (for testing only -- in the paper, alpha=1)
	local_IE_weighting = gamma*c*LI_to_E_factor
	global_IE_weighting = gamma*(1-c*alpha_var)*GI_to_E_factor

	for t = 1:Int(sim_length/0.1)
		# uncomment for progress printing
		# if( (floor(Int, t*100.0))%10000==0 )
		# 	@info "Simulating timestep" t
		# end
		LI_H = zeros(NLI)
		if Int(stim_start/0.1) <= t <= Int(stim_end_time/0.1)
			if graphical
				E_H = zeros(NE)
				E_H[findall(isequal(stim_target), E_patterns)] .= 0.2
			else
				E_H = E_patterns[stim_target,:].*(0.2)
			end
		else
			E_H = zeros(NE)
		end

		# neural dyanmics
		E_current += (-E_current + E_to_E_factor*(J_EE*E_rate) .- local_IE_weighting*(J_LIE'*LI_rate) .- global_IE_weighting*(J_GIE'*GI_rate) + E_H)/tau_exc
		E_rate = phi.(E_current) .+ abs.(rand(Normal(0.0, 0.00015), NE))
		LI_current += (-LI_current .+ E_to_LI_factor*(J_ELI'*E_rate) + LI_H)/tau_inh_local
		LI_rate = psi.(LI_current)
		GI_current += (-GI_current .+ E_to_GI_factor*(J_EGI'*E_rate))/tau_inh_global
		GI_rate = psi.(GI_current)

		E_rates[t, :] = E_rate
		LI_rates[t, :] = LI_rate
		GI_rates[t, :] = GI_rate
	end

	return E_rates, LI_rates, GI_rates
end

# main function for a single simulation

function main(gamma, c, alpha_var, rebalanced, grinasty, graphical, graph_name, sim_length=250, p=100::Int, NE=4_000::Int, NLI=500::Int, NGI=500::Int, a=0.5::Float64, f=0.01::Float64, target_dishib=false)
	seed = round(Int, time())
	@info "Random" seed
	Random.seed!(seed)

	@info "Creating patterns"
	if graphical
		if graph_name == "triK5"
			g = generate_triK5_graph()
		elseif graph_name == "karate"
			g = generate_karate_club_graph()
		elseif graph_name == "tutte"
			g = generate_tutte_graph()
		elseif graph_name == "gridworld"
			g = generate_gridworld_graph()
		end

		p = length(vertices(g))
		E_patterns_length = Int(f*NE)
		E_patterns = vcat(fill.(collect(1:p),E_patterns_length)...)

		I_assembly_length = Int(f*NLI)
		I_assemblies = vcat(fill.(collect(1:p),I_assembly_length)...)

		@info "Creating connection matrices"
		if grinasty
			J_EE = create_XX_graphical_grinasty_connection_matrix(NE, E_patterns, g)
		else
			J_EE = create_XX_graphical_amit_connection_matrix(NE, E_patterns, g, a)
		end
		J_ELI, J_LIE = create_graphical_local_EI_IE_matrices(E_patterns, I_assemblies, NE, NLI)
	else
		E_patterns = zeros(p, NE)
		E_pattern_length = Int(f*NE)
		for E_pattern = 1:p
			E_patterns[E_pattern, randperm(NE)[1:E_pattern_length]] .= 1
		end
		I_assemblies = zeros(p, NLI)
		I_assembly_length = Int(f*NLI)
		for I_assembly = 1:p
			I_assemblies[I_assembly, randperm(NLI)[1:I_assembly_length]] .= 1
		end

		@info "Creating connection matrices"
		if grinasty
			J_EE = create_XX_grinasty_connection_matrix(E_patterns)
		else
			J_EE = create_XX_amit_connection_matrix(E_patterns, a)
		end
		J_ELI, J_LIE = create_local_EI_IE_matrices(E_patterns, I_assemblies)
		g = nothing
	end

	J_GIE = rand(Categorical([(1-GI_to_E_conn_prob), GI_to_E_conn_prob]), NGI, NE) .- 1
	J_EGI = rand(Categorical([(1-E_to_GI_conn_prob), E_to_GI_conn_prob]), NE, NGI) .- 1

	if rebalanced
		# reweight the global and local I such that it is proportionally more on each E neuron which receives proprtionally more E
		sum_of_inputs = sum(J_EE, dims=1)
		mean_nonzero = sum(J_EE)/length(sum_of_inputs[sum_of_inputs.>0])
		weighted_sum_of_inputs = sum_of_inputs/mean_nonzero
		weighted_sum_of_inputs = dropdims(weighted_sum_of_inputs, dims=1)
		J_LIE = weighted_sum_of_inputs'.*J_LIE
		J_GIE = weighted_sum_of_inputs'.*J_GIE
	end

	@info "Starting simulation"
	E_rates, LI_rates, GI_rates = simulate(graphical, graph_name, grinasty, gamma, c, alpha_var, J_EE, J_ELI, J_LIE, J_GIE, J_EGI, E_patterns, I_assemblies, NE, NLI, NGI, sim_length, 1, f, target_dishib)
	@info "Simulation complete"
	return g, E_rates, LI_rates, GI_rates, J_EE, J_ELI, J_LIE, E_patterns, I_assemblies, seed
end

# plotting functions

function make_simple_figure()
	stim_neurons = findall(>(0), E_patterns[1,:])
	neighbour1 = findall(>(0), E_patterns[2,:])
	neighbour2 = findall(>(0), E_patterns[3,:])
	neighbour3 = findall(>(0), E_patterns[4,:])
	other_neurons = findall(<(1), E_patterns[1,:])

	stimed_rates = E_rates[:,stim_neurons]
	neighbour1_rates = E_rates[:,neighbour1]
	neighbour2_rates = E_rates[:,neighbour2]
	neighbour3_rates = E_rates[:,neighbour3]
	other_rates = E_rates[:,other_neurons]

	mean_stim_rates = mean(stimed_rates, dims=2)
	std_stim_rates = std(stimed_rates, dims=2)

	mean_neighbour1_rates = mean(neighbour1_rates, dims=2)
	std_neighbour1_rates = std(neighbour1_rates, dims=2)

	mean_neighbour2_rates = mean(neighbour2_rates, dims=2)
	std_neighbour2_rates = std(neighbour2_rates, dims=2)

	mean_neighbour3_rates = mean(neighbour3_rates, dims=2)
	std_neighbour3_rates = std(neighbour3_rates, dims=2)

	mean_other_rates = mean(other_rates, dims=2)
	std_other_rates = mean(other_rates, dims=2)

	plot([mean_stim_rates mean_stim_rates], grid=true, ylims=[0, 0.09], label=["Stimulated pattern" ""],
	    fillrange=[mean_stim_rates+std_stim_rates mean_stim_rates-std_stim_rates], fillalpha=0.3, c=:green)
	plot!([mean_neighbour1_rates mean_neighbour1_rates], grid=true, ylims=[0, 0.09], label=["1st neighbour" ""],
		fillrange=[mean_neighbour1_rates+std_neighbour1_rates mean_neighbour1_rates-std_neighbour1_rates], fillalpha=0.3, c=:orange4)
	plot!([mean_neighbour2_rates mean_neighbour2_rates], grid=true, ylims=[0, 0.09], label=["2nd neighbour" ""],
		fillrange=[mean_neighbour2_rates+std_neighbour2_rates mean_neighbour2_rates-std_neighbour2_rates], fillalpha=0.3, c=:orange3)
	plot!([mean_neighbour3_rates mean_neighbour3_rates], grid=true, ylims=[0, 0.09], label=["3rd neighbour" ""],
		fillrange=[mean_neighbour3_rates+std_neighbour3_rates mean_neighbour3_rates-std_neighbour3_rates], fillalpha=0.3, c=:orange2)
	plot!([mean_other_rates mean_other_rates], grid=true, ylims=[0, 0.09], label=["Other Patterns" ""],
	    fillrange=[mean_other_rates+std_other_rates mean_other_rates-std_other_rates], fillalpha=0.3, c=:orange,
	    xlabel="Time (ms)", ylabel=string("Spiking rate (mean ± S.D.)"), legend=:bottomright, fmt=:pdf)
end

function make_simple_inhib_figure()
	stim_neurons = findall(isequal(1), I_assemblies)
	neighbour1 = findall(isequal(2), I_assemblies)
	neighbour2 = findall(isequal(3), I_assemblies)
	neighbour3 = findall(isequal(4), I_assemblies)
	other_neurons = findall(>(4), I_assemblies)

	stimed_rates = LI_rates[:,stim_neurons]
	neighbour1_rates = LI_rates[:,neighbour1]
	neighbour2_rates = LI_rates[:,neighbour2]
	neighbour3_rates = LI_rates[:,neighbour3]
	other_rates = LI_rates[:,other_neurons]

	mean_stim_rates = mean(stimed_rates, dims=2)
	std_stim_rates = std(stimed_rates, dims=2)

	mean_neighbour1_rates = mean(neighbour1_rates, dims=2)
	std_neighbour1_rates = std(neighbour1_rates, dims=2)

	mean_neighbour2_rates = mean(neighbour2_rates, dims=2)
	std_neighbour2_rates = std(neighbour2_rates, dims=2)

	mean_neighbour3_rates = mean(neighbour3_rates, dims=2)
	std_neighbour3_rates = std(neighbour3_rates, dims=2)

	mean_other_rates = mean(other_rates, dims=2)
	std_other_rates = mean(other_rates, dims=2)

	plot([mean_stim_rates mean_stim_rates], grid=true, ylims=[0, 0.09], label=["Stimulated pattern" ""],
	    fillrange=[mean_stim_rates+std_stim_rates mean_stim_rates-std_stim_rates], fillalpha=0.3, c=:green)
	plot!([mean_neighbour1_rates mean_neighbour1_rates], grid=true, ylims=[0, 0.09], label=["1st neighbour" ""],
		fillrange=[mean_neighbour1_rates+std_neighbour1_rates mean_neighbour1_rates-std_neighbour1_rates], fillalpha=0.3, c=:orange4)
	plot!([mean_neighbour2_rates mean_neighbour2_rates], grid=true, ylims=[0, 0.09], label=["2nd neighbour" ""],
		fillrange=[mean_neighbour2_rates+std_neighbour2_rates mean_neighbour2_rates-std_neighbour2_rates], fillalpha=0.3, c=:orange3)
	plot!([mean_neighbour3_rates mean_neighbour3_rates], grid=true, ylims=[0, 0.09], label=["3rd neighbour" ""],
		fillrange=[mean_neighbour3_rates+std_neighbour3_rates mean_neighbour3_rates-std_neighbour3_rates], fillalpha=0.3, c=:orange2)
	plot!([mean_other_rates mean_other_rates], grid=true, ylims=[0, 0.09], label=["Other Patterns" ""],
	    fillrange=[mean_other_rates+std_other_rates mean_other_rates-std_other_rates], fillalpha=0.3, c=:orange,
	    xlabel="Time (ms)", ylabel=string("Spiking rate (mean ± S.D.)"), legend=:bottomright, fmt=:pdf)
end

function make_graph_figure(save_string, savefig_bool, graph_name)
	n_vertices = length(vertices(g))
	vertex_activities = zeros(n_vertices)
	for v in vertices(g)
		vertex_activities[v] = mean(E_rates[end-20:end,(1+(v-1)*40):(v*40)])
	end

	nodelabel = 1:nv(g)
	membership = vertex_activities./0.09
	nodecolor = cgrad(:grays, rev=true)
	# membership color
	nodefillc = nodecolor[membership.*256]

	if (graph_name == "triK5")
		locs_x = zeros(length(collect(vertices(g))))
		locs_y = zeros(length(collect(vertices(g))))
		locs_y[1:5] .+= 0.5
		locs_y[6:15] .-= 0.5
		locs_x[6:10] .-= 0.5
		locs_x[11:15] .+= 0.5

		# community A
		locs_x[2] += 0.2
		locs_x[3] -= 0.2
		locs_x[4] += 0.15
		locs_x[5] -= 0.15
		locs_y[1] += 0.8
		locs_y[2:3] .+= 0.5
		locs_y[4:5] .+= 0.1

		# community B
		locs_x[9] += 0.2
		locs_x[8] -= 0.2
		locs_x[10] += 0.15
		locs_x[7] -= 0.15
		locs_y[6] += 0.8
		locs_y[9] += 0.5
		locs_y[8] += 0.5
		locs_y[10] += 0.1
		locs_y[7] += 0.1

		# community C
		locs_x[12] += 0.2
		locs_x[15] -= 0.2
		locs_x[13] += 0.15
		locs_x[14] -= 0.15
		locs_y[11] += 0.8
		locs_y[12] += 0.5
		locs_y[15] += 0.5
		locs_y[13] += 0.1
		locs_y[14] += 0.1

		draw(PDF(save_string), gplot(g, locs_x, locs_y, nodefillc=nodefillc))
	elseif (graph_name == "karate")
		karate_plot_xy = readdlm("karate_plot_xy.txt", '\t', Float64)
		karate_plot_x = karate_plot_xy[1,:]
		karate_plot_y = karate_plot_xy[2,:]
		draw(PDF(save_string), gplot(g, karate_plot_x, karate_plot_y, nodefillc=nodefillc))
	elseif (graph_name == "tutte")
		tutte_plot_xy = readdlm("tutte_plot_xy.txt", '\t', Float64)
		tutte_plot_x = tutte_plot_xy[1,:]
		tutte_plot_y = tutte_plot_xy[2,:]
		draw(PDF(save_string), gplot(g, tutte_plot_x, tutte_plot_y, nodefillc=nodefillc))
	elseif (graph_name == "gridworld")
		locs_x = zeros(100)
		locs_y = zeros(100)
		for v in 1:25
			locs_x[v] = mod(v-1, 5)/4
			locs_y[v] = floor((v/5)+1)/4
			if (v % 5) == 0
				locs_y[v] -= 0.25
			end
		end
		locs_x[1:25] .-= 1
		locs_y[1:25] .-= 1

		for v in 26:50
			locs_x[v] = mod(v-1, 5)/4
			locs_y[v] = floor((v/5)+1)/4
			if (v % 5) == 0
				locs_y[v] -= 0.25
			end
		end
		locs_x[26:50] .-= 1
		locs_y[26:50] .-= 0.5

		for v in 51:75
			locs_x[v] = mod(v-1, 5)/4
			locs_y[v] = floor((v/5)+1)/4
			if (v % 5) == 0
				locs_y[v] -= 0.25
			end
		end
		locs_x[51:75] .+= 0.5
		locs_y[51:75] .-= 3.5

		for v in 76:100
			locs_x[v] = mod(v-1, 5)/4
			locs_y[v] = floor((v/5)+1)/4
			if (v % 5) == 0
				locs_y[v] -= 0.25
			end
		end
		locs_x[76:100] .+= 0.5
		locs_y[76:100] .-= 3

		draw(PDF(save_string), gplot(g, locs_x, locs_y, nodefillc=nodefillc))
	else
		if savefig_bool
			draw(PDF(save_string), gplot(g, nodefillc=nodefillc, nodelabel=nodelabel))
		else
			gplot(g, layout=circular_layout, nodefillc=nodefillc, nodelabel=nodelabel)
		end
	end
end

# example simulations and plots

simlen = 500
c = 0.7
f = 0.02
p = 50
g, E_rates, LI_rates, GI_rates, J_EE, J_ELI, J_LIE, E_patterns, I_assemblies, seed = main(1, c, 1.0, true, true, false, "", simlen, p, 4_000, 500, 500, 0, f, false)
make_simple_figure()

savefig_bool = true
simlen = 500
graph_name = "gridworld"
c = 0.525
g, E_rates, LI_rates, GI_rates, J_EE, J_ELI, J_LIE, E_patterns, I_assemblies, seed = main(1, c, 1.0, true, true, true, graph_name, simlen, 100, 4_000, 500, 500, 0, 0.01, false)
save_string = string("seed",seed,"_c",c,"_graph",graph_name,"_simlen",simlen,"_grinasty_normd E rates on vertices.pdf")
plot(GI_rates,legend=false, xlabel="Time (ms)", ylabel=string("Spiking rate"), color=:grey, ylims=[0,0.35])
savefig(string("seed",seed,"_global I rates.pdf"))
plot(LI_rates,legend=false, xlabel="Time (ms)", ylabel=string("Spiking rate"), color=:grey, ylims=[0,0.065])
savefig(string("seed",seed,"_local I rates.pdf"))
make_graph_figure(save_string, savefig_bool, graph_name)

# main correlations simulations function

function main_correlations(gamma, c, alpha_var, rebalanced, grinasty, graphical, graph_name, sim_length=250, p=100::Int, NE=4_000::Int, NLI=500::Int, NGI=500::Int, a=0.5::Float64, f=0.01::Float64, target_dishib=false)
	seed = round(Int, time())
	@info "Random" seed
	Random.seed!(seed)

	@info "Creating patterns"
	if graphical
		if graph_name == "triK5"
			g = generate_triK5_graph()
		elseif graph_name == "karate"
			g = generate_karate_club_graph()
		elseif graph_name == "tutte"
			g = generate_tutte_graph()
		elseif graph_name == "gridworld"
			g = generate_gridworld_graph()
		end

		p = length(vertices(g))
		E_patterns_length = Int(f*NE)
		E_patterns = vcat(fill.(collect(1:p),E_patterns_length)...)

		I_assembly_length = Int(f*NLI)
		I_assemblies = vcat(fill.(collect(1:p),I_assembly_length)...)

		@info "Creating connection matrices"
		if grinasty
			J_EE = create_XX_graphical_grinasty_connection_matrix(NE, E_patterns, g)
		else
			J_EE = create_XX_graphical_amit_connection_matrix(NE, E_patterns, g, a)
		end
		J_ELI, J_LIE = create_graphical_local_EI_IE_matrices(E_patterns, I_assemblies, NE, NLI)
	else
		E_patterns = zeros(p, NE)
		E_pattern_length = Int(f*NE)
		for E_pattern = 1:p
			E_patterns[E_pattern, randperm(NE)[1:E_pattern_length]] .= 1
		end
		I_assemblies = zeros(p, NLI)
		I_assembly_length = Int(f*NLI)
		for I_assembly = 1:p
			I_assemblies[I_assembly, randperm(NLI)[1:I_assembly_length]] .= 1
		end

		@info "Creating connection matrices"
		if grinasty
			J_EE = create_XX_grinasty_connection_matrix(E_patterns)
		else
			J_EE = create_XX_amit_connection_matrix(E_patterns, a)
		end
		J_ELI, J_LIE = create_local_EI_IE_matrices(E_patterns, I_assemblies)
		g = nothing
	end

	J_GIE = rand(Categorical([(1-GI_to_E_conn_prob), GI_to_E_conn_prob]), NGI, NE) .- 1
	J_EGI = rand(Categorical([(1-E_to_GI_conn_prob), E_to_GI_conn_prob]), NE, NGI) .- 1

	if rebalanced
		# reweight the global and local I such that it is proportional for each E neuron with its incoming E
		sum_of_inputs = sum(J_EE, dims=1)
		mean_nonzero = sum(J_EE)/length(sum_of_inputs[sum_of_inputs.>0])
		weighted_sum_of_inputs = sum_of_inputs/mean_nonzero
		weighted_sum_of_inputs = dropdims(weighted_sum_of_inputs, dims=1)
		J_LIE = weighted_sum_of_inputs'.*J_LIE
		J_GIE = weighted_sum_of_inputs'.*J_GIE
	end

	@info "Starting simulations"
	max_rep = 1
	rep_attractors = zeros(max_rep, p, NE)
	for rep = 1:max_rep
		@info "Starting simulation repetition" rep
		for stim_id = 1:p
			E_rates, LI_rates, GI_rates = E_rates, LI_rates, GI_rates = simulate(graphical, graph_name, grinasty, gamma, c, alpha_var, J_EE, J_ELI, J_LIE, J_GIE, J_EGI, E_patterns, I_assemblies, NE, NLI, NGI, sim_length, stim_id, f, target_dishib)
			rep_attractors[rep,stim_id,:] = mean(E_rates[end-Int(19/0.1):end,:],dims=1)
		end
	end
	attractors = dropdims(mean(rep_attractors,dims=1),dims=1)

	@info "Removing 'non-selective' cells"
	maximum_FRs_in_attractors_all_neurons = maximum(attractors,dims=1)
	histogram([maximum_FRs_in_attractors_all_neurons...], legend=false, xlabel="Maximum firing rate in attractor state of any attractor", ylabel="Number of neurons", fmt=:pdf)
	if graphical
		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_histogram_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_histogram_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_histogram_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_histogram_unbalanced.pdf"))
			end
		end
	else
		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_histogram_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_histogram_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_histogram_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_maximum_FRs_in_attractors_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_histogram_unbalanced.pdf"))
			end
		end
	end
	sig_units_refs = findall(>(0.02), maximum_FRs_in_attractors_all_neurons)
	sig_units_IDs = LinearIndices(maximum_FRs_in_attractors_all_neurons)[sig_units_refs]
	sig_units_FRs = maximum_FRs_in_attractors_all_neurons[sig_units_IDs]
	sig_units_attractors = attractors[:,sig_units_IDs]

	@info "Calculating correlations"
	# all neurons
	correlations_array = zeros(p,p)
	attractor_array_analysed = attractors
	N_s = length(attractor_array_analysed)/p
	for mu = 1:p
		for ve = 1:p
			V_mu_bar = 1/N_s * sum(attractor_array_analysed[mu,:]) # "mean"
			V_mu_bar_sq = 1/N_s * sum(attractor_array_analysed[mu,:].^2) - V_mu_bar^2 # "variance"
			V_ve_bar = 1/N_s * sum(attractor_array_analysed[ve,:])
			V_ve_bar_sq = 1/N_s * sum(attractor_array_analysed[ve,:].^2) - V_ve_bar^2
			COV_mu_ve = 1/N_s * sum(attractor_array_analysed[mu,:].*attractor_array_analysed[ve,:]) - V_mu_bar*V_ve_bar # "covariance"
			C_mu_ve = COV_mu_ve / sqrt(V_mu_bar_sq*V_ve_bar_sq) # "correlation"
			correlations_array[mu,ve] = C_mu_ve
		end
	end
	writedlm(string("seed",seed,"correlations_array.csv"), correlations_array, ", ")

	if graphical
		heatgray = cgrad(:greys, rev=true)
		if (graph_name == "karate") && (c == 0.1)
			clustering_result = kmeans(correlations_array,5)
			idx = sortperm(assignments(clustering_result))
			members = correlations_array[:,idx]
			members = members[idx,:]
			plot(heatmap(members, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		elseif (graph_name == "karate") && (c == 0.525)
			clustering_result = kmeans(correlations_array,3)
			idx = sortperm(assignments(clustering_result))
			members = correlations_array[:,idx]
			members = members[idx,:]
			plot(heatmap(members, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		elseif (graph_name == "tutte") && (c == 0.1)
			clustering_result = kmeans(correlations_array,3)
			idx = sortperm(assignments(clustering_result))
			members = correlations_array[:,idx]
			members = members[idx,:]
			plot(heatmap(members, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		elseif (graph_name == "tutte") && (c == 0.525)
			clustering_result = kmeans(correlations_array,4)
			idx = sortperm(assignments(clustering_result))
			members = correlations_array[:,idx]
			members = members[idx,:]
			plot(heatmap(members, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		else
			plot(heatmap(correlations_array, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		end

		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_unbalanced.pdf"))
			end
		end

		Q = 0
		g_labels, _ = label_propagation(g)
		for mu in 1:p
			for v in 1:p
				if (mu != v)
					lpa = g_labels[mu]==g_labels[v] ? 1 : -1
					Q += lpa*correlations_array[mu, v]
				end
			end
		end
		Q /= ((p*p)-p)
		open(string("seed",seed,"_Q.txt"),"a") do io
		   println(io,"Q ",Q)
		end

	else
		attractor_correlations = zeros(p,p)
		for attractor = 1:p
			attractor_correlations[attractor,:] = diag(circshift(correlations_array, attractor-1))
		end
		plot(collect(0:Int(p/2)-1),attractor_correlations[1:Int(p/2),:],color="grey",linewidth=1)
		plot!(collect(1:Int(p/2)),reverse(attractor_correlations[Int(p/2)+1:p,:],dims=1),color="grey",linewidth=1)
		plot!(collect(0:Int(p/2)),mean(attractor_correlations[1:Int(p/2)+1,:]',dims=1)',color="black",linewidth=4,legend=false,
				xlabel="Attractor distance",ylabel="Attractor correlation (all neurons)")
		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_all_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_unbalanced.pdf"))
			end
		end
	end
	# 'selective' neurons
	correlations_array = zeros(p,p)
	attractor_array_analysed = sig_units_attractors
	N_s = length(attractor_array_analysed)/p
	for mu = 1:p
		for ve = 1:p
			V_mu_bar = 1/N_s * sum(attractor_array_analysed[mu,:]) # "mean"
			V_mu_bar_sq = 1/N_s * sum(attractor_array_analysed[mu,:].^2) - V_mu_bar^2 # "variance"
			V_ve_bar = 1/N_s * sum(attractor_array_analysed[ve,:])
			V_ve_bar_sq = 1/N_s * sum(attractor_array_analysed[ve,:].^2) - V_ve_bar^2
			COV_mu_ve = 1/N_s * sum(attractor_array_analysed[mu,:].*attractor_array_analysed[ve,:]) - V_mu_bar*V_ve_bar # "covariance"
			C_mu_ve = COV_mu_ve / sqrt(V_mu_bar_sq*V_ve_bar_sq) # "correlation"
			correlations_array[mu,ve] = C_mu_ve
		end
	end
	writedlm(string("seed",seed,"correlations_array_selective.csv"), correlations_array, ", ")

	if graphical
		heatgray = cgrad(:greys, rev=true)
		plot(heatmap(correlations_array, xlabel="Vertex",ylabel="Vertex",color = heatgray, clim=(-1,1)))
		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_graph",graph_name,"_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_unbalanced.pdf"))
			end
		end

		Q_selective = 0
		g_labels, _ = label_propagation(g)
		for mu in 1:p
			for v in 1:p
				if (mu != v)
					lpa = g_labels[mu]==g_labels[v] ? 1 : -1
					Q_selective += lpa*correlations_array[mu, v]
				end
			end
		end
		Q_selective /= ((p*p)-p)
		open(string("seed",seed,"_Q_selective.txt"),"a") do io
		   println(io,"Q_selective ",Q_selective)
		end
	else
		attractor_correlations = zeros(p,p)
		for attractor = 1:p
			attractor_correlations[attractor,:] = diag(circshift(correlations_array, attractor-1))
		end
		plot(collect(0:Int(p/2)-1),attractor_correlations[1:Int(p/2),:],color="grey",linewidth=1)
		plot!(collect(1:Int(p/2)),reverse(attractor_correlations[Int(p/2)+1:p,:],dims=1),color="grey",linewidth=1)
		plot!(collect(0:Int(p/2)),mean(attractor_correlations[1:Int(p/2)+1,:]',dims=1)',color="black",linewidth=4,legend=false,
				xlabel="Attractor distance",ylabel="Attractor correlation (all neurons)")
		if rebalanced
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_rebalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_rebalanced.pdf"))
			end
		else
			if grinasty
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_grinasty_unbalanced.pdf"))
			else
				savefig(string("seed",seed,"_attractor_correlations_over_distance_selective_neurons_c",c,"_gamma",gamma,"_alpha",alpha_var,"_NE",NE,"_NGI",NGI,"_NLI",NLI,"_a",a,"_unbalanced.pdf"))
			end
		end
	end
end

# examples (warning: compute time is linear in the number of vertices of the memory graph)

main_correlations(1, 0.1, 1.0, true, true, true, "triK5", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.525, 1.0, true, true, true, "triK5", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.1, 1.0, true, true, true, "tutte", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.525, 1.0, true, true, true, "tutte", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.1, 1.0, true, true, true, "karate", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.525, 1.0, true, true, true, "karate", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.1, 1.0, true, true, true, "gridworld", 500, 0, 4000, 500, 500, 0.0, 0.01, false)
main_correlations(1, 0.525, 1.0, true, true, true, "gridworld", 500, 0, 4000, 500, 500, 0.0, 0.01, false)

max_stars_repo_name: tfburns/BurnsHagaFukai

max_issues_repo_name: tfburns/BurnsHagaFukai

max_forks_repo_name: tfburns/BurnsHagaFukai


text:
## Sorting

issorted(x::RLEVector) = issorted(x.runvalues)

reverse(x::RLEVector) = RLEVector(reverse(x.runvalues), cumsum(reverse(widths(x))))

function reverse!(x::RLEVector)
    reverse!(x.runvalues)
    x.runends[:] = cumsum(reverse(widths(x)))
    return x
end

function permute_runs(x::RLEVector, indices)
    RLEVector(rvalue(x)[indices], cumsum(widths(x)[indices]))
end

function sort(x::RLEVector{T1,T2}) where {T1,T2}
    ord = sortperm(x.runvalues)
    rle = RLEVector{T1,T2}(x.runvalues[ord], cumsum(widths(x)[ord]))
    return (rle)
end

function sort!(x::RLEVector)
    ord = sortperm(x.runvalues)
    x.runvalues[:] = x.runvalues[ord]
    x.runends[:] = cumsum(widths(x)[ord])
    return (x)
end

function sortperm(x::RLEVector)
    ord = sortperm(x.runvalues)
    RLEVector(ord, cumsum(widths(x)[ord]))
end


text:
"""
    This package contains both, the script that will produce results
    and the common tools used for produced.
    The module itself only contains the common tools
"""
module Chemostat_Rath2017

    using ProjAssistant
    @gen_top_proj

    import Chemostat
    const Ch = Chemostat
    import Chemostat.MetNets
    import Chemostat.MetEP
    import Chemostat.MetLP

    include("RathData/RathData.jl")
    const Rd = RathData

    include("Human1/Human1.jl")
    const H1 = Human1
    const HG = H1.HumanGEM

    include("Utils/Utils.jl")

    function __init__()
        @create_proj_dirs
    end

end # module


text:
#=
Function for basic statistically-relevant properties, e.g. mutual information
fit errors, etc.
=#
using Statistics
export meanstd, standard
meanstd(x) = mean(x), std(x)
standard(x) = (x .- mean(x)) ./ std(x)

#########################################################################
# Basic statistical properties
#########################################################################
export yearly_variability
using Dates: AbstractDateTime

"""
    yearly_variability(t, x) → dates, res
Decompose some vector `x` to its yearly variability with respect to time vector `t`,
i.e. to all the values of `x` that correspond to a unique day and month.
Return `dates`, which are the unique dates (tuples (day, month)) and `res`,
where `res[i]` are all the values of `x` that correspond to date `dates[i]`.
"""
function yearly_variability(t, xs)
    time = Array(t)
    dates = unique(daymonth.(time))
    res = [Float64[] for u in dates]
    for (t, x) in zip(time, xs)
        i = findfirst(isequal(daymonth(t)), dates)
        push!(res[i], x)
    end
    return dates, res
end

using Statistics
export periodic_correlation
"""
	periodic_correlation(x, y, lags, demean=true)
Calculate the correlation function of periodic signals both of which have period
their length (e.g. functions of longitude).
"""
function periodic_correlation(x,y,τ, demean=true)
	if length(unique(x)) == length(unique(y)) == 1
		# Dividing by `std` if `x` is constant leads to `Inf` results.
		return ones(eltype(x), length(τ))
	end
    if demean
        x = x .- mean(x)
        y = y .- mean(y)
    end
    N = length(x)
    r = [sum(x[n]*y[mod1(n+m, N)] for n in 1:N)/N for m in τ]
    r ./ (std(x)*std(y)) # follow normalization of `StatsBase.crosscor`
end

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
default_version() = VersionNumber(VERSION.major)

"""
    Template(; kwargs...) -> Template

Records common information used to generate a package. If you don't wish to manually
create a template, you can use [`interactive_template`](@ref) instead.

# Keyword Arguments
* `user::AbstractString=""`: GitHub (or other code hosting service) username. If left
  unset, it will take the the global git config's value (`github.user`). If that is not
  set, an `ArgumentError` is thrown. **This is case-sensitive for some plugins, so take
  care to enter it correctly.**
* `host::AbstractString="github.com"`: URL to the code hosting service where your package
  will reside. Note that while hosts other than GitHub won't cause errors, they are not
  officially supported and they will cause certain plugins will produce incorrect output.
* `license::AbstractString="MIT"`: Name of the package license. If an empty string is
  given, no license is created. [`available_licenses`](@ref) can be used to list all
  available licenses, and [`show_license`](@ref) can be used to print out a particular
  license's text.
* `authors::Union{AbstractString, Vector{<:AbstractString}}=""`: Names that appear on the
  license. Supply a string for one author or an array for multiple. Similarly to `user`,
  it will take the value of of the global git config's value if it is left unset.
* `dir::AbstractString=$(replace(Pkg.devdir(), homedir() => "~"))`: Directory in which the
  package will go. Relative paths are converted to absolute ones at template creation time.
* `julia_version::VersionNumber=$(default_version())`: Minimum allowed Julia version.
* `ssh::Bool=false`: Whether or not to use SSH for the git remote. If `false` HTTPS will be used.
* `dev::Bool=true`: Whether or not to `Pkg.develop` generated packages.
* `manifest::Bool=false`: Whether or not to commit the `Manifest.toml`.
* `plugins::Vector{<:Plugin}=Plugin[]`: A list of `Plugin`s that the package will include.
"""
struct Template
    user::String
    host::String
    license::String
    authors::String
    dir::String
    julia_version::VersionNumber
    ssh::Bool
    dev::Bool
    manifest::Bool
    plugins::Dict{DataType, <:Plugin}

    function Template(;
        user::AbstractString="",
        host::AbstractString="https://github.com",
        license::AbstractString="MIT",
        authors::Union{AbstractString, Vector{<:AbstractString}}="",
        dir::AbstractString=Pkg.devdir(),
        julia_version::VersionNumber=default_version(),
        ssh::Bool=false,
        dev::Bool=true,
        manifest::Bool=false,
        plugins::Vector{<:Plugin}=Plugin[],
        git::Bool=true,
    )
        # Check for required Git options for package generation
        # (you can't commit to a repository without them).
        git && isempty(LibGit2.getconfig("user.name", "")) && missingopt("user.name")
        git && isempty(LibGit2.getconfig("user.email", "")) && missingopt("user.email")

        # If no username was set, look for one in the global git config.
        # Note: This is one of a few GitHub specifics (maybe we could use the host value).
        if isempty(user)
            user = LibGit2.getconfig("github.user", "")
        end
        if isempty(user)
            throw(ArgumentError("No GitHub username found, set one with user=username"))
        end

        host = URI(startswith(host, "https://") ? host : "https://$host").host

        if !isempty(license) && !isfile(joinpath(LICENSE_DIR, license))
            throw(ArgumentError("License '$license' is not available"))
        end

        # If no author was set, look for one in the global git config.
        if isempty(authors)
            authors = LibGit2.getconfig("user.name", "")
            email = LibGit2.getconfig("user.email", "")
            isempty(email) || (authors *= " <$email>")
        elseif authors isa Vector
            authors = join(authors, ", ")
        end

        dir = abspath(expanduser(dir))

        plugin_dict = Dict{DataType, Plugin}(typeof(p) => p for p in plugins)
        if (length(plugins) != length(plugin_dict))
            @warn "Plugin list contained duplicates, only the last of each type was kept"
        end

        new(user, host, license, authors, dir, julia_version, ssh, dev, manifest, plugin_dict)
    end
end

function Base.show(io::IO, t::Template)
    maybe(s::String) = isempty(s) ? "None" : s
    spc = "  "

    println(io, "Template:")
    println(io, spc, "→ User: ", maybe(t.user))
    println(io, spc, "→ Host: ", maybe(t.host))

    print(io, spc, "→ License: ")
    if isempty(t.license)
        println(io, "None")
    else
        println(io, t.license, " ($(t.authors) ", year(today()), ")")
    end

    println(io, spc, "→ Package directory: ", replace(maybe(t.dir), homedir() => "~"))
    println(io, spc, "→ Minimum Julia version: v", version_floor(t.julia_version))
    println(io, spc, "→ SSH remote: ", t.ssh ? "Yes" : "No")
    println(io, spc, "→ Add packages to main environment: ", t.dev ? "Yes" : "No")
    println(io, spc, "→ Commit Manifest.toml: ", t.manifest ? "Yes" : "No")

    print(io, spc, "→ Plugins:")
    if isempty(t.plugins)
        print(io, " None")
    else
        for plugin in sort(collect(values(t.plugins)); by=string)
            println(io)
            buf = IOBuffer()
            show(buf, plugin)
            print(io, spc^2, "• ")
            print(io, join(split(String(take!(buf)), "\n"), "\n$(spc^2)"))
        end
    end
end

"""
    interactive_template(; fast::Bool=false) -> Template

Interactively create a [`Template`](@ref). If `fast` is set, defaults will be assumed for
all values except username and plugins.
"""
function interactive_template(; git::Bool=true, fast::Bool=false)
    @info "Default values are shown in [brackets]"
    # Getting the leaf types in a separate thread eliminates an awkward wait after
    # "Select plugins" is printed.
    plugin_types = @async leaves(Plugin)
    kwargs = Dict{Symbol, Any}()

    default_user = LibGit2.getconfig("github.user", "")
    print("Username [", isempty(default_user) ? "REQUIRED" : default_user, "]: ")
    user = readline()
    kwargs[:user] = if !isempty(user)
        user
    elseif !isempty(default_user)
        default_user
    else
        throw(ArgumentError("Username is required"))
    end

    kwargs[:host] = if fast || !git
        "https://github.com"  # If Git isn't enabled, this value never gets used.
    else
        default_host = "github.com"
        print("Code hosting service [$default_host]: ")
        host = readline()
        isempty(host) ? default_host : host
    end

    kwargs[:license] = if fast
        "MIT"
    else
        println("License:")
        io = IOBuffer()
        available_licenses(io)
        licenses = ["" => "", collect(LICENSES)...]
        menu = RadioMenu(String["None", split(String(take!(io)), "\n")...])
        # If the user breaks out of the menu with Ctrl-c, the result is -1, the absolute
        # value of which correponds to no license.
        first(licenses[abs(request(menu))])
    end

    # We don't need to ask for authors if there is no license,
    # because the license is the only place that they matter.
    kwargs[:authors] = if fast || isempty(kwargs[:license])
        LibGit2.getconfig("user.name", "")
    else
        default_authors = LibGit2.getconfig("user.name", "")
        default_str = isempty(default_authors) ? "None" : default_authors
        print("Package author(s) [$default_str]: ")
        authors = readline()
        isempty(authors) ? default_authors : authors
    end

    kwargs[:dir] = if fast
        Pkg.devdir()
    else
        default_dir = Pkg.devdir()
        print("Path to package directory [$default_dir]: ")
        dir = readline()
        isempty(dir) ? default_dir : dir
    end

    kwargs[:julia_version] = if fast
        VERSION
    else
        default_julia_version = VERSION
        print("Minimum Julia version [", version_floor(default_julia_version), "]: ")
        julia_version = readline()
        isempty(julia_version) ? default_julia_version : VersionNumber(julia_version)
    end

    kwargs[:ssh] = if fast || !git
        false
    else
        print("Set remote to SSH? [no]: ")
        uppercase(readline()) in ["Y", "YES", "T", "TRUE"]
    end

    kwargs[:dev] = if fast
        true
    else
        print("Add packages to main environment? [yes]: ")
        uppercase(readline()) in ["", "Y", "YES", "T", "TRUE"]
    end

    kwargs[:manifest] = if fast
        false
    else
        print("Commit Manifest.toml? [no]: ")
        uppercase(readline()) in ["Y", "YES", "T", "TRUE"]
    end

    println("Plugins:")
    # Only include plugin types which have an `interactive` method.
    plugin_types = filter(t -> hasmethod(interactive, (Type{t},)), fetch(plugin_types))
    type_names = map(t -> split(string(t), ".")[end], plugin_types)
    menu = MultiSelectMenu(String.(type_names); pagesize=length(type_names))
    selected = collect(request(menu))
    kwargs[:plugins] = Vector{Plugin}(map(interactive, getindex(plugin_types, selected)))

    return Template(; git=git, kwargs...)
end

leaves(T::Type)::Vector{DataType} = isconcretetype(T) ? [T] : vcat(leaves.(subtypes(T))...)

missingopt(name) = @warn "Git config option '$name' missing, package generation will fail unless you supply a GitConfig"


text:
using Random, LinearAlgebra
using Stheno: FiniteKernel, FiniteCrossKernel, LhsFiniteCrossKernel, RhsFiniteCrossKernel,
    LhsFiniteZeroCrossKernel, RhsFiniteZeroCrossKernel

@testset "indexing" begin

let
    # Set up some GPs.
    rng, N, N′, D, gpc = MersenneTwister(123456), 4, 5, 2, GPC()
    X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
    μ1, μ2 = ConstantMean(1), ConstantMean(-1.5)
    k1, k2 = EQ(), Linear(0.5)
    f1, f2 = GP(μ1, k1, gpc), GP(μ2, k2, gpc)
    f3 = f1(X)

    # Check mean and marginal covariance under indexing.
    @test mean(f3) == FiniteMean(μ1, X)
    @test kernel(f3) == FiniteKernel(k1, X)
    @test kernel(f3, f1) == LhsFiniteCrossKernel(kernel(f1), X)
    @test kernel(f1, f3) == RhsFiniteCrossKernel(kernel(f1), X)
    @test kernel(f3, f2) == LhsFiniteZeroCrossKernel(X)
    @test kernel(f2, f3) == RhsFiniteZeroCrossKernel(X)

    # Check that kernel types are correct.
    f4 = f1(X′)
    @test kernel(f4) isa FiniteKernel
    @test kernel(f4, f4) isa FiniteKernel
    @test kernel(f4, f3) isa FiniteCrossKernel
    @test kernel(f3, f4) isa FiniteCrossKernel
    @test kernel(f4, f1) isa LhsFiniteCrossKernel
    @test kernel(f1, f4) isa RhsFiniteCrossKernel
    @test kernel(f2, f4) isa RhsFiniteZeroCrossKernel
    @test kernel(f4, f2) isa LhsFiniteZeroCrossKernel

    # These tests are currently unnecessary as nested indexing has been removed.
#     # Check that nested indexing works as expected.
#     f5 = f3(1:2)
#     @test typeof(kernel(f5)) <: Finite
#     @test typeof(kernel(f4, f5)) <: Finite
#     @test typeof(kernel(f5, f4)) <: Finite
#     @test typeof(kernel(f5, f3)) <: Finite
#     @test typeof(kernel(f3, f5)) <: Finite
#     @test typeof(kernel(f5, f2)) <: LhsFinite
#     @test typeof(kernel(f2, f5)) <: RhsFinite
#     @test typeof(kernel(f5, f1)) <: LhsFinite
#     @test typeof(kernel(f1, f5)) <: RhsFinite

#     # Check that the kernels evaluate correctly.
#     id5, id4, id3 = eachindex(f5), eachindex(f4), eachindex(f3)
#     @test kernel(f5).(id5, id5') == kernel(f3).(id5, id5')
#     @test kernel(f5, f4).(id5, id4') == kernel(f1).(x[id5], y[id4]')
#     @test kernel(f4, f5).(id4, id5') == kernel(f5, f4).(id5, id4')'
#     @test kernel(f5, f3).(id5, id3') == kernel(f1).(x[id5], x[id3]')
#     @test kernel(f3, f5).(id3, id5') == kernel(f5, f3).(id5, id3')'
end

let
    # Set up some GPs.
    rng, N, N′, D, gpc = MersenneTwister(123456), 4, 5, 2, GPC()
    X, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
    μ1, μ2 = ConstantMean(1.0), ConstantMean(-1.5)
    k1, k2 = EQ(), Linear(0.5)
    f1, f2 = GP(μ1, k1, gpc), GP(μ2, k2, gpc)
    f3, f4 = f1(X), f2(X′)
    fj = BlockGP([f1, f2])
    fjXX′ = fj(BlockData([X, X′]))

    @test fjXX′ isa BlockGP

    Σ11, Σ12, Σ22 = xcov(f3, f3), xcov(f3, f4), xcov(f4, f4)
    Σ_manual = BlockMatrix(reshape([Σ11, Σ12', Σ12, Σ22], 2, 2))

    # Test that indexing a BlockGP is consistent with separate indexing.
    @test mean_vec(fjXX′) == vcat(mean_vec(f3), mean_vec(f4))
    @test cov(fjXX′) == Σ_manual
    @test xcov(fjXX′, fjXX′) == cov(fjXX′)
    @test xcov(fjXX′, f3) == BlockMatrix(reshape([xcov(f3, f3), xcov(f4, f3)], 2, 1))

    # Check that broadcasting functionality works appropriately.
    @test cov(BlockGP([f1, f2])(X)) == cov(BlockGP([f1(X), f2(X)]))
    @test cov(BlockGP([f1, f2])(X)(1:N-1)) == cov(BlockGP([f1(X), f2(X)])(1:N-1))
end

end # @testset indexing


text:
@testset "Query" begin
    @testset "is_bidirected" begin
        g = rand_graph(10, 20, bidirected=true, graph_type=GRAPH_T)
        @test is_bidirected(g)

        g = rand_graph(10, 20, bidirected=false, graph_type=GRAPH_T)
        @test !is_bidirected(g)
    end

    @testset "has_multi_edges" begin
        if GRAPH_T == :coo
            s = [1, 1, 2, 3]
            t = [2, 2, 2, 4]
            g = GNNGraph(s, t, graph_type=GRAPH_T)
            @test has_multi_edges(g)

            s = [1, 2, 2, 3]
            t = [2, 1, 2, 4]
            g = GNNGraph(s, t, graph_type=GRAPH_T)
            @test !has_multi_edges(g)
        end
    end

    @testset "has_self_loops" begin
        s = [1, 1, 2, 3]
        t = [2, 2, 2, 4]
        g = GNNGraph(s, t, graph_type=GRAPH_T)
        @test has_self_loops(g)

        s = [1, 1, 2, 3]
        t = [2, 2, 3, 4]
        g = GNNGraph(s, t, graph_type=GRAPH_T)
        @test !has_self_loops(g)
    end

    @testset "degree" begin
        @testset "unweighted" begin
            s = [1, 1, 2, 3]
            t = [2, 2, 2, 4]
            g = GNNGraph(s, t, graph_type=GRAPH_T)

            @test degree(g) == degree(g; dir=:out) == [2, 1, 1, 0] # default is outdegree
            @test degree(g; dir=:in) == [0, 3, 0, 1]
            @test degree(g; dir=:both) == [2, 4, 1, 1]
            @test eltype(degree(g, Float32)) == Float32

            if TEST_GPU
                g_gpu = g |> gpu
                d = degree(g)
                d_gpu = degree(g_gpu)
                @test d_gpu isa CuVector{Int}
                @test Array(d_gpu) == d
            end
        end

        @testset "weighted" begin
            # weighted degree
            s = [1, 1, 2, 3]
            t = [2, 2, 2, 4]
            eweight = [0.1, 2.1, 1.2, 1]
            g = GNNGraph((s, t, eweight), graph_type=GRAPH_T)
            @test degree(g) ==  [2.2, 1.2, 1.0, 0.0]
            @test degree(g, edge_weight=nothing) ==  degree(g)
            d = degree(g, edge_weight=false)
            if GRAPH_T == :coo
                @test d ==  [2, 1, 1, 0]
            else
                # Adjacency matrix representation cannot disambiguate multiple edges
                # and edge weights
                @test d ==  [1, 1, 1, 0]
            end
            @test eltype(d) <: Integer
            if GRAPH_T == :coo
                @test degree(g, edge_weight=2*eweight) == [4.4, 2.4, 2.0, 0.0]
            end

            if TEST_GPU
                g_gpu = g |> gpu
                d = degree(g)
                d_gpu = degree(g_gpu)
                @test d_gpu isa CuVector{Float32}
                @test Array(d_gpu) ≈ d
            end
        end
    end

    @testset "laplacian_matrix" begin
        g = rand_graph(10, 30, graph_type=GRAPH_T)
        A = adjacency_matrix(g)
        D = Diagonal(vec(sum(A, dims=2)))
        L = laplacian_matrix(g)
        @test eltype(L) == eltype(g)
        @test L ≈ D - A
    end

    @testset "adjacency_matrix" begin
        a = sprand(5, 5, 0.5)
        g = GNNGraph(a, graph_type=GRAPH_T)
        A = adjacency_matrix(g, Float32)
        @test a ≈ A
        @test eltype(A) == Float32
    end
end


text:
@symbol_func function K_DV(cur_reactor::AbstractReactor)
  cur_K = _K_DV(cur_reactor)

  cur_K *= K_n(cur_reactor) ^ 2

  cur_K /= cur_reactor.max_q_DV

  cur_K ^= (1/3.2)

  cur_K
end

@symbol_func function _K_DV(cur_reactor::AbstractReactor)
  cur_K = 27.75e-3

  cur_K /= cur_reactor.epsilon ^ 1.2

  cur_K *= K_P(cur_reactor)

  cur_K /= ( 1 + cur_reactor.kappa_95 ^ 2 ) ^ 0.6

  cur_K
end


text:
struct Room
    region::CartesianIndices{2,Tuple{UnitRange{Int},UnitRange{Int}}}
end

Room(origin, height, width) = Room(CartesianIndices((origin.I[1] : origin.I[1] + height - 1, origin.I[2] : origin.I[2] + width - 1)))

get_origin(region::CartesianIndices{2}) = region[1, 1]
get_origin(room::Room) = get_origin(room.region)

get_interior(room::Room) = CartesianIndices((room.region.indices[1].start + 1 : room.region.indices[1].stop - 1,
                                             room.region.indices[2].start + 1 : room.region.indices[2].stop - 1))

function is_intersecting(room1::Room, room2::Room)
    intersection = intersect(get_interior(room1), room2.region)
    length(intersection) > 0 ? true : false
end

place_room!(env::AbstractGridWorld, room::Room) = place_room!(get_world(env), room)

function place_room!(world::GridWorldBase, room::Room)
    top = room.region.indices[1].start
    bottom = room.region.indices[1].stop
    left = room.region.indices[2].start
    right = room.region.indices[2].stop

    world[WALL, top, left:right] .= true
    world[WALL, bottom, left:right] .= true
    world[WALL, top:bottom, left] .= true
    world[WALL, top:bottom, right] .= true
end

include("empty_room.jl")
include("grid_rooms.jl")
include("sequential_rooms.jl")
include("maze.jl")
include("go_to_target.jl")
include("door_key.jl")
include("collect_gems.jl")
include("dynamic_obstacles.jl")
include("sokoban/sokoban.jl")
include("snake.jl")
include("catcher.jl")
include("transport.jl")


text:
"""
```julia
PlutoUI.wrapped(render_function::Function)
```

Wrap an existing widget inside HTML content.

`render_function` is a function that you write yourself, take a look at the examples below.

# Examples

## 🐶

We use the [`do` syntax](https://docs.julialang.org/en/v1/manual/functions/#Do-Block-Syntax-for-Function-Arguments) to write our `render_function`. The `Child` function is wrapped around the input that we want to capture.

```julia
DogsSlider(range) = wrapped() do Child
	md""\"
	### How many dogs do you have?

	I have \$(
		Child(Slider(range))
	) dogs.
	""\"
end

@bind num_dogs DogsSlider(5:10)

num_dogs == 5 # (initially)
```

> The output looks like:
>
> ![screenshot of running the code above inside Pluto](https://user-images.githubusercontent.com/6933510/148256466-4c525a00-de86-40c4-a6c4-388b5be3dc1c.png)

"""
wrapped(f::Function) = transformed_value(only, combine(f))


text:
using OhMyJulia
using JLD2
using FileIO

const ϵ = .1

struct FFM{T<:Real}
    w::Vector{T}
    V::Matrix{Vector{T}}
end

function FFM{T}(n::Int, nf::Int, k::Int=4) where T<:Real
    FFM{T}(zeros(T, n+1), [rand(T, k) - .5 for i in 1:n, j in 1:nf])
end

struct Trainer{T<:Real} # RMSProp
    gw::Vector{T}
    gV::Matrix{Vector{T}}
    vw::Vector{T}
    vV::Matrix{Vector{T}}
end

function Trainer(m::FFM{T}) where T
    t = Trainer{T}(deepcopy(m.w), deepcopy(m.V), deepcopy(m.w), deepcopy(m.V))
    reset!(t, true)
    t
end

function reset!(t::Trainer, acc=false)
    t.gw .= 0
    for v in t.gV
        v .= 0
    end

    acc || return

    t.vw .= 0
    for v in t.vV
        v .= 0
    end
end

function predict(m, X)
    n = length(X)

    s = car(m.w) + mapreduce(x->m.w[car(x) + 1], +, X)
    for k1 in 1:n, k2 in k1+1:n
        i, fi = X[k1]
        j, fj = X[k2]

        s += m.V[i, fj]' * m.V[j, fi]
    end

    s
end

function learn!(m, t, X, y, clip=Inf)
    Δ = max(min(predict(m, X) - y, clip), -clip)

    t.gw[1] += Δ
    for (i, fi) in X
        t.gw[i+1] += Δ
    end

    n = length(X)
    for k1 in 1:n, k2 in k1+1:n
        i, fi = X[k1]
        j, fj = X[k2]

        t.gV[i, fj] .+= Δ .* m.V[j, fi]
        t.gV[j, fi] .+= Δ .* m.V[i, fj]
    end
end

function step!(m, t, lr=.1, L2=1e-4, γ = .99)
    @. t.vw = γ * t.vw + t.gw ^ 2
    for (vv, gv) in zip(t.vV, t.gV)
        @. vv = γ * vv + gv ^ 2
    end

    m.w .*= 1 - L2
    @. m.w -= lr * t.gw / sqrt(t.vw + ϵ)
    for (v, vv, gv) in zip(m.V, t.vV, t.gV)
        v .*= 1 - L2
        @. v -= lr * gv / sqrt(vv + ϵ)
    end

    reset!(t)
end

function loss(m, X, y)
    .5(predict(m, X) - y) ^ 2
end

function train(train_file, test_file)
    trainX, trainy, n1, nf1 = read_libsvm(train_file)
    testX,  testy,  n2, nf2 = read_libsvm(test_file)

    n, nf = max(n1, n2), max(nf1, nf2)
    m = FFM{f64}(n, nf)
    t = Trainer(m)

    L = mapreduce(x->loss(m, car(x), cadr(x)), +, zip(testX, testy))
    println("epoch: 0, loss: $L")

    best = (nothing, Inf)
    for epoch = 1:50
        for i in 1:length(trainy)
            learn!(m, t, trainX[i], trainy[i])
            i % 10 == 0 && step!(m, t)
        end

        L = mapreduce(x->loss(m, car(x), cadr(x)), +, zip(testX, testy))
        best = L < cadr(best) ? (deepcopy(m), L) : best
        println("epoch: $epoch, loss: $L")
    end

    save("models.jld2", Dict("last" => m, "best" => car(best)))
end

function read_libsvm(file)
    Xs, ys = [], []
    n, nf = Ref(0), Ref(0)
    for line in eachline(split, file)
        push!(ys, parse(Int, car(line)))
        push!(Xs, map(cdr(line)) do x
            field, feature = split(x, ':') # file encoding starts with 0
            field, feature = parse.(Int, (field, feature)) .+ 1
            n[], nf[] = max(n[], feature), max(nf[], field)
            feature, field
        end)
    end
    Xs, ys, n[], nf[]
end


text:
const PARALLEL_DOC = """
    @parallel kernelcall

!!! note "Advanced"
        @parallel ranges kernelcall
        @parallel nblocks nthreads kernelcall
        @parallel ranges nblocks nthreads kernelcall
        @parallel (...) kwargs... kernelcall

Declare the `kernelcall` parallel. The kernel will automatically be called as required by the package for parallelization selected with [`@init_parallel_kernel`](@ref). Synchronizes at the end of the call (if a stream is given via keyword arguments, then it synchronizes only this stream).

# Arguments
- `kernelcall`: a call to a kernel that is declared parallel.
!!! note "Advanced optional arguments"
    - `ranges::Tuple{UnitRange{},UnitRange{},UnitRange{}} | Tuple{UnitRange{},UnitRange{}} | Tuple{UnitRange{}} | UnitRange{}`: the ranges of indices in each dimension for which computations must be performed.
    - `nblocks::Tuple{Integer,Integer,Integer}`: the number of blocks to be used if the package CUDA was selected with [`@init_parallel_kernel`](@ref).
    - `nthreads::Tuple{Integer,Integer,Integer}`: the number of threads to be used if the package CUDA was selected with [`@init_parallel_kernel`](@ref).
    - `kwargs...`: keyword arguments to be passed further to CUDA (ignored for Threads).

!!! note "Performance note"
    Kernel launch parameters are automatically defined with heuristics, where not defined with optional kernel arguments. For CUDA `nthreads` is whenever reasonable set to (32,8,1) and `nblocks` accordingly to ensure that enough threads are launched.

See also: [`@init_parallel_kernel`](@ref)
"""
@doc PARALLEL_DOC
macro parallel(args...) check_initialized(); checkargs_parallel(args...); esc(parallel(args...)); end

##
const PARALLEL_INDICES_DOC = """
    @parallel_indices indices kernel

Declare the `kernel` parallel and generate the given parallel `indices` inside the `kernel` using the package for parallelization selected with [`@init_parallel_kernel`](@ref).
"""
@doc PARALLEL_INDICES_DOC
macro parallel_indices(args...) check_initialized(); checkargs_parallel_indices(args...); esc(parallel_indices(args...)); end

##
const PARALLEL_ASYNC_DOC = """
@parallel_async kernelcall

!!! note "Advanced"
        @parallel_async ranges kernelcall
        @parallel_async nblocks nthreads kernelcall
        @parallel_async ranges nblocks nthreads kernelcall
        @parallel_async (...) kwargs... kernelcall

Declare the `kernelcall` parallel as with [`@parallel`](@ref) (see [`@parallel`](@ref) for more information); deactivates however automatic synchronization at the end of the call. Use [`@synchronize`](@ref) for synchronizing.

!!! note "Performance note"
    @parallel_async falls currently back to running synchronously if the package Threads was selected with [`@init_parallel_kernel`](@ref).

See also: [`@synchronize`](@ref), [`@parallel`](@ref)
"""
@doc PARALLEL_ASYNC_DOC
macro parallel_async(args...) check_initialized(); checkargs_parallel(args...); esc(parallel_async(args...)); end

##
const SYNCHRONIZE_DOC = """
    @synchronize()

Synchronize the GPU/CPU.

See also: [`@parallel_async`](@ref)
"""
@doc SYNCHRONIZE_DOC
macro synchronize(args...) check_initialized(); esc(synchronize(args...)); end

## MACROS FORCING PACKAGE, IGNORING INITIALIZATION

macro parallel_cuda(args...)            check_initialized(); checkargs_parallel(args...); esc(parallel(args...; package=PKG_CUDA)); end
macro parallel_threads(args...)         check_initialized(); checkargs_parallel(args...); esc(parallel(args...; package=PKG_THREADS)); end
macro parallel_indices_cuda(args...)    check_initialized(); checkargs_parallel_indices(args...); esc(parallel_indices(args...; package=PKG_CUDA)); end
macro parallel_indices_threads(args...) check_initialized(); checkargs_parallel_indices(args...); esc(parallel_indices(args...; package=PKG_THREADS)); end
macro parallel_async_cuda(args...)      check_initialized(); checkargs_parallel(args...); esc(parallel_async(args...; package=PKG_CUDA)); end
macro parallel_async_threads(args...)   check_initialized(); checkargs_parallel(args...); esc(parallel_async(args...; package=PKG_THREADS)); end
macro synchronize_cuda(args...)         check_initialized(); esc(synchronize(args...; package=PKG_CUDA)); end
macro synchronize_threads(args...)      check_initialized(); esc(synchronize(args...; package=PKG_THREADS)); end

## ARGUMENT CHECKS

function checkargs_parallel(args...)
    if isempty(args) @ArgumentError("arguments missing.") end
    if !is_call(args[end]) @ArgumentError("the last argument must be a kernel call (obtained: $(args[end])).") end
    posargs, = split_args(args)
    if length(posargs) > 3 @ArgumentError("too many positional arguments.") end
end

function checkargs_parallel_indices(args...)
    if (length(args) != 2) @ArgumentError("wrong number of arguments.") end
    if !is_function(args[end]) @ArgumentError("the last argument must be a function definition (obtained: $(args[end])).") end
end

## GATEWAY FUNCTIONS

parallel_async(args::Union{Symbol,Expr}...; package::Symbol=get_package()) = parallel(args...; package=package, async=true)

function parallel(args::Union{Symbol,Expr}...; package::Symbol=get_package(), async::Bool=false)
    posargs, kwargs, kernelarg = split_args(args)
    if     (package == PKG_CUDA)    parallel_call_cuda(posargs..., kernelarg, kwargs, async)
    elseif (package == PKG_THREADS) parallel_call_threads(posargs..., kernelarg, async) # Ignore keyword args as they are not for the threads case (noted in doc).
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

function parallel_indices(args::Union{Symbol,Expr}...; package::Symbol=get_package(), async::Bool=false)
    numbertype = get_numbertype()
    parallel_kernel(package, numbertype, args...)
end

function synchronize(args::Union{Symbol,Expr}...; package::Symbol=get_package())
    if     (package == PKG_CUDA)    synchronize_cuda(args...)
    elseif (package == PKG_THREADS) synchronize_threads(args...)
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

## @PARALLEL KERNEL FUNCTIONS

function parallel_kernel(package::Symbol, numbertype::DataType, indices::Union{Symbol,Expr}, kernel::Expr)
    if (!isa(indices,Symbol) && !isa(indices.head,Symbol)) @ArgumentError("@parallel_indices: argument 'indices' must be a tuple of indices or a single index (e.g. (ix, iy, iz) or (ix, iy) or ix ).") end
    if isa(indices,Expr) indices = indices.args else indices = [indices] end
    body = get_body(kernel)
    body = remove_return(body)
    if (package == PKG_CUDA)
        kernel = substitute(kernel, :(Data.Array), :(Data.DeviceArray))
    end
    kernel = push_to_signature!(kernel, :($RANGES_VARNAME::$RANGES_TYPE))
    if (package == PKG_CUDA)
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[1])::$INT_CUDA))
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[2])::$INT_CUDA))
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[3])::$INT_CUDA))
    end
    ranges = [:($RANGES_VARNAME[1]), :($RANGES_VARNAME[2]), :($RANGES_VARNAME[3])]
    if (package == PKG_CUDA)
        body = add_threadids(indices, ranges, body)
        body = literaltypes(numbertype, body)
    elseif (package == PKG_THREADS)
        body = add_loop(indices, ranges, body)
        body = literaltypes(numbertype, body)
    else
        @ArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
    body = add_return(body)
    set_body!(kernel, body)
    return kernel
end

## @PARALLEL CALL FUNCTIONS

function parallel_call_cuda(ranges::Union{Symbol,Expr}, nblocks::Union{Symbol,Expr}, nthreads::Union{Symbol,Expr}, kernelcall::Expr, kwargs::Array, async::Bool)
    ranges = :(ParallelStencil.ParallelKernel.promote_ranges($ranges))
    push!(kernelcall.args, ranges)
    push!(kernelcall.args, :($INT_CUDA(length($ranges[1]))))
    push!(kernelcall.args, :($INT_CUDA(length($ranges[2]))))
    push!(kernelcall.args, :($INT_CUDA(length($ranges[3]))))
    synccall = async ? :(begin end) : create_synccall_cuda(kwargs)
    return :( CUDA.@cuda blocks=$nblocks threads=$nthreads $(kwargs...) $kernelcall; $synccall )
end

function parallel_call_cuda(nblocks::Union{Symbol,Expr}, nthreads::Union{Symbol,Expr}, kernelcall::Expr, kwargs::Array, async::Bool)
    maxsize = :( $nblocks .* $nthreads )
    ranges  = :(ParallelStencil.ParallelKernel.compute_ranges($maxsize))
    parallel_call_cuda(ranges, nblocks, nthreads, kernelcall, kwargs, async)
end

function parallel_call_cuda(ranges::Union{Symbol,Expr}, kernelcall::Expr, kwargs::Array, async::Bool)
    maxsize  = :(length.(ParallelStencil.ParallelKernel.promote_ranges($ranges)))
    nthreads = :( ParallelStencil.ParallelKernel.compute_nthreads($maxsize) )
    nblocks  = :( ParallelStencil.ParallelKernel.compute_nblocks($maxsize, $nthreads) )
    parallel_call_cuda(ranges, nblocks, nthreads, kernelcall, kwargs, async)
end

function parallel_call_cuda(kernelcall::Expr, kwargs::Array, async::Bool)
    ranges = :( ParallelStencil.ParallelKernel.get_ranges($(kernelcall.args[2:end]...)) )
    parallel_call_cuda(ranges, kernelcall, kwargs, async)
end

function parallel_call_threads(ranges::Union{Symbol,Expr}, kernelcall::Expr, async::Bool)
    ranges = :(ParallelStencil.ParallelKernel.promote_ranges($ranges))
    push!(kernelcall.args, ranges)
    if async
        return kernelcall # NOTE: This cannot be used currently as there is no obvious solution how to sync in this case.
    else
        return kernelcall
    end
end

function parallel_call_threads(ranges::Union{Symbol,Expr}, nblocks::Union{Symbol,Expr}, nthreads::Union{Symbol,Expr}, kernelcall::Expr, async::Bool)
    parallel_call_threads(ranges, kernelcall, async)
end

function parallel_call_threads(nblocks::Union{Symbol,Expr}, nthreads::Union{Symbol,Expr}, kernelcall::Expr, async::Bool)
    maxsize = :( $nblocks .* $nthreads )
    ranges  = :(ParallelStencil.ParallelKernel.compute_ranges($maxsize))
    parallel_call_threads(ranges, kernelcall, async)
end

function parallel_call_threads(kernelcall::Expr, async::Bool)
    ranges = :( ParallelStencil.ParallelKernel.get_ranges($(kernelcall.args[2:end]...)) )
    parallel_call_threads(ranges, kernelcall, async)
end

## @SYNCHRONIZE FUNCTIONS

synchronize_cuda(args::Union{Symbol,Expr}...) = :(CUDA.synchronize($(args...)))
synchronize_threads(args::Union{Symbol,Expr}...) = :(begin end)

## MACROS AND FUNCTIONS TO MODIFY THE TYPE OF LITERALS

macro literaltypes(type, expr) check_literaltype(type); esc(literaltypes(@eval(__module__, $type), expr)) end
macro literaltypes(type1, type2, expr) check_literaltype(type1); check_literaltype(type2); esc(literaltypes(@eval(__module__, $type1), @eval(__module__, $type2), expr)) end

function literaltypes(type::DataType, expr::Expr)
    args = expr.args
    head = expr.head
    for i=1:length(args)
        if type <: Integer && typeof(args[i]) <: Integer && head == :call && args[1] in OPERATORS # NOTE: only integers in operator calls are modified (not in other function calls as e.g. size).
            literal = type(args[i])
            args[i] = :($literal)
        elseif type <: AbstractFloat && typeof(args[i]) <: AbstractFloat
            literal = type(args[i])
            args[i] = :($literal)
        elseif typeof(args[i]) == Expr
            args[i] = literaltypes(type, args[i])
        end
    end
    return expr
end

function literaltypes(type1::DataType, type2::DataType, expr::Expr)
    if !((type1 <: Integer && type2 <: AbstractFloat) || (type1 <: AbstractFloat && type2 <: Integer))
        @ModuleInternalError("incoherent type arguments: one must be a subtype of Integer and the other a subtype of AbstractFloat.")
    end
    return literaltypes(type2, literaltypes(type1,expr))
end

## FUNCTIONS TO ADD THREAD-IDS / LOOPS IN KERNELS

function add_threadids(indices::Array, ranges::Array, block::Expr)
    if !(length(ranges)==3) @ModuleInternalError("ranges must be an Array or Tuple of size 3.") end # E.g. (5:28,5:28,1:1) in 2D. Note that for simplicity for the users and developers, 1D and 2D problems are always expressed like 3D problems...
    rangelength_x, rangelength_y, rangelength_z = RANGELENGTHS_VARNAMES
    tx, ty, tz = THREADIDS_VARNAMES
    ndims = length(indices)
    if ndims == 1
        ix, = indices
        range_x, = ranges
        quote
            $tx = (CUDA.blockIdx().x-1) * CUDA.blockDim().x + CUDA.threadIdx().x;  # thread ID, dimension x
            if ($tx > $rangelength_x) return; end
            $ix = $range_x[$tx]                                                    # index, dimension x
            $block
        end
    elseif ndims == 2
        ix, iy = indices
        range_x, range_y = ranges
        quote
            $tx = (CUDA.blockIdx().x-1) * CUDA.blockDim().x + CUDA.threadIdx().x;  # thread ID, dimension x
            $ty = (CUDA.blockIdx().y-1) * CUDA.blockDim().y + CUDA.threadIdx().y;  # thread ID, dimension y
            if ($tx > $rangelength_x) return; end
            if ($ty > $rangelength_y) return; end
            $ix = $range_x[$tx]                                                    # index, dimension x
            $iy = $range_y[$ty]                                                    # index, dimension y
            $block
        end
    elseif ndims == 3
        ix, iy, iz = indices
        range_x, range_y, range_z = ranges
        quote
            $tx = (CUDA.blockIdx().x-1) * CUDA.blockDim().x + CUDA.threadIdx().x;  # thread ID, dimension x
            $ty = (CUDA.blockIdx().y-1) * CUDA.blockDim().y + CUDA.threadIdx().y;  # thread ID, dimension y
            $tz = (CUDA.blockIdx().z-1) * CUDA.blockDim().z + CUDA.threadIdx().z;  # thread ID, dimension z
            if ($tx > $rangelength_x) return; end
            if ($ty > $rangelength_y) return; end
            if ($tz > $rangelength_z) return; end
            $ix = $range_x[$tx]                                                    # index, dimension x
            $iy = $range_y[$ty]                                                    # index, dimension y
            $iz = $range_z[$tz]                                                    # index, dimension z
            $block
        end
    end
end

function add_loop(indices::Array, ranges::Array, block::Expr)
    if !(length(ranges)==3) @ModuleInternalError("ranges must be an Array or Tuple of size 3") end # E.g. (5:28,5:28,1:1) in 2D.
    ndims = length(indices)
    if ndims == 1
        quote
            Base.Threads.@threads for $(indices[1]) in $(ranges[1])
                $block
            end
        end
    elseif ndims == 2
        quote
            Base.Threads.@threads for $(indices[2]) in $(ranges[2])
                for $(indices[1]) in $(ranges[1])
                    $block
                end
            end
        end
    elseif ndims == 3
        quote
            Base.Threads.@threads for $(indices[3]) in $(ranges[3])
                for $(indices[2]) in $(ranges[2]), $(indices[1]) in $(ranges[1])
                    $block
                end
            end
        end
    end
end

## FUNCTIONS TO GET RANGES, PROMOTE THEM TO 3D, AND COMPUTE NTHREADS, NBLOCKS AND RANGES

promote_ranges(ranges::RANGES_TYPE_1D)          = (ranges,    1:1, 1:1)
promote_ranges(ranges::RANGES_TYPE_1D_TUPLE)    = (ranges..., 1:1, 1:1)
promote_ranges(ranges::RANGES_TYPE_2D)          = (ranges..., 1:1)
promote_ranges(ranges::RANGES_TYPE)             = ranges
promote_ranges(ranges)                          = @ModuleInternalError("ranges must be a Tuple of UnitRange of size 1, 2 or 3 (obtained: $ranges; its type is: $(typeof(ranges))).")

promote_maxsize(maxsize::MAXSIZE_TYPE_1D)       = (maxsize,    1, 1)
promote_maxsize(maxsize::MAXSIZE_TYPE_1D_TUPLE) = (maxsize..., 1, 1)
promote_maxsize(maxsize::MAXSIZE_TYPE_2D)       = (maxsize..., 1)
promote_maxsize(maxsize::MAXSIZE_TYPE)          = maxsize
promote_maxsize(maxsize)                        = @ModuleInternalError("maxsize must be a Tuple of Integer of size 1, 2 or 3 (obtained: $maxsize; its type is: $(typeof(maxsize))).")

function get_ranges(args...)
    complex_args = [x for x in args if !(isa(x, AbstractArray) || isa(x, Number))]
    if length(complex_args) > 0 @ArgumentError("automatic detection of ranges not possible in @parallel call: some kernel arguments are neither arrays nor scalars. Specify ranges or nthreads and nblocks manually.") end
    arrays = [x for x in args if isa(x, AbstractArray)];
    if (length(arrays) == 0) @ArgumentError("automatic detection of ranges not possible: no arrays detected in kernel arguments. Specify ranges or nthreads and nblocks manually.") end
    maxsize = ( maximum([size(A,1) for A in arrays]), maximum([size(A,2) for A in arrays]), maximum([size(A,3) for A in arrays]) )
    return (1:max(maxsize[1],1), 1:max(maxsize[2],1), 1:max(maxsize[3],1))
end

function compute_ranges(maxsize)
    maxsize = promote_maxsize(maxsize)
    return (1:maxsize[1], 1:maxsize[2], 1:maxsize[3])
end

function compute_nthreads(maxsize) # This is a heuristic, which results in (32,8,1) threads, except if maxsize[1] < 32 or maxsize[2] < 8.
    maxsize = promote_maxsize(maxsize)
    nthreads_x = min(32, maxsize[1])
    nthreads_y = min(ceil(Int,NTHREADS_MAX/nthreads_x), maxsize[2])
    nthreads_z = min(ceil(Int,NTHREADS_MAX/(nthreads_x*nthreads_y)), maxsize[3])
    return (nthreads_x, nthreads_y , nthreads_z)
end

function compute_nblocks(maxsize, nthreads)
    maxsize = promote_maxsize(maxsize)
    if !(isa(nthreads, Union{AbstractArray,Tuple}) && length(nthreads)==3) @ArgumentError("nthreads must be an Array or Tuple of size 3 (obtained: $nthreads; its type is: $(typeof(nthreads))).") end
    return ceil.(Int, maxsize./nthreads)
end

## FUNCTIONS TO CREATE SYNCHRONIZATION CALLS

function create_synccall_cuda(kwargs::Array)
    kwarg_stream = [x for x in kwargs if x.args[1]==:stream]
    if length(kwarg_stream) == 0
        synchronize_cuda()
    elseif length(kwarg_stream) == 1
        stream = kwarg_stream[1].args[2]
        synchronize_cuda(stream)
    else
        @KeywordArgumentError("there can only be one keyword argument stream in a @parallel call.")
    end
end

function create_synccall_threads(kwargs::Array)
    synchronize_threads()
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

UnionVector.eval(quote

FlatBuffers.@with_kw mutable struct Attacker
    sword_attack_damage::Int32 = 0
end
FlatBuffers.@ALIGN(Attacker, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:Attacker} = [
    0x00000004
]

Attacker(buf::AbstractVector{UInt8}) = FlatBuffers.read(Attacker, buf)
Attacker(io::IO) = FlatBuffers.deserialize(io, Attacker)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
"""
    periods(tdf::TimeDataFrame)
returns the periods of the TimeDataFrame
"""
periods(tdf::TimeDataFrame) = getfield(tdf, :periods)

"""
    firstperiod(tdf::TimeDataFrame)
returns the first period of the TimeDataFrame
"""
firstperiod(tdf::TimeDataFrame) = periods(tdf)[1]

"""
    lasttperiod(tdf::TimeDataFrame)
returns the last period of the TimeDataFrame
"""
lastperiod(tdf::TimeDataFrame) = periods(tdf)[end]

"""
    dataframe(tdf::TimeDataFrame)
returns the DataFrame inside the TimeDataFrame
"""
dataframe(tdf::TimeDataFrame) = getfield(tdf, :data)

"""
    iscontinuous(tdf::TimeDataFrame)
returns whether the TimeDataFrame is continuous
"""
iscontinuous(tdf::TimeDataFrame) = getfield(tdf, :continuous)


text:
using Sim, Printf
import Sim.init!

abstract type PState end
struct Idle <: PState end
struct Wait <: PState end
struct Unalert <: PState end

abstract type PEvent end
struct Start <: PEvent end
struct Serve <: PEvent end
struct Return <: PEvent end
struct Miss <: PEvent end

mutable struct Player
    name::AbstractString
    opp::Union{Number,Player}
    state::PState
    accuracy::Float64
    attentiveness::Float64
    score::Int64

    Player(name, acc, att) = new(name, 0, Idle(), acc, att, 0)
end

const dist = 3 # distance for ball to fly [m]
const vs   = 10 # serve velocity [m/s]
const vr   = 20 # return velocity [m/s]

rd(s::Float64) = randn()*s + 1

function init!(p::Player, opp::Player)
    p.opp = opp
    if rand() ≤ p.attentiveness
        p.state = Wait()
    else
        p.state = Unalert()
    end
end

function serve(p::Player)
    ts = 3 + dist*rd(0.15)/(vs*rd(0.25))
    if (rand() ≤ p.accuracy) && (p.state == Wait())
        event!(𝐶, :(step!($(p.opp), Serve())), after, ts)
        @printf("%.2f: %s serves %s\n", tau()+ts, p.name, p.opp.name)
    else
        event!(𝐶, :(step!($(p.opp), Miss())), after, ts)
        @printf("%.2f: %s serves and misses %s\n", tau()+ts, p.name, p.opp.name)
    end
    if rand() ≥ p.attentiveness
        p.state = Unalert()
    end
end

function ret(p::Player)
    tr = dist*rd(0.15)/(vr*rd(0.25))
    if rand() ≤ p.accuracy
        event!(𝐶, :(step!($(p.opp), Return())), after, tr)
        @printf("%.2f: %s returns %s\n", tau()+tr, p.name, p.opp.name)
    else
        event!(𝐶, :(step!($(p.opp), Miss())), after, tr)
        @printf("%.2f: %s returns and misses %s\n", tau()+tr, p.name, p.opp.name)
    end
    if rand() ≥ p.attentiveness
        p.state = Unalert()
    end
end

"default transition for players"
step!(p::Player, q::PState, σ::PEvent) =
        println("undefined transition for $(p.name), $q, $σ")

"player p gets a start command"
step!(p::Player, ::Union{Wait, Unalert}, ::Start) = serve(p)

"player p is waiting and gets served or returned"
step!(p::Player, ::Wait, ::Union{Serve, Return}) = ret(p)

"player p is unalert and gets served or returned"
function step!(p::Player, ::Unalert, ::Union{Serve, Return})
    @printf("%.2f: %s looses ball\n", tau(), p.name)
    p.opp.score += 1
    p.state = Wait()
    serve(p)
end

"player p is waiting or unalert and gets missed"
function step!(p::Player, ::Union{Wait, Unalert}, ::Miss)
    p.score += 1
    p.state = Wait()
    serve(p)
end

"simplified `step!` call"
step!(p::Player, σ::PEvent) = step!(p, p.state, σ)

ping = Player("Ping", 0.90, 0.90)
pong = Player("Pong", 0.90, 0.90)
init!(ping, pong)
init!(pong, ping)
step!(ping, Start())

run!(𝐶, 30)
println("Ping scored $(ping.score)")
println("Pong scored $(pong.score)")

reset!(𝐶)