text:
using BenchmarkTools, RData, MixedModels, StatsModels

const SUITE = BenchmarkGroup()

const dat = Dict(Symbol(k)=>v for (k,v) in load(joinpath(dirname(pathof(MixedModels)), "..", "test", "dat.rda")));

const mods = Dict{Symbol,Vector{Expr}}(
    :Alfalfa => [:(1+A*B+(1|G)), :(1+A+B+(1|G))],
    :Animal => [:(1+(1|G)+(1|H))],
    :Arabidopsis => [],             # glmm and rename variables
    :Assay => [:(1+A+B*C+(1|G)+(1|H))],
    :AvgDailyGain => [:(1+A*U+(1|G)), :(1+A+U+(1|G))],
    :BIB => [:(1+A*U+(1|G)), :(1+A+U+(1|G))],
    :Bond => [:(1+A+(1|G))],
    :Chem97 => [:(1+(1|G)+(1|H)),:(1+U+(1|G)+(1|H))],
    :Contraception => [],           # glmm and rename variables
    :Cultivation => [:(1+A*B+(1|G)), :(1+A+B+(1|G)), :(1+A+(1|G))],
    :Demand => [:(1+U+V+W+X+(1|G)+(1|H))],
    :Dyestuff => [:(1+(1|G))],
    :Dyestuff2 => [:(1+(1|G))],
    :Early => [:(1+U+U&A+(1+U|G))],
    :Exam => [:(1+A*U+B+(1|G)), :(1+A+B+U+(1|G))],
    :Gasoline => [:(1+U+(1|G))],
    :Gcsemv => [:(1+A+(1|G))],      # variables must be renamed
    :Genetics => [:(1+A+(1|G)+(1|H))],
    :HR => [:(1+A*U+V+(1+U|G))],
    :Hsb82 => [:(1+A+B+C+U+(1|G))],
    :IncBlk => [:(1+A+U+V+W+Z+(1|G))],
    :InstEval => [:(1+A+(1|G)+(1|H)+(1|I)),:(1+A*I+(1|G)+(1|H))],
    :KKL => [],                    # variables must be renamed
    :KWDYZ => [],                  # variables must be renamed
    :Mississippi => [:(1+A+(1|G))],
    :Mmmec => [],                  # glmm (and offset) and variables renamed
    :Multilocation => [:(1+A+(0+A|G)+(1|H))],
    :Oxboys => [:(1+U+(1+U|G))],
    :PBIB => [:(1+A+(1|G))],
    :Pastes => [:(1+(1|G)+(1|H))],
    :Penicillin => [:(1+(1|G)+(1|H))],
    :Pixel => [:(1+U+V+(1+U|G)+(1|H))],  # variables must be renamed
    :Poems => [:(1+U+V+W+(1|G)+(1|H)+(1|I))],
    :Rail => [:(1+(1|G))],
    :SIMS => [:(1+U+(1+U|G))],
    :ScotsSec => [:(1+A+U+V+(1|G)+(1|H))],
    :Semi2 => [:(1+A+(1|G)+(1|H))],
    :Semiconductor => [:(1+A*B+(1|G))],
    :Socatt => [],                 # variables must be renamed - binomial glmm?
    :TeachingII => [:(1+A+T+U+V+W+X+Z+(1|G))],
    :VerbAgg => [:(1+A+B+C+U+(1|G)+(1|H))], # Bernoulli glmm and rename variables
    :Weights => [:(1+A*U+(1+U|G))],
    :WWheat => [:(1+U+(1+U|G))],
    :bdf => [],                   # rename variables and look up model
    :bs10 => [:(1+U+V+W+((1+U+V+W)|G)+((1+U+V+W)|H))],
    :cake => [:(1+A*B+(1|G))],
    :cbpp => [:(1+A+(1|G))],      # Binomial glmm, create and rename variables
    :d3 => [:(1+U+(1|G)+(1|H)+(1|I)), :(1+U+(1+U|G)+(1+U|H)+(1+U|I))],
    :dialectNL => [:(1+A+T+U+V+W+X+(1|G)+(1|H)+(1|I))],
    :egsingle => [:(1+A+U+V+(1|G)+(1|H))],
    :epilepsy => [],              # unknown origin
    :ergoStool => [:(1+A+(1|G))],
    :gb12 => [:(1+S+T+U+V+W+X+Z+((1+S+U+W)|G)+((1+S+T+V)|H))],
    :grouseticks => [],           # rename variables, glmm needs formula
    :guImmun => [],               # rename variables, glmm needs formula
    :guPrenat => [],              # rename variables, glmm needs formula
    :kb07 => [:(1+S+T+U+V+W+X+Z+((1+S+T+U+V+W+X+Z)|G)+((1+S+T+U+V+W+X+Z)|H)),
              :(1+S+T+U+V+W+X+Z+(1|G)+((0+S)|G)+((0+T)|G)+((0+U)|G)+((0+V)|G)+((0+W)|G)+
              ((0+X)|G)+((0+Z)|G)+(1|H)+((0+S)|H)+((0+T)|H)+((0+U)|H)+((0+V)|H)+
              ((0+W)|H)+((0+X)|H)+((0+Z)|H))],
    :ml1m => [:(1+(1|G)+(1|H))],
    :paulsim => [:(1+S+T+U+(1|H)+(1|G))],  # names of H and G should be reversed
    :sleepstudy => [:(1+U+(1+U|G)), :(1+U+(1|G)+(0+U|G))],
    :s3bbx => [],                 # probably drop this one
    :star => []                   # not sure it is worthwhile working with these data
    );

fitbobyqa(rhs::Expr, dsname::Symbol) = fit(LinearMixedModel, @eval(@formula(Y ~ $rhs)), dat[dsname])
compactstr(ds,rhs) = replace(string(ds, ':', rhs), ' ' => "")

