text:
########################
## Lift to Type Level ##
########################

abstract type TExpression end
struct TOperation{F,A1,A2} <: TExpression end

TExpression(x) = convert(TExpression, x)
Base.convert(::Type{TExpression}, v::Variable) = v.name
Base.convert(::Type{TExpression}, x::Constant) = x.value
function Base.convert(::Type{TExpression}, op::Operation)
    if length(op.args) == 1
        TOperation{op.func,TExpression(op.args[1]),Nothing}()
    else
        TOperation{op.func,TExpression(op.args[1]),TExpression(op.args[2])}()
    end
end

Base.convert(::Type{Expression}, ::TOperation{F,A1,Nothing}) where {F,A1} =
    Operation(F, Expression[Expression(A1)])
Base.convert(::Type{Expression}, ::TOperation{F,A1,A2}) where {F,A1,A2} =
    Operation(F, Expression[Expression(A1), Expression(A2)])

#############
## TSystem ##
#############
struct TSystem{TE,V,P} end

type_level(sys::System) = typeof(TSystem(sys.expressions, sys.variables, sys.parameters))
function TSystem(
    exprs::Vector{<:Expression},
    var_order::AbstractVector{<:Variable},
    param_order::AbstractVector{<:Variable} = Variable[],
)
    V = tuple((var.name for var in var_order)...)
    if !isempty(param_order)
        P = tuple((var.name for var in param_order)...)
    else
        P = Nothing
    end
    TE = tuple(convert.(TExpression, exprs)...)
    TSystem{TE,V,P}()
end

Base.show(io::IO, ::Type{T}) where {T<:TSystem} = show_info(io, T)
function Base.show(io::IO, TS::TSystem)
    show_info(io, typeof(TS))
    print(io, "()")
end
function show_info(io::IO, ::Type{TSystem{TE,V,P}}) where {TE,V,P}
    n = length(TE)
    m = length(V)
    mp = (P == Nothing ? 0 : length(P))
    print(io, "TSystem{$n,$m,$mp,#$(hash(TE))}")
end

interpret(TS::TSystem) = interpret(typeof(TS))
function interpret(::Type{TSystem{TE,V,P}}) where {TE,V,P}
    exprs = [Expression(e) for e in TE]
    vars = [Variable(e) for e in V]
    if P == Nothing
        params = Variable[]
    else
        params = [Variable(v) for v in P]
    end
    System(exprs, vars, params)
end

###############
## THomotopy ##
###############

struct THomotopy{TE,V,T,P} end

type_level(sys::Homotopy) =
    typeof(THomotopy(sys.expressions, sys.variables, sys.t, sys.parameters))
function THomotopy(
    exprs::Vector{<:Expression},
    var_order::AbstractVector{<:Variable},
    t::Variable,
    param_order::AbstractVector{<:Variable} = Variable[],
)
    TE = tuple(convert.(TExpression, exprs)...)
    V = tuple((var.name for var in var_order)...)
    T = t.name
    P = isempty(param_order) ? Nothing : tuple((var.name for var in param_order)...)
    THomotopy{TE,V,T,P}()
end

Base.show(io::IO, ::Type{T}) where {T<:THomotopy} = show_info(io, T)
function Base.show(io::IO, TS::THomotopy)
    show_info(io, typeof(TS))
    print(io, "()")
end
function show_info(io::IO, ::Type{THomotopy{TE,V,T,P}}) where {TE,V,T,P}
    n = length(TE)
    m = length(V)
    mp = (P == Nothing ? 0 : length(P))
    print(io, "THomotopy{$n,$m,$mp,#$(hash(TE))}")
end

interpret(TS::THomotopy) = interpret(typeof(TS))
function interpret(::Type{THomotopy{TE,V,T,P}}) where {TE,V,T,P}
    exprs = [Expression(e) for e in TE]
    vars = [Variable(e) for e in V]
    t = Variable(T)
    params = P == Nothing ? Variable[] : [Variable(v) for v in P]
    Homotopy(exprs, vars, t, params)
end

###########################
## CODEGEN + COMPILATION ##
###########################

struct InstructionId
    name::Symbol
end
name(id::InstructionId) = id.name

Base.:(==)(x::InstructionId, y::InstructionId) = x.name == y.name
Base.show(io::IO, v::InstructionId) = print(io, v.name)
Base.convert(::Type{Expr}, x::InstructionId) = x.name

const SimpleExpression = Union{InstructionId,Constant,Variable}

struct Instruction
    func::Union{Nothing,Symbol} # module, function name
    args::Vector{SimpleExpression}
end

Instruction(f::Union{Nothing,Symbol}, a::SimpleExpression) = Instruction(f, [a])
Instruction(x::SimpleExpression) = Instruction(nothing, x)

Base.hash(I::Instruction, h::UInt) = foldr(hash, I.args, init = hash(I.func, h))
Base.:(==)(a::Instruction, b::Instruction) = a.func == b.func && a.args == b.args

function Base.convert(::Type{Expr}, op::Instruction)
    if op.func === nothing && length(op.args) == 1
        convert(Expr, op.args[1])
    else
        Expr(:call, op.func, convert.(Expr, op.args)...)
    end
end
Base.show(io::IO, op::Instruction) = print(io, convert(Expr, op))

struct InstructionList
    instructions::OrderedDict{Instruction,InstructionId}
    var::Symbol
    n::Base.RefValue{Int}
end

function InstructionList(; var::Symbol = :ι, n::Base.RefValue{Int} = Ref(0))
    instructions = OrderedDict{Instruction,InstructionId}()
    InstructionList(instructions, var, n)
end

function Base.push!(v::InstructionList, i::Instruction; id_var::Symbol = v.var)
    if haskey(v.instructions, i)
        v.instructions[i]
    else
        id = InstructionId(Symbol(id_var, v.n[] += 1))
        push!(v.instructions, i => id)
        id
    end
end

Base.push!(list::InstructionList, x::SimpleExpression) = push!(list, Instruction(x))

function Base.push!(list::InstructionList, op::Operation)
    if op.func == :^ &&
       op.args[2] isa Constant && op.args[2].value isa Integer && op.args[2].value > 1
        n = op.args[2].value
        if n == 2
            i = _push!(list, op.args[1])
            push!(list, Instruction(:*, [i, i]))
        else
            n₁ = div(n, 2) + 1 # make sure n₁ > n₂
            n₂ = n - n₁
            id1 = push!(list, Operation(:^, op.args[1], Constant(n₁)))
            if n₂ == 1
                push!(list, Instruction(:*, [id1, _push!(list, op.args[1])]))
            else
                id2 = push!(list, Operation(:^, op.args[1], Constant(n₂)))
                push!(list, Instruction(:*, [id1, id2]))
            end
        end
    else
        push!(list, Instruction(op.func, map(arg -> _push!(list, arg), op.args)))
    end
end
_push!(list::InstructionList, op::Operation) = push!(list, op)
_push!(list::InstructionList, op::Expression) = op

Base.length(v::InstructionList) = length(v.instructions)
Base.iterate(v::InstructionList) = iterate(v.instructions)
Base.iterate(v::InstructionList, state) = iterate(v.instructions, state)
Base.eltype(v::Type{InstructionList}) = eltype(v.instructions)

function Base.show(io::IO, ::MIME"text/plain", list::InstructionList)
    println(io, "InstructionList:")
    for (instr, id) in list
        println(io, convert(Expr, id), " = ", convert(Expr, instr))
    end
end

function Base.convert(::Type{Expr}, list::InstructionList)
    Expr(:block, map(list) do (instr, id)
        :($(convert(Expr, id)) = $(convert(Expr, instr)))
    end...)
end

## CODEGEN
struct Compiled{T1<:Union{TSystem,THomotopy},T2<:Union{System,Homotopy}}
    obj::T2
end
Compiled(obj::T) where {T<:Union{System,Homotopy}} = Compiled{type_level(obj),T}(obj)

Base.size(C::Compiled) = size(C.obj)
Base.size(C::Compiled, i::Integer) = size(C.obj, i)
Base.length(C::Compiled) = length(C.obj)

const CompiledSystem{T<:TSystem} = Compiled{T,System}
const CompiledHomotopy{T<:THomotopy} = Compiled{T,Homotopy}

compile(F::System) = Compiled(F)
compile(H::Homotopy) = Compiled(H)

interpreted(C::Compiled) = C.obj

#################
## evaluations ##
#################
const TExpr = Union{TSystem,THomotopy}

make_indexing(F::System) = make_indexing(F.variables, nothing, F.parameters)
make_indexing(H::Homotopy) = make_indexing(H.variables, H.t, H.parameters)
function make_indexing(vars, t, params)
    lhs = Any[convert(Expr, v) for v in vars]
    rhs = Any[:(x[$i]) for i = 1:length(vars)]
    if t !== nothing
        push!(lhs, convert(Expr, t))
        push!(rhs, :t)
    end
    if params !== nothing
        append!(lhs, convert.(Expr, params))
        for i = 1:length(params)
            push!(rhs, :(p[$i]))
        end
    end
    :($(Expr(:tuple, lhs...)) = $(Expr(:tuple, rhs...)))
end

interpret(T::Type{<:TSystem}) = System(T)
interpret(T::Type{<:THomotopy}) = Homotopy(T)

function jacobian!(U::AbstractMatrix, S::Compiled, args...)
    evaluate_jacobian!(nothing, U, S, args...)
    U
end
function evaluate_jacobian(S::Compiled, args...)
    evaluate(S, args...), jacobian(S, args...)
end

function _evaluate_impl(::Type{T}) where {T<:TExpr}
    I = interpret(T)
    quote
        let $(make_indexing(I))
            $(begin
                list = InstructionList()
                ids = [convert(Expr, push!(list, fi)) for fi in I.expressions]
                quote
                    $(convert(Expr, list))
                    @SVector $(Expr(:vect, ids...))
                end
            end)
        end
    end
end

@generated function evaluate(F::CompiledSystem{T}, x::AbstractVector, p = nothing) where {T}
    _evaluate_impl(T)
end
(F::CompiledSystem)(x, p = nothing) = evaluate(F, x, p)
@generated function evaluate(
    F::CompiledHomotopy{T},
    x::AbstractVector,
    t,
    p = nothing,
) where {T}
    _evaluate_impl(T)
end
(H::CompiledHomotopy)(x, t, p = nothing) = evaluate(H, x, t, p)

function _evaluate!_impl(::Type{T}) where {T<:TExpr}
    u = gensym(:u)
    I = interpret(T)
    quote
        let $u = u
            let $(make_indexing(I))
                $(begin
                    list = InstructionList()
                    ids = [convert(Expr, push!(list, fi)) for fi in I.expressions]
                    quote
                        $(convert(Expr, list))
                        $(map(i -> :($u[$i] = $(ids[i])), 1:length(ids))...)
                    end
                end)
            end
        end
        u
    end
end

@generated function evaluate!(
    u::AbstractVector,
    F::CompiledSystem{T},
    x::AbstractVector,
    p = nothing,
) where {T}
    _evaluate!_impl(T)
end
@generated function evaluate!(
    u::AbstractVector,
    F::CompiledHomotopy{T},
    x::AbstractVector,
    t,
    p = nothing,
) where {T}
    _evaluate!_impl(T)
end

function _evaluate_jacobian!_impl(::Type{T}) where {T<:TExpr}
    u, U = gensym(:u), gensym(:U)
    I = interpret(T)
    n, m = size(I)
    quote
        let ($u, $U) = (u, U)
            let $(make_indexing(I))
                $(begin
                    list = InstructionList()
                    eval_ids = Vector{Any}(undef, n)
                    jac_ids = Matrix{Any}(undef, n, m)
                    for i = 1:n
                        fᵢ = I.expressions[i]
                        eval_ids[i] = convert(Expr, push!(list, fᵢ))
                        for j = 1:m
                            jac_ids[i, j] = convert(Expr, push!(list, I.jacobian[i, j]))
                        end
                    end
                    quote
                        $(convert(Expr, list))
                        if !($u isa Nothing)
                            $(map(i -> :($u[$i] = $(eval_ids[i])), 1:n)...)
                        end
                        $(vec([:($U[$i, $j] = $(jac_ids[i, j])) for i = 1:n, j = 1:m])...)
                        nothing
                    end
                end)
            end
        end
    end