SUITE["simplescalar"] = BenchmarkGroup(["single", "simple", "scalar"])
for ds in [:Alfalfa, :AvgDailyGain, :BIB, :Bond, :cake, :Cultivation, :Dyestuff,
    :Dyestuff2, :ergoStool, :Exam, :Gasoline, :Hsb82, :IncBlk, :Mississippi,
    :PBIB, :Rail, :Semiconductor, :TeachingII]
    for rhs in mods[ds]
        SUITE["simplescalar"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa($(QuoteNode(rhs)), $(QuoteNode(ds)))
    end
end

SUITE["singlevector"] = BenchmarkGroup(["single", "vector"])
for ds in [:Early, :HR, :Oxboys, :SIMS, :sleepstudy, :Weights, :WWheat]
    for rhs in mods[ds]
        SUITE["singlevector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa($(QuoteNode(rhs)), $(QuoteNode(ds)))
    end
end

SUITE["nested"] = BenchmarkGroup(["multiple", "nested", "scalar"])
for ds in [:Animal, :Chem97, :Genetics, :Pastes, :Semi2]
    for rhs in mods[ds]
        SUITE["nested"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa($(QuoteNode(rhs)), $(QuoteNode(ds)))
    end
end

SUITE["crossed"] = BenchmarkGroup(["multiple", "crossed", "scalar"])
for ds in [:Assay, :Demand, :InstEval, :Penicillin, :ScotsSec,
           :dialectNL, :egsingle, :ml1m, :paulsim]
    for rhs in mods[ds]
        SUITE["crossed"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa($(QuoteNode(rhs)), $(QuoteNode(ds)))
    end
end

SUITE["crossedvector"] = BenchmarkGroup(["multiple", "crossed", "vector"])
for ds in [:bs10, :d3, :gb12, :kb07]
    for rhs in mods[ds]
        SUITE["crossedvector"][compactstr(ds, rhs)] = @benchmarkable fitbobyqa($(QuoteNode(rhs)), $(QuoteNode(ds)))
    end
end


text:
@symbol_func function G_V(cur_reactor::AbstractReactor)
  cur_G = safe_symbol(cur_reactor, :R_0) ^ 2

  cur_G -= ( R_CS(cur_reactor) + d(cur_reactor) ) ^ 2

  cur_G *= cur_reactor.T_bar ^ (3/2)

  cur_G *= cur_reactor.pi

  cur_G /= 1e7

  cur_G
end

ext: jl


text:
instances=[
    "afiro",
    "adlittle",
    "sc50a",
    "sc50b",
    "sc105",
    "sc205"
]
truefalse = [
    # true,
    false,
]
phaseone_methods = [
    Simplex.PhaseOne.ARTIFICIAL,
    Simplex.PhaseOne.CPLEX,
]
pivot_rules = [
    Simplex.Bland,
    Simplex.Steepest,
    Simplex.Dantzig,
]

function run_netlib_instance(netlib, use_gpu, method, pivot)::Simplex.Status
    c, xl, xu, bl, bu, A = Simplex.GLPK_read_mps(netlib)
    lp =  MatOI.LPForm{Float64, typeof(A), typeof(c)}(
        MOI.MIN_SENSE, c, A, bl, bu, xl, xu
    )
    Simplex.summary(lp)
    status = Simplex.run(lp,
        gpu = use_gpu,
        phase_one_method = method,
        pivot_rule = pivot)
    return status
end

@testset "netlib instances" begin
    for use_gpu in truefalse
        arch = use_gpu ? "GPU" : "CPU"
        @testset "$arch" begin
            for i in instances
                @testset "$(i)" begin
                    for method in phaseone_methods
                        @testset "$(Symbol(method))" begin
                            for pivot in pivot_rules
                                @testset "$(Symbol(pivot))" begin
                                    netlib = "../netlib/$i.mps"
                                    @test Simplex.Optimal == run_netlib_instance(netlib, use_gpu, method, pivot)
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 79789782-8d97-11eb-112c-fd8c2c294527
using Pkg, DrWatson

# ╔═╡ dd81850e-8d97-11eb-0cff-8bf4f38f0660
begin
	using Distributions
	using StatsPlots
	using StatsBase
	using LaTeXStrings
	using CSV
	using DataFrames
	using LinearAlgebra
	using Random
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ a3565770-8e31-11eb-1635-3b842d541ab6
md" ## Clip-11-01-08s.jl"

# ╔═╡ dd81af5c-8d97-11eb-10e4-73daebcd8875
begin
	df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame)
	df[!, "treatment"] = 1 .+ df.prosoc_left + 2 * df.condition
end;

# ╔═╡ 691258d8-8dab-11eb-3ec4-8b2059af0832
unique(df.treatment)

# ╔═╡ dd82254a-8d97-11eb-3f3d-59b7642304f5
stan11_1_1 = "
parameters {
 	real a;
}
model {
	a ~ normal(0, 10.0);
}
";

# ╔═╡ dd8c8ce4-8d97-11eb-1249-03382bfe3feb
begin
	m11_1_1s = SampleModel("m11.1.1s", stan11_1_1)
	rc11_1_1s = stan_sample(m11_1_1s)
end;

# ╔═╡ dd96b924-8d97-11eb-0d47-ffa72b39a875
if success(rc11_1_1s)
	post11_1_1s_df = read_samples(m11_1_1s, :dataframe)
	PRECIS(post11_1_1s_df)
end

# ╔═╡ 8436e374-8e3a-11eb-0e48-d10c66f248a4
if success(rc11_1_1s)
	p1 = logistic.(post11_1_1s_df.a)
	mean(p1)
end

# ╔═╡ dce2c76c-8e3b-11eb-1b40-a739dda46d0b
stan11_1_2 = "
parameters {
 	real a;
}
model {
	a ~ normal(0, 1.5);
}
";

# ╔═╡ 61c85042-8daa-11eb-392e-45b1d7adb978
begin
	m11_1_2s = SampleModel("m11.1.2s", stan11_1_2)
	rc11_1_2s = stan_sample(m11_1_2s)
end;

# ╔═╡ dc1a5f76-8e3a-11eb-1a9d-4374c4fadbc2
if success(rc11_1_2s)
	post11_1_2s_df = read_samples(m11_1_2s, :dataframe)
	PRECIS(post11_1_2s_df)
end

# ╔═╡ dc1a8a32-8e3a-11eb-0173-11ea33425350
if success(rc11_1_2s)
	p2 = logistic.(post11_1_2s_df.a)
end

# ╔═╡ dc1b3d88-8e3a-11eb-3167-0726ac404587
if success(rc11_1_1s) && success(rc11_1_2s)
	density(p1, lab="p1")
	density!(p2, lab="p2")
end

# ╔═╡ 83a05b7a-8e30-11eb-23f4-cdf5fa972dd4
stan11_2 = "
data {
  int N;
  int K;
  int<lower=1, upper=4> treatment[N];
  //int<lower=0,upper=1> pulled_left[N];
}
parameters {
  real a;
  vector[K] b;
}
model {
	a ~ normal(0, 1.5);
    b ~ normal(0, 0.5);
 	//pulled_left ~ binomial_logit(a + b[treatment]);
}
";

# ╔═╡ dcd3f5b0-8e30-11eb-1a78-451a3437b6d6
begin
	m11_2s = SampleModel("m11.2s", stan11_2)
	data3 = (N = size(df, 1), K = length(unique(df.treatment)),
		pulled_left = df.pulled_left, treatment = df.treatment)
	rc11_2s = stan_sample(m11_2s; data=data3)
end;

# ╔═╡ dcd4235c-8e30-11eb-2130-a5027eb07957
if success(rc11_2s)
	post11_2s_df = read_samples(m11_2s, :dataframe)
	PRECIS(post11_2s_df)
end

# ╔═╡ dcd4dc46-8e30-11eb-33e4-67633257f5df
begin
	nt11_2s = read_samples(m11_2s, :namedtuple)
	p3 = [logistic.(nt11_2s.a .+ nt11_2s.b[i]) for i in 1:4]
	mean.([p3[i] for i in 1:4])
end

# ╔═╡ dce32f32-8e30-11eb-0db6-670aa2e9e017
mean(abs.(p3[1] - p3[2]))

# ╔═╡ dce3ebc0-8e30-11eb-0c16-53439160c3b9
begin
	density(p1, ylims=(0, 10), lab="p1")
	density!(p2, lab="p2")
	density!(abs.(p3[1] - p3[2]), lab="p3")
end

# ╔═╡ b6a6c8e6-8e31-11eb-3276-2d3f9297fe70
md" ## End of clip-11-01-08s.jl"

# ╔═╡ Cell order:
# ╟─a3565770-8e31-11eb-1635-3b842d541ab6
# ╠═79789782-8d97-11eb-112c-fd8c2c294527
# ╠═dd81850e-8d97-11eb-0cff-8bf4f38f0660
# ╠═dd81af5c-8d97-11eb-10e4-73daebcd8875
# ╠═691258d8-8dab-11eb-3ec4-8b2059af0832
# ╠═dd82254a-8d97-11eb-3f3d-59b7642304f5
# ╠═dd8c8ce4-8d97-11eb-1249-03382bfe3feb
# ╠═dd96b924-8d97-11eb-0d47-ffa72b39a875
# ╠═8436e374-8e3a-11eb-0e48-d10c66f248a4
# ╠═dce2c76c-8e3b-11eb-1b40-a739dda46d0b
# ╠═61c85042-8daa-11eb-392e-45b1d7adb978
# ╠═dc1a5f76-8e3a-11eb-1a9d-4374c4fadbc2
# ╠═dc1a8a32-8e3a-11eb-0173-11ea33425350
# ╠═dc1b3d88-8e3a-11eb-3167-0726ac404587
# ╠═83a05b7a-8e30-11eb-23f4-cdf5fa972dd4
# ╠═dcd3f5b0-8e30-11eb-1a78-451a3437b6d6
# ╠═dcd4235c-8e30-11eb-2130-a5027eb07957
# ╠═dcd4dc46-8e30-11eb-33e4-67633257f5df
# ╠═dce32f32-8e30-11eb-0db6-670aa2e9e017
# ╠═dce3ebc0-8e30-11eb-0c16-53439160c3b9
# ╟─b6a6c8e6-8e31-11eb-3276-2d3f9297fe70


text:
module PreSaveSanitizers

using ..Ahorn, Maple
using Gtk, Gtk.ShortNames

function sortRoomNames(side::Maple.Side, config::Ahorn.Config)
    map = side.map
    sortRoomNames = get(config, "sort_rooms_on_save", true)

    if sortRoomNames
        sort!(map.rooms, by=r -> r.name)
        Ahorn.updateTreeView!(Ahorn.roomList, Ahorn.getTreeData(map), row -> row[1] == Ahorn.loadedState.roomName, updateByReplacement=true)
    end
end

function warnMissingDecals(side::Maple.Side, config::Ahorn.Config)
    map = side.map
    warnMissing = get(config, "warn_missing_decals_on_save", true)

    if warnMissing
        Ahorn.loadChangedExternalSprites!()

        atlas = Ahorn.getAtlas("Gameplay")
        textures = Ahorn.spritesToDecalTextures(atlas)

        missingDecals = Tuple{Maple.Room, Array{Maple.Decal, 1}, Maple.Decal}[]
        existingDecalNames = Dict{String, Bool}()

        for room in map.rooms
            for decals in Array{Maple.Decal, 1}[room.fgDecals, room.bgDecals]
                for decal in decals
                    if get(existingDecalNames, decal.texture, false)
                        continue
                    end

                    texture = Ahorn.fixTexturePath(decal.texture)

                    if !(texture in textures || splitext(texture)[1] in textures)
                        push!(missingDecals, (room, decals, decal))

                    else
                        existingDecalNames[decal.texture] = true
                    end
                end
            end
        end

        missingCount = length(missingDecals)
        plural = missingCount == 1 ? "" : "s"

        if missingCount > 0
            uniqueNames = unique(String[decal.texture for (room, decals, decal) in missingDecals])

            dialogText = "You have $missingCount decal$plural with no longer existing texture$plural.\n" *
                "These can potentially crash Celeste when their room is loaded.\n" *
                "Do you want to automatically delete these decals?\n\n" *
                join(uniqueNames, "\n")
            confirmed = ask_dialog(dialogText, Ahorn.window)

            if confirmed
                Ahorn.History.addSnapshot!(Ahorn.History.MapSnapshot("Deleting missing decals", map))

                for (room, decals, decal) in missingDecals
                    filter!(d -> !(d.texture in uniqueNames), decals)
                end
            end
        end
    end
end

function warnDuplicateIds(side::Maple.Side, config::Ahorn.Config)
    map = side.map
    warnDuplicate = get(config, "warn_duplicate_entity_ids_on_save", true)

    if warnDuplicate
        idGroupDict = Ahorn.EntityIds.groupByIds(map)
        duplicateIndices = Int[]
        duplicateCount = 0

        for (id, targets) in idGroupDict
            if length(targets) > 1
                push!(duplicateIndices, id)
                duplicateCount += length(targets)
            end
        end

        if duplicateCount > 0
            plural1 = duplicateCount == 1 ? "" : "s"
            plural2 = duplicateCount == 1 ? "y" : "ies"
            dialogText = "You have $duplicateCount entit$plural2/trigger$plural1 with duplicate ID$plural1.\n" *
                "This can cause entities or triggers to get disabled when they shouldn't be, or other weird behaviour in Celeste.\n" *
                "Do you want to automatically reassign IDs?"
            confirmed = ask_dialog(dialogText, Ahorn.window)

            if confirmed
                Ahorn.History.addSnapshot!(Ahorn.History.MapSnapshot("Reassigning entity/trigger IDs", map))

                success, reassigned = Ahorn.EntityIds.attemptIdFixing(map, idGroupDict)

                if !success
                    dialogText = "Some entities/triggers were not safe to reassign IDs to.\n" *
                        "Do you want to reassign IDs anyways? This might cause entities that depend on raw IDs to need reconfiguring."
                    confirmed = ask_dialog(dialogText, Ahorn.window)

                    if confirmed
                        Ahorn.EntityIds.attemptIdFixing(map, idGroupDict, true)
                    end
                end
            end
        end
    end
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
"""
    ComplexD64

A complex value formed or Double64 real and imaginary parts
Also a constructor for that.
"""

struct ComplexD64 <: Number
    hilo::Tuple{Complex{Float64}, Complex{Float64}}
end

ComplexD64(x::ComplexD64) = x
ComplexD64(x::Double64) = ComplexD64((x, ZeroD64))
ComplexD64(x::Float64) = ComplexD64((Double64(x), ZeroD64)))

function ComplexD64(hi::Complex{Float64}, lo::Complex{Float64})
    rehi, relo = two_sum(real(hi), real(lo))
    imhi, imlo = two_sum(imag(hi), imag(lo))
    hi = ComplexF64(rehi, imhi)
    lo = ComplexF64(relo, imlo)
    ComplexD64((hi, lo))
end

function ComplexD64(real::Double64, imag::Double64)
    rehi, relo = HiLo(real)
    imhi, imlo = HiLo(imag)
     hi = ComplexF64(rehi, imhi)
    lo = ComplexF64(relo, imlo)
    ComplexD64((hi, lo))
end

Hi(x::ComplexD64) = x.hilo[1]
Lo(x::ComplexD64) = x.hilo[2]
HiLo(x::ComplexD64) = x.hilo

Base.real(x::ComplexD64) = Double64(real(hi(x)), real(lo(x)))
Base.imag(x::ComplexD64) = Double64(imag(hi(x)), imag(lo(x)))

ComplexD64(real::Double64, imag::Double64) =
    ComplexD64(Complex{Float64}(Hi(real), Hi(imag)), Complex(Float64(Lo(real), Lo(imag)))

const CplxD64 = Complex{Double64}

Hi(x::ComplexD64) = ComplexF64(Hi(x.re), Hi(x.im))
Lo(x::ComplexD64) = ComplexF64(Lo(x.re), Lo(x.im))
HiLo(x::ComplexD64) = Hi(x), Lo(x)

"""
    FloatComplexD64 <: Union
Union of Double64 and ComplexD64 types.
"""
const FloatComplexD64 = Union{Double64, ComplexD64}

Base.real(x::Double64) = x
Base.imag(x::Double64) = Double64((0.0, 0.0))


text:
(::Type{CA})(::UndefInitializer) where {CA<:Caggregate} = init(bitstype(CA), Val(sizeof(CA)))
@generated function init(::Type{CA}, ::Val{size}) where {CA<:Caggregate, size}
	init = ((zero(UInt8) for _ in 1:size)...,)
	return :(CA(undef, $(init)))
end

function (::Type{CA})(init::Union{NTuple{N, UInt8} where {N}, CA, Cbinding{CA}, Cptrs{CA}, UndefInitializer} = undef; kwargs...) where {CA<:Caggregate}
	result =
		init isa UndefInitializer ? bitstype(CA)(undef) :
		init isa Tuple ? bitstype(CA)(undef, init) :
		init isa Cbinding ? init[] :
		init isa Cptrs ? init[] :
		init

	for kwarg in kwargs
		result = initproperty(result, kwarg...)
	end

	return result
end

Base.read(io::IO, ::Type{CA}) where {CA<:Caggregate} = CA(((read(io, UInt8) for _ in 1:sizeof(CA))...,))

Base.sizeof(::Type{CA}) where {CA<:Caggregate} = invoke(sizeof, Tuple{Any}, bitstype(CA))
Base.zero(::Type{CA}) where {CA<:Caggregate} = CA()
Base.convert(::Type{CA}, nt::NamedTuple) where {CA<:Caggregate} = CA(; nt...)

Base.fieldcount(::Type{CA}) where {CA<:Caggregates} = length(fieldnames(CA))
Base.fieldtypes(::Type{CA}) where {CA<:Caggregates} = ((Base.tuple_type_head(t) for (n, t) in fields(CA))...,)
Base.fieldnames(::Type{CA}) where {CA<:Caggregates} = ((n for (n, t) in fields(CA))...,)
Base.propertynames(::CA) where {CA<:Caggregates} = fieldnames(CA)

Base.getproperty(ca::CA, sym::Symbol) where {CA<:Caggregates} = field(CA, sym) === nothing ? error("type $(CA) has no field $(sym)") : load(ca, field(CA, sym))
Base.setproperty!(ca::CA, sym::Symbol, val) where {CA<:Caggregate} = error("cannot set value of field $(sym) in type $(CA)")
initproperty(ca::CA, sym::Symbol, val) where {CA<:Caggregate} = field(CA, sym) === nothing ? error("type $(CA) has no field $(sym)") : store!(ca, field(CA, sym), val)

function Base.show(io::IO, m::MIME"text/plain", ::Type{C}) where {C<:Union{Caggregates, Cenums, Carrays, Cptrs}}
	if isbitstype(C)
		# TODO: pretty-print types into c"Module.libmodule.struct S" format
		Base.show_datatype(io, C)
	else
		invoke(show, Tuple{IO, MIME"text/plain", Type}, io, m, C)
	end
end

function Base.show(io::IO, ca::Caggregate)
	show(io, typeof(ca))
	print(io, '(')
	if get(io, :compact, false)
		print(io, "...")
	else
		for (ind, name) in enumerate(fieldnames(typeof(ca)))
			ind > 1 && print(io, ", ")
			print(io, name, "=")
			show(io, getproperty(ca, name))
		end
	end
	print(io, ')')
end

function Base.show(io::IO, m::MIME"text/plain", ca::Caggregate)
	show(io, typeof(ca))
	println(io, " (", sizeof(ca), "-byte ", ca isa Cstruct ? "struct)" : "union)")

	numfields = length(fieldnames(typeof(ca)))
	maxwidth = 0
	for (ind, name) in enumerate(fieldnames(typeof(ca)))
		if ind < numfields && getproperty(ca, name) isa Union{Caggregates, Cenums}
			maxwidth = 0
			break
		end
		maxwidth = max(maxwidth, length(String(name)))
	end

	indent = get(io, :indent, 1)
	prefix = repeat("  ", indent)
	sub = IOContext(io, :indent => indent+1)
	for (ind, name) in enumerate(fieldnames(typeof(ca)))
		field = getproperty(ca, name)
		print(sub, prefix, rpad(String(name), maxwidth), " = ")
		field isa Union{Caggregates, Cenums} ? show(sub, m, field) : show(sub, field)
		ind < numfields && println(sub)
	end
end


text:
using FFTW: plan_fft!
using Statistics: mean, cov, var
using Distributed: workers, remotecall, fetch, myid

function linreg(x,y)
  b = cov(x,y)/var(x)
  a = mean(y) - b*mean(x)
  a,b
end

function detrendTS!(outar::AbstractMatrix,xin::AbstractVector{T}) where T
    x=collect(T,eachindex(xin))
    a,b=linreg(x,xin)
    for i in eachindex(xin)
        outar[i,1]=a+b*x[i]
    end
end

function tscale2ind(b::Float64,l::Int)
    i=round(Int,l/b+1)
    return i
end
mirror(i,l)=l-i+2

"""
    filterTSFFT(c)

Filter each time series using a Fourier filter and return the decomposed series
in 4 time windows (Trend, Long-Term Variability, Annual Cycle, Fast Oscillations)

**Input Axis** `Time`axis

**Output Axes** `Time`axis, `Scale`axis

"""
function filterTSFFT(c;kwargs...)
  indims = InDims("time",filter=AnyMissing())
  outdims = OutDims("time",(c,p)->CategoricalAxis("Scale",["Trend", "Long-Term Variability", "Annual Cycle", "Fast Oscillations"]))
  ntime = length(getAxis("Time",c))
  plans = map(workers()) do id
    remotecall(id,ntime,eltype(c)) do nt,et
      testar = zeros(Complex{Base.nonmissingtype(et)},nt)
      fftplan = plan_fft!(testar)
      ifftplan = inv(fftplan)
      fftplan, ifftplan
    end
  end
  plandict = Dict(zip(workers(), plans))
  mapCube(filterTSFFT,c,getNpY(c),plandict;indims=indims,outdims=outdims,kwargs...)
end

function filterTSFFT(outar::AbstractMatrix,y::AbstractVector, annfreq::Number,
  plandict; nharm::Int=3)

  fftplan, ifftplan = fetch(plandict[myid()])

    detrendTS!(outar,y)
    l        = length(y)

    fy       = Complex{Base.nonmissingtype(eltype(y))}[y[i]-outar[i,1] for i=1:l]
    fftplan * fy #Apply fftplan on fy inplace changes fy
    fyout    = similar(fy)
    czero    = zero(eltype(fy))

    #Remove annual cycle
    fill!(fyout,zero(eltype(fyout)))

    for jharm = 1:nharm
        iup=tscale2ind(annfreq*1.1/jharm,l)

        idown=tscale2ind(annfreq*0.9/jharm,l)

        for i=iup:idown
            fyout[i]  = fy[i]
            i2        = mirror(i,l)
            fyout[i2] = fy[i2]
            fy[i]     = czero
            fy[i2]    = czero
        end
    end

    ifftplan * fyout #Apply ifftplan on fyout inplace, changes fyout
    for i=1:l
        outar[i,3]=real(fyout[i])
    end

    iup   = tscale2ind(annfreq*1.1,l)
    idown = tscale2ind(annfreq*0.9,l)
    #Now split the remaining parts
    fill!(fyout,czero)
    for i=2:iup-1
        i2        = mirror(i,l)
        fyout[i]  = fy[i]
        fyout[i2] = fy[i2]
        fy[i]     = czero
        fy[i2]    = czero
    end
    ifftplan * fyout
    ifftplan * fy

    for i=1:l
        outar[i,2]=real(fyout[i])
        outar[i,4]=real(fy[i])
    end
end


text:
"""
    start_parameter_table(model; parameter_table)

Return a vector of starting values taken from `parameter_table`.
"""
function start_parameter_table end

# splice model and loss functions
function start_parameter_table(model::AbstractSemSingle; kwargs...)
    return start_parameter_table(
        model.observed,
        model.imply,
        model.optimizer,
        model.loss.functions...;
        kwargs...)
end

# RAM(Symbolic)
function start_parameter_table(observed, imply, optimizer, args...; kwargs...)
    return start_parameter_table(
        ram_matrices(imply);
        kwargs...)
end

function start_parameter_table(ram_matrices::RAMMatrices; parameter_table::ParameterTable, kwargs...)

    start_val = zeros(0)

    for identifier_ram in ram_matrices.parameters
        found = false
        for (i, identifier_table) in enumerate(parameter_table.identifier)
            if identifier_ram == identifier_table
                push!(start_val, parameter_table.start[i])
                found = true
                break
            end
        end
        if !found throw(ErrorException("At least one parameter could not be found in the parameter table.")) end
    end

    return start_val

end


text:
"""
Transient heating of a rectangular steel block suddenly subject to surface temperature of 100° C.

The initial temperature of the block is 20° C. The temperature at the center
climbs from 20° C to approximately 90° C in 1200 seconds.
"""

module steel_block
using FinEtools
using FinEtoolsHeatDiff
using FinEtoolsHeatDiff.AlgoHeatDiffModule
import LinearAlgebra: cholesky, mul!
using Gnuplot

function steel_block_m()
	# Thermal conductivity
	thermal_conductivity =  [i==j ? 44.0*phun("W/K/m") : zero(FFlt) for i=1:3, j=1:3]; # conductivity matrix
	# Specific heat: the specific heat is per unit volume
	rho = 8000*phun("kg/m^3")
	specific_heat = 470.0*phun("J/kg/K")*rho
	theta = 1.0; # generalized trapezoidal method parameter
	dt = 1.0; # time step
	tend = 1200; # length of the time interval
	T0(x, y, z) = 20.0;# Initial distribution of temperature
	Ts = 100.0 # Prescribed surface temperature
	dimensions = (0.2, 0.15, 0.3).*phun("m")
	nels = (10, 10, 10)

	tolerance =min(dimensions...)/max(nels...)/100;

	fens,fes = H8block(dimensions..., nels...)
	@show lcenter  = selectnode(fens; nearestto=[0.0 0.0 0.0],  inflate = tolerance)

	geom = NodalField(fens.xyz)
	Temp = NodalField(zeros(size(fens.xyz,1),1))

	# Initial temperature
	for i in 1:count(fens)
		Temp.values[i] = T0(fens.xyz[i, :]...)
	end

	lx  = selectnode(fens; box=[dimensions[1] dimensions[1] -Inf Inf -Inf Inf],  inflate = tolerance)
	ly  = selectnode(fens; box=[-Inf Inf dimensions[2] dimensions[2] -Inf Inf],  inflate = tolerance)
	lz  = selectnode(fens; box=[-Inf Inf -Inf Inf dimensions[3] dimensions[3]],  inflate = tolerance)
	for (i, l) in enumerate([lx, ly, lz])
		setebc!(Temp, l, true, 1, Ts)
	end
	applyebc!(Temp)
    numberdofs!(Temp)

    centerdof = Temp.dofnums[lcenter]

    material = MatHeatDiff(thermal_conductivity, specific_heat)
    femm = FEMMHeatDiff(IntegDomain(fes, GaussRule(3, 3)), material)

    K = conductivity(femm, geom, Temp)
    C = capacity(femm, geom, Temp)
    L = nzebcloadsconductivity(femm, geom, Temp)

    A = cholesky((1/dt*C+theta*K))
    B = (1/dt*C-(1-theta)*K)

    Tn = gathersysvec(Temp)
    Tn1 = deepcopy(Tn)

    tstart = time()
    t = 0.0
    ts = Float64[t]
    @show Tn1[centerdof, 1][1]
    center_T = Float64[Tn1[centerdof, 1][1]];
    while true # Time stepping
    	t = t+dt; # Advance the current time
    	Tn .= Tn1 # current temperature system vector

    	# Solve for the new temperature  vector
    	rhs = B*Tn + L; # complete right-hand side (loads)
    	Tn1 .= A\rhs; # compute next temperature
    	push!(ts, t)
    	push!(center_T, Tn1[centerdof, 1][1])
    	if (t == tend) # Have we reached the end?  If so jump out.
	    	break;
	    end
	    if (t+dt > tend) # Adjust the last time step so that we exactly reach tend
	    	dt = tend-t;
	    end
	end
	println("Time = $(time() - tstart)")

	@show minimum(center_T), maximum(center_T)
	@gp "set terminal windows 0 " :-
	@gp  :- vec(ts) vec(center_T) "lw 2 lc rgb 'red' with lines title 'Temperature at the corner' "
	@gp  :- "set xlabel 'Time [s]'"
	@gp  :- "set ylabel 'Temperature [degrees C]'"

    true
end

end # module steel_block

using .steel_block
steel_block.steel_block_m()

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
@testset "test_demosystems" begin

# Just some very simple tests of systemdepot
@test size(DemoSystems.woodberry()) == (2,2)
@test size(DemoSystems.fotd()) == (1,1)
@test freqresp(DemoSystems.fotd(τ=2, T=5), [0.0])[:] == [1.0]

Random.seed!(10)

sys = ssrand(1,proper=true)
@test size(sys) == (1,1)
@test ControlSystems.nstates(sys) == 2
@test iszero(sys.D)

@test all(isstable, ssrand(1,1,2,stable=true) for _=1:100)
@test all(isstable, ssrand(1,1,2,stable=true,Ts=0.1) for _=1:100)
@test all(isstable, ssrand(1,1,2,stable=true,Ts=10) for _=1:100)

sys = ssrand(2,2,5,proper=false,stable=true)
@test size(sys) == (2,2)
@test ControlSystems.nstates(sys) == 5
@test isstable(sys)
@test !iszero(sys.D)

Random.seed!(20)
sys = ssrand(2,2,5,stable=false)
@test !isstable(sys)

sys = ssrand(2,2,5,proper=false,stable=true, Ts=0.5)
@test size(sys) == (2,2)
@test ControlSystems.nstates(sys) == 5
@test isstable(sys)
@test !iszero(sys.D)

sys = ssrand(2,2,5,proper=false,stable=false, Ts=0.5)
@test !isstable(sys)

end


text:
# This file is an attempt to walk two named tuples simultaneously. It's very
# rough, maybe not even usable in its current state. It's not yet clear whether
# this approach will have any advantages over others like `LazyMerge`.

function merge(x::Nested, y::Nested)
    _merge_ordered(x,y)
end

function _merge_ordered(x::NamedTuple{Kx,Vx}, y::NamedTuple{Ky,Vy}) where {Kx,Vx,Ky,Vy}

end

function walk2(f, xs::Tx, ys::Tx) where {Tx<:Tuple, Ty<:Tuple}
    (new_x, new_xs) = untuple(xs)
    (new_y, new_ys) = untuple(ys)
    return walk2(f, new_x, new_y, new_xs, new_ys, ())
end

function walk2(f, x::Nothing, y, xs::Tuple{}, ys, acc)
    new_acc = f(y, acc)
    return foldl(f, ys; init=new_acc)
end

function walk2(f, x, y::Nothing, xs, ys::Tuple{}, acc)
    new_acc = f(x, acc)
    return foldl(f, xs; init=new_acc)``
end

function walk2(f, x::Nothing, y::Nothing, xs::Tuple{}, ys::Tuple{}, acc)
    return acc
end

function walk2(f, x, y, xs, ys, acc)
    new_x = x
    new_xs = xs
    new_y = y
    new_ys = ys

    if x==y
        new_acc = f(x, acc)
        (new_x, new_xs) = untuple(xs)
        (new_y, new_ys) = untuple(ys)
    elseif x < y
        new_acc = f(x, acc)
        (new_x, new_xs) = untuple(xs)
    else
        new_acc = f(y, acc)
        (new_y, new_ys) = untuple(ys)
    end

    return walk2(f, new_x, new_y, new_xs, new_ys, new_acc)
end

function untuple(x::Tuple{})
    return (nothing, ())
end

function untuple(x::NTuple{N,T}) where {N,T}
    f(h,t...) = (h,t)
    return f(x...)
end

using BangBang

@inline function slice_namedtuple(nt::NamedTuple{K,V}, r::UnitRange) where {K,V}
    @boundscheck begin
        n = length(K)
        if r.start < 1
            throw(BoundsError(nt, 0))
        elseif r.stop > n
            throw(BoundsError(nt, n+1))
        end
    end

    ks = @inbounds keys(nt)[r]
    vs = @inbounds values(nt)[r]
    return namedtuple(ks)(vs)
end

@generated function keymerge(x::NamedTuple{Nx}, y::NamedTuple{Ny}) where {Nx, Ny}
    Tuple(sort(union(Nx, Ny)))
end

@inline function mergent(xs, ys)
    m = length(xs)
    n = length(ys)

    xkeys = keys(xs)
    xvals = values(xs)

    ykeys = keys(ys)
    yvals = values(ys)

    i = 1
    j = 1

    ks = keymerge(xs, ys)

    vs = ()
    while true
        # We're done with xs
        if i > m
            vs = append!!(vs, @inbounds yvals[j:n])
            break
        end

        # We're done with ys
        if j > n
            vs = append!!(vs, @inbounds xvals[i:m])
            break
        end

        xkey = @inbounds xkeys[i]
        ykey = @inbounds ykeys[j]

        if xkey == ykey
            xval = @inbounds xvals[i]
            yval = @inbounds yvals[j]

            vs = push!!(vs, mergent(xval, yval))

            i += 1
            j += 1
        elseif xkey < ykey
            xval = @inbounds xvals[i]
            vs = push!!(vs, xval)
            i += 1
        else
            yval = @inbounds yvals[i]
            vs = push!!(vs, yval)
            j += 1
        end
    end

    return namedtuple(ks)(vs)
end


text:
using JSON
using ZMQ

const ctx = Context()
sock = Socket(ctx, REP)
ZMQ.bind(sock, "tcp://*:7788")

function triple(x)
	sleep(10)
	3x
end

while true
	println("Server running.")

	msg = JSON.parse(bytestring(ZMQ.recv(sock)))

	result =""
	try
		result = triple(float(msg["value"]))
	catch
		result = "julia function returned an error"
	end
	println(result)

	ZMQ.send(sock, JSON.json({"result"=>result}))
end


text:
# export jlmie_Isff, jlmie_Isff_EM_n

function jlmie_Isff_EM_n(nmat,radius,lbd0,nenv,theta,phi,n)
# calculate far-feild intensity of scattered field
# from n-th order Electric and Magnetic multipoles separately
    m = nmat./nenv
    lbd = lbd0./nenv  # effective wavelength in the environment
    k = 2*pi./lbd  # free space wavenumber in the environment
    x = k.*radius
    mu = cos.(theta)
    an,bn,_,_ = jlmie_abcd(m,x,n)
    pin,taun = jlmie_pt(mu,n)
    S1EDn = (2*n+1)/(n*(n+1)).*(an.*pin)
    S2EDn = (2*n+1)/(n*(n+1)).*(an.*taun)
    S1MDn = (2*n+1)/(n*(n+1)).*(bn.*taun)
    S2MDn = (2*n+1)/(n*(n+1)).*(bn.*pin)
    # settings
    E0 = 1  # amplitude of incident electric field
    r = 1  # distance from the sphere to an observation point (1 m)
    # main
    Est_EDn = E0 * exp.(1im.*k*r)./(-1im.*k*r) .* cos.(phi).*S2EDn
    Esp_EDn = -1*E0 * exp.(1im.*k*r)./(-1im.*k*r) .* sin.(phi).*S1EDn
    Est_MDn = E0 * exp.(1im.*k*r)./(-1im.*k*r) .* cos.(phi).*S2MDn
    Esp_MDn = -1*E0 * exp.(1im.*k*r)./(-1im.*k*r) .* sin.(phi).*S1MDn
    Is_EDn = abs.(Est_EDn).^2 .+ abs.(Esp_EDn).^2
    Is_MDn = abs.(Est_MDn).^2 .+ abs.(Esp_MDn).^2
    return Is_EDn, Is_MDn
end

function jlmie_Isff(nmat,radius,lbd0,nenv,theta,phi,nmax::Int=-1)
# calculate far-field intensity of scattered light
# (as a total from n = 1 to nmax)
    if nmax == -1
        nmax = largenmax(jlmie_mx(nmat,radius,lbd0,nenv)[2])
    end

    Est,Esp = jlmie_Esff_n(nmat,radius,lbd0,nenv,theta,phi,1)
    if nmax > 1
        for n = 2:nmax
            Est_n,Esp_n = jlmie_Esff_n(nmat,radius,lbd0,nenv,theta,phi,n)
            Est = Est .+ Est_n
            Esp = Esp .+ Esp_n
        end
    end
    Ist = abs.(Est).^2
    Isp = abs.(Esp).^2
    Is = Ist .+ Isp
    return Is,Ist,Isp
end

function jlmie_Esff_n(nmat,radius,lbd0,nenv,theta,phi,n::Int)
# calculate scattered electric far-fields E_s,theta and E_s,phi
# resulting from n-th order resonances
# See Sec. 4.4.4 of B&H for derivation
# Input:
# nmat: n_material; refractive index of the material of the sphere
# radius: radius of the sphere
# lbd0: target vacuum wavelength (range)
# nenv: n_environement; refractive index of the environment
#       around the sphere
    # preparation
    m = nmat./nenv
    lbd = lbd0./nenv  # effective wavelength in the environment
    k = 2*pi./lbd  # free space wavenumber in the environment
    x = k.*radius
    mu = cos.(theta)
    S1n,S2n = jlmie_S12(m,x,mu,n)

    # settings
    E0 = 1  # amplitude of incident electric field
    r = 1  # distance from the sphere to an observation point (1 m)

    # main
    Est_n = E0 * exp.(1im.*k*r)./(-1im.*k*r) .* cos.(phi).*S2n
    Esp_n = -1*E0 * exp.(1im.*k*r)./(-1im.*k*r) .* sin.(phi).*S1n
    return Est_n,Esp_n
end

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
# %% Scheduler tests
@testset "Scheduler tests" begin
@testset "Simple Schedulers" begin
    N = 1000

    # by_id
    model = ABM(Agent0; scheduler = by_id)
    for i in 1:N
        add_agent!(model)
    end
    @test sort!(collect(keys(model.agents))) == 1:N
    @test model.scheduler(model) == 1:N

    # fastest
    Random.seed!(12)
    model = ABM(Agent0; scheduler = fastest)
    for i in 1:N
        add_agent!(model)
    end
    @test sort!(collect(model.scheduler(model))) == 1:N

    # random
    Random.seed!(12)
    model = ABM(Agent0; scheduler = random_activation)
    for i in 1:N
        add_agent!(model)
    end
    @test model.scheduler(model)[1:3] == [913, 522, 637] # reproducibility test

    # partial
    Random.seed!(12)
    model = ABM(Agent0; scheduler = partial_activation(0.1))
    for i in 1:N
        add_agent!(model)
    end

    a = model.scheduler(model)
    @test length(a) < N
    @test a[1] == 74 # reproducibility test

    # by property
    model = ABM(Agent2; scheduler = property_activation(:weight))
    for i in 1:N
        add_agent!(model, rand() / rand())
    end

    Random.seed!(12)
    a = model.scheduler(model)

    ids = collect(keys(model.agents))
    properties = [model.agents[id].weight for id in ids]

    @test ids[sortperm(properties)] == a
end

@testset "Union Types" begin
    @test Agents.union_types(Union{Agent0}) == (Agent0,)
    @test Agents.union_types(Union{Agent0,Agent1}) == (Agent0, Agent1)
    @test Agents.union_types(Union{Agent0,Agent1,Agent2,Agent3}) ==
          (Agent0, Agent1, Agent2, Agent3)
    #Union types are not order preserving
    @test Agents.union_types(Union{Agent1,Agent0}) == (Agent0, Agent1)
    @test Agents.union_types(Union{Agent1,Agent3,Agent2,Agent0}) ==
          (Agent0, Agent1, Agent2, Agent3)
end

# Mixed model
function init_mixed_model(choices = [3, 3, 3, 3]; scheduler = fastest)
    model = ABM(Union{Agent0,Agent1,Agent2,Agent3}, scheduler = scheduler, warn = false)
    atypes = (Agent0,Agent1,Agent2,Agent3)
    id = 1
    for i in 1:choices[1]
        a0 = Agent0(id)
        add_agent!(a0, model)
        id += 1
    end
    for i in 1:choices[2]
        a1 = Agent1(id, (0, 0))
        add_agent!(a1, model)
        id +=1
    end
    for i in 1:choices[3]
        a2 = Agent2(id, 5.0)
        add_agent!(a2, model)
        id += 1
    end
    for i in 1:choices[4]
        a3 = Agent3(id, (0, 0), 5.0)
        add_agent!(a3, model)
        id += 1
    end
    return model
end

@testset "Mixed Scheduler" begin
    # standard scheduler
    model = init_mixed_model()
    @test sort!(collect(model.scheduler(model))) == 1:12

    # shuffling types scheduler
    Random.seed!(12)
    model = init_mixed_model(scheduler = by_type(true, false))
    s1 = model.scheduler(model)
    s2 = model.scheduler(model)
    @test unique([typeof(model[id]) for id in s1]) != unique([typeof(model[id]) for id in s2])
    @test count(model[id] isa Agent2 for id in model.scheduler(model)) == 3
    c = begin
        x = 0; s = model.scheduler(model)
        x += count(a -> a == Agent0, typeof(model[id]) for id in s)
        x += count(a -> a == Agent1, typeof(model[id]) for id in s)
        x += count(a -> a == Agent2, typeof(model[id]) for id in s)
        x += count(a -> a == Agent3, typeof(model[id]) for id in s)
    end
    @test c == 12

    # NOT shuffling types scheduler
    Random.seed!(12)
    model = init_mixed_model(scheduler = by_type(false, false))
    s1 = model.scheduler(model)
    s2 = model.scheduler(model)
    @test unique([typeof(model[id]) for id in s1]) == unique([typeof(model[id]) for id in s2])

    # Not shuffling types, but shuffling agents
    Random.seed!(12)
    model = init_mixed_model(scheduler = by_type(false, true))
    s1 = model.scheduler(model)
    s2 = model.scheduler(model)
    @test [typeof(model[id]) for id in s1] == [typeof(model[id]) for id in s2]
    # here we actually check whether agents of same type are shuffled
    @test model[s1[1]].id ≠ model[s2[1]] || model[s1[2]].id ≠ model[s2[2]]

    # Explicit order of types scheduling
    Random.seed!(12)
    model =
        ABM(Union{Agent1,Agent0}, scheduler = by_type((Agent1, Agent0), true), warn = false)
    for id in 1:3
        a1 = Agent1(id, (0, 0))
        add_agent!(a1, model)
    end
    for id in 4:6
        a0 = Agent0(id)
        add_agent!(a0, model)
    end
    s = model.scheduler(model)
    @test [typeof(model[id]) for id in s] ==
          [Agent1, Agent1, Agent1, Agent0, Agent0, Agent0]
    @test all(x -> x < 4, s[1:3])
    @test all(x -> x > 3, s[4:6])
end

@testset "Scheduler as struct" begin
    mutable struct MyScheduler
        n::Int # step number
        w::Float64
    end
    function (ms::MyScheduler)(model::ABM)
        ms.n += 1 # increment internal counter by 1 for each step
        if ms.n < 5
            return keys(model.agents) # order doesn't matter in this case
        else
            ids = collect(allids(model))
            # filter all ids whose agents have `w` less than some amount
            filter!(id -> model[id].weight > ms.w, ids)
            return ids
        end
    end

    model = ABM(Agent2;
        properties = Dict{Int, Bool}(),
        scheduler = MyScheduler(0, 5.0)
    )
    for w in 1.0:10.0
        add_agent!(model, w)
    end
    for i in 1:10
        ids = sort!(collect(model.scheduler(model)))
        if i < 5
            @test ids == 1:10
        else
            @test ids == 6:10
        end
    end
end
end


text:
1:10 # 1:10 with *both ends* inclusive!

typeof(1:10) # UnitRange{Int64}

collect(1:10) # 10-element Array{Int64,1}: 1, ..., 10


text:
# Autogenerated wrapper script for XML2_jll for x86_64-w64-mingw32
export libxml2, xmlcatalog, xmllint

using Zlib_jll
using Libiconv_jll
JLLWrappers.@generate_wrapper_header("XML2")
JLLWrappers.@declare_library_product(libxml2, "libxml2-2.dll")
JLLWrappers.@declare_executable_product(xmlcatalog)
JLLWrappers.@declare_executable_product(xmllint)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Libiconv_jll)
    JLLWrappers.@init_library_product(
        libxml2,
        "bin\\libxml2-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        xmlcatalog,
        "bin\\xmlcatalog.exe",
    )

    JLLWrappers.@init_executable_product(
        xmllint,
        "bin\\xmllint.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module VLKeggSDK

# include -
include("Include.jl")

# type exports -
export KEGGReaction
export KEGGOrgansim
export KEGGSequence
export KEGGCompound
export KEGGPathway

# method exports -
export get_ec_number_for_gene
export get_ec_numbers_for_pathway
export get_genes_in_organism_pathway
export get_pathways_for_organism
export get_kegg_organism_codes
export get_sequence_for_gene
export get_sequence_for_protein
export get_compound_records
export get_compound_records_for_reaction
export get_reactions_for_ec_number
export get_reaction_for_rn_number
export build_metabolite_keggid_matching_table

export check
export extract_db_file_section
export extract_metabolite_symbols
export extract_stoichiometric_dictionary
export extract_atom_dictionary
export build_reaction_table_from_vff_file

export http_get_call_with_url

end # module


text:
include("learners/learners.jl")
include("policies/policies.jl")
include("approximators/approximators.jl")
include("explorers/explorers.jl")
include("trajectories/trajectories.jl")
include("preprocessors.jl")
include("agents/agents.jl")


text:
"""
    mutable struct SnapshotHeader( [ **Fields ])

Contains the data of the `HEAD` block of a Gadget snapshot.

# Fields
| Name                                 | Meaning                                                                                |
| :-------------------------------     | :------------------------------------------------------------------------------------- |
| `npart::Vector{Int32}`               | an array of particle numbers per type in this snapshot                                 |
| `massarr::Vector{Float64}`           | an array of particle masses per type in this snapshot - if zero: MASS block present    |
| `time::Float64`                      | time / scale factor of the simulation                                                  |
| `z::Float64`                         | redshift of the simulation                                                             |
| `flag_sfr::Int32`                    | 1 if simulation was run with star formation, else 0                                    |
| `flag_feedback::Int32`               | 1 if simulation was run with stellar feedback, else 0                                  |
| `nall::Vector{UInt32}`               | total number of particles in the simulation                                            |
| `flag_cooling::Int32`                | 1 if simulation was run with cooling, else 0                                           |
| `num_files::Int32`                   | number of snapshots over which the simulation is distributed                           |
| `omega_0::Float64`                   | Omega matter                                                                           |
| `boxsize::Float64`                   | total size of the simulation box                                                       |
| `omega_l::Float64`                   | Omega dark enery                                                                       |
| `h0::Float64`                        | little h                                                                               |
| `flag_stellarage::Int32`             | 1 if simulation was run with stellar age, else 0                                       |
| `flag_metals::Int32`                 | 1 if simulation was run with metals, else 0                                            |
| `npartTotalHighWord::Vector{UInt32}` | weird                                                                                |
| `flag_entropy_instead_u::Int32`      | 1 if snapshot U field contains entropy instead of internal energy, else 0              |
| `flag_doubleprecision::Int32`        | 1 if snapshot is in double precision, else 0                                           |
| `flag_ic_info::Int32`                | 1 if initial snapshot file contains an info block, else 0                              |
| `lpt_scalingfactor::Float32`         | factor to use second order ic generation                                               |
| `fill::Vector{Int32}`                | the HEAD block needs to be filled with zeros to have a size of 256 bytes               |
"""
mutable struct SnapshotHeader <: AbstractGadgetHeader
    npart::Vector{Int32}                # an array of particle numbers per type in this snapshot
    massarr::Vector{Float64}            # an array of particle masses per type in this snapshot - if zero: MASS block present
    time::Float64                       # time / scale factor of the simulation
    z::Float64                          # redshift of the simulation
    flag_sfr::Int32                     # 1 if simulation was run with star formation, else 0
    flag_feedback::Int32                # 1 if simulation was run with stellar feedback, else 0
    nall::Vector{UInt32}                # total number of particles in the simulation
    flag_cooling::Int32                 # 1 if simulation was run with cooling, else 0
    num_files::Int32                    # number of snapshots over which the simulation is distributed
    boxsize::Float64                    # total size of the simulation box
    omega_0::Float64                    # Omega matter
    omega_l::Float64                    # Omega dark enery
    h0::Float64                         # little h
    flag_stellarage::Int32              # 1 if simulation was run with stellar age, else 0
    flag_metals::Int32                  # 1 if simulation was run with metals, else 0
    npartTotalHighWord::Vector{UInt32}  # weird
    flag_entropy_instead_u::Int32       # 1 if snapshot U field contains entropy instead of internal energy, else 0
    flag_doubleprecision::Int32         # 1 if snapshot is in double precision, else 0
    flag_ic_info::Int32                 # 1 if initial snapshot file contains an info block, else 0
    lpt_scalingfactor::Float32          # factor to use second order ic generation
    fill::Vector{Int32}                 # the HEAD block needs to be filled with zeros to have a size of 256 bytes

    function SnapshotHeader(npart::Vector{Int32}=Int32.([0,0,0,0,0,0]),
           massarr::Vector{Float64}=zeros(6),
           time::Float64=0.,
           z::Float64=0.,
           flag_sfr::Int32=Int32(0),
           flag_feedback::Int32=Int32(0),
           nall::Vector{UInt32}=UInt32.([0,0,0,0,0,0]),
           flag_cooling::Int32=Int32(0),
           num_files::Int32=Int32(0),
           boxsize::Float64=0.,
           omega_0::Float64=0.,
           omega_l::Float64=0.,
           h0::Float64=0.,
           flag_stellarage::Int32=Int32(0),
           flag_metals::Int32=Int32(0),
           npartTotalHighWord::Vector{UInt32}=UInt32.([0,0,0,0,0,0]),
           flag_entropy_instead_u::Int32=Int32(0),
           flag_doubleprecision::Int32=Int32(0),
           flag_ic_info::Int32=Int32(0),
           lpt_scalingfactor::Float32=Float32(0.),
           fill::Vector{Int32}=Int32.(zeros(12)))

          new(npart,
              massarr,
              time,
              z,
              flag_sfr,
              flag_feedback,
              nall,
              flag_cooling,
              num_files,
              boxsize,
              omega_0,
              omega_l,
              h0,
              flag_stellarage,
              flag_metals,
              npartTotalHighWord,
              flag_entropy_instead_u,
              flag_doubleprecision,
              flag_ic_info,
              lpt_scalingfactor,
              fill)
    end
end

"""
    head_to_obj(filename::String)

Returns the header of a snapshot as a `SnapshotHeader` object.
"""
function head_to_obj(filename)

    h = SnapshotHeader()

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize == 8
        swap = false
        snap_format = 2
    elseif blocksize == 256
        swap = false
        snap_format = 1
    else
        blocksize = bswap(blocksize)
        if blocksize == 8
            swap = true
            snap_format = 2
        elseif blocksize == 256
            swap = true
            snap_format = 1
        else
            error("incorrect file format encountered when reading header of $filename")
        end
    end

    #println("Reading snapshot format: ", snap_format)

    if snap_format == 2
        seek(f, 16)
        skip_line = read(f, Int32)
    end

    h.npart = read!(f, Array{Int32,1}(undef,6))
    h.massarr = read!(f, Array{Float64,1}(undef,6))
    h.time = read(f, Float64)
    h.z = read(f, Float64)
    h.flag_sfr = read(f, Int32)
    h.flag_feedback = read(f, Int32)
    h.nall = read!(f, Array{UInt32,1}(undef,6))
    h.flag_cooling = read(f, Int32)
    h.num_files = read(f, Int32)
    h.boxsize = read(f, Float64)
    h.omega_0 = read(f, Float64)
    h.omega_l = read(f, Float64)
    h.h0 = read(f, Float64)
    h.flag_stellarage = read(f, Int32)
    h.flag_metals = read(f, Int32)
    h.npartTotalHighWord = read!(f, Array{UInt32,1}(undef,6))
    h.flag_entropy_instead_u = read(f, Int32)
    h.flag_doubleprecision = read(f, Int32)
    h.flag_ic_info = read(f, Int32)
    h.lpt_scalingfactor = read(f, Float32)

    close(f)

    return h
end

"""
    head_to_dict(filename::String)

Returns the header of a snapshot as a dictionary.
"""
function head_to_dict(filename::String)

        header = Dict()

        f = open(filename)
        blocksize = read(f, Int32)

        if blocksize[1] == 8
            swap = 0
            snap_format = 2
        elseif blocksize[1] == 256
            swap = 0
            snap_format = 1
        else
            blocksize[1] = bswap(blocksize[1])
            if blocksize[1] == 8
                swap = 1
                snap_format = 2
            elseif blocksize[1] == 256
                swap = 1
                snap_format = 1
            else
                println("incorrect file format encountered when reading header of", filename)
            end
        end

        if snap_format == 2
            seek(f, 16)
            skip_line = read(f, Int32)
        end

        header["snap_format"] = snap_format
        header["PartTypes"] = ["PartType0", "PartType1", "PartType2",
                               "PartType3", "PartType4", "PartType5"]
        header["npart"] = read!(f, Array{Int32,1}(undef,6))
        header["massarr"] = read!(f, Array{Float64,1}(undef,6))
        header["time"] = read(f, Float64)
        header["redshift"] = read(f, Float64)
        header["flag_sfr"] = read(f, Int32)
        header["flag_feedback"] = read(f, Int32)
        header["nall"] = read!(f, Array{UInt32,1}(undef,6))
        header["flag_cooling"] = read(f, Int32)
        header["num_files"] = read(f, Int32)
        header["boxsize"] = read(f, Float64)
        header["omega_m"] = read(f, Float64)
        header["omega_l"] = read(f, Float64)
        header["hubble"] = read(f, Float64)
        header["flag_stellarage"] = read(f, Int32)
        header["flag_metals"] = read(f, Int32)
        header["npartTotalHighWord"] = read!(f, Array{UInt32,1}(undef,6))
        header["flag_entropy_instead_u"] = read(f, Int32)
        header["flag_doubleprecision"] = read(f, Int32)
        header["flag_ic_info"] = read(f, Int32)
        header["lpt_scalingfactor"] = read(f, Float32)

        close(f)

        return header

end

"""
    read_header(filename::String)

Reads the header of a snapshot file or file base (without .0, .1, etc.)
and returns a SnapshotHeader object.

See also: [`head_to_obj`](@ref)
"""
function read_header(filename::String)
    filename = select_file(filename, 0)
    return head_to_obj(filename)
end


text:
# functions to log data

@doc """
### Utils.sharedFaceLogging

  This function writes files for the function that evalutes the shared face
  integrals.  This function should only be called in debug mode.

  If opts["writeqface"] is true,  writes the q values at the face to a file
  q_sharedface_i_myrank.dat, where i is the peer process index (not MPI rank)
  and myrank is the MPI rank of the current process.

  If opts["write_fluxface"] is true, it writes eqn.flux_sharedface to a file
  flux_sharedface_i_myrank.dat

  Inputs:
    mesh
    sbp
    eqn: and AbstractSolutionData
    opts: options dictonary
    data: SharedFaceData object
    qL_arr: an array holding solution values at face nodes on the
            left side of the interface,  numDofPerNode x numfacenodes x
            numsharedfaces on this partition boundary
    qR_arr: solution values at face nodes on the right side of the interface.
            same shape as qL_arr

  also, the eqn.flux_shared is used to write the face flux.  It is the same
  shape as qL_arr

  Aliasing restrictions: qL_arr, qR_arr, and eqn.flux_sharedface must not alias.
"""->
function sharedFaceLogging(mesh, sbp, eqn::AbstractSolutionData{Tsol},
                           opts, data::SharedFaceData, qL_arr, qR_arr) where Tsol

  if opts["writeqface"]
    myrank = data.myrank
    i = data.peeridx
    tmp_arr = zeros(Tsol, mesh.numDofPerNode, 2, mesh.numNodesPerFace, mesh.peer_face_counts[i])
    for j = 1:mesh.peer_face_counts[i]
      for k=1:mesh.numNodesPerFace
        tmp_arr[:, 1, k, j] = qL_arr[:, k, j]
        tmp_arr[:, 2, k, j] = qR_arr[:, k, j]
      end
    end
    println(eqn.params.f, "q_sharedface $i = \n", tmp_arr)
    fname = string("qsharedface_", i, "_", myrank, ".dat")
    writedlm(fname, tmp_arr)

  end  # end if

  if opts["write_fluxface"]
    for i=1:mesh.npeers
      fname = string("fluxsharedface_", i, "_", myrank, ".dat")
      writedlm(fname, eqn.flux_sharedface[i])
    end
  end

  flush(eqn.params.f)
  return nothing
end


text:
module OOPStyle

using Requires

# ref: The article posted by yatra9.
# https://qiita.com/yatra9/items/0a1a9a5ba19e9efe08c0#%E6%BC%94%E7%AE%97%E5%AD%90%E3%81%AE%E3%82%AA%E3%83%BC%E3%83%90%E3%83%BC%E3%83%AD%E3%83%BC%E3%83%89%E3%81%8Cjulia%E7%95%8C%E3%81%B8%E4%B8%8E%E3%81%88%E3%82%8B%E5%BD%B1%E9%9F%BF
Base.getproperty(obj, s::Symbol) = isdefined(obj, s) ? getfield(obj, s) : eval(quote (args...) ->  $(s)($obj, args...) end)

# Base.getproperty(obj::PyObject, s::Symbol) = isdefined(obj, s) ? getfield(obj, s) : (obj)[s]
function __init__()
    @require PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
        using PyCall
        Base.getproperty(obj::PyCall.PyObject, s::Symbol) = isdefined(obj, s) ? getfield(obj, s) : (obj)[s]
    end
    @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" isdefined(Main, :PyCall) || begin
        using PyPlot
        Base.getproperty(obj::PyPlot.PyObject, s::Symbol) = isdefined(obj, s) ? getfield(obj, s) : (obj)[s]
    end
end

end # module


text:
ns = 5
nc = 3
na = 2
y = randn(ns)
z = randn(nc)
x = randn(na)

struct MyType{T} <: AbstractVector{T}
    data :: Vector{T}
end
Base.similar(A::MyType{T}) where {T} = MyType{T}(similar(A.data))
Base.similar(A::MyType{T},::Type{S}) where {T,S} = MyType(similar(A.data,S))
Base.size(A::MyType) = size(A.data)
Base.getindex(A::MyType, i::Int) = getindex(A.data,i)
Base.setindex!(A::MyType, v, i::Int) = setindex!(A.data,v,i)
Base.IndexStyle(::MyType) = IndexLinear()

Base.BroadcastStyle(::Type{<:MyType}) = Broadcast.ArrayStyle{MyType}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MyType}},::Type{T}) where {T}
    similar(find_mt(bc),T)