end

@generated function evaluate_jacobian!(
    u::Union{Nothing,AbstractVector},
    U::AbstractMatrix,
    F::CompiledSystem{T},
    x::AbstractVector,
    p::Union{Nothing,AbstractVector} = nothing,
) where {T}
    _evaluate_jacobian!_impl(T)
end
@generated function evaluate_jacobian!(
    u::Union{Nothing,AbstractVector},
    U::AbstractMatrix,
    F::CompiledHomotopy{T},
    x::AbstractVector,
    t,
    p::Union{Nothing,AbstractVector} = nothing,
) where {T}
    _evaluate_jacobian!_impl(T)
end

function _jacobian_impl(::Type{T}) where {T<:TExpr}
    I = interpret(T)
    n, m = size(I)
    quote
        let $(make_indexing(I))
            $(begin
                list = InstructionList()
                jac_ids = map(e -> convert(Expr, push!(list, e)), I.jacobian)
                quote
                    $(convert(Expr, list))
                    @SMatrix $(Expr(:vcat, (Expr(:row, jac_ids[i, :]...) for i = 1:n)...))
                end
            end)
        end
    end
end

@generated function jacobian(F::CompiledSystem{T}, x, p = nothing) where {T}
    _jacobian_impl(T)
end

@generated function jacobian(F::CompiledHomotopy{T}, x, t, p = nothing) where {T}
    _jacobian_impl(T)
end

function _dt_impl(::Type{T}) where {T<:THomotopy}
    I = interpret(T)
    quote
        let $(make_indexing(I))
            $(begin
                list = InstructionList()
                ids = map(e -> convert(Expr, push!(list, e)), I.dt)
                quote
                    $(convert(Expr, list))
                    @SVector $(Expr(:vect, ids...))
                end
            end)
        end
    end
end

@generated function dt(H::CompiledHomotopy{T}, x::AbstractVector, t, p = nothing) where {T}
    _dt_impl(T)
end

function _dt!_impl(::Type{T}) where {T<:THomotopy}
    u = gensym(:u)
    I = interpret(T)
    quote
        let $u = u
            let $(make_indexing(I))
                $(begin
                    list = InstructionList()
                    ids = map(e -> convert(Expr, push!(list, e)), I.dt)
                    quote
                        $(convert(Expr, list))
                        $((:($u[$i] = $(ids[i])) for i = 1:length(ids))...)
                    end
                end)
            end
        end
        u
    end
end

@generated function dt!(
    u::AbstractVector,
    H::CompiledHomotopy{T},
    x::AbstractVector,
    t,
    p = nothing,
) where {T}
    _dt!_impl(T)
end

####

function _dt_jacobian!_impl(::Type{T}) where {T<:THomotopy}
    u, U = gensym(:u), gensym(:U)
    I = interpret(T)
    n, m = size(I)
    quote
        let ($u, $U) = (u, U)
            let $(make_indexing(I))
                $(begin
                    list = InstructionList()
                    dt_ids = Vector{Any}(undef, n)
                    jac_ids = Matrix{Any}(undef, n, m)
                    for i = 1:n
                        fᵢ = I.expressions[i]
                        dt_ids[i] = convert(Expr, push!(list, I.dt[i]))
                        for j = 1:m
                            jac_ids[i, j] = convert(Expr, push!(list, I.jacobian[i, j]))
                        end
                    end
                    quote
                        $(convert(Expr, list))
                        if !($u isa Nothing)
                            $(map(i -> :($u[$i] = $(dt_ids[i])), 1:n)...)
                        end
                        $(vec([:($U[$i, $j] = $(jac_ids[i, j])) for i = 1:n, j = 1:m])...)
                        nothing
                    end
                end)
            end
        end
    end
end

@generated function dt_jacobian!(
    u::AbstractVector,
    U::AbstractMatrix,
    H::CompiledHomotopy{T},
    x::AbstractVector,
    t,
    p = nothing,
) where {T}
    _dt_jacobian!_impl(T)
end

function dt_jacobian(H::CompiledHomotopy{T}, x::AbstractVector, t, p = nothing) where {T}
    dt(H, x, t, p), jacobian(H, x, t, p)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ ca5a9969-dc3d-4154-a2fd-c0057c5650bf
using PlutoUI

# ╔═╡ d9300c0c-f329-49ce-a1eb-fd9a52923709
md"### Defining Genetic Algortihm Structs"

# ╔═╡ ccac7604-7006-41c7-bfea-ccf25ef2cdb5
mutable struct population
	genes
end

# ╔═╡ 312b5890-b172-4c97-bb96-26e092b9d4b4
mutable struct genes
	genePool
	fitness
end

# ╔═╡ 873e0df3-6b7d-441b-95d6-ee09eb5d325c
md"# -----------------------------------------"

# ╔═╡ 647f2f4c-55c0-444c-a804-b54b6c96a588
function maximum(optG)
	max = length(optG)
	return max
end

# ╔═╡ 9bafd5f1-935a-434e-8a30-3415e0eed5a6
function populate(input,size)
	G = Array{genes}(undef, input)
	defaultGenes = rand([1,0],size)
	defaultFitness = 0
	for i = 1:input
		G[i] = genes(defaultGenes,defaultFitness)
	end
	Population = population(G)
end

# ╔═╡ 59ae2a0a-0c3f-47e7-9abd-bdff524dd018
md"### // Fitness Function"

# ╔═╡ 54e2cec7-539a-4e4d-a7d2-901d7598601c
function fitness(gene,optG)
	size = length(gene.genePool)
	gene.fitness = 0
	for i = 1:size
		if gene.genePool[i] == optG[i]
			gene.fitness+=1
		end
	end
	return gene.fitness
end

# ╔═╡ d3612443-07d5-41b4-849d-ec6ac1d053e8
function getFit(P)
	gSize = size(P.genes[1].genePool)
	f = genes(zeros(gSize),0)
	l = length(P.genes)
	for i = 1:l

			if f.fitness <= P.genes[i].fitness
				f = P.genes[i]
			end
	end
	return f
end

# ╔═╡ 377a605a-ecd4-4896-a2d3-895a7e1f4230
md"### // Selection Probability"

# ╔═╡ d2efa0ab-63f6-4908-83cf-1e46900eb914
function selection(P)
	poolSize = 5
	size = length(P.genes)
	gSize = length(P.genes[1].genePool)

	pool = populate(poolSize,size)

	for i = 1:poolSize
		random = rand(1:size)
		pool.genes[i] = P.genes[random]
	end

	f = getFit(pool)
	return f
end

# ╔═╡ 15feb3d0-4b5e-4f45-8e29-18984a885d3e
md"### // Crossover"

# ╔═╡ 22aaf6b4-aa7d-4844-b4a5-745b1edaf9fb
function crossover(g1,g2)
	rate = 0.5

	size=length(g1.genePool)

	gSize = length(g2.genePool)

	offSpring = genes(zeros(Int64,gSize),0)

	for i = 1:size
		if rand() <= rate
			offSpring = genes(g1.genePool,g1.fitness)
		else
			offSpring = genes(g2.genePool,g2.fitness)
		end
	end
	return offSpring
end

# ╔═╡ 299598be-8188-4451-a824-4067dd374c14
md"### // Mutation Probability"

# ╔═╡ 44480dcc-ab4f-44a5-b014-faba689620bf
function mutation(offSpring)
	mRate::Float64 = 0.5
	gSize =length(offSpring.genePool)

	for i = 1:gSize
		if rand() < mRate
			offSpring.genePool[i] = rand([1,0])
		end
	end

end

# ╔═╡ 7d9044ff-dfa9-472b-ba39-ff6fd70ffa31
md"### // Evolution Function"

# ╔═╡ fa5f3319-31be-4a7a-b26f-fab51a9a92f0
function evolution(P)

	size = length(P.genes)
	 y = population(P.genes)

	for i in 1:size
		g1 = selection(P)
        g2 = selection(P)
        oS = crossover(g1, g2)
        y.genes[i] = oS

	end

	s2 = length(y.genes)
	for i = 1:s2

	mutation(y.genes[i])
	end

    return y

end

# ╔═╡ 47cd8b1f-9690-4353-b902-33703ad9d892
md"### // Best Solution Function"

# ╔═╡ 0e4c481d-03ce-41df-b544-9371fc681954
md"# Test GA"

# ╔═╡ 52247fdf-548c-40ae-8c06-381d2359154a
x = populate(50,10)

# ╔═╡ 6997d095-b88a-4465-a67f-54807e1c9dd4
optimalGene = [1,1,1,1,1,1,1,1,1,1]

# ╔═╡ 3559305f-1825-412b-b377-863afad9b820
function best(P)
	gen = 0;

while(fitness(getFit(P),optimalGene) < maximum(optimalGene))
  gen += 1

  println("Generation: ", gen," Fitness: ",getFit(P).fitness )
  evolved = evolution(P)
	end
	best = getFit(P).genePool
	println("\nBest Solution in Generation ", gen)
	println("\nGene:", best)
end

# ╔═╡ 8d6a968f-5d37-4033-a8b1-1a201563a700
getFit(x)

# ╔═╡ b80e64b1-bf8d-402e-b542-70d2e29d64e2
with_terminal() do
best(x)
end

# ╔═╡ Cell order:
# ╠═ca5a9969-dc3d-4154-a2fd-c0057c5650bf
# ╟─d3e85cf5-c1d9-427f-b75b-db74298ebf30
# ╟─58370141-22fe-4c79-bffc-534bba496ea8
# ╟─b5ef8b65-3ada-4b70-a905-efa2978d4a6c
# ╠═d9300c0c-f329-49ce-a1eb-fd9a52923709
# ╠═ccac7604-7006-41c7-bfea-ccf25ef2cdb5
# ╠═312b5890-b172-4c97-bb96-26e092b9d4b4
# ╟─873e0df3-6b7d-441b-95d6-ee09eb5d325c
# ╠═647f2f4c-55c0-444c-a804-b54b6c96a588
# ╠═9bafd5f1-935a-434e-8a30-3415e0eed5a6
# ╟─59ae2a0a-0c3f-47e7-9abd-bdff524dd018
# ╠═54e2cec7-539a-4e4d-a7d2-901d7598601c
# ╠═d3612443-07d5-41b4-849d-ec6ac1d053e8
# ╟─377a605a-ecd4-4896-a2d3-895a7e1f4230
# ╠═d2efa0ab-63f6-4908-83cf-1e46900eb914
# ╟─15feb3d0-4b5e-4f45-8e29-18984a885d3e
# ╠═22aaf6b4-aa7d-4844-b4a5-745b1edaf9fb
# ╟─299598be-8188-4451-a824-4067dd374c14
# ╠═44480dcc-ab4f-44a5-b014-faba689620bf
# ╟─7d9044ff-dfa9-472b-ba39-ff6fd70ffa31
# ╠═fa5f3319-31be-4a7a-b26f-fab51a9a92f0
# ╟─47cd8b1f-9690-4353-b902-33703ad9d892
# ╠═3559305f-1825-412b-b377-863afad9b820
# ╟─0e4c481d-03ce-41df-b544-9371fc681954
# ╠═52247fdf-548c-40ae-8c06-381d2359154a
# ╠═6997d095-b88a-4465-a67f-54807e1c9dd4
# ╠═8d6a968f-5d37-4033-a8b1-1a201563a700
# ╠═b80e64b1-bf8d-402e-b542-70d2e29d64e2

ext: jl

max_stars_repo_name: LeonardoTPetrus/geneticalgorithm

max_issues_repo_name: LeonardoTPetrus/geneticalgorithm

max_forks_repo_name: LeonardoTPetrus/geneticalgorithm


text:
# ....................................................... FORTRAN-format object .............................................................................

"""
    FORTRAN_format

Object to hold a FORTRAN format specifier decomposed into its fields.
Accepted datatype specifiers are:  Aw,  Iw,  Fw.d,  Ew.d,  Dw.d.
Accepted output formating specifiers are: Aw,  Iw.m,  Bw.m,  Ow.m,  Zw.m,  Fw.d,  Ew.dEe,  ENw.d,  ESw.d,  Gw.dEe,  Dw.dEe.
Notation: 'w' - width, 'm' (optional) - minimum number of digits, 'd' - number of digits to right of decimal,
'e' - number of digits in exponent; 'N'/'S' (optional) - switch for engineering/scientific formating of the 'E' type.

The fields are:
* `Type::String`: primary FORTRAN datatype
* `TypeChar::Char`: primary FORTRAN datatype character
* `EngSci::Union{Char,Nothing}`: secundary datatype character (N for engineering, S for scientific)
* `width::Int`: width of numeric field
* `nmin::Int`: minimum number of digits displayed
* `ndec::Int`: number of digits to right of decimal
* `nexp::Int`: number of digits in exponent
"""
struct FORTRAN_format

    Type::String
    TypeChar::Char
    EngSci::Union{Char,Nothing}
    width::Int
    nmin::Int
    ndec::Int
    nexp::Int

end

# ................... cast format or datatype string 'format' into FORTRAN_format object ...........................

"""
    cast_FORTRAN_format(format::String)

Decompose the format specifier 'format' into its fields and cast this into the FORTRAN_format object.
Allowed format specifiers are of the types: Aw, Iw.m, Bw.m, Ow.m, Zw.m, Fw.d, Ew.dEe, ENw.d, ESw.d, Gw.dEe, Dw.dEe,
with: 'w' - width, 'm' (optional) - minimum number of digits, 'd' - number of digits to right of decimal,
'e' - number of digits in exponent; 'N'/'S' (optional) - switch for engineering/scientific formating of the 'E' type.
#### Examples:
```
t = cast_FORTRAN_format("I10")
FORTRAN_format("Iw", 'I', nothing, 10, 0, 0, 0)

t = cast_FORTRAN_format("I10.12")
FORTRAN_format("Iw.m", 'I', nothing, 10, 12, 0, 0)

t = cast_FORTRAN_format("E10.5E3")
FORTRAN_format("Ew.dEe", 'E', nothing, 10, 0, 5, 3)

t.Type, t.TypeChar, t.EngSci, t.width, t.nmin, t.ndec, t.nexp
("Ew.dEe", 'E', nothing, 10, 0, 5, 3)
```
"""
function cast_FORTRAN_format(str::String)

    s = strip(str,['\'',' ']); w = m = d = e = 0

    strErr = "FitsError: $s: not a valid FORTRAN print format "
    accept = ['A','I','B','O','Z','F','E','G','D','N','S','.','1','2','3','4','5','6','7','8','9','0']

    sum([s[i] ∉ accept for i ∈ eachindex(s)]) == 0 ? n=length(s) : return error(strErr * "(unknown type character)")
    n > 1 ?  X = s[1] : return error(strErr * " (width field not specified)")
    s[2] ∈ ['N','S'] ? (ns = s[2]; i = 3) : (ns = nothing; i = 2)
    rst = s[i:end]; length(rst) < 1 && return error(strErr * "(decimal side not specified)")

    if !occursin('.', rst)
        sum(.!isnumeric.(collect(rst))) > 0 && return error(strErr)
        X ∈ ['A','I','B','O','Z'] ? (t = X * "w"; w = parse(Int,rst)) : 0
        X ∈ ['F','E','G','D'] && return error(strErr * "(decimal field not specified)")
    else
        X ∈ ['A'] && return error(strErr * "(decimal point incompatible with $X type)")
        spl = split(rst,'.'); length(spl) > 2 && return error(strErr * "(two decimal points not allowed)")
        lhs = spl[1]
        rhs = spl[2]
        length(lhs) < 1 && return error(strErr * "(width field not specified)")
        length(rhs) < 1 && return error(strErr * "(decimal side not specified)")
        sum(.!isnumeric.(collect(lhs))) > 0 && return error(strErr * "(width field not numeric)")

        if !occursin('E',rhs)
            sum(.!isnumeric.(collect(rhs))) > 0 && return error(strErr * "(decimal field not numeric)")
            x = isnothing(ns) ? "" : ns
            X ∈ ['E','G','D'] ? (t = X * x * "w.d"; w = parse(Int,lhs); d = parse(Int,rhs)) : 0
        else
            X ∈ ['F','I','B','O','Z'] && return error(strErr * "(exponent incompatible with $X type)")
            splE = split(rhs,'E'); length(splE) > 2 && return error(strErr * "(unexpected E character)")
            rhs1 = splE[1]
            rhs2 = splE[2]
            length(rhs1) < 1 && return error(strErr * "(decimal field not specified)")
            length(rhs2) < 1 && return error(strErr * "(exponent field not specified)")
            sum(.!isnumeric.(collect(rhs1))) > 0 && return error(strErr * "(decimal field not numeric)")
            sum(.!isnumeric.(collect(rhs2))) > 0 && return error(strErr * "(exponent field not numeric)")
            X ∈ ['E','G','D'] ? (t = X * "w.dEe"; w = parse(Int,lhs); d = parse(Int,rhs1); e = parse(Int,rhs2)) : 0
        end
        X ∈ ['I','B','O','Z'] ? (t = X * "w.m"; w = parse(Int,lhs); m = parse(Int,rhs)) : 0
        X ∈ ['F'] ? (t = X * "w.d"; w = parse(Int,lhs); d = parse(Int,rhs)) : 0
    end

    return FORTRAN_format(t,X,ns,w,m,d,e)

end

"""
    cast_FORTRAN_datatype(format::String)

Decompose the format specifier 'format' into its fields and cast this into the FORTRAN_format object.
Allowed format specifiers are of the types: Aw, Iw, Fw.d, Ew.d, Dw.d,
where: 'w' - width, 'd' - number of digits to right of decimal point.
#### Examples:
```
f = cast_FORTRAN_datatype("I10")
FORTRAN_format("Iw", 'I', nothing, 10, 0, 0, 0)

f = cast_FORTRAN_datatypet("F10.4")
FORTRAN_format("Fw.d", 'F', nothing, 10, 0, 4, 0)

f = cast_FORTRAN_datatype("E10.5")
FORTRAN_format("Ew.d", 'E', nothing, 10, 0, 5, 0)

t.Type, t.TypeChar, t.EngSci, t.width, t.nmin, t.ndec, t.nexp
("Ew.d", 'E', nothing, 10, 0, 5, 0)
```
"""
function cast_FORTRAN_datatype(str::String)

    s = strip(str,['\'',' ']); w = d = 0

    strErr = "FitsError: $s: not a valid FORTRAN datatype "
    accept = ['A','I','F','E','D','.','1','2','3','4','5','6','7','8','9','0']

    sum([s[i] ∉ accept for i ∈ eachindex(s)]) == 0 ? n=length(s) : return error(strErr * "(unknown type character)")
    n > 1 ?  X = s[1] : return error(strErr * " (width field not specified)")
    rst = s[2:end]; length(rst) < 1 && return error(strErr * "(decimal side not specified)")

    if !occursin('.', rst)
        sum(.!isnumeric.(collect(rst))) > 0 && return error(strErr)
        X ∈ ['A','I'] ? (t = X * "w"; w = parse(Int,rst)) : 0
        X ∈ ['F','E','D'] && return error(strErr * "(decimal field not specified)")
    else
        X ∈ ['A','I'] && return error(strErr * "(decimal point incompatible with $X type)")
        spl = split(rst,'.'); length(spl) > 2 && return error(strErr * "(two decimal points not allowed)")
        lhs = spl[1]
        rhs = spl[2]
        length(lhs) < 1 && return error(strErr * "(width field not specified)")
        length(rhs) < 1 && return error(strErr * "(decimal side not specified)")
        sum(.!isnumeric.(collect(lhs))) > 0 && return error(strErr * "(width field not numeric)")

        sum(.!isnumeric.(collect(rhs))) > 0 && return error(strErr * "(decimal field not numeric)")
        X ∈ ['E','D'] ? (t = X * "w.d"; w = parse(Int,lhs); d = parse(Int,rhs)) : 0
        X ∈ ['F'] ? (t = X * "w.d"; w = parse(Int,lhs); d = parse(Int,rhs)) : 0
    end

    return FORTRAN_format(t,X,nothing,w,0,d,0)

end


text:
module Cincinnati

export greeting, bye

include("greetings.jl")
include("types.jl")

end # module

max_stars_repo_name: JuliaTagBot/Cincinnati

max_issues_repo_name: JuliaTagBot/Cincinnati

max_forks_repo_name: JuliaTagBot/Cincinnati


text:
using Test, PositionVelocityTime, GNSSDecoder
import PositionVelocityTime: GNSSDecoderState, GPSData, GPSL1Constants, ECEF, SatelliteState

include("test_data.jl")
include("sat_position.jl")
include("user_position.jl")
include("position_velocity_time.jl")
include("sv_time.jl")


text:
const USE_xpu = false
using ParallelStencil
using ParallelStencil.FiniteDifferences1D
@static if USE_xpu
    @init_parallel_stencil(CUDA, Float64, 1)
else
    @init_parallel_stencil(Threads, Float64, 1)
end
using Plots, Printf

macro Hx(ix)
    esc(:(H[$ix+2] - H[$ix]))
end

"""
    compute_U!(u2, u, H, g, dt_dx, size_u1_2)

ParallelStencil kernel that computes the Lax-Friedrichs for the momentum in x-direction of the shallow water
equations.
"""
@parallel_indices (ix) function compute_U!(u2, u, H, g, dt_dx, size_u1_2)
    if (ix <= size_u1_2)
        u2[ix+1] = u[ix+1] + 1 / 2 * dt_dx * (-1 / 2 * (u[ix+2] * u[ix+2] - u[ix] * u[ix]) - g * @Hx(ix))
    end
    return
end

"""
	compute_H!(H2, H, dt_dx, u, size_H1_2)

ParallelStencil kernel that computes a time step of the shallow water 1D
process.
"""
@parallel_indices (ix) function compute_H!(H2, H, dt_dx, u, size_H1_2)
    if (ix <= size_H1_2)
        H2[ix+1] = 1 / 2 * (H[ix+2] + H[ix]) - 1 / 2 * dt_dx * ((H[ix+2] * u[ix+2]) - (H[ix] * u[ix])) #continuity with Lax-Friedrichs Method
    end
    return
end