end

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MyType}})
    similar(find_mt(bc))
end

find_mt(bc::Base.Broadcast.Broadcasted) = find_mt(bc.args)
find_mt(args::Tuple) = find_mt(find_mt(args[1]), Base.tail(args))
find_mt(x) = x
find_mt(::Tuple{}) = nothing
find_mt(a::MyType, rest) = a
find_mt(::Any, rest) = find_mt(rest)

@testset "Solution structure" begin

  u = solvector(state=y,constraint=z)
  @test state(u) === y && constraint(u) === z
  @test mainvector(u) === ArrayPartition(y,z)

  u = solvector(state=y,constraint=z,aux_state=x)
  @test state(u) === y
  @test constraint(u) === z
  @test aux_state(u) === x

  u = solvector(state=y)
  @test state(u) === y
  @test constraint(u) == empty(y)
  @test aux_state(u) == nothing
  @test mainvector(u) == ArrayPartition(y,empty(y))

  e = ConstrainedSystems._empty(y)
  @test isempty(e)

  u2 = solvector(state=MyType(zero(y)),constraint=zero(z),aux_state=zero(x))
  u2p = u2 .+ 1
  @test typeof(state(u2p)) <: MyType
  @test typeof(constraint(u2p)) <: typeof(z)
  @test typeof(aux_state(u2p)) <: typeof(x)
  @test state(u2p) == MyType(fill(1.0,size(y)))
  @test constraint(u2p) == fill(1.0,size(z))
  @test aux_state(u2p) == fill(1.0,size(x))

  u2m = u2p .* 2
  @test typeof(state(u2m)) <: MyType
  @test state(u2m) == MyType(fill(2.0,size(y)))
  @test constraint(u2m) == fill(2.0,size(z))
  @test aux_state(u2m) == fill(2.0,size(x))

  u2m = similar(u2p)
  u2m .= u2p .* 2 .- 3 .+ u2p ./ 2
  @test state(u2m) == MyType(fill(-0.5,size(y)))
  @test constraint(u2m) == fill(-0.5,size(z))
  @test aux_state(u2m) == fill(-0.5,size(x))