"""
    shallow_water_1D_xpu(;nx=1024, do_visu=false)

1D shallow water equations solver for an instantaneous dam break.
The Lax-Friedrichs Method was applied to the continuity equation.
Geometry (length of 40 meters) and initial conditions (half of domain
have initial water level of 5 meters, other half is dry) match
BASEMENT version 2.8 Test Case H_1 "Dam break in a closed channel

# Arguments
    - `nx`: number of discretised cells.
    - `do_visu`: if true, each physical time step will be ploted.

# Return values
    - `H`: The solution arrray (Water surface height in m).
"""
@views function shallow_water_1D_xpu(;
    # Numerics
    nx = 1024,
    # Visualisation
    do_visu = false)
    # Physics
    Lx = 40.0
    g = 9.81
    u_max = 6 #from review of results with very small timesteps
    ttot = 20.0
    # Numerics
    nout = 100
    # Derived numerics
    dx = Lx / nx
    dt = dx / u_max / 10 #stability condition of Lax-Friedrichs Method, further divided by 10 for stability of initially discontinuous wave front
    #println(dt)
    nt = cld(ttot, dt)
    xc = LinRange(dx / 2, Lx - dx / 2, nx)
    dt_dx = dt / dx
    # Array initialisation
    H = @zeros(nx)
    len = round(Int64(nx / 2))
    H[1:len] = 5.0 .+ @zeros(len)
    H2 = copy(H)
    u = @zeros(nx)
    u2 = copy(u)
    size_H1_2, size_u1_2 = size(H, 1) - 2, size(u, 1) - 2
    # Preparation of visualisation
    if do_visu
        ENV["GKSwstype"] = "nul"
        if isdir("plots/part-2/temp") == false
            mkdir("plots/part-2/temp")
        end
        loadpath = "./plots/part-2/temp/"
        anim = Animation(loadpath, String[])
        println("Animation directory: $(anim.dir)")
    end
    t_tic = 0.0
    niter = 0
    # Time loop
    for it = 0:nt
        if (it == 11)
            t_tic = Base.time()
            niter = 0
        end
        @parallel compute_H!(H2, H, dt_dx, u, size_H1_2)
        @parallel compute_U!(u2, u, H, g, dt_dx, size_u1_2)
        H, H2 = H2, H # pointer swap
        u, u2 = u2, u # pointer swap

        # Boundary conditions
        H[1] = H[2]
        H[end] = H[end-1]
        u[1] = -u[2]
        u[end] = -u[end-1]
        niter += 1
        if do_visu && (it % nout == 0)
            plot(xc, Array(H), xlims = (xc[1], xc[end]), ylims = (0, 10),
                xlabel = "Lx (m)", ylabel = "water surface elevation (m)", label = "h",
                title = "time = $(round(it*dt, sigdigits=3)) s, stability: $(round(maximum(abs.(u))*dt/dx,sigdigits=3)) /1",
                linewidth = :1.0, framestyle = :box)
            plot!(xc, u, label = "u", linewidth = :1.0)
            frame(anim)
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (2 * 2) / 1e9 * nx * sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it = t_toc / niter                      # Execution time per iteration [s]
    T_eff = A_eff / t_it                       # Effective memory throughput [GB/s]
    @printf("Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits = 3), niter)
    if do_visu
        gif(anim, "plots/part-2/shallow_water_1D_LF.gif", fps = 15)
    end
    return Array(H)
end

# shallow_water_1D_xpu(;do_visu=true)


text:
module MatchCore
import MLStyle.toolz: bind
using MLStyle.toolz.List
using MLStyle.toolz
using MLStyle.Err
using MLStyle.Render

# a token to denote matching failure
export Failed, failed
struct Failed end
const failed = Failed()

struct err_spec
    location :: LineNumberNode
    msg      :: String
end

struct config
    loc       :: LineNumberNode
    errs      :: linkedlist
end

# lens of config.
# Variables named in snake_case are marked as internal use
loc(conf      :: config)               = conf.loc
errs(conf     :: config)               = conf.errs
set_loc(loc   :: LineNumberNode)       = (conf) -> config(loc, conf.errs)
set_errs(errs :: linkedlist)           = (conf) -> config(conf.loc, errs)

# some useful utils:
err!(msg :: String) =
    bind(getBy $ loc)  do loc
    bind(getBy $ errs) do errs
    lens = set_errs $ cons(err_spec(loc, msg), errs)
    putBy(lens)
    end
    end

check_do(f, check_fn, err_fn) = expr ->
     if !check_fn(expr)
         err_fn(expr)
     else
         f(expr)
     end

# TODO: for friendly error report
recog_err(expr) :: String = "syntax error"

check_syntax(f, predicate) = begin
    check_fn = expr -> predicate(expr)
    err_fn   = expr -> bind(err! ∘ recog_err $ expr) do _
        return! $ nothing
    end
    check_do(f, check_fn, err_fn)
end

const init_state = config(LineNumberNode(1), nil())

# implementation of qualified pattern matching

export Qualifier
Qualifier = Function

export internal, invasive, share_with, share_through

internal = (my_mod, umod) -> my_mod === umod
invasive = (my_mod, umod) -> true
share_with(ms::Set{Module}) = (_, umod) -> umod in ms
share_through(symbol, value) = (_, umod) -> isdefined(umod, symbol) && getfield(umod, symbol) === true

export qualifier_test
function qualifier_test(qualifiers :: Set{Qualifier}, use_mod, def_mod)
    any(qualifiers) do q
        q(def_mod, use_mod)
    end
end

export PDesc
struct PDesc
    # for gapp patterns: (spec_vars :: Vector{Any}, gapobject :: Object, tl :: Vector{AST}) -> bool
    # for app patterns: (appobject :: Object, tl :: Vector{AST}) -> bool
    # for general    : (case:AST) -> bool

    predicate  :: Function

    # for gapp:
    #    (to_match : Symbol,
    #     forall :: Vector{Any},
    #     spec_vars :: Vector{Any},
    #     appobject :: Object,
    #     args :: Vector{AST},
    #     mod :: Module) -> (AST -> AST)
    # for app: (to_match  : Symbol, appobject: Object, args::Vector{AST}, mod::Module) -> (AST -> AST)
    # general: (to_match  : Symbol, case:AST, mod :: Module) -> (AST -> AST)
    rewrite    :: Function

    qualifiers :: Set{Qualifier}
end

PDesc(;predicate, rewrite, qualifiers) =
    PDesc(predicate, rewrite, qualifiers)

# A bug occurred when key is of Module type.
# Tentatively we use associate list(vector?).

const PATTERNS = Vector{Tuple{Module, PDesc}}()

export register_pattern

function register_pattern(pdesc :: PDesc, defmod :: Module)
    push!(PATTERNS, (defmod, pdesc))
end

# a simple example to define pattern `1`:
# tp = pattern_descriptor(
#       x -> x === 1,
#       (s, c, m) -> quote $s == 1 end,
#       Set([invasive])
#      )
# registerPattern(tp, MatchCore)

function get_pattern(case, use_mod :: Module)
    for (def_mod, desc) in PATTERNS
        if qualifier_test(desc.qualifiers, use_mod, def_mod) && desc.predicate(case)
           return desc.rewrite
        end
    end
end

is_head_eq(s :: Symbol) = (e::Expr) -> e.head == s

function collect_cases(expr :: Expr) :: State
    expr |>
    check_syntax(is_head_eq(:block))       do block
    bind(forM(collect_case, block.args))   do cases
    return! $ filter(a -> a !== nothing, cases)
    end
    end
end

function collect_case(expr :: LineNumberNode) :: State
    bind(putBy ∘ set_loc $ expr) do _
    return! $ nothing
    end
end

function collect_case(expr :: Expr) :: State
    expr |>
    check_syntax(is_head_eq(:call))            do expr
    expr.args |>
    check_syntax(args ->
                length(args) == 3 &&
                args[1]      == :(=>))      do (_, case, body)
    bind(getBy $ loc)                       do loc
    return! $ (loc, case, body)
    end
    end
    end
end

const INTERNAL_COUNTER = Dict{Module, Int}()

function remove_module_patterns(mod :: Module)
    delete!(INTERNAL_COUNTER, mod)
end

function get_name_of_module(m::Module) :: String
    string(m)
end

# allocate names for anonymous temporary variables.
export mangle
function mangle(mod::Module)
    get!(INTERNAL_COUNTER, mod) do
       0
    end |> id -> begin
    INTERNAL_COUNTER[mod] = id + 1
    mod_name = get_name_of_module(mod)
    Symbol("$mod_name $id")
    end

end

function match_impl(target, cbl, mod)
    # cbl: case body list
    # cbl = (fst ∘ (flip $ runState $ init_state) ∘ collectCases) $ cbl
    bind(collect_cases(cbl))                          do cbl
    # cbl :: [(LineNumberNodem, Expr, Expr)]
    tag_sym = mangle(mod)
    mk_match_body(target, tag_sym, cbl, mod)
    end
end

throw_from(errs) = begin
    # TODO: pretty print
    s = string(errs)
    throw(SyntaxError("$s"))
end

# the form:
# @match begin
#     ...
# end
export @match
macro match(target, cbl)
   (a, s) = runState $ match_impl(target, cbl, __module__) $ init_state
   if isempty(s.errs)
       esc $ a
   else
       throw_from(s.errs)
   end
end

function mk_match_body(target, tag_sym, cbl, mod)
    bind(getBy $ loc) do loc # start 1
    final =
        @format [loc, throw, InternalException] quote
            loc
            throw(InternalException("Non-exhaustive pattern found!"))
        end
    result = mangle(mod)
    cbl = collect(cbl)
    main_logic =
       foldr(cbl, init=final) do (loc, case, body), last # start 2
           expr = mk_pattern(tag_sym, case, mod)(body)
           @format [
               result,
               expr,
               loc,
               failed,
               last
           ] quote
              loc
              result = expr
              result === failed ? last : result
           end
       end  # end 2
    return! $ @format [tag_sym, target, main_logic] quote
       let tag_sym = target
           main_logic
       end
    end
    end # end 1
end

export mk_pattern
function mk_pattern(tag_sym :: Symbol, case :: Any, mod :: Module)
    rewrite = get_pattern(case, mod)
    if rewrite !== nothing
        return rewrite(tag_sym, case, mod)
    end
    case = string(case)
    throw $ PatternUnsolvedException("invalid usage or unknown case $case")
end

end # module end


text:
struct PDESystem <: ModelingToolkit.AbstractSystem
  eqs
  bcs
  domain
  indvars
  depvars

  @add_kwonly function PDESystem(eqs, bcs, domain, indvars, depvars)
      new(eqs, bcs, domain, indvars, depvars)
  end
end

Base.getproperty(x::PDESystem, sym::Symbol) = getfield(x, sym)

Base.summary(prob::PDESystem) = string(nameof(typeof(prob)))
function Base.show(io::IO, sys::PDESystem)
  println(io,summary(sys))
  println(io,"eqs: ", sys.eqs)
  println(io,"bcs: ", sys.bcs)
  println(io,"domain: ", sys.domain)
  println(io,"depvars: ", sys.depvars)
  println(io,"indvars: ", sys.indvars)
  return nothing
end


text:
using JuMP
using SeaPearl
using MathOptInterface

const MOI = MathOptInterface

struct GraphColoringVariableSelection  <: SeaPearl.AbstractVariableSelection{true}
    sortedPermutation
    degrees
end
function (vs::GraphColoringVariableSelection)(model::SeaPearl.CPModel)
    sortedPermutation, degrees = vs.sortedPermutation, vs.degrees
    maxDegree = 0
    toReturn = nothing
    for i in sortedPermutation
        if !SeaPearl.isbound(model.variables[string(i)])
            if isnothing(toReturn)
                toReturn = model.variables[string(i)]
                maxDegree = degrees[i]
            end
            if degrees[i] < maxDegree
                return toReturn
            end

            if length(model.variables[string(i)].domain) < length(toReturn.domain)
                toReturn = model.variables[string(i)]
            end
        end
    end
    return toReturn
end

struct KnapsackVariableSelection <: SeaPearl.AbstractVariableSelection{true} end
function (::KnapsackVariableSelection)(model::SeaPearl.CPModel)
    i = 1
    while SeaPearl.isbound(model.variables[string(i)])
        i += 1
    end
    return model.variables[string(i)]
end

@testset "Using SeaPearl with JuMP" begin
    @testset "GraphColoring model" begin
        model = Model(SeaPearl.Optimizer)

        numberOfVertices = 4

        @variable(model, 1 <= x[1:numberOfVertices] <= numberOfVertices)

        degrees = zeros(Int, numberOfVertices)

        @constraint(model, [x[1], x[2]] in SeaPearl.NotEqualSet())
        @constraint(model, [x[2], x[3]] in SeaPearl.NotEqualSet())
        @constraint(model, [x[2], x[4]] in SeaPearl.NotEqualSet())

        degrees[1] = 1
        degrees[2] = 3
        degrees[3] = 1
        degrees[4] = 1

        @variable(model, 1 <= y <= numberOfVertices)
        @constraint(model, x[1] <= y)
        @constraint(model, x[2] <= y)
        @constraint(model, x[3] <= y)
        @constraint(model, x[4] <= y)
        @objective(model, Min, y)

        sortedPermutation = sortperm(degrees; rev=true)

        # define the heuristic used for variable selection
        variableheuristic = GraphColoringVariableSelection(sortedPermutation, degrees)
        MOI.set(model, SeaPearl.MOIVariableSelectionAttribute(), variableheuristic)

        # Define the heuristic used for value selection
        numberOfSteps = 0
        valueheuristic = SeaPearl.BasicHeuristic((x; cpmodel=nothing) -> (numberOfSteps += 1; minimum(x.domain)))
        MOI.set(model, SeaPearl.MOIValueSelectionAttribute(), valueheuristic)

        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())

        @test numberOfSteps == 2
        @test status == MOI.OPTIMAL
        @test has_values(model)
        @test value.(x) == [1, 2, 1, 1] || value.(x) == [2, 1, 2, 2] #TODO: See Issue #43
        @test value(y) == 2
        println(model)
        println(status)
        println()
        println(value.(x))
        println(value(y))
    end

    @testset "Knapsack with JuMP" begin
        n = 4
        capacity = 11

        model = Model(SeaPearl.Optimizer)

        @variable(model, 0 <= x[1:n] <= 1)
        @constraint(model, 4 * x[1] + 5 * x[2] + 8 * x[3] + 3 * x[4] <= capacity)

        @expression(model, val_sum, 8 * x[1] + 10 * x[2] + 15 * x[3] + 4 * x[4])
        @objective(model, Min, -val_sum)

        variableheuristic = KnapsackVariableSelection()
        MOI.set(model, SeaPearl.MOIVariableSelectionAttribute(), variableheuristic)

        # Define the heuristic used for value selection
        numberOfSteps = 0
        valueheuristic = SeaPearl.BasicHeuristic((x; cpmodel=nothing) -> (numberOfSteps += 1; maximum(x.domain)))
        MOI.set(model, SeaPearl.MOIValueSelectionAttribute(), valueheuristic)

        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())

        @test numberOfSteps == 3
        @test status == MOI.OPTIMAL
        @test has_values(model)
        @test value.(x) == [0, 0, 1, 1]
        @test value(val_sum) == 19
        @test value(4 * x[1] + 5 * x[2] + 8 * x[3] + 3 * x[4]) == 11

        println(model)
        println(status)
        println(has_values(model))
        println(value.(x))
    end
end


text:
module PandorasBoxLaserNoteBlock

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/laserNoteBlock" LaserNoteBlock(x::Integer, y::Integer, pitch::Int=69, volume::Number=1.0, sound::String="game_03_deskbell_again", direction::String="Horizontal", atPlayer::Bool=false)

const directions = String["Horizontal", "Vertical"]

const placements = Ahorn.PlacementDict(
    "Laser Note Block ($direction) (Pandora's Box)" => Ahorn.EntityPlacement(
        LaserNoteBlock,
        "point",
        Dict{String, Any}(
            "direction" => direction
        )
    ) for direction in directions
)

Ahorn.editingOptions(entity::LaserNoteBlock) = Dict{String, Any}(
    "direction" => directions
)

const sprites = Dict{String, String}(
    "Horizontal" => "objects/pandorasBox/laser/noteblock/noteblock_horizontal",
    "Vertical" => "objects/pandorasBox/laser/noteblock/noteblock_vertical"
)

function Ahorn.selection(entity::LaserNoteBlock)
    x, y = Ahorn.position(entity)
    direction = get(entity, "direction", "Horizontal")
    sprite = get(sprites, direction, sprites["Horizontal"])

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=0.5)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::LaserNoteBlock, room::Maple.Room)
    direction = get(entity, "direction", "Horizontal")
    sprite = get(sprites, direction, sprites["Horizontal"])

    Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=0.5)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
# Test functions for checking Stirling1 and Stirling2

"""
Common code for the two Stirling matrix functions.
"""
function _matrix_maker(n::Int, f::Function)
    if n < 0
        throw(DomainError())
    end

    M = zeros(BigInt, n + 1, n + 1)
    for i = 0:n
        for j = 0:n
            M[i+1, j+1] = f(i, j)
        end
    end
    return M
end

"""
`Stirling1matrix(n)` creates an `n+1`-by-`n+1` matrix
of Stirling numbers of the first kind (from `0,0` to `n,n`).
"""
function Stirling1matrix(n::Int)
    return _matrix_maker(n, Stirling1)
end

"""
`Stirling2matrix(n)` creates an `n+1`-by-`n+1` matrix
of Stirling numbers of the second kind (from `0,0` to `n,n`).
"""
function Stirling2matrix(n::Int)
    return _matrix_maker(n, Stirling2)
end


text:
using DONEs
using Test
using Distributions
using LinearAlgebra
using Random
Random.seed!(1)

@testset "DONEs.jl" begin
    # Write your own tests here.

    @testset "RFE / Random Fourier Expansion" begin

        # case n > 1
        D = 10
        b = 2π*rand(D)
        c = randn(D)
        P = zeros(Float64,D,D) + Diagonal(ones(D))

        λ = 1E-6
        variable_offset = false
        h = zeros(Float64,D)
        offset = 0.

        @testset "Construction" for n in (1,4)
            Ω = randn(D,n)
            @test RFE(D,n,Ω,b,c,P,λ,variable_offset,h,offset) isa RFE
            @test RFE(D,n,0.1) isa RFE
            @test RFE(D,n,Distributions.Normal(0.,1.),Distributions.Uniform(0.,2π),λ=λ,variable_offset=variable_offset) isa RFE
        end

        @testset "Input validation" begin
            n = 4
            Ω = randn(D,n)

            @test_throws DimensionMismatch RFE(D-1,n,Ω,b,c,P,λ,variable_offset,h,offset)
            @test_throws DimensionMismatch RFE(D,n-1,Ω,b,c,P,λ,variable_offset,h,offset)
            @test_throws DimensionMismatch RFE(D,n,randn(D-1,n),b,c,P,λ,variable_offset,h,offset)
            @test_throws DimensionMismatch RFE(D,n,Ω,randn(D-1),c,P,λ,variable_offset,h,offset)
            @test_throws DimensionMismatch RFE(D,n,Ω,b,randn(D-1),P,λ,variable_offset,h,offset)
            @test_throws DimensionMismatch RFE(D,n,Ω,b,c,randn(D-1,D),λ,variable_offset,h,offset)
            @test_throws AssertionError RFE(D,n,Ω,b,c,P,-1.,variable_offset,h,offset)
            @test_throws AssertionError RFE(D,n,Ω,b,c,Diagonal(vcat(-1.,ones(D-1))) |> Matrix,λ,variable_offset,h,offset)
            @test_throws AssertionError RFE(D,n,Ω,b,c,randn(D,D),λ,variable_offset,h,offset)
        end

        @testset "Updating and evaluating the RFE" begin
            f(x) = x^2
            n = 1
            rfe = RFE(D,n,0.1)
            P = rfe.P |> copy
            N = 40
            points = LinRange(-1,1,N)
            @test begin
                try
                    for x in points
                        DONEs.updateRFE!(rfe,[x],f(x))
                    end
                catch
                    error()
                end
                rfe.P != P
            end
            # approximation power
            @test begin
                e = [DONEs.evaluateRFE(rfe,[x]) - f(x) for x in points]
                norm(e) < 3E-2
            end
            # downdating
            @test begin
                try
                    for x in LinRange(-1,1,5)
                        DONEs.downdateRFE!(rfe,[x],f(x))
                    end
                catch
                    error()
                end
                rfe.P != P
            end
        end
    end

    @testset "DONE functions, n=$n, sliding window=$sliding_window" for n in (1,2), sliding_window in (true,false)
        D = 30
        lb = -ones(n)
        ub = ones(n)
        current_optimal_x = 0.1*rand(Float64,n)
        sliding_window_length = 100
        σ_surrogate_exploration = 0.1
        σ_function_exploration = 0.1
        f(x) = dot(x,x)
        @testset "Construction" begin
            rfe = RFE(D,n,0.1)
            @test DONE(rfe, lb, ub, σ_surrogate_exploration, σ_function_exploration,sliding_window=sliding_window, sliding_window_length=sliding_window_length ) isa DONE
        end

        # @testset "Input validation" begin
        # end

        @testset "Adding measurements" begin
            rfe = RFE(D,n,0.1)
            done = DONE(rfe, lb, ub, σ_surrogate_exploration, σ_function_exploration,sliding_window=sliding_window, sliding_window_length=sliding_window_length )
            N = 200
            points = [2*rand(n) .- 1.0 for i in 1:N]
            P = rfe.P |> copy
            @test begin
                try
                    for x in points
                        add_measurement!(done,x,f(x))
                    end
                catch
                    error()
                end
                rfe.P != P
            end
            @test done.iteration == N

            # approximation power
            @test begin
                points = [2*rand(n) .- 1.0 for i in 1:50]
                e = [DONEs.evaluateRFE(rfe,x) - f(x) for x in points]
                norm(e) < 1E-2
            end

            @testset "Optimization with NLopt of the surrogate function" begin
                @test begin
                    try
                        update_optimal_input!(done)
                    catch
                        error()
                    end
                    norm(done.current_optimal_x - zeros(n)) < 1E-2
                end

            end
        end

    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test
using SuiteSparse

if Base.USE_GPL_LIBS
   include("umfpack.jl")
   include("cholmod.jl")
   include("spqr.jl")
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
function lab_15(r::Robot)
    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    for i ∈ (Nord, Ost, Sud, West)
        go_to_border_and_return_path!(r, i; markers = true)
    end
    go_by_path!(r, path)

end


text:
include("1.jl")
include("11.jl")
include("14.jl")
function mark_allzicper(r::Robot)
    angleofsquare1111(r::Robot)
    putmarker!(r)
    snake11per(r)
    if (isborder(r,Ost))
        putmarkholeper(r,West)
    else
        putmarkholeper(r,inverse(West))
    end
end

function snake11per(r::Robot)
    while !isborder(r,Nord)
        if (isborder(r,Ost))
            putmarkholeper(r,West)
        else
            putmarkholeper(r,inverse(West))
        end
        if (isborder(r,Ost) || isborder(r,West))
            disexit(r,Nord)
            putmarker!(r)
        end
    end

end
function angleofsquare1111(r::Robot)
    while ((!isborder(r,Sud)) || (!isborder(r,West)))
        if (!isborder(r,Sud))
            disexit(r,Sud)
        end
        if (!isborder(r,West))
            disexit(r,West)
        end
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
using JuliaVariables
using MLStyle
include("ngg/ngg.jl")
include("utils.jl")
include("closure.jl")
using .NGG
include("func_arg_decs.jl")

"""`ex` should be a scoped expression
"""
function closure_conv(top::Any, ex::Any)
    function conv(ex::Expr)
        @when Expr(:scoped, scope, inner) = ex begin
            block = Any[]
            for var in scope.bounds
                if var.is_mutable && var.is_shared
                    name = var.name
                    if var.name in scope.bound_inits
                        push!(block, :($name = Core.Box($name)))
                    else
                        push!(block, :($name = Core.Box()))
                    end
                end
            end
            push!(block, conv(inner))
            Expr(:block, block...)
            @when Expr(:function, head, inner && Expr(:scoped, scope, _)) = ex

            freenames = Symbol[f.name for f in scope.freevars]
            # If the evaluation module is a symbol and not in arguments
            if top isa Symbol && all(scope.bounds) do e
                e.name !== top
            end
                push!(freenames, top)
            end

            head = conv(head)
            fh = func_header(head)
            lambda_n = Symbol(:function)
            name = fh.name === unset ? lambda_n : fh.name
            fh = @with fh.args = FuncArg[map(func_arg, freenames)..., fh.args...]

            if fh.fresh !== unset || fh.ret !== unset
                error("GG doesn't support type parameters or return type annotations.")
            end

            args = of_args(fh.args)
            kwargs = of_args(fh.kwargs)
            inner = conv(inner)

            fn = mkngg(Symbol(name), args, kwargs, inner)
            if !isempty(freenames)
                closure_vars = Expr(:tuple, freenames...)
                fn = quote
                    let freevars = $closure_vars
                        $Closure{$fn,Base.typeof(freevars)}(freevars)
                    end
                end
            end

            if name !== lambda_n
                fn = Expr(:block, :($(fh.name) = $fn))
            end
            fn
            @when Expr(hd, args...) = ex
            Expr(hd, map(conv, args)...)
        end
    end

    function conv(s::Var)
        name = s.name
        s.is_global && return :($top.$name)
        s.is_mutable && s.is_shared && return begin
            :($name.contents)
        end
        name
    end
    conv(s) = s

    conv(ex.args[2])