end

struct MyType1{F}
  x :: F
end

struct MyType2{F}
  x :: F
end

@testset "Solution structure with user type" begin

  a = MyType1(x)
  b = MyType2(y)
  @test isempty(ConstrainedSystems._empty(a))
  u = solvector(state=a)
  u = solvector(state=a,constraint=b)
  u = solvector(state=a,constraint=b,aux_state=a)

end

@testset "Function structure" begin

  u = solvector(state=y,constraint=z,aux_state=x)
  du = similar(u)
  fill!(du,0.0)

  function state_r1!(dy,y,p,t)
    fill!(dy,1.0)
  end

  function aux_r1!(dy,y,p,t)
    fill!(dy,2.0)
  end

  r1! = ArrayPartition((state_r1!,aux_r1!))

  r1_c! = ConstrainedSystems._complete_r1(r1!,Val(true),_func_cache=du)
  r1_c!(du,u,nothing,0.0)

  @test state(du) == fill(1.0,length(y))
  @test constraint(du) == fill(0.0,length(z))
  @test aux_state(du) == fill(2.0,length(x))

  function r2!(dz,p,t)
    fill!(dz,3.0)
  end

  fill!(du,0.0)
  r2_c! = ConstrainedSystems._complete_r2(r2!,Val(true),_func_cache=du)
  r2_c!(du,u,nothing,0.0)
  @test state(du) == fill(0.0,length(y))
  @test constraint(du) == fill(3.0,length(z))
  @test aux_state(du) == fill(0.0,length(x))

  B1 = ones(Float64,length(y),length(z))
  B2 = ones(Float64,length(z),length(y))

  function B1!(dy,z,p)
    dy .= p[1]*z
  end
  function B2!(dz,y,p)
    dz .= p[2]*y
  end

  fill!(du,0.0)
  fill!(u,1.0)
  p = [B1,B2]
  B1_c! = ConstrainedSystems._complete_B1(B1!,Val(true),_func_cache=du)
  B1_c!(du,u,p,0.0)
  @test state(du) == fill(-float(length(z)),length(y))
  @test constraint(du) == fill(0.0,length(z))
  @test aux_state(du) == fill(0.0,length(x))

  fill!(du,0.0)
  fill!(u,1.0)
  p = [B1,B2]
  B2_c! = ConstrainedSystems._complete_B2(B2!,Val(true),_func_cache=du)
  B2_c!(du,u,p,0.0)
  @test state(du) == fill(0.0,length(y))
  @test constraint(du) == fill(-float(length(y)),length(z))
  @test aux_state(du) == fill(0.0,length(x))