end

function _get_body(::RuntimeFn{Args,Kwargs,Body}) where {Args,Kwargs,Body}
    Body
end

function _get_body(ex)
    error(ex)
end

struct UnderGlobal
    mod::Any
    ex::Any
end

macro under_global(m, ex)
    esc(:($UnderGlobal($m, $ex)))
end

function gg(compmod::Module, runmod::Any, source::Union{Nothing,LineNumberNode}, ex)
    (head, body) = @match ex begin
        Expr(:(=), head, body) => (head, body)
        Expr(:function, head, body) => (head, body)
        Expr(:->, head, body) => (head, body)
        _ => error("Malformed generated function at $source.")
    end

    fh = func_header(head)
    locals = Any[]
    if fh.args !== unset
        for arg in fh.args
            push!(locals, arg.name)
        end
    end
    if fh.kwargs !== unset
        for arg in fh.kwargs
            push!(locals, arg.name)
        end
    end
    if fh.fresh !== unset
        for name in map(extract_tvar, fh.fresh)
            push!(locals, name)
        end
    end
    if fh.name !== unset
        push!(locals, fh.name)
    end

    pseudo_head = Expr(:tuple, locals...)

    genbody = quote
        let body = $body
            if body isa $UnderGlobal
                ast = Base.macroexpand($compmod, body.ex)
                fake_ast = Base.Expr(:function, $(QuoteNode(pseudo_head)), ast)
                fake_ast = $simplify_ex(fake_ast)
                fake_ast = $solve!(fake_ast)
                fake_fn = $closure_conv(body.mod, fake_ast)
                $from_type($_get_body(fake_fn))
            else
                ast = Base.macroexpand($compmod, body)
                fake_ast = Base.Expr(:function, $(QuoteNode(pseudo_head)), ast)
                fake_ast = $simplify_ex(fake_ast)
                fake_ast = $solve!(fake_ast)
                fake_fn = $closure_conv($(QuoteNode(runmod)), fake_ast)
                $from_type($_get_body(fake_fn))
            end
        end
    end

    generator = Expr(:function, head, genbody)
    Expr(:macrocall, Symbol("@generated"), source, generator)
end

macro gg(ex)
    ex = gg(__module__, __module__, __source__, ex)
    esc(ex)
end


text:
module TestEltypeArrays

using Test
using ModifiedArrays
using ModifiedArrays: EltypeMod

@testset "Modifying element type" begin
    a = [1, 2, 3]
    A = modify(a, EltypeMod{Float64}())
    @test eltype(A) == Float64
    @test eltype(parent(a)) == Int64
    @test A[1] === 1.0
    A[2] = 4.0
    @test A[2] === 4.0
    @test parent(A)[2] === 4
    @test_throws InexactError A[2] = 4.5
    @test parent(A)[2] === 4
end

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using AbstractAlgebra
AbstractAlgebra.PrettyPrinting.set_html_as_latex(true)

eqstr(x) = "\\ \\text{is " * string(x) * "}"
eqstr(x::SetElem) = "=" * sprint(x) do io, x; show(io, "text/latex", x) end

function dispeq(lhs, rhs)
    latexstr = "\$\\displaystyle $lhs " * eqstr(rhs) * "\$"
    display("text/html", latexstr)
end

function euclidean(f, g, stop = 1)
    while true
        r = mod(f, g)
        iszero(r) && return g
        degree(r) ≤ stop && return r
        f, g = g, r
    end
end

function rootdeg1(f)
    C = coefficients(f) |> collect
    -first(C)//last(C)
end

function commonroot(f, g)
    m = euclidean(f, g, 1)
    rootdeg1(m)
end

function revpoly(z, x, f)
    m = degree(f)
    C = coefficients(f)
    sum(a*z^(i-1)*x^(m-i+1) for (i, a) in enumerate(C))
end

function val(h::FracElem, θ, L)
    n = L(numerator(h)(θ))
    d = L(denominator(h)(θ))
    n//d
end

function testval(h::FracElem, θ, β, L)
    n = L(numerator(h)(θ))
    d = L(denominator(h)(θ))
    n == β*d
end

function calcallresults(z, x, f, g)
    f_plus = (-1)^degree(f)*f(z - x)
    F_plus = resultant(f_plus, g)
    β_plus = commonroot(g, f_plus)

    f_mult = (-1)^degree(f)*revpoly(z, x, f)
    F_mult = resultant(f_mult, g)
    β_mult = commonroot(g, f_mult)

    K = FractionField(base_ring(z))
    B1, α = K["α"]
    L1 = ResidueField(B1, numerator(f(z))(α))
    B2, β = L1["β"]
    L = ResidueField(B2, numerator(g(z))(β))

    val_F_plus = val(F_plus, α + β, L)
    test_β_plus = testval(β_plus, α + β, β, L)
    val_F_mult = val(F_mult, α * β, L)
    test_β_mult = testval(β_mult, α * β, β, L)

    F_plus, β_plus, F_mult, β_mult, val_F_plus, test_β_plus, val_F_mult, test_β_mult
end

function dispallresults(z, x, f, g)
    F_plus, β_plus, F_mult, β_mult, val_F_plus, test_β_plus, val_F_mult, test_β_mult =
        @time calcallresults(z, x, f, g)
    flush(stdout)
    dispeq("F_\\alpha(x)", f)
    dispeq("F_\\beta(x)", g)
    dispeq("R_{\\alpha + \\beta}(z)", F_plus)
    dispeq("R_{\\alpha + \\beta}(\\alpha + \\beta)", val_F_plus)
    dispeq("\\beta_\\mathrm{plus}(z)", β_plus)
    dispeq("\\beta_\\mathrm{plus}(\\alpha + \\beta) = \\beta", test_β_plus)
    dispeq("R_{\\alpha\\beta}(z)", F_mult)
    dispeq("R_{\\alpha\\beta}(\\alpha\\beta)", val_F_mult)
    dispeq("\\beta_\\mathrm{mult}(z)", β_mult)
    dispeq("\\beta_\\mathrm{mult}(\\alpha\\beta) = \\beta", test_β_mult)
end

# %%
R, (a, b, c, p, q, r, s) = ZZ["a", "b", "c", "p", "q", "r", "s"]
K = FractionField(R)
Rz, z = R["z"]
Kz = FractionField(Rz)
Rx, x = Kz["x"]

# %%
f = x^2 - a
g = x^2 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^3 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^6 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^7 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^6 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^7 - p
dispallresults(z, x, f, g)

# %%
f = x^4 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
f = x^4 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^5 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^3 + p*x + q
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^4 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^5 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = a*x^2 + b*x + c
g = p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^2 + b*x + c
g = x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^2 + b*x + c
g = x^3 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^2 + a*x + b
g = x^4 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = x^2 + a*x + b
g = x^5 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = x^2 + a*x + b
g = x^6 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = x^2 + a*x + b
g = x^7 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^3 + p*x + q
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^4 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^5 + p*x^3 + q*x^2 + r*x + s
dispallresults(z, x, f, g)

# %%
f = x^3 + a*x + b
g = x^3 + p*x + q
dispallresults(z, x, f, g)

# %%
f = x^3 + a*x + b
g = x^4 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module DimensionalArrayTraits

export Order, Ordered, Unordered

export Grid, IndependentGrid, UnknownGrid

export AbstractAllignedGrid, AllignedGrid

export AbstractCategoricalGrid, CategoricalGrid,

export DependentGrid, TransformedGrid, LookupGrid

export Sampling, Single, Mean

export Locus, Start, Center, End

abstract type Order end

"""
Trait container for dimension and array ordering in AllignedGrid.

The default is `Ordered(Forward()`, `Forward())`

All combinations of forward and reverse order for data and indices seem to occurr
in real datasets, as strange as that seems. We cover these possibilities by specifying
the order of both explicitly.

Knowing the order of indices is important for using methods like `searchsortedfirst()`
to find indices in sorted lists of values. Knowing the order of the data is then
required to map to the actual indices. It's also used to plot the data later - which
always happens in smallest to largest order.

Base also defines Forward and Reverse, but they seem overly complicated for our purposes.
"""
struct Ordered{D,A} <: Order
    index::D
    array::A
end
Ordered() = Ordered(Forward(), Forward())

"""
Trait indicating that the array or dimension has no order.
"""
struct Unordered <: Order end

indexorder(order::Ordered) = order.index
indexorder(order::Unordered) = Unordered()
arrayorder(order::Ordered) = order.array
arrayorder(order::Unordered) = Unordered()

"""
Trait indicating that the array or dimension is in the normal forward order.
"""
struct Forward <: Order end

"""
Trait indicating that the array or dimension is in the reverse order.
Selector lookup or plotting will be reversed.
"""
struct Reverse <: Order end

Base.reverse(::Reverse) = Forward()
Base.reverse(::Forward) = Reverse()
Base.reverse(o::Ordered) = Ordered(revese(indexorder(o)), revese(arrayorder(o)))
Base.reverse(o::Unordered) = Unordered()

"""
Indicates wether the cell value is specific to the locus point
or is related to the whole the span.

The span may contain a value if the distance between locii if known.
This will often be identical to the distance between any two sequential
cell values, but may be distinct due to rounding errors in a vector index,
or context-dependent spans such as `Month`.
"""
abstract type Sampling end

struct Mean <: Sampling end

struct Single <: Sampling end

"""
Indicate the position of index values in grid cells.

This is frequently `Start` for time series, but may be `Center`
for spatial data.
"""
abstract type Locus end

"""
Indicates dimensions that are defined by their center coordinates/time/position.
"""
struct Center <: Locus end

"""
Indicates dimensions that are defined by their start coordinates/time/position.
"""
struct Start <: Locus end

"""
Indicates dimensions that are defined by their end coordinates/time/position.
"""
struct End <: Locus end

"""
Traits describing the grid type of a dimension.
"""
abstract type Grid end

dims(g::Grid) = nothing

"""
Fallback grid type
"""
struct UnknownGrid <: Grid end

"""
Traits describing a grid dimension that is independent of other grid dimensions.
"""
abstract type IndependentGrid{O} <: Grid end

abstract type AbstractAllignedGrid{O} <: IndependentGrid{O} end

"""
Trait describing a grid aligned with a dimension, independent of
other dimensions.

## Fields
- `order`: `Order` trait indicating array and index order
- `locus`: `Locus` trait indicating the position of the indexed point within the cell span
- `sampling`: `Sampling` trait indicating wether the grid cells are single samples or means
- `span`: the size of a grid step, such as 1u"km" or `Month(1)`
"""
struct AllignedGrid{O<:Order,L<:Locus,Sa<:Sampling,Sp} <: AbstractAllignedGrid{O}
    order::O
    locus::L
    sampling::Sa
    span::Sp
end
AllignedGrid(; order=Ordered(), locus=Center(), sampling=Mean(), span=nothing) =
    AllignedGrid(order, locus, sampling, span)

order(g::AllignedGrid) = g.order
span(g::AllignedGrid) = g.span
locus(g::AllignedGrid) = g.locus
sampling(g::AllignedGrid) = g.sampling

abstract type AbstractCategoricalGrid{O} <: IndependentGrid{O} end

"""
Traits describing a dimension where the values are categories.

## Fields
- `order`: `Order` trait indicating array and index order
"""
struct CategoricalGrid{O<:Order} <: AbstractCategoricalGrid{O}
    order::O
end
CategoricalGrid(;order=Ordered()) = CategoricalGrid(order)

order(g::CategoricalGrid) = g.order

"""
Traits describing a grid dimension that is dependent on other grid dimensions.

Indexing into a dependent dimension must provide all other dependent dimensions.
"""
abstract type DependentGrid <: Grid end

"""
Grid type using an affine transformation to convert dimension from
`dim(grid)` to `dims(array)`.

## Fields
- `dims`: a tuple containing dimenension types or symbols matching the order
          needed by the transform function.
- `sampling`: a `Sampling` trait indicating wether the grid cells are sampled points or means
"""
struct TransformedGrid{D,Sa<:Sampling} <: DependentGrid
    dims::D
    sampling::Sa