end

@testset "ConstrainedODEFunction" begin

  ns = 5
  nc = 2
  na = 3
  y = ones(Float64,ns)
  z = ones(Float64,nc)
  x = ones(Float64,na)

  B1 = Array{Float64}(undef,ns,nc)
  B2 = Array{Float64}(undef,nc,ns)
  B1 .= 1:ns
  B2 .= transpose(B1)
  p = [B1,B2]

  ode_rhs!(dy,y,p,t) = dy .= 1.01*y
  constraint_force!(dy,z,p) = dy .= p[1]*z
  constraint_rhs!(dz,p,t) = dz .= 1.0
  constraint_op!(dz,y,p) = dz .= p[2]*y

  u₀ = solvector(state=y,constraint=z)
  du = zero(u₀)

  f = ConstrainedODEFunction(ode_rhs!,constraint_rhs!,
                              constraint_force!,constraint_op!,_func_cache=u₀)
  f(du,u₀,p,0.0)
  @test state(du) == 1.01*y .- B1*z
  @test constraint(du) == 1.0 .- B2*y

  ode_rhs(y,p,t) = 1.01*y
  constraint_force(z,p) = p[1]*z
  constraint_rhs(p,t) = fill(1.0,nc)
  constraint_op(y,p) = p[2]*y

  f = ConstrainedODEFunction(ode_rhs,constraint_rhs,
                            constraint_force,constraint_op)
  du = f(u₀,p,0.0)
  @test state(du) == 1.01*y .- B1*z
  @test constraint(du) == 1.0 .- B2*y

  L = 2*I
  f = ConstrainedODEFunction(ode_rhs!,constraint_rhs!,
                              constraint_force!,constraint_op!,L,_func_cache=u₀)

  f(du,u₀,p,0.0)
  @test state(du) == 1.01*y .- B1*z .+ L*y
  @test state(du) ≈ [1.01,-0.99,-2.99,-4.99,-6.99] atol=1e-12
  @test constraint(du) == 1.0 .- B2*y
  @test constraint(du) ≈ [-14.0, -14.0] atol=1e-12

  f = ConstrainedODEFunction(ode_rhs,constraint_rhs,
                            constraint_force,constraint_op,L)

  du = f(u₀,p,0.0)
  @test state(du) == 1.01*y .- B1*z .+ L*y
  @test state(du) ≈ [1.01,-0.99,-2.99,-4.99,-6.99] atol=1e-12
  @test constraint(du) == 1.0 .- B2*y
  @test constraint(du) ≈ [-14.0, -14.0] atol=1e-12

  # Unconstrained system
  u₀ = solvector(state=y)
  du = zero(u₀)
  f = ConstrainedODEFunction(ode_rhs!,L,_func_cache=u₀)
  f(du,u₀,p,0.0)
  @test state(du) == 1.01*y .+ L*y
  @test state(du) ≈ fill(3.01,ns) atol=1e-12
  @test constraint(du) == empty(y)

  # Test that unconstrained systems create a saddle system that works properly
  # Should only invert the upper left operator
  u = deepcopy(u₀)
  S = SaddleSystem(2*I,f,p,p,deepcopy(du),Direct)
  u .= S\u₀
  @test state(u) ≈ fill(0.5,ns) atol=1e-12
  @test constraint(u) == empty(y)

  f = ConstrainedODEFunction(ode_rhs,L)
  du = f(u₀,p,0.0)
  @test state(du) == 1.01*y .+ L*y
  @test state(du) ≈ fill(3.01,ns) atol=1e-12
  @test constraint(du) == empty(y)

  # Test that unconstrained systems create a saddle system that works properly
  # Should only invert the upper left operator
  u = deepcopy(u₀)
  S = SaddleSystem(2*I,f,p,p,deepcopy(du),Direct)
  u .= S\u₀
  @test state(u) ≈ fill(0.5,ns) atol=1e-12
  @test constraint(u) == empty(y)