end
TransformedGrid(dims=(), sampling=Mean()) = TransformedGrid(dims, sampling)

dims(g::TransformedGrid) = g.dims
sampling(g::TransformedGrid) = g.sampling

"""
Grid type that uses an array lookup to convert dimension from
`dim(grid)` to `dims(array)`.

## Fields
- `dims`: a tuple containing dimenension types or symbols matching the order
          needed to index the lookup matrix.
- `sampling`: a `Sampling` trait indicating wether the grid cells are sampled points or means
"""
struct LookupGrid{D,Sa<:Sampling} <: DependentGrid
    dims::D
    sampling::Sa
end
LookupGrid(dims=(), sampling=Mean()) = LookupGrid(dims, sampling)

dims(g::LookupGrid) = g.dims
sampling(g::LookupGrid) = g.sampling

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## efficient value-based hashing of integers ##

function hash_integer(n::Integer, h::UInt)
    h ⊻= hash_uint((n % UInt) ⊻ h)
    n = abs(n)
    n >>>= sizeof(UInt) << 3
    while n != 0
        h ⊻= hash_uint((n % UInt) ⊻ h)
        n >>>= sizeof(UInt) << 3
    end
    return h
end

function hash_integer(n::BigInt, h::UInt)
    s = n.size
    s == 0 && return hash_integer(0, h)
    p = convert(Ptr{UInt}, n.d)
    b = unsafe_load(p)
    h ⊻= hash_uint(ifelse(s < 0, -b, b) ⊻ h)
    for k = 2:abs(s)
        h ⊻= hash_uint(unsafe_load(p, k) ⊻ h)
    end
    return h
end

## generic hashing for rational values ##

function hash(x::Real, h::UInt)
    # decompose x as num*2^pow/den
    num, pow, den = decompose(x)

    # handle special values
    num == 0 && den == 0 && return hash(NaN, h)
    num == 0 && return hash(ifelse(den > 0, 0.0, -0.0), h)
    den == 0 && return hash(ifelse(num > 0, Inf, -Inf), h)

    # normalize decomposition
    if den < 0
        num = -num
        den = -den
    end
    z = trailing_zeros(num)
    if z != 0
        num >>= z
        pow += z
    end
    z = trailing_zeros(den)
    if z != 0
        den >>= z
        pow -= z
    end

    # handle values representable as Int64, UInt64, Float64
    if den == 1
        left = ndigits0z(num,2) + pow
        right = trailing_zeros(num) + pow
        if -1074 <= right
            if 0 <= right && left <= 64
                left <= 63                     && return hash(Int64(num) << Int(pow), h)
                signbit(num) == signbit(den)   && return hash(UInt64(num) << Int(pow), h)
            end # typemin(Int64) handled by Float64 case
            left <= 1024 && left - right <= 53 && return hash(ldexp(Float64(num),pow), h)
        end
    end

    # handle generic rational values
    h = hash_integer(den, h)
    h = hash_integer(pow, h)
    h = hash_integer(num, h)
    return h
end

#=
`decompose(x)`: non-canonical decomposition of rational values as `num*2^pow/den`.

The decompose function is the point where rational-valued numeric types that support
hashing hook into the hashing protocol. `decompose(x)` should return three integer
values `num, pow, den`, such that the value of `x` is mathematically equal to

    num*2^pow/den

The decomposition need not be canonical in the sense that it just needs to be *some*
way to express `x` in this form, not any particular way – with the restriction that
`num` and `den` may not share any odd common factors. They may, however, have powers
of two in common – the generic hashing code will normalize those as necessary.

Special values:

 - `x` is zero: `num` should be zero and `den` should have the same sign as `x`
 - `x` is infinite: `den` should be zero and `num` should have the same sign as `x`
 - `x` is not a number: `num` and `den` should both be zero
=#

decompose(x::Integer) = x, 0, 1
decompose(x::Rational) = numerator(x), 0, denominator(x)

function decompose(x::Float16)::NTuple{3,Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt16, x)
    s = (n & 0x03ff) % Int16
    e = (n & 0x7c00 >> 10) % Int
    s |= Int16(e != 0) << 10
    d = ifelse(signbit(x), -1, 1)
    s, e - 25 + (e == 0), d
end

function decompose(x::Float32)::NTuple{3,Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt32, x)
    s = (n & 0x007fffff) % Int32
    e = (n & 0x7f800000 >> 23) % Int
    s |= Int32(e != 0) << 23
    d = ifelse(signbit(x), -1, 1)
    s, e - 150 + (e == 0), d
end

function decompose(x::Float64)::Tuple{Int64, Int, Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt64, x)
    s = (n & 0x000fffffffffffff) % Int64
    e = (n & 0x7ff0000000000000 >> 52) % Int
    s |= Int64(e != 0) << 52
    d = ifelse(signbit(x), -1, 1)
    s, e - 1075 + (e == 0), d
end

function decompose(x::BigFloat)::Tuple{BigInt, Int, Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return x.sign, 0, 0
    x == 0 && return 0, 0, x.sign
    s = BigInt()
    s.size = cld(x.prec, 8*sizeof(GMP.Limb)) # limbs
    b = s.size * sizeof(GMP.Limb)            # bytes
    ccall((:__gmpz_realloc2, :libgmp), Cvoid, (Ref{BigInt}, Culong), s, 8b) # bits
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), s.d, x.d, b) # bytes
    s, x.exp - 8b, x.sign
end

## streamlined hashing for smallish rational types ##

function hash(x::Rational{<:BitInteger64}, h::UInt)
    num, den = Base.numerator(x), Base.denominator(x)
    den == 1 && return hash(num, h)
    den == 0 && return hash(ifelse(num > 0, Inf, -Inf), h)
    if isodd(den)
        pow = trailing_zeros(num)
        num >>= pow
    else
        pow = trailing_zeros(den)
        den >>= pow
        pow = -pow
        if den == 1 && abs(num) < 9007199254740992
            return hash(ldexp(Float64(num),pow))
        end
    end
    h = hash_integer(den, h)
    h = hash_integer(pow, h)
    h = hash_integer(num, h)
    return h
end

## hashing Float16s ##

hash(x::Float16, h::UInt) = hash(Float64(x), h)

## hashing strings ##

const memhash = UInt === UInt64 ? :memhash_seed : :memhash32_seed
const memhash_seed = UInt === UInt64 ? 0x71e729fd56419c81 : 0x56419c81

function hash(s::Union{String,SubString{String}}, h::UInt)
    h += memhash_seed
    # note: use pointer(s) here (see #6058).
    ccall(memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), pointer(s), sizeof(s), h % UInt32) + h
end
hash(s::AbstractString, h::UInt) = hash(String(s), h)

max_stars_repo_name: ashaywalke/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
#! format: off

abstract type AbstractStorageFormulation <: AbstractDeviceFormulation end
abstract type AbstractEnergyManagement  <: AbstractStorageFormulation end
struct BookKeeping <: AbstractStorageFormulation end
struct BookKeepingwReservation <: AbstractStorageFormulation end
struct EndOfPeriodEnergyTarget <: AbstractEnergyManagement end

get_variable_sign(_, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = NaN
########################### ActivePowerInVariable, Storage #################################

get_variable_binary(::ActivePowerInVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ActivePowerInVariable, ::Type{<:PSY.Storage}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = nothing
get_variable_sign(::ActivePowerInVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = -1.0

########################### ActivePowerOutVariable, Storage #################################

get_variable_binary(::ActivePowerOutVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ActivePowerOutVariable, ::Type{<:PSY.Storage}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = nothing
get_variable_sign(::ActivePowerOutVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0

############## ReactivePowerVariable, Storage ####################
get_variable_sign(::PowerSimulations.ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0
get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_expression_name(::ReactivePowerVariable, ::Type{<:PSY.Storage}) = :nodal_balance_reactive

############## EnergyVariable, Storage ####################

get_variable_binary(::EnergyVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_lower_bound(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_initial_value(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_initial_energy(d)

############## ReserveVariable, Storage ####################

get_variable_binary(::ReserveVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = true

#! format: on

################################## output power constraints#################################

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:AbstractStorageFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            limits_func = x -> PSY.get_output_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:AbstractStorageFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            limits_func = x -> PSY.get_input_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:BookKeepingwReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_output_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:BookKeepingwReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.Storage}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_input_active_power_limits(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function add_constraints!(
    optimization_container::OptimizationContainer,
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_reactive_power_limits(d)
        constraint_infos[ix] = DeviceRangeConstraintInfo(name, limits)
    end

    device_range!(
        optimization_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(RangeConstraint, ReactivePowerVariable, St),
            make_variable_name(ReactivePowerVariable, St),
        ),
    )
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    formulation::AbstractStorageFormulation,
) where {St <: PSY.Storage}
    storage_energy_initial_condition!(optimization_container, devices, formulation)
    return
end

######################### Initialize Functions for Storage #################################
# TODO: This IC needs a cache for Simulation over long periods of tim
function storage_energy_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.Storage, D <: AbstractStorageFormulation}
    key = ICKey(EnergyLevel, T)
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        EnergyVariable(),
        key,
        _make_initial_condition_energy,
        _get_variable_initial_value,
        StoredEnergy,
    )

    return
end

############################ Energy Capacity Constraints####################################

function energy_capacity_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    model::DeviceModel{St, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_state_of_charge_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end

    device_range!(
        optimization_container,
        RangeConstraintSpecInternal(
            constraint_infos,
            make_constraint_name(ENERGY_CAPACITY, St),
            make_variable_name(ENERGY, St),
        ),
    )
    return
end

############################ book keeping constraints ######################################

function make_efficiency_data(
    devices::IS.FlattenIteratorWrapper{St},
) where {St <: PSY.Storage}
    names = Vector{String}(undef, length(devices))
    in_out = Vector{InOut}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        names[ix] = PSY.get_name(d)
        in_out[ix] = PSY.get_efficiency(d)
    end

    return names, in_out
end

function energy_balance_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    ::Type{D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    efficiency_data = make_efficiency_data(devices)
    key = ICKey(EnergyLevel, St)

    if !has_initial_conditions(optimization_container.initial_conditions, key)
        throw(
            IS.DataFormatError(
                "Initial Conditions for $(St) Energy Constraints not in the model",
            ),
        )
    end

    energy_balance(
        optimization_container,
        get_initial_conditions(optimization_container, ICKey(EnergyLevel, St)),
        efficiency_data,
        make_constraint_name(ENERGY_LIMIT, St),
        (
            make_variable_name(ACTIVE_POWER_IN, St),
            make_variable_name(ACTIVE_POWER_OUT, St),
            make_variable_name(ENERGY, St),
        ),
    )
    return
end

function energy_target_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    ::Type{EndOfPeriodEnergyTarget},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {St <: PSY.BatteryEMS, S <: PM.AbstractPowerModel}
    constraint_infos_target =
        Vector{DeviceEnergyTargetConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        constraint_info_target = DeviceEnergyTargetConstraintInfo(
            PSY.get_name(d),
            PSY.get_rating(d),
            PSY.get_storage_target(d),
        )
        constraint_infos_target[ix] = constraint_info_target
    end

    energy_soft_target(
        optimization_container,
        constraint_infos_target,
        make_constraint_name(ENERGY_TARGET, St),
        (make_variable_name(ENERGY, St), make_variable_name(ENERGY_TARGET_SLACK, St)),
    )

    return
end

############################ Energy Management constraints ######################################

function AddCostSpec(
    ::Type{PSY.BatteryEMS},
    ::Type{EndOfPeriodEnergyTarget},
    optimization_container::OptimizationContainer,
)
    variable_cost_func = x -> -PSY.get_energy_value(x) + PSY.get_penalty_cost(x)
    return AddCostSpec(;
        variable_type = EnergyVariable,
        component_type = PSY.BatteryEMS,
        variable_cost = variable_cost_func,
        multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    )
end

function add_to_cost!(
    optimization_container::OptimizationContainer,
    spec::AddCostSpec,
    cost_data::Float64,
    component::T,
) where {T <: PSY.Storage}
    component_name = PSY.get_name(component)
    time_steps = model_time_steps(optimization_container)

    linear_gen_cost!(
        optimization_container,
        make_variable_name(spec.variable_type, spec.component_type),
        component_name,
        cost_data,
        time_steps[end],
    )
    return
end

"""
Add variables to the OptimizationContainer for a Storage device.
"""
function cost_function!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward} = nothing,
) where {T <: PSY.Storage, U <: AbstractStorageFormulation}
    for d in devices
        spec = AddCostSpec(T, U, optimization_container)
        @debug T, spec
        add_to_cost!(optimization_container, spec, spec.variable_cost(d), d)
    end
    return
end


text:
@testset "Simple DSOS/SDSOS example with $(factory.constructor)" for factory in sdp_factories
    function sdsos_simple(cone)
        @testset "with $cone" begin
            model = SOSModel(factory)
            @variable(model, α)
            @objective(model, Max, α)
            @polyvar x y
            @constraint(model, x^2 + α*x*y + y^2 in SDSOSCone())
            optimize!(model)
            @test JuMP.value(α) ≈ 2.0 rtol=1e-5
            @test JuMP.objective_value(model) ≈ 2.0 rtol=1e-5
        end
    end
    sdsos_simple(DSOSCone())
    sdsos_simple(SDSOSCone())
    sdsos_simple(SOSCone())
end


text:
function coctailsort(a::Vector)
    b = copy(a)
    isordered = false
    lo, hi = 1, length(b)
    while !isordered && hi > lo
        isordered = true
        for i in lo+1:hi
            if b[i] < b[i-1]
                b[i-1], b[i] = b[i], b[i-1]
                isordered = false
            end
        end
        hi -= 1
        if isordered || hi ≤ lo break end
        for i in hi:-1:lo+1
            if b[i-1] > b[i]
                b[i-1], b[i] = b[i], b[i-1]
                isordered = false
            end
        end
        lo += 1
    end
    return b
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", cocktailsort(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Gallium

using POMDPs
using ARDESPOT
using ProfileView
using POMDPModels

pomdp = BabyPOMDP()
bounds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = DESPOTSolver(epsilon_0=0.0,
                      bounds=bounds,
                      rng=MersenneTwister(4),
                      random_source=MersenneSource(500, 1, MersenneTwister(4))
                     )
p = solve(solver, pomdp)
b0 = initial_state_distribution(pomdp)
@enter ARDESPOT.build_despot(p, b0)


text:
module MatFac

using Flux, Functors, Zygote, ChainRules, ChainRulesCore, CUDA, BSON

include("batch_iter.jl")
include("util.jl")
include("viewable.jl")
include("noise_models.jl")
include("model.jl")
include("fit.jl")

end

max_stars_repo_name: dpmerrell/BatchMatFac

max_issues_repo_name: dpmerrell/BatchMatFac

max_forks_repo_name: dpmerrell/BatchMatFac


text:
include(joinpath("..", "biped.jl"))

# test kinematics
q = rand(nq)

norm(kinematics_1(model, q, body = :torso, mode = :ee)
    - [q[1] - model.l_torso * sin(q[3]); q[2] + model.l_torso * cos(q[3])])
norm(kinematics_1(model, q, body = :torso, mode = :com)
    - [q[1] - model.d_torso * sin(q[3]); q[2] + model.d_torso * cos(q[3])])
norm(kinematics_1(model, q, body = :thigh_1, mode = :ee)
    - [q[1] + model.l_thigh1 * sin(q[4]); q[2] - model.l_thigh1 * cos(q[4])])
norm(kinematics_1(model, q, body = :thigh_1, mode = :com)
    - [q[1] + model.d_thigh1 * sin(q[4]); q[2] - model.d_thigh1 * cos(q[4])])
norm(kinematics_1(model, q, body = :thigh_2, mode = :ee)
    - [q[1] + model.l_thigh2 * sin(q[6]); q[2] - model.l_thigh2 * cos(q[6])])
norm(kinematics_1(model, q, body = :thigh_2, mode = :com)
    - [q[1] + model.d_thigh2 * sin(q[6]); q[2] - model.d_thigh2 * cos(q[6])])

norm(kinematics_2(model, q, body = :calf_1, mode = :ee)
    - [q[1] + model.l_thigh1 * sin(q[4]) + model.l_calf1 * sin(q[5]);
    q[2] - model.l_thigh1 * cos(q[4]) - model.l_calf1 * cos(q[5])])
norm(kinematics_2(model, q, body = :calf_1, mode = :com)
    - [q[1] + model.l_thigh1 * sin(q[4]) + model.d_calf1 * sin(q[5]);
    q[2] - model.l_thigh1 * cos(q[4]) - model.d_calf1 * cos(q[5])])
norm(kinematics_2(model, q, body = :calf_2, mode = :ee)
    - [q[1] + model.l_thigh2 * sin(q[6]) + model.l_calf2 * sin(q[7]);
    q[2] - model.l_thigh2 * cos(q[6]) - model.l_calf2 * cos(q[7])])
norm(kinematics_2(model, q, body = :calf_2, mode = :com)
    - [q[1] + model.l_thigh2 * sin(q[6]) + model.d_calf2 * sin(q[7]);
    q[2] - model.l_thigh2 * cos(q[6]) - model.d_calf2 * cos(q[7])])

k1(z) = kinematics_1(model, z, body = :torso, mode = :ee)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :torso, mode = :ee))

k1(z) = kinematics_1(model, z, body = :torso, mode = :com)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :torso, mode = :com))

k1(z) = kinematics_1(model, z, body = :thigh_1, mode = :ee)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_1, mode = :ee))

k1(z) = kinematics_1(model, z, body = :thigh_1, mode = :com)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_1, mode = :com))

k1(z) = kinematics_1(model, z, body = :thigh_2, mode = :ee)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_2, mode = :ee))

k1(z) = kinematics_1(model, z, body = :thigh_2, mode = :com)
norm(ForwardDiff.jacobian(k1,q) - jacobian_1(model, q, body = :thigh_2, mode = :com))

k2(z) = kinematics_2(model, z, body = :calf_1, mode = :ee)
norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_1, mode = :ee))

k2(z) = kinematics_2(model, z, body = :calf_1, mode = :com)
norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_1, mode = :com))

k2(z) = kinematics_2(model, z, body = :calf_2, mode = :ee)
norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_2, mode = :ee))

k2(z) = kinematics_2(model, z, body = :calf_2, mode = :com)
norm(ForwardDiff.jacobian(k2,q) - jacobian_2(model, q, body = :calf_2, mode = :com))

q̇ = rand(nq)
lagrangian(model, q, q̇)
dLdq(model, q, q̇)
dLdq̇(model, q, q̇)

eigen(M_func(model, q))

tmp_q(z) = dLdq̇(model, z, q̇)
tmp_q̇(z) = dLdq̇(model, q, z)
norm(ForwardDiff.jacobian(tmp_q̇,q̇) - M_func(model, q))


text:
function device_model!(x,
                      output_ode::Vector{Float64},
                      voltage_r,
                      voltage_i,
                      current_r,
                      current_i,
                      ix_range::UnitRange{Int64},
                      ode_range::UnitRange{Int64},
                      controls,
                      device::DynG,
                      sys::DynamicSystem) where {DynG <: DynGenerator}
       #Obtain local device states
    n_states = total_device_states(device)
    device_states = @view x[ix_range]
       #Obtain references
    sys_Sbase = get_sys_base(sys)
    sys_f = get_sys_f(sys)

    #Update Voltage data
    device.inner_vars[VR_gen_var] = voltage_r[1]
    device.inner_vars[VI_gen_var] = voltage_i[1]

       #Obtain ODEs and Mechanical Power for Turbine Governor
    mdl_tg_ode!(device_states,
                view(output_ode, ode_range),
                device)
       #Obtain ODEs for AVR
    mdl_pss_ode!(device_states,
                 view(output_ode, ode_range),
                 device)
       #Obtain ODEs for AVR
    mdl_avr_ode!(device_states,
                view(output_ode, ode_range),
                device)
                   #Obtain ODEs for Machine
    mdl_machine_ode!(device_states,
                    view(output_ode, ode_range),
                    current_r,
                    current_i,
                    sys_Sbase,
                    sys_f,
                    device)
                       #Obtain ODEs for Shaft
    mdl_shaft_ode!(device_states,
                   view(output_ode, ode_range),
                   sys_f,
                   device)

    return

end

function device_model!(voltage_r,
                       voltage_i,
                       current_r,
                       current_i,
                       device::StaticSource,
                       sys::DynamicSystem)

    mdl_source!(voltage_r, voltage_i, current_r, current_i, device, sys)

    return
end

function device_model!(voltage_r,
                       voltage_i,
                       current_r,
                       current_i,
                       device::PSY.PowerLoad,
                       sys::DynamicSystem)

    mdl_Zload!(voltage_r, voltage_i, current_r, current_i, device, sys)

    return
end

function device_model!(x,
                        output_ode::Vector{Float64},
                        voltage_r,
                        voltage_i,
                        current_r,
                        current_i,
                        ix_range::UnitRange{Int64},
                        ode_range::UnitRange{Int64},
                        controls,
                        device::DynI,
                        sys::DynamicSystem) where {DynI <: DynInverter}

        #Obtain local device states
        n_states = total_device_states(device)
        device_states = @view x[ix_range]

        #Obtain references
        Sbase = get_sys_base(sys)
        sys_f = get_sys_f(sys)

        #Update Voltage data
        device.inner_vars[VR_inv_var] = voltage_r[1]
        device.inner_vars[VI_inv_var] = voltage_i[1]
        #device.inner_vars[Vh_var] = sqrt(voltage_r[1]^2 + voltage_i[1]^2)

        #Update V_ref
        V_ref = get_inverter_Vref(device)
        device.inner_vars[v_control_var] = V_ref

        #Obtain ODES for DC side
        mdl_DCside_ode!(device)

        #Obtain ODEs for OuterLoop
        mdl_outer_ode!(device_states,
                       view(output_ode, ode_range),
                       sys_f,
                       device)

#Obtain ODEs for PLL
mdl_freq_estimator_ode!(device_states,
                        view(output_ode, ode_range),
                        sys_f,
                        device)

#Obtain inner controller ODEs and modulation commands
mdl_VScontrol_ode!(device_states,
                   view(output_ode, ode_range),
                   device)

#Obtain converter relations
mdl_converter_ode!(device)

#Obtain ODEs for output filter
mdl_filter_ode!(device_states,
                        view(output_ode, ode_range),
                        current_r,
                        current_i,
                        Sbase,
                        sys_f,
                        device)

        return
end


text:
# Author: AJ Tyler
# Date: 24/01/22

# This script calls the compare and classify function from the main module.

include("main.jl")
using .main
compare("DianaConfigs")
classify("DianaConfigs")


text:
using TakingImmutabilitySeriously
using Test

@testset "TakingImmutabilitySeriously.jl" begin
    # Write your own tests here.
end


text:
@testset "csr_to_sparse" begin
    ei = [1;2;3]
    ej = [3;4;4]
    ev = [8;9;10]
    (rp,ci,ai,m) = sparse_to_csr(ei,ej,ev)
    (nzi,nzj,nzv) = csr_to_sparse(rp,ci,ai)
    A = sparse(nzi,nzj,nzv,length(rp)-1,maximum(ci))

    # more tests added here
    # clique to sparse test
    rp = collect(1:5:26)
    #ci = vec(reshape(repmat(1:5,5,1)',25,1))
    ci = copy(vec(repeat(1:5,5,1)'))
    ai = ones(Int64,25)
    A = csr_to_sparse_matrix(rp,ci,ai,5,5)
    @test Matrix(A) == ones(5,5)

    A = csr_to_sparse_matrix(rp,ci,ai)
    @test Matrix(A) == ones(5,5)

    # 100 random trials
    for t = 1:100
        A = sprand(100,80,0.01)
        (rp,ci,ai) = sparse_to_csr(A)
        A2 = csr_to_sparse_matrix(rp,ci,ai,100,80)
        @test A == A2
    end
end