end


text:
export fraction, arclength, intersection

function Arc(from::Point, through::Point, to::Point)
    circle = Circle(from, through, to)
    start_angle = angle(from_to(circle.center, from))
    end_angle = angle(from_to(circle.center, to))

    # correct for the angle points being incorrectly transformed on the unit circle
    is_cw = rad(signed_angle_to(through → from, through → to)) >= 0
    if is_cw && end_angle > start_angle
        end_angle -= rad(2π)
    elseif !is_cw && end_angle < start_angle
        end_angle += rad(2π)
    end

    Arc(circle.center, circle.radius, start_angle, end_angle)
end

function Arc(from::Point, to::Point, radiusfraction::Real)
    heightpoint = between(from, to, 0.5) + rotate(from → to, deg(90)) * radiusfraction / 2
    arc = Arc(from, heightpoint, to)
end

function fraction(a::Arc, f::Real)
    angle = a.start_angle + (a.end_angle - a.start_angle) * f
    a.center + Point(angle) * a.radius
end

function fractionangle(a::Arc, f::Real)
    angle = a.start_angle + (a.end_angle - a.start_angle) * f
end

needed_attributes(::Type{Arc}) = (Visible, Linewidth, Stroke, Linestyle)

function arclength(a::Arc)
    a.radius * abs(rad(a.end_angle - a.start_angle))
end

function lengthen(a::Arc, ang::Angle)
    Arc(a.center, a.radius, a.start_angle, a.end_angle + ang)
end

function between(v, a, b)
    if b > a
        a <= v <= b
    else
        b <= v <= a
    end
end

function intersection(a::Arc, l::Line)
    c = Circle(a.center, a.radius)

    points = intersection(c, l)
    angles = angle.(points .- a.center)

    [p for (an, p) in zip(angles, points) if between(an, a.start_angle, a.end_angle)]
end


text:
using MCMCBenchmarks, Test, Random
Random.seed!(7841)
path = pathof(MCMCBenchmarks)
@testset "LBA Generative Model " begin
    include(joinpath(path, "../../Models/LBA/LBA_Models.jl"))
    include(joinpath(path, "../../Models/LBA/LinearBallisticAccumulator.jl"))
    v = [1.0, 1.5, 2.0]
    A = .8
    k = .2
    tau = .3
    Nd = 10^5
    data = simulateLBA(;Nd=Nd, v=v, A=A, k=k, tau=tau)
    @test typeof(data) == NamedTuple{(:choice, :rt, :N, :Nc),Tuple{Array{Int64, 1},Array{Float64, 1},Int64, Int64}}
    p1 = mean(data.choice .== 1)
    @test p1 ≈ .18 atol=2e-2
    p2 = mean(data.choice .== 2)
    @test p2 ≈ .32 atol=2e-2
end

@testset "Gaussian Generative Model " begin
    include(joinpath(path, "../../Models/Gaussian/Gaussian_Models.jl"))
    mu = 0
    sigma = 1
    Nd = 10^5
    data = simulateGaussian(;Nd=Nd, μ=mu, σ=sigma)
    @test typeof(data) == NamedTuple{(:y, :N),Tuple{Array{Float64, 1}, Int64}}
    @test mu ≈ mean(data.y) atol=1e-1
    @test sigma ≈ std(data.y) atol=1e-1
end

@testset "Hierarchical Poisson Generative Model " begin
    include(joinpath(path,
    "../../Models/Hierarchical_Poisson/Hierarhical_Poisson_Models.jl"))
    Nd = 1
    Ns = 10
    a0 = 1.0
    a1 = .5
    a0_sig = .3
    data = simulatePoisson(;Nd=Nd, Ns=Ns, a0=a0, a1=a1, a0_sig=a0_sig)
    @test typeof(data) == NamedTuple{(:y, :x, :idx, :N, :Ns),Tuple{Array{Int64,1},Array{Float64,1},Array{Int64,1},Int64,Int64}}
end

@testset "Signal Detection Generative Model " begin
    include(joinpath(path,"../../Models/SDT/SDT_Functions.jl"))
    include(joinpath(path,"../../Models/SDT/SDT.jl"))
    Nd = 10^4
    d = 2.0
    c = 0.0
    data = simulateSDT(;Nd=Nd, d=d, c=c)
    @test typeof(data) == NamedTuple{(:hits, :fas, :Nd),Tuple{Int64,Int64,Int64}}
    hr′ = data.hits/Nd
    far′ = data.fas/Nd
    hr = cdf(Normal(0, 1), d/2-c)
    far = cdf(Normal(0, 1), -d/2-c)
    @test hr′ ≈ hr atol = .02
    @test far′ ≈ far atol = .02
end

@testset "Linear Regression Generative Model " begin
    include(joinpath(path,
    "../../Models/Linear_Regression/Linear_Regression_Models.jl"))
    Nd=10
    Nc = 2
    β0=1.
    β=fill(.5, Nc)
    σ=1
    data = simulateRegression(;Nd=Nd, Nc=Nc, β=β, β0=β0)
    @test typeof(data) == NamedTuple{(:x, :y, :Nd, :Nc),Tuple{Array{Float64,2},Array{Float64,1},Int64,Int64}}
end


text:
if isdefined(@__MODULE__, :context_variable)
    context_variable += 2
end


text:
module TestVariables

using ModiaLang
@usingModiaPlot
using Measurements
using StaticArrays
using Test

@testset "... Test variable declarations" begin

v1 = input | output | flow | Var(min=0, nominal=10, init=5)
@test isequal(v1, Var(input = true, output = true, flow = true, min = 0, nominal = 10, init = 5))

v2 = parameter | 5 | Var(info="v2 info")
@test isequal(v2, Var(parameter = true, value = 5, info = "v2 info"))

v3 = "string value" | parameter | info"v3 info"
@test v3 == Var(parameter = true, value = "string value", info = "v3 info")

v4 = parameter | 5.5u"m/s" | Var(min=0)
@test isequal(v4, Var(parameter = true, value = 5.5u"m*s^-1", min = 0))

v5 = parameter | Var(5, min=0)
@test isequal(v5, Var(parameter = true, min = 0, value = 5))

Pin = Model( v = potential | Var(nominal=10), i = flow )
@test isequal(Pin, Model(v = Var(potential = true, nominal = 10), i = Var(flow = true)))

v6 = Var(2.0)
@test isequal(v6, Var(value = 2.0))

v7 = Var(2.0 ± 0.1)
@test isequal(v7, Var(value = 2.0 ± 0.1))

v8 = Var(SVector(1.0, 2.0, 3.0))
@test isequal(v8, Var(value = [1.0, 2.0, 3.0]))

v9 = v4 | Var(parameter=false)
@test isequal(v9, Var(parameter = false, value = 5.5u"m*s^-1", min = 0))

v10 = v4 | Var(parameter=nothing)
@test isequal(v10, Var(value = 5.5u"m*s^-1", min = 0))

v11 = v4 | interval(1,10)
@test isequal(v11, Var(parameter = true, value = 5.5u"m*s^-1", min = 1, max = 10))

TestVar1 = Model(
    p = parameter | -1 | info"A parameter",
    s = parameter | "string value" | info"A string parameter", # Not used
    x = Var(init=0),
    y = Var(start=0), # Not used
    equations = :[
        der(x) = p*x + 1
    ]
)

# @showModel TestVar1

model = @instantiateModel(TestVar1, logCode=true)
simulate!(model, stopTime= 1, requiredFinalStates = [1-exp(-1)])
plot(model, "x")

simulate!(model, stopTime= 1, merge=Map(p=-2, x=2), requiredFinalStates = [2-1.5*(1-exp(-2))])
plot(model, "x")

end

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
export Ellipse,Circle,Plate,NACA4

"""
    Ellipse(a,b,n) <: Body

Construct an elliptical body with semi-major axis `a` and semi-minor axis `b`,
with `n` points distributed on the body perimeter.
"""
mutable struct Ellipse{N} <: Body{N}
  a :: Float64
  b :: Float64
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function Ellipse(a::Float64,b::Float64,N::Int)
    x̃ = zeros(N)
    ỹ = zeros(N)
    θ = range(0,stop=2π,length=N+1)
    @. x̃ = a*cos(θ[1:N])
    @. ỹ = b*sin(θ[1:N])

    Ellipse{N}(a,b,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)
end

"""
    Circle(a,n) <: Body

Construct a circular body with radius `a`
and with `n` points distributed on the body perimeter.
"""
Circle(a::Float64,N::Int) = Ellipse(a,a,N)

function Base.show(io::IO, body::Ellipse{N}) where {N}
    if body.a == body.b
      println(io, "Circular body with $N points and radius $(body.a)")
    else
      println(io, "Elliptical body with $N points and semi-axes ($(body.a),$(body.b))")
    end
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

"""
    Plate(length,thick,n,[λ=1.0]) <: Body

Construct a flat plate with length `length` and thickness `thick`,
with `n` points distributed on the body perimeter.

The optional parameter `λ` distributes the points differently. Values between `0.0`
and `1.0` are accepted.

The constructor `Plate(length,n,[λ=1.0])` creates a plate of zero thickness.
"""
mutable struct Plate{N} <: Body{N}
  len :: Float64
  thick :: Float64
  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function Plate(len::Float64,N::Int;λ::Float64=1.0)

    # set up points on plate
    #x = [[len*(-0.5 + 1.0*(i-1)/(N-1)),0.0] for i=1:N]

    Δϕ = π/(N-1)
    Jϕa = [sqrt(sin(ϕ)^2+λ^2*cos(ϕ)^2) for ϕ in range(π-Δϕ/2,stop=Δϕ/2,length=N-1)]
    Jϕ = len*Jϕa/Δϕ/sum(Jϕa)
    x̃ = -0.5*len .+ Δϕ*cumsum([0.0; Jϕ])
    ỹ = zero(x̃)

    Plate{N}(len,0.0,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)

end

function Plate(len::Float64,thick::Float64,N::Int;λ::Float64=1.0)
    # input N is the number of panels on one side only

    # set up points on flat sides
    Δϕ = π/N
    Jϕa = [sqrt(sin(ϕ)^2+λ^2*cos(ϕ)^2) for ϕ in range(π-Δϕ/2,stop=Δϕ/2,length=N)]
    Jϕ = len*Jϕa/Δϕ/sum(Jϕa)
    xtopface = -0.5*len .+ Δϕ*cumsum([0.0; Jϕ])
    xtop = 0.5*(xtopface[1:N] + xtopface[2:N+1])

    Δsₑ = Δϕ*Jϕ[1]
    Nₑ = 2*floor(Int,0.25*π*thick/Δsₑ)
    xedgeface = [0.5*len + 0.5*thick*cos(ϕ) for ϕ in range(π/2,stop=-π/2,length=Nₑ+1)]
    yedgeface = [          0.5*thick*sin(ϕ) for ϕ in range(π/2,stop=-π/2,length=Nₑ+1)]
    xedge = 0.5*(xedgeface[1:Nₑ]+xedgeface[2:Nₑ+1])
    yedge = 0.5*(yedgeface[1:Nₑ]+yedgeface[2:Nₑ+1])

    x̃ = Float64[]
    ỹ = Float64[]
    for xi in xtop
      push!(x̃,xi)
      push!(ỹ,0.5*thick)
    end
    for i = 1:Nₑ
      push!(x̃,xedge[i])
      push!(ỹ,yedge[i])
    end
    for xi in reverse(xtop,dims=1)
      push!(x̃,xi)
      push!(ỹ,-0.5*thick)
    end
    for i = Nₑ:-1:1
      push!(x̃,-xedge[i])
      push!(ỹ,yedge[i])
    end

    Plate{length(x̃)}(len,thick,(0.0,0.0),0.0,x̃,ỹ,x̃,ỹ)

end

function Base.show(io::IO, body::Plate{N}) where {N}
    println(io, "Plate with $N points and length $(body.len) and thickness $(body.thick)")
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end

"""
    NACA4(cam,pos,thick[;np=20][,Xc=(0.0,0.0)][,len=1.0]) <: Body{N}

Generates points in the shape of a NACA 4-digit airfoil of chord length 1. The
relative camber is specified by `cam`, the position of
maximum camber (as fraction of chord) by `pos`, and the relative thickness
by `thick`.

The optional parameter `np` specifies the number of points on the upper
or lower surface. The optional parameter `Zc` specifies the mean position of
the vertices (which is set to the origin by default). The optional parameter
`len` specifies the chord length.

# Example

```jldoctest
julia> w = Bodies.NACA4(0.0,0.0,0.12);
```
"""
mutable struct NACA4{N} <: Body{N}
  len :: Float64
  camber :: Float64
  pos :: Float64
  thick :: Float64

  cent :: Tuple{Float64,Float64}
  α :: Float64

  x̃ :: Vector{Float64}
  ỹ :: Vector{Float64}

  x :: Vector{Float64}
  y :: Vector{Float64}

end

function NACA4(cam::Number,pos::Number,t::Number;np=20,Xc=(0.0,0.0),len=1.0)

# Here, cam is the fractional camber, pos is the fractional chordwise position
# of max camber, and t is the fractional thickness.

npan = 2*np-2

# Trailing edge bunching
an = 1.5
anp = an+1
x = zeros(np)

θ = zeros(size(x))
yc = zeros(size(x))

for j = 1:np
    frac = Float64((j-1)/(np-1))
    x[j] = 1 - anp*frac*(1-frac)^an-(1-frac)^anp;
    if x[j] < pos
        yc[j] = cam/pos^2*(2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam/pos*(1-x[j]/pos))
        end
    else
        yc[j] = cam/(1-pos)^2*((1-2*pos)+2*pos*x[j]-x[j]^2)
        if pos > 0
            θ[j] = atan(2*cam*pos/(1-pos)^2*(1-x[j]/pos))
        end
    end
end

xu = zeros(size(x))
yu = xu
xl = xu
yl = yu

yt = t/0.20*(0.29690*sqrt.(x)-0.12600*x-0.35160*x.^2+0.28430*x.^3-0.10150*x.^4)

xu = x-yt.*sin.(θ)
yu = yc+yt.*cos.(θ)

xl = x+yt.*sin.(θ)
yl = yc-yt.*cos.(θ)

rt = 1.1019*t^2;
θ0 = 0
if pos > 0
    θ0 = atan(2*cam/pos)
end
# Center of leading edge radius
xrc = rt*cos(θ0)
yrc = rt*sin(θ0)
θle = collect(0:π/50:2π)
xlec = xrc .+ rt*cos.(θle)
ylec = yrc .+ rt*sin.(θle)

# Assemble data
coords = [xu yu xl yl x yc]
cole = [xlec ylec]

# Close the trailing edge
xpanold = [0.5*(xl[np]+xu[np]); reverse(xl[2:np-1],dims=1); xu[1:np-1]]
ypanold = [0.5*(yl[np]+yu[np]); reverse(yl[2:np-1],dims=1); yu[1:np-1]]

xpan = zeros(npan)
ypan = zeros(npan)
for ipan = 1:npan
    if ipan < npan
        xpan1 = xpanold[ipan]
        ypan1 = ypanold[ipan]
        xpan2 = xpanold[ipan+1]
        ypan2 = ypanold[ipan+1]
    else
        xpan1 = xpanold[npan]
        ypan1 = ypanold[npan]
        xpan2 = xpanold[1]
        ypan2 = ypanold[1]
    end
    xpan[ipan] = 0.5*(xpan1+xpan2)
    ypan[ipan] = 0.5*(ypan1+ypan2)
end
w = ComplexF64[1;reverse(xpan,dims=1)+im*reverse(ypan,dims=1)]*len
w .-= mean(w)

x̃ = real.(w)
ỹ = imag.(w)

NACA4{length(x̃)}(len,cam,pos,t,Xc,0.0,x̃,ỹ,x̃,ỹ)

end

function Base.show(io::IO, body::NACA4{N}) where {N}
    println(io, "NACA 4-digit airfoil with $N points and length $(body.len) and thickness $(body.thick)")
    println(io, "   Current position: ($(body.cent[1]),$(body.cent[2]))")
    println(io, "   Current angle (rad): $(body.α)")
end


text:
#=
# A set of time steppers
=#

mutable struct Stepper
    type::String;       # The constant for the stepper type
    implicit::Bool;     # implicit or explicit
    Nsteps::Int;        # number of steps
    dt::Float64;        # step size
    cfl::Float64;       # CFL number
    stages::Int;        # how many stages
    a::Array{Float64};  # for RK steppers
    b::Array{Float64};  # for RK steppers
    c::Array{Float64};  # for RK steppers

    Stepper(t, c) = new(t, false, 0, 0, c, 0, [], [], []);
end

function init_stepper(x, stepper)
    if stepper.type == EULER_EXPLICIT
        stepper.implicit = false;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.25;
        end
        stepper.dt = stepper.cfl*dxmin*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == EULER_IMPLICIT
        stepper.implicit = true;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == PECE
        stepper.implicit = false;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 1;
        end
        stepper.dt = stepper.cfl*dxmin*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == CRANK_NICHOLSON
        stepper.implicit = true;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.25;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == LSRK4
        stepper.implicit = false;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 5;
        stepper.a = [0; -0.41789047449985195; -1.192151694642677; -1.6977846924715279; -1.5141834442571558];
        stepper.b = [0.14965902199922912; 0.37921031299962726;  0.8229550293869817; 0.6994504559491221; 0.15305724796815198];
        stepper.c = [0.0; 0.14965902199922912; 0.37040095736420475; 0.6222557631344432; 0.9582821306746903];

        return stepper;

    elseif stepper.type == RK4
        stepper.implicit = false;
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 4;
        stepper.a = [   0   0   0   0;
                        0.5 0   0   0;
                        0   0.5 0   0;
                        0   0   1   0];
        stepper.b = [1/6; 1/3; 1/3; 1/6];
        stepper.c = [0.0; 0.5; 0.5; 1];

        return stepper;
    end
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
"""
Computes BFS tree of a Graph from a given source. It returns the levels of each node with respect to source starting from level 1.

Input:
-----
* G::SparseMatrixCSC
* root::Int

Output:
------
* bfs_labels::Vector{Int32}

"""
function bfs(G::SparseMatrixCSC, root::Int)
	rows = G.colptr - 1
    cols = G.rowval - 1
    nodes = length(rows) - 1
    edges = length(cols)
    rows = map(Int32, rows)
    cols = map(Int32, cols)
    bfs_label = Array(Int32, nodes)
	ccall((:bfs, gunrock_lib), Void, (Ptr{Cint}, Int, Int, Ptr{Cint}, Ptr{Cint}, Int), bfs_label, nodes, edges, rows, cols, root - 1)
    bfs_label += 1
	bfs_label
end


text:
@testset "GaussMeter" begin
  @testset "SimulatedGaussMeter" begin
    include("SimulatedGaussMeterTest.jl")
  end

  if "lakeshore" in ARGS
    @testset "LakeShore" begin
      include("LakeShoreTest.jl")
    end
  end

  if "lakeshoref71" in ARGS
    @testset "LakeShoreF71" begin
      include("LakeShoreF71Test.jl")
    end
  end
end


text:
using CRRao
using Documenter

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

makedocs(;
    modules=[CRRao],
    authors="xKDR Forum",
    repo="https://github.com/xKDR/CRRao.jl/blob/{commit}{path}#{line}",
    sitename="CRRao.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xKDR.github.io/CRRao.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md"
    ],
    strict = :doctest
)

deploydocs(;
    repo="github.com/xKDR/CRRao.jl",
    target = "build",
    devbranch = "main"
)


text:
"""
    SimpleSDMPredictor(::Type{EarthEnv}, ::Type{HabitatHeterogeneity}, layer::Integer=1; resolution::Int64=25, kwargs...)

| Variable | Explanation              | Measure                                                          | Value range                         | Relationship to heterogeneity |
| -------- | ------------------------ | ---------------------------------------------------------------- | ----------------------------------- | ----------------------------- |
| 1        | Coefficient of variation | Normalized dispersion of EVI                                     | >=0                                 | Positive                      |
| 2        | Evenness                 | Evenness of EVI                                                  | >=0; <=1                            | Positive                      |
| 3        | Range                    | Range of EVI                                                     | >=0                                 | Positive                      |
| 4        | Shannon                  | Diversity of EVI                                                 | >=0; <=ln(max # of different EVI)   | Positive                      |
| 5        | Simpson                  | Diversity of EVI                                                 | >=0; <=1-1/(max # of different EVI) | Positive                      |
| 6        | Standard deviation       | Dispersion of EVI                                                | >=0                                 | Positive                      |
| 7        | Contrast                 | Exponentially weighted difference in EVI between adjacent pixels | >=0                                 | Positive                      |
| 8        | Correlation              | Linear dependency of EVI on adjacent pixels                      | >=-1; <=1                           | Nonlinear                     |
| 9        | Dissimilarity            | Difference in EVI between adjacent pixels                        | >=0                                 | Positive                      |
| 10       | Entropy                  | Disorderliness of EVI                                            | >=0                                 | Positive                      |
| 11       | Homogeneity              | Similarity of EVI between adjacent pixels                        | >=0; <=1                            | Negative                      |
| 12       | Maximum                  | Dominance of EVI combinations between adjacent pixels            | >=0; <=1                            | Negative                      |
| 13       | Uniformity               | Orderliness of EVI                                               | >=0; <=1                            | Negative                      |
| 14       | Variance                 | Dispersion of EVI combinations between adjacent pixels           | >=0                                 | Positive                      |

These data are released under a CC-BY-NC license to Tuanmu & Jetz.
"""
function SimpleSDMPredictor(::Type{EarthEnv}, ::Type{HabitatHeterogeneity}, layer::Integer=1; resolution::Int64=25, kwargs...)
    @assert resolution in [1, 5, 25]
    file = _get_raster(EarthEnv, HabitatHeterogeneity, layer, resolution)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end

function SimpleSDMPredictor(::Type{EarthEnv}, ::Type{HabitatHeterogeneity}, layers::AbstractArray; kwargs...)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(EarthEnv, HabitatHeterogeneity, l; kwargs...) for l in layers]
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

for (func,typ) in ((:test_csrsm,:CudaSparseMatrixCSR),
                   (:test_cscsm,:CudaSparseMatrixCSC))
    @eval begin
        function $func(elty)
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = $typ(sparse(A))
            info = CUSPARSE.sm_analysis('N','U',d_A,'O')
            d_Y = CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test_approx_eq(Y,h_Y)
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test_approx_eq(y,h_y)
            d_X = CudaArray(rand(elty,n,n))
            @test_throws(DimensionMismatch, CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O'))
            A = sparse(rand(elty,m,n))
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.sm_analysis('T','U',d_A,'O'))
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end
    end
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_csrsm(elty)
    test_cscsm(elty)
    println("sm took ", toq(), " for ", elty)
end


text:
#ENV["LD_LIBRARY_PATH"]=joinpath(@__DIR__,"shared")

module CppHello
    using CxxWrap
    @wrapmodule(joinpath(@__DIR__,"libhello"))

    function __init__()
        @initcxx
    end
end

@show CppHello.greet()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function demorgan(a::Bool, b::Bool)::Bool
    return (a && b) == !(!(a) || !(b))
end

@assert(demorgan(true, true))
@assert(demorgan(true, false))
@assert(demorgan(false, true))
@assert(demorgan(false, false))
println("OK")


text:
abstract type Logger end

struct NoLogger <:Logger end

struct ThermoLogger <:Logger
    freq::Int64
    Step::Bool
    Temp::Bool
    Energy::Bool
    Momentum::Bool
    index::Array{Int64,1}
    Write::Bool
    File::String
    Info::Array{Array{Float64,1}}
    function ThermoLogger(freq::Int64,
                            ; Energy=false,
                            Momentum=false,
                            Write=false,
                            File="" )
                            Step=true
                            Temp=true
                            index=[1]
                            info=Array{Array{Float64,1},1}(undef,0)
                            new(freq::Int64,
                                Step::Bool,
                                Temp::Bool,
                                Energy::Bool,
                                Momentum::Bool,
                                index::Array{Int64,1},
                                Write::Bool,
                                File::String,
                                info::Array{Array{Float64,1}}
                            )
    end
end

function init(logger::Logger,input::Dict,start_step::Int64)
    #init!()
    thermofreq=input["thermo information"]["freq"]

    @info "Log Thermo Information Every: $(thermofreq)"

    L=2
    energy_=false
    momentum_=false
    write_=false
    file_=""
    for i in 1:length(keys(input["thermo information"]))
        if "energy" in keys(input["thermo information"])
            if input["thermo information"]["energy"]
                energy_=true
                L+=1
            end
        end
        if "momentum" in keys(input["thermo information"])
            if input["thermo information"]["momentum"]
                momentum_=true
                L+=1
            end
        end
        if "write" in keys(input["thermo information"])
            write_=input["thermo information"]["write"]
        end
        if "file" in keys(input["thermo information"])
            file_=input["thermo information"]["file"]
        end
    end

    thermologger=ThermoLogger(thermofreq,Energy=energy_,Momentum=momentum_,Write=write_,File=file_)

    if thermofreq!=1
        push!(thermologger.Info,zeros(L))
    end
    for stepi in 1:input["steps"]

        if stepi%thermofreq==0.0
            push!(thermologger.Info,zeros(L))
        end
    end
    return thermologger
end