text:
@inline Base.zero(::Mask{W,U}) where {W,U} = Mask{W}(zero(U))

@inline extract_data(m::Mask) = m.u
@inline Base.:(&)(m1::Mask{W}, m2::Mask{W}) where {W} = Mask{W}(m1.u & m2.u)
@inline Base.:(&)(m::Mask{W}, u::UIntTypes) where {W} = Mask{W}(m.u & u)
@inline Base.:(&)(u::UIntTypes, m::Mask{W}) where {W} = Mask{W}(u & m.u)

@inline Base.:(&)(m::Mask{W}, b::Bool) where {W} = Mask{W}(b ? m.u : zero(m.u))
@inline Base.:(&)(b::Bool, m::Mask{W}) where {W} = Mask{W}(b ? m.u : zero(m.u))

@inline Base.:(|)(m1::Mask{W}, m2::Mask{W}) where {W} = Mask{W}(m1.u | m2.u)
@inline Base.:(|)(m::Mask{W}, u::UIntTypes) where {W} = Mask{W}(m.u | u)
@inline Base.:(|)(u::UIntTypes, m::Mask{W}) where {W} = Mask{W}(u | m.u)

@inline Base.:(|)(m::Mask{W,U}, b::Bool) where {W,U} = b ? max_mask(Mask{W,U}) : m
@inline Base.:(|)(b::Bool, m::Mask{W,U}) where {W,U} = b ? max_mask(Mask{W,U}) : m
@inline Base.:(|)(m::Mask{16,UInt16}, b::Bool) where {W} = Mask{W}(b ? 0xffff : m.u)
@inline Base.:(|)(b::Bool, m::Mask{16,UInt16}) where {W} = Mask{W}(b ? 0xffff : m.u)
@inline Base.:(|)(m::Mask{8,UInt8}, b::Bool) where {W} = Mask{W}(b ? 0xff : m.u)
@inline Base.:(|)(b::Bool, m::Mask{8,UInt8}) where {W} = Mask{W}(b ? 0xff : m.u)
@inline Base.:(|)(m::Mask{4,UInt8}, b::Bool) where {W} = Mask{W}(b ? 0x0f : m.u)
@inline Base.:(|)(b::Bool, m::Mask{4,UInt8}) where {W} = Mask{W}(b ? 0x0f : m.u)
@inline Base.:(|)(m::Mask{2,UInt8}, b::Bool) where {W} = Mask{W}(b ? 0x03 : m.u)
@inline Base.:(|)(b::Bool, m::Mask{2,UInt8}) where {W} = Mask{W}(b ? 0x03 : m.u)

@inline Base.:(⊻)(m1::Mask{W}, m2::Mask{W}) where {W} = Mask{W}(m1.u ⊻ m2.u)
@inline Base.:(⊻)(m::Mask{W}, u::UIntTypes) where {W} = Mask{W}(m.u ⊻ u)
@inline Base.:(⊻)(u::UIntTypes, m::Mask{W}) where {W} = Mask{W}(u ⊻ m.u)

@inline Base.:(⊻)(m::Mask{W}, b::Bool) where {W} = Mask{W}(b ? ~m.u : m.u)
@inline Base.:(⊻)(b::Bool, m::Mask{W}) where {W} = Mask{W}(b ? ~m.u : m.u)

@inline Base.:(<<)(m::Mask{W}, i) where {W} = Mask{W}(m.u << i)
@inline Base.:(>>)(m::Mask{W}, i) where {W} = Mask{W}(m.u >> i)
@inline Base.:(>>>)(m::Mask{W}, i) where {W} = Mask{W}(m.u >>> i)

@inline Base.:(~)(m::Mask{W}) where {W} = Mask{W}( ~m.u )
@inline Base.:(!)(m::Mask{W}) where {W} = Mask{W}( ~m.u )

@inline Base.:(==)(m1::Mask{W}, m2::Mask{W}) where {W} = m1.u == m2.u
@inline Base.:(==)(m::Mask{W}, u::UIntTypes) where {W} = m.u == u
@inline Base.:(==)(u::UIntTypes, m::Mask{W}) where {W} = u == m.u
@inline Base.:(!=)(m1::Mask{W}, m2::Mask{W}) where {W} = m1.u != m2.u
@inline Base.:(!=)(m::Mask{W}, u::UIntTypes) where {W} = m.u != u
@inline Base.:(!=)(u::UIntTypes, m::Mask{W}) where {W} = u != m.u

@inline Base.count_ones(m::Mask) = count_ones(m.u)

@inline Base.:(+)(m::Mask, i::Integer) = i + count_ones(m.u)
@inline Base.:(+)(i::Integer, m::Mask) = i + count_ones(m.u)

function mask_type(W)
    if W <= 8
        return UInt8
    elseif W <= 16
        return UInt16
    elseif W <= 32
        return UInt32
    elseif W <= 64
        return UInt64
    else#if W <= 128
        return UInt128
    end
end
mask_type(::Val{4}) = UInt8
mask_type(::Val{8}) = UInt8
mask_type(::Val{16}) = UInt16
mask_type(::Val{32}) = UInt32
mask_type(::Val{64}) = UInt64

@generated function mask_type(::Type{T}, ::Val{P}) where {T,P}
    mask_type(pick_vector_width(P, T))
end
@generated function mask_type(::Type{T}) where {T}
    mask_type(pick_vector_width(T))
end

@generated function max_mask(::Type{T}) where {T}
    W = pick_vector_width(T)
    U = mask_type(W)
    Mask{W,U}(one(U)<<W - one(U))
end
@generated max_mask(::Type{Mask{W,U}}) where {W,U} = Mask{W,U}(one(U)<<W - one(U))

@generated function mask(::Type{T}, l::Integer) where {T}
    M = mask_type(T)
    W = pick_vector_width(T)
    # tup = Expr(:tuple, [Base.unsafe_trunc(M, 1 << w - 1) for w in 0:W]...)
    quote
        $(Expr(:meta,:inline))
        # @inbounds $tup[rem+1]
        # rem = valrem(Val{$W}(), l - 1) + 1
        # Mask{$W,$M}(one($M) << (rem & $(typemax(M))) - $(one(M)))
        rem = valrem(Val{$W}(), (l % $M) - one($M)) + one($M)
        Mask{$W,$M}($(typemax(M)) >>> ($(M(8sizeof(M))) - rem))
    end
end

@generated function mask(::Val{W}, l::Integer) where {W}
    M = mask_type(W)
#    W = pick_vector_width(T)
    # tup = Expr(:tuple, [Base.unsafe_trunc(M, 1 << w - 1) for w in 0:W]...)
    quote
        $(Expr(:meta,:inline))
        # @inbounds $tup[rem+1]
        # rem = valrem(Val{$W}(), l % $M)
        # Mask{$W,$M}($(typemax(M)) >>> ($(M(8sizeof(M))) - rem))
        # rem = valrem(Val{$W}(), l - 1) + 1
        rem = valrem(Val{$W}(), (l % $M) - one($M)) + one($M)
        Mask{$W,$M}($(typemax(M)) >>> ($(M(8sizeof(M))) - rem))
        # Mask{$W,$M}(one($M) << rem)
        # Mask{$W,$M}(one($M) << (rem) - $(one(M)))
    end
end
@generated mask(::Val{W}, ::Static{L}) where {W, L} = mask(Val(W), L)

unstable_mask(W, rem) = mask(Val(W), rem)

@generated function masktable(::Val{W}, rem::Integer) where {W}
    masks = Expr(:tuple)
    for w ∈ 0:W-1
        push!(masks.args, extract_data(unstable_mask(W, w == 0 ? W : w)))
    end
    Expr(
        :block,
        Expr(:meta,:inline),
        Expr(:call, Expr(:curly, :Mask, W), Expr(
            :macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)),
            Expr(:call, :getindex, masks, Expr(:call, :+, 1, Expr(:call, :valrem, Expr(:call, Expr(:curly, W)), :rem)))
        ))
    )
end

@inline tomask(m::Unsigned) = m
@inline tomask(m::Mask) = m
@generated function tomask(v::Vec{W,Bool}) where {W}
    usize = W > 8 ? nextpow2(W) : 8
    utyp = "i$(usize)"
    U = mask_type(W)
    instrs = String[]
    push!(instrs, "%bitvec = trunc <$W x i8> %0 to <$W x i1>")
    if usize == W
        push!(instrs, "%mask = bitcast <$W x i1> %bitvec to i$(W)")
    else
        push!(instrs, "%maskshort = bitcast <$W x i1> %bitvec to i$(W)")
        push!(instrs, "%mask = zext i$(W) %maskshort to i$(usize)")
    end
    push!(instrs, "ret i$(usize) %mask")
    quote
        $(Expr(:meta, :inline))
        Mask{$W}(Base.llvmcall(
            $(join(instrs, "\n")), $U, Tuple{Vec{$W,Bool}}, v
        ))
    end
end
@inline tomask(v::AbstractStructVec{<:Any,Bool}) = tomask(extract_data(v))

@inline getindexzerobased(m::Mask, i) = (m.u >>> i) % Bool
@inline function Base.getindex(m::Mask{W}, i::Integer) where {W}
    @boundscheck i > W && throw(BoundsError(m, i))
    getindexzerobased(m, i - 1)
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{2})
    Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >> 1
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{4})
    Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >> 2
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{8})
    Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), i.i >> 3
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{16})
    Base.unsafe_convert(Ptr{UInt16}, ptr.ptr), i.i >> 4
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{32})
    Base.unsafe_convert(Ptr{UInt32}, ptr.ptr), i.i >> 5
end
@inline function ptr_index(ptr::AbstractBitPointer, i::_MM{64})
    Base.unsafe_convert(Ptr{UInt64}, ptr.ptr), i.i >> 6
end

@inline function bitload(ptr::AbstractBitPointer, i::_MM{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    Mask{W}(vload(ptr, ind))
end
@inline bitload(ptr::AbstractBitPointer, i, ::Union{UIntTypes,Mask}) = bitload(ptr, i)
@inline bitload(ptr::AbstractBitPointer, i::Integer) = getindexzerobased(bitload(ptr, _MM{8}(i)), i & 7)

# @inline function vstore!(ptr::AbstractBitPointer, m::Mask{8}, i::Integer)
    # vstore!(Base.unsafe_convert(Ptr{UInt8}, ptr.ptr), (m.u % Bool), i)
# end
@inline function bitstore!(ptr::AbstractBitPointer, m::Mask{W}, i::_MM{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    vstore!(ptr, m.u, ind)
end
@inline function bitstore!(ptr::AbstractBitPointer, m::Mask{W}, i::_MM{W}, mask::Mask{W}) where {W}
    ptr, ind = ptr_index(ptr, i)
    vstore!(ptr, m.u, ind)
end

@generated function bitstore!(
    ptr::Ptr{T}, v::Mask{W,U}, mask::Mask{W,U}
) where {W,T,U<:UIntTypes}
    @assert isa(Aligned, Bool)
    ptyp = JuliaPointerType
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$W"
    decls = String[]
    instrs = String[]
    align = sizeof(U)
    push!(instrs, "%ptr = inttoptr $ptyp %0 to <$W x i1>*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_input %2 to <$W x i1>")
    else
        push!(instrs, "%vtrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%masktrunc = trunc $mtyp_input %2 to $mtyp_trunc")
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$W x i1>")
    end
    push!(decls,
        "declare void @llvm.masked.store.v$(W)i1(<$W x i1>, <$W x i1>*, i32, <$W x i1>)"
    )
    push!(instrs,
        "call void @llvm.masked.store.v$(W)i1(<$W x i1> %v, <$W x i1>* %ptr, i32 $align, <$W x i1> %mask)"
    )
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, $U, $U},
            ptr, v.u, mask.u)
    end
end
@generated function bitstore!(
    ptr::Ptr{T}, v::Mask{W,U}, ind::I, mask::Mask{W,U}
) where {W,T,I<:Integer,U<:UIntTypes}
    @assert isa(Aligned, Bool)
    ptyp = JuliaPointerType
    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$W"
    decls = String[]
    instrs = String[]
    align = sizeof(U)
    if mtyp_input == mtyp_trunc
        push!(instrs, "%ptr = inttoptr $ptyp %0 to <$W x i1>*")
        push!(instrs, "%offsetptr = getelementptr inbounds <$W x i1>, <$W x i1>* %ptr, i$(8sizeof(I)) %2")
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_input %3 to <$W x i1>")
    else
        push!(instrs, "%ptr = inttoptr $ptyp %0 to i$(8align)*")
        push!(instrs, "%tempptr = getelementptr inbounds i$(8align),  i$(8align)* %ptr, i$(8sizeof(I)) %2")
        push!(instrs, "%offsetptr = bitcast i$(8align) %tempptr to <$W x i1>*")
        push!(instrs, "%vtrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%masktrunc = trunc $mtyp_input %3 to $mtyp_trunc")
        push!(instrs, "%v = bitcast $mtyp_input %1 to <$W x i1>")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$W x i1>")
    end
    push!(decls,
        "declare void @llvm.masked.store.v$(W)i1(<$W x i1>, <$W x i1>*, i32, <$W x i1>)"
    )
    push!(instrs,
        "call void @llvm.masked.store.v$(W)i1(<$W x i1> %v, <$W x i1>* %offsetptr, i32 $align, <$W x i1> %mask)"
    )
    push!(instrs, "ret void")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Cvoid, Tuple{Ptr{$T}, $U, $I, $U},
            ptr, v.u, ind, mask.u)
    end
end

@inline vload(ptr::AbstractBitPointer, i::Tuple) = bitload(ptr, offset(ptr, staticm1(i)))
@inline vload(ptr::AbstractBitPointer, i::Tuple, ::Mask) = vload(ptr, i)
@inline function vload(bptr::PackedStridedBitPointer{1}, (i,j)::Tuple{_MM{W},<:Integer}) where {W}
    j -= 1
    s = bptr.strides[1]
    shift = (s * j) & (W - 1)
    U = mask_type(Val{W}())
    UW = widen(U)
    ptr, ind = ptr_index(bptr, _MM{W}(i.i - 1 + j*s))
    u = vload(Base.unsafe_convert(Ptr{UW}, gep(ptr, ind)))
    # @show ind, shift, u
    Mask{W}((u >>> shift) % U)
end
# @inline function vstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, (i,j)::Tuple{_MM{W},<:Integer}) where {W}
#     j -= 1
#     s = bptr.strides[1]
#     shift = (s * j) & (W - 1)
#     U = mask_type(Val{W}())
#     UW = widen(U)
#     ptr, ind = ptr_index(bptr, _MM{W}(i.i - 1 + j*s))
#     um = ((vload(Base.unsafe_convert(Ptr{UW}, gep(ptr, ind))) )
#     u = (v.u % UW) << shift
#     # @show ind, shift, u
#     vstore!(Base.unsafe_convert(Ptr{UW}, gep(ptr, ind)), u | um)

# end
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::SVec{W,Bool}, i::Tuple{_MM{W},<:Integer}) where {W} = vstore!(bptr, tomask(v), i)
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{_MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, v, i)
# @inline vstore!(bptr::PackedStridedBitPointer{1}, v::SVec{W,Bool}, i::Tuple{_MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, tomask(v), i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{_MM{W},<:Integer}) where {W} = vstore!(bptr, v, i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::Mask{W}, i::Tuple{_MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, v, i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::SVec{W,Bool}, i::Tuple{_MM{W},<:Integer}) where {W} = vstore!(bptr, tomask(v), i)
# @inline vnoaliasstore!(bptr::PackedStridedBitPointer{1}, v::SVec{W,Bool}, i::Tuple{_MM{W},<:Integer}, ::AbstractMask) where {W} = vstore!(bptr, tomask(v), i)

@inline vstore!(ptr::AbstractBitPointer, v::Mask, i::Tuple) = bitstore!(ptr, v, offset(ptr, staticm1(i)))
@inline vstore!(ptr::AbstractBitPointer, v::Mask, i::Tuple, u::AbstractMask) = bitstore!(ptr, v, offset(ptr, staticm1(i)), tomask(u))
@inline vnoaliasstore!(ptr::AbstractBitPointer, v::Mask, i::Tuple) = bitstore!(ptr, v, offset(ptr, staticm1(i)))
@inline vnoaliasstore!(ptr::AbstractBitPointer, v::Mask, i::Tuple, u::AbstractMask) = bitstore!(ptr, v, offset(ptr, staticm1(i)), tomask(u))
@inline vstore!(ptr::AbstractBitPointer, v::SVec{<:Any,Bool}, i::Tuple) = bitstore!(ptr, tomask(v), offset(ptr, staticm1(i)))
@inline vstore!(ptr::AbstractBitPointer, v::SVec{<:Any,Bool}, i::Tuple, u::AbstractMask) = bitstore!(ptr, tomask(v), offset(ptr, staticm1(i)), tomask(u))
@inline vnoaliasstore!(ptr::AbstractBitPointer, v::SVec{<:Any,Bool}, i::Tuple) = bitstore!(ptr, tomask(v), offset(ptr, staticm1(i)))
@inline vnoaliasstore!(ptr::AbstractBitPointer, v::SVec{<:Any,Bool}, i::Tuple, u::AbstractMask) = bitstore!(ptr, tomask(v), offset(ptr, staticm1(i)), tomask(u))


text:
# This is a 1D simple test
include("../src/BTG_main.jl")
include("../src/validation.jl")
include("../src/kernels/kernel.jl")

function test_fun_1d(x)
    return sin(x) + 1.1
end

# training points
x_train = [-4, -3, -2, 0, 1, 2, 3, 4]
y_train = test_fun_1d.(x_train)
n_train = size(x_train,1)
training_data = [x_train, y_train]

# Test data
x_test = range(-5, stop = 5, step = 0.23)
y_test = test_fun_1d.(x_test)
testing_data = [x_test, y_test];

# model setting
kernel = kernel_BTG
nonlinfun = BoxCox

# cross validation
# print("STARTING CROSS VALIDATION \n")
# cross_validation(training_data, kernel, nonlinfun)

# call the BTG model
print("STARTING PREDICTION \n")
BTG_main(training_data, testing_data, kernel, nonlinfun)
println("FINISHING PREDICTION")

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
	callout()

	receives a String message
	produces a sound representation of the message through speakers
"""
function callout(mess::String)
	# padsp is necessary for me to play through pulseaudio
	# ignorestatus is useful since flite may perform the audio but does not exit cleanly with 0
	cmd = ignorestatus(`padsp flite "$mess"`)
	run(cmd);
end
#
"""
	checkmem()

	checks memory remaining as percentage (headroom)
	returns decimal fraction in range 0.0:1.0
"""
function checkmem()
	# -b returns values in bytes
	cmd = `free -b`
	lins = readlines(cmd)
	lina = split(lins[2]," ",keepempty=false)
	# available divided by total
	return parse(Int,lina[7])/parse(Int,lina[2])
end
"""
	checkcpu()

	reports back the temperature of the CPU
"""
function checkcpu()
	# function to check state of cpu temperature
	# assumes that Julia is using all cores equally
	cmd = `sensors`
	# alternatively get the json output from sensors
	sens = readlines(cmd)
	# examine only those lines related to cores
	cs = sens[contains.(sens,"Core")]
	# create a vector to store the individual core temperatures
	ts = Vector{Real}()
	# reg expression to extract real values
	R = r"\d+\.\d+"
	ncores = length(cs)
	for i in 1:ncores
		f = findall(R,cs[i])
		t = cs[i][f[1]]
		push!(ts,parse(Float32,t))
	end
	# return the average of the core temperatures
	return sum(ts)/ncores
end


text:
module m7test13a
using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNonlinearModule: totalLagrangean2current!
using LinearAlgebra: norm
using BenchmarkTools
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3

    m1 = FinEtoolsDeforNonlinear.MatDeforNeohookeanModule.MatDeforNeohookean(mr, E, nu)
    m2 = FinEtoolsDeforNonlinear.MatDeforNeohookeanADModule.MatDeforNeohookeanAD(mr, E, nu)
    # @show m
    update! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.update!
    tangentmoduli! = FinEtoolsDeforNonlinear.MatDeforNonlinearModule.tangentmoduli!

    stress = fill(0.0, 6)
    D = fill(0.0, 6, 6)
    output = FFlt[]
    statev = FFlt[]
    tn = 0.0
    dtn = 0.0
    loc = [0.0 0.0 0.0]
    label = 0
    quantity=:nothing
    Fn1 = [1.1 0 0; 0 1.0 0; 0 0 1.0]
    Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]

    update!(m1, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[863901.0097711232, 349914.99578510894, 349914.99578510894, 0.0, 0.0, 0.0]) / E < 1.0e-7
    tangentmoduli!(m1, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @btime $tangentmoduli!($m1, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    D1 = deepcopy(D)

    update!(m2, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # @show stress
    @test norm(stress-[863901.0097711232, 349914.99578510894, 349914.99578510894, 0.0, 0.0, 0.0] ) / E < 1.0e-7
    tangentmoduli!(m2, D, statev, Fn1, Fn, tn, dtn, loc, label)
    @btime $tangentmoduli!($m2, $D, $statev, $Fn1, $Fn, $tn, $dtn, $loc, $label)
    D2 = deepcopy(D)
    @test norm(D1-D2) / E < 1.0e-7

@btime totalLagrangean2current!($D, $D2, rand(3, 3))

    # Fn1 = [1.0 0 0; 0 1.01 0; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[40184.6915460777, 95648.94230769236, 40184.69154607771, 0.0, 0.0, 0.0] ) / E < 1.0e-7

    # Fn1 = [1.0 0 0; 0 1.0 0.01; 0 0 1.0]
    # Fn = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    # update!(m, statev, stress, output, Fn1, Fn, tn, dtn, loc, label, quantity)
    # # @show stress
    # @test norm(stress-[201.92307692305477, 740.4317307692086, 471.1538461537944, 0.0, 0.0, 26927.78846153846]) / E < 1.0e-7
end
end
using .m7test13a
m7test13a.test()


text:
using Pkg; Pkg.activate(joinpath(@__DIR__, ".."))
using CSV, DataFrames

const DATADIR = joinpath(@__DIR__, "..", "data")
d1 = CSV.read(joinpath(DATADIR, "order_95966_data.txt"),
              DataFrame, missingstring = "-") |> dropmissing
d2 = CSV.read(joinpath(DATADIR, "order_95967_data.txt"),
              DataFrame, missingstring = "-") |> dropmissing
d = vcat(d1, d2)
rename!(d, :fu3010h1 => :wind_peak,
           :prestah0 => :pressure,
           :tre200h0 => :temperature,
           :rre150h0 => :precipitation,
           :sre000h0 => :sunshine_duration,
           :fu3010h0 => :wind_mean,
           :dkl010h0 => :wind_direction)
gs = [rename!(DataFrame(g), [n => "$(g.stn[1])_$n"
                             for n in names(g)[3:end]]...)[:, 2:end]
      for g in groupby(d, :stn)]
df = innerjoin(gs..., on = :time) |> dropmissing
train = df.time .< 2018000000;
CSV.write(joinpath(DATADIR, "weather2015-2018.csv"), df[train, :])
CSV.write(joinpath(DATADIR, "weather2019-2020.csv"), df[(!).(train), :])

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
using LinearAlgebra
using SparseArrays
using CUDA
using Printf

include("bicgstab.jl")

function square_preconditioned(n :: Int=10)
        A   = ones(n, n) + (n-1) * Matrix(I, n, n)
        b   = 10.0 * [1:n;]
        M⁻¹ = 1/n * Matrix(I, n, n)
        return A, b, M⁻¹
end

A, b, N = square_preconditioned(100)
x = similar(b); r = similar(b)
spA = sparse(A)
spN = sparse(N)
y, iter, status = bicgstab(spA, b, spN, x)
@show y
@show iter
@show status

cuspA = CUDA.CUSPARSE.CuSparseMatrixCSR(spA)
cuspN = CUDA.CUSPARSE.CuSparseMatrixCSR(spN)
cux = CuVector(x)
cub = CuVector(b)
cuy, iter, status = bicgstab(cuspA, cub, cuspN, cux)
@show cuy
@show iter
@show status


text:
# Get all state switing indices
function get_all_indices(state_time_series)
	tᵦ = []
	# Initialise first index
	push!(tᵦ,1)

	# Push body indices
	for i in 2:size(state_time_series,1)
		if state_time_series[i-1] == 0 && state_time_series[i] == 1
			push!(tᵦ,i)
		elseif state_time_series[i-1] == 1 && state_time_series[i] == 0
			push!(tᵦ,i)
		end
	end
	push!(tᵦ,size(state_time_series,1))
	return tᵦ
end

# Get indices per each state
function get_indices(state_time_series,t⃗)

	indss = get_all_indices(state_time_series)[1:end-1]

	# Separate times
	# If first state is membrane
	if state_time_series[1] == 0
		memᵢ = indss[1:2:end]
		basᵢ = indss[2:2:end]
	# If first state is basal
	else
		memᵢ = indss[2:2:end]
		basᵢ = indss[1:2:end]
	end

	return telegraph_points(memᵢ, basᵢ)
end

# Get associated times
function get_times(gi::telegraph_points,t⃗,s₀)
	t̄ = size(t⃗,1)
	tᵢ = sort(vcat(gi.i₀,gi.i₁,t̄))
	Δts = diff(t⃗[tᵢ])
	sΔts = size(Δts,1)

	# Turn times in each type
	if mod(sΔts,2) == 0
		rΔts = reshape(Δts,2,sΔts÷2)'
	elseif mod(sΔts,2) == 1
		push!(Δts,1)
		rΔts = reshape(Δts,2,(sΔts÷2+1))'
	else
		# Do nothing
		error("This is weird, this is supposeed to be divisible by 2, instead it is $(sΔts)")
	end

	# Which is s = 0 and which is s = 1
	if s₀ == 0
		return telegraph_time_dist(rΔts[:,1], rΔts[1:end-1,2])
	elseif s₀ == 1
		return telegraph_time_dist(rΔts[1:end-1,2],rΔts[:,1])

	else
		error("This is weird, this is supposeed to be divisible by 2, instead it is $(s₀)")
	end
end


text:
__precompile__()

module CoordinateTransformations

using StaticArrays

using Rotations
export RotMatrix, Quat, SpQuat, AngleAxis, RodriguesVec,
       RotX, RotY, RotZ,
       RotXY, RotYX, RotZX, RotXZ, RotYZ, RotZY,
       RotXYX, RotYXY, RotZXZ, RotXZX, RotYZY, RotZYZ,
       RotXYZ, RotYXZ, RotZXY, RotXZY, RotYZX, RotZYX

# Core methods
export compose, ∘, transform_deriv, transform_deriv_params
export Transformation, IdentityTransformation

# 2D coordinate systems and their transformations
export Polar
export PolarFromCartesian, CartesianFromPolar

# 3D coordinate systems and their transformations
export Spherical, Cylindrical
export SphericalFromCartesian, CartesianFromSpherical,
       CylindricalFromCartesian, CartesianFromCylindrical,
       CylindricalFromSpherical, SphericalFromCylindrical

# Common transformations
export AbstractAffineMap
export AffineMap, LinearMap, Translation

include("core.jl")
include("coordinatesystems.jl")
include("affine.jl")

# Deprecations
export transform
Base.@deprecate_binding AbstractTransformation Transformation
Base.@deprecate transform(transformation::Transformation, x) transformation(x)

end # module


text:
## Diagonal Preconditioners

struct DiagonalPreconditioner{D}
    diag::D
end
struct InvDiagonalPreconditioner{D}
    diag::D
end

Base.eltype(A::Union{DiagonalPreconditioner,InvDiagonalPreconditioner}) = eltype(A.diag)

function LinearAlgebra.ldiv!(A::DiagonalPreconditioner, x)
    x .= x ./ A.diag
end

function LinearAlgebra.ldiv!(y, A::DiagonalPreconditioner, x)
    y .= x ./ A.diag
end

#=
function LinearAlgebra.ldiv!(y::Matrix, A::DiagonalPreconditioner, b::Matrix)
    @inbounds @simd for j ∈ 1:size(y, 2)
        for i ∈ 1:length(A.diag)
            y[i,j] = b[i,j] / A.diag[i]
        end
    end
    return y
end
=#

function LinearAlgebra.ldiv!(A::InvDiagonalPreconditioner, x)
    x .= x .* A.diag
end

function LinearAlgebra.ldiv!(y, A::InvDiagonalPreconditioner, x)
    y .= x .* A.diag
end

#=
function LinearAlgebra.ldiv!(y::Matrix, A::InvDiagonalPreconditioner, b::Matrix)
    @inbounds @simd for j ∈ 1:size(y, 2)
        for i ∈ 1:length(A.diag)
            y[i,j] = b[i,j] * A.diag[i]
        end
    end
    return y
end
=#

LinearAlgebra.mul!(y, A::DiagonalPreconditioner, x) = LinearAlgebra.ldiv!(y, InvDiagonalPreconditioner(A.diag), x)
LinearAlgebra.mul!(y, A::InvDiagonalPreconditioner, x) = LinearAlgebra.ldiv!(y, DiagonalPreconditioner(A.diag), x)

## Compose Preconditioner

struct ComposePreconditioner{Ti,To}
    inner::Ti
    outer::To
end

Base.eltype(A::ComposePreconditioner) = promote_type(eltype(A.inner), eltype(A.outer))

function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
    @unpack inner, outer = A

    ldiv!(inner, x)
    ldiv!(outer, x)
end

function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
    @unpack inner, outer = A

    ldiv!(y, inner, x)
    ldiv!(outer, y)
end

struct InvPreconditioner{T}
    P::T
end

Base.eltype(A::InvPreconditioner) = Base.eltype(A.P)
LinearAlgebra.mul!(y, A::InvPreconditioner, x) = ldiv!(y, A.P, x)


text:
#=
  A/B testing, given lists of successes for A and B

  Testing somewhat larger instance (750 entries)

  True pA: 0.5
  True pB: 0.4

  Not too bad, but I'm not happy with the rhat values for pB.
  Here's the result for SMC():
  """
  Summary Statistics
     parameters      mean       std   naive_se      mcse       ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64   Float64   Float64       Float64

             pA    0.0626    0.0122     0.0001    0.0010   48.0229    1.0040        0.6051
             pB    0.0369    0.0086     0.0001    0.0007   41.4276    1.0830        0.5220
           diff    0.0257    0.0154     0.0002    0.0013   42.5037    1.0447        0.5355
  aGreaterThanB    0.9450    0.2280     0.0023    0.0192   49.0300    1.0013        0.6178

  Quantiles
     parameters      2.5%     25.0%     50.0%     75.0%     97.5%
         Symbol   Float64   Float64   Float64   Float64   Float64

             pA    0.0383    0.0558    0.0635    0.0760    0.0900
             pB    0.0238    0.0311    0.0387    0.0399    0.0575
           diff   -0.0039    0.0171    0.0248    0.0316    0.0562
  aGreaterThanB    0.0000    1.0000    1.0000    1.0000    1.0000

  Credible interval for diff with mass 0.9: (0.002490..0.049147)
  """

  Cf ~/pymc3/ab_test.py
     ~/blog/ab_testing2.blog
     ~/webppl/ab_testing2.wppl

=#

using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function ab_testing2(as,bs)
    n = length(as)

    pA ~ Beta(1,1) # true prob: 0.05
    pB ~ Beta(1,1) # true prob: 0.04

    for i in 1:n
        as[i] ~ flip(pA)
        bs[i] ~ flip(pB)
    end

    diff ~ Dirac(pA - pB)
    aGreaterThanB ~ Dirac(diff > 0)

end

# Generated by ./gener_data.pl
# True prob 0.05
# random_binomial(750,1,0.05);
as = [0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,
      0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
      0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0];

# True prob = 0.04
# random_binomial(750,1,0.04);
bs = [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0];

model = ab_testing2(as,bs)

num_chns = 4
# chns = sample(model, Prior(), 1000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(5), 10_000)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, SGLD(), 10_000)

display(chns)
# display(plot(chns))

credible_interval(chns,"diff",0.9)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

function aster_parse_nodes(section; strip_characters=true)
    nodes = Dict{Any, Vector{Float64}}()
    has_started = false
    for line in split(section, '\n')
        m = collect((string(m.match) for m in eachmatch(r"[\w.-]+", line)))
        if length(m) == 1
            if (m[1] == "COOR_2D") || (m[1] == "COOR_3D")
                has_started = true
                continue
            end
            if m[1] == "FINSF"
                break
            end
        end
        if !has_started # we have not found COOR_2D or COOR_3D yet.
            continue
        end
        to(T, x) = map(x -> parse(T, x), x)
        if length(m) == 4
            if strip_characters
                nodes[parse_node_id(m[1])] = to(Float64, m[2:end])
            else
                nodes[nid] = to(Float64, m[2:end])
            end
        end
    end
    return nodes
end

""" Code Aster binary file (.med). """
mutable struct MEDFile
    data :: Dict
end

function MEDFile(fn::String)
    return MEDFile(h5read(fn, "/"))
end

function get_mesh_names(med::MEDFile)
    return sort(collect(keys(med.data["FAS"])))
end

""" Convert vector of Int8 to ASCII string. """
function to_ascii(data::Vector{Int8})
    return ascii(unsafe_string(pointer(convert(Vector{UInt8}, data))))
end

function get_mesh(med::MEDFile, mesh_name::String)
    if !haskey(med.data["FAS"], mesh_name)
        @warn("Mesh $mesh_name not found from med file.")
        meshes = get_mesh_names(med)
        all_meshes = join(meshes, ", ")
        @warn("Available meshes: $all_meshes")
        error("Mesh $mesh_name not found.")
    end
    return med.data["FAS"][mesh_name]
end

""" Return node sets from med file.

Notes
-----
One node set id can have multiple names.

"""
function get_node_sets(med::MEDFile, mesh_name::String)::Dict{Int64, Vector{String}}
    mesh = get_mesh(med, mesh_name)
    node_sets = Dict{Int64, Vector{String}}(0 => ["OTHER"])
    if !haskey(mesh, "NOEUD")
        return node_sets
    end
    for (k, v) in mesh["NOEUD"]
        nset_id = parse(Int, split(k, "_")[2])
        node_sets[nset_id] = collect(to_ascii(d) for d in v["GRO"]["NOM"])
    end
    return node_sets
end

"""
    get_element_sets(med, mesh_name)

Return element sets from med file. Return type is a dictionary, where the key is
the element set id number (integer) and value is a vector of strings, containing
human-readable name for element set.

# Notes

One element set id can have multiple names.

"""
function get_element_sets(med::MEDFile, mesh_name::String)::Dict{Int64, Vector{String}}
    mesh = get_mesh(med, mesh_name)
    element_sets = Dict{Int64, Vector{String}}()
    if !haskey(mesh, "ELEME")
        return element_sets
    end
    elset_keys = sort(collect(keys(mesh["ELEME"])))
    for (i, k) in enumerate(elset_keys)
        v = mesh["ELEME"][k]
        if startswith(k, "FAM") # exported from Salome
            elset_id = parse(Int, split(k, '_')[2])
        else # exported from Gmsh
            elset_id = -i
        end
        if !isempty(v)
            element_sets[elset_id] = map(strip, collect(to_ascii(d) for d in v["GRO"]["NOM"]))
        else
            element_sets[elset_id] = [""]
        end
    end
    return element_sets
end

function get_nodes(med::MEDFile, nsets::Dict{Int, Vector{String}}, mesh_name::String)
    increments = keys(med.data["ENS_MAA"][mesh_name])
    @assert length(increments) == 1
    increment = first(increments)
    nodes = med.data["ENS_MAA"][mesh_name][increment]["NOE"]
    nset_ids = nodes["FAM"]
    nnodes = length(nset_ids)
    node_ids = get(nodes, "NUM", collect(1:nnodes))
    node_coords = nodes["COO"]
    dim = round(Int, length(node_coords)/nnodes)
    node_coords = reshape(node_coords, nnodes, dim)'
    d = Dict{Int64}{Tuple{Vector{String}, Vector{Float64}}}()
    for i=1:nnodes
        nset = nsets[nset_ids[i]]
        d[node_ids[i]] = (nset, node_coords[:, i])
    end
    return d
end

function get_connectivity(med::MEDFile, elsets::Dict{Int64, Vector{String}}, mesh_name::String)
    if !haskey(elsets, 0)
        elsets[0] = ["OTHER"]
    end
    increments = keys(med.data["ENS_MAA"][mesh_name])
    @assert length(increments) == 1
    increment = first(increments)
    all_elements = med.data["ENS_MAA"][mesh_name][increment]["MAI"]
    d = Dict{Int64, Tuple{Symbol, Vector{String}, Vector{Int64}}}()
    for eltype in keys(all_elements)
        elements = all_elements[eltype]
        elset_ids = elements["FAM"]
        nelements = length(elset_ids)
        element_ids = get(elements, "NUM", collect(1:nelements))
        element_connectivity = elements["NOD"]
        element_dim = round(Int, length(element_connectivity)/nelements)
        element_connectivity = reshape(element_connectivity, nelements, element_dim)'
        for i=1:nelements
            eltype = Symbol(eltype)
            elco = element_connectivity[:, i]
            elset = elsets[elset_ids[i]]
            d[element_ids[i]] = (eltype, elset, elco)
        end
    end
    return d
end

"""
    aster_read_mesh(filename, mesh_name=nothing)

Parse code aster .med file and return mesh data in a dictionary.

Dictionary contains additional dictionaries `nodes`, `node_sets`, `elements`,
`element_sets`, `element_types`, `surface_sets` and `surface_types`.

If mesh file contains several meshes, one must provide the mesh name as
additional argument or expcetion will be thrown.

"""
function aster_read_mesh(filename::String, mesh_name=nothing)
    isfile(filename) || error("Cannot read mesh file $filename: file not found.")
    aster_read_mesh_(MEDFile(filename), mesh_name)
end

function aster_read_mesh_(med::MEDFile, mesh_name=nothing)
    mesh_names = get_mesh_names(med)
    all_meshes = join(mesh_names, ", ")
    if mesh_name == nothing
        length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes")
        mesh_name = mesh_names[1]
    else
        mesh_name in mesh_names || error("Mesh $mesh_name not found from mesh file $fn. Available meshes: $all_meshes")
    end

    mesh = Dict{String, Dict}()
    mesh["nodes"] = Dict{Int, Vector{Float64}}()
    mesh["node_sets"] = Dict{String, Vector{Int}}()
    mesh["elements"] = Dict{Int, Vector{Int}}()
    mesh["element_types"] = Dict{Int, Symbol}()
    mesh["element_sets"] = Dict{String, Vector{Int}}()
    mesh["surface_sets"] = Dict{String, Vector{Tuple{Int, Symbol}}}()
    mesh["surface_types"] = Dict{String, Symbol}()

    elsets = get_element_sets(med, mesh_name)
    nsets = get_node_sets(med, mesh_name)
    for (nid, (nset_, coords)) in get_nodes(med, nsets, mesh_name)
        mesh["nodes"][nid] = coords
        for nset in nset_
            if !haskey(mesh["node_sets"], nset)
                mesh["node_sets"][nset] = []
            end
            push!(mesh["node_sets"][nset], nid)
        end
    end
    for (elid, (eltyp, elset_, elcon)) in get_connectivity(med, elsets, mesh_name)
        mesh["elements"][elid] = elcon
        mesh["element_types"][elid] = eltyp
        for elset in elset_
            if !haskey(mesh["element_sets"], elset)
                mesh["element_sets"][elset] = []
            end
            push!(mesh["element_sets"][elset], elid)
        end
    end
    return mesh
end


text:
# BBI - Shared Parts of bigWig and bigBed
# =======================================
#
# Detailed description of bigWig and bigBed file formats are found in the paper and its supplement:
# https://doi.org/10.1093/bioinformatics/btq351
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2922891/bin/supp_btq351_bbiSuppFINAL.doc
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module BBI

import BioCore: Mem
import Libz

include("header.jl")
include("summary.jl")
include("btree.jl")
include("rtree.jl")
include("section.jl")
include("zoom.jl")
include("compress.jl")

end # module


text:
using GLPK
using JuMP
using LinearAlgebra

# N-Queens
N = 8

model = Model(GLPK.Optimizer);

@variable(model, x[i=1:N, j=1:N], Bin)

# There must be exactly one queen in a given row/column
for i=1:N
    @constraint(model, sum(x[i, :]) == 1)
    @constraint(model, sum(x[:, i]) == 1)
end

# There can only be one queen on any given diagonal
for i in -(N-1):(N-1)
    @constraint(model, sum(diag(x,i)) <= 1)
    @constraint(model, sum(diag(reverse(x,dims=1), i)) <=1)
end

optimize!(model)

solution = convert.(Int,value.(x))


text:
const _GW_PS_CMD_OUT_STRIP_PATTERN = ['\t', ' ', '\n']
_get_proc_lstart(pid) = try
    pid = string(pid)
    pid == "-1" && return ""
    out = read(`ps -p $(pid) -o lstart=`, String)
    return string(strip(out, _GW_PS_CMD_OUT_STRIP_PATTERN))
    catch; return ""
end

_try_getpid(proc) = try
    return getpid(proc)
    catch; return -1
end

_force_kill(pid) = try
    pid = string(pid)
    pid == "-1" && return false
    run(`kill -SIGKILL $(pid)`; wait = true)
    return true
    catch; return false
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using ConeProgramDiff
using Documenter

makedocs(;
    modules=[ConeProgramDiff],
    authors="Chandler Squires, Theo Diamandis",
    repo="https://github.com/tjdiamandis/ConeProgramDiff.jl/blob/{commit}{path}#L{line}",
    sitename="ConeProgramDiff.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tjdiamandis.github.io/ConeProgramDiff.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/tjdiamandis/ConeProgramDiff.jl",
)


text:
module SignedDistanceFunction
include("./sdistance.jl")
using .Sdistance # exportされたmethodのみ使える
export signedDistance2D
end

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
using OrdinaryDiffEq, Test, RecursiveArrayTools, DiffEqDevTools, Statistics

u0 = fill(0., 2)
v0 = ones(2)
function f1_harmonic(dv,v,u,p,t)
  dv .= -u
end
function f2_harmonic(du,v,u,p,t)
  du .= v
end
function harmonic_analytic(y0, p, x)
  v0, u0 = y0
  ArrayPartition(-u0*sin(x) + v0*cos(x),u0*cos(x) + v0*sin(x))
end
ff_harmonic = DynamicalODEFunction(f1_harmonic,f2_harmonic;analytic=harmonic_analytic)
prob = DynamicalODEProblem(ff_harmonic,v0,u0,(0.0,5.0))

sol = solve(prob,SymplecticEuler(),dt=1/2)
sol_verlet = solve(prob,VelocityVerlet(),dt=1/100)
sol_ruth3 = solve(prob,Ruth3(),dt=1/100)

interp_time = 0:0.001:5
interp = sol(0.5)
interps = sol(interp_time)

sol_tsit5 = solve(prob,Tsit5())

prob = SecondOrderODEProblem(f1_harmonic,v0,u0,(0.0,5.0))

sol2 = solve(prob,SymplecticEuler(),dt=1/2)
sol2_verlet = solve(prob,VelocityVerlet(),dt=1/100)
sol2_ruth3 = solve(prob,Ruth3(),dt=1/100)

sol2_verlet(0.1)

@test sol[end][1] == sol2[end][1]
@test sol_verlet[end][1] == sol2_verlet[end][1]
@test sol_ruth3[end][1] == sol2_ruth3[end][1]
@test sol[end][3] == sol2[end][3]
@test sol_verlet[end][3] == sol2_verlet[end][3]
@test sol_ruth3[end][3] == sol2_ruth3[end][3]

prob = DynamicalODEProblem(ff_harmonic,v0,u0,(0.0,5.0))
println("Convergence tests")

dts = 1 .//2 .^(6:-1:3)
# Symplectic Euler
sim = test_convergence(dts,prob,SymplecticEuler(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 1 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 1 rtol = 1e-1
# Verlet
sim = test_convergence(dts,prob,VelocityVerlet(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
# Test that position converges faster for Verlet
position_error = :final => [mean(sim[i].u[2].x[1] - sim[i].u_analytic[2].x[1]) for i in 1:length(sim)]
@test first(DiffEqDevTools.calc𝒪estimates(position_error).second) ≈ 4.0 rtol=1e-1

# 2nd Order Tableaus
sim = test_convergence(dts,prob,VerletLeapfrog(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
sim = test_convergence(dts,prob,PseudoVerletLeapfrog(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
sim = test_convergence(dts,prob,McAte2(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1

# Ruth
sim = test_convergence(dts,prob,Ruth3(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
sim = test_convergence(dts,prob,McAte3(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
sim = test_convergence(dts,prob,CandyRoz4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,McAte4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,McAte42(),dense_errors=true)
@test_broken sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test_broken sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,CalvoSanz4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

dts = 1 .//2 .^(4:-1:0)
sim = test_convergence(dts,prob,McAte5(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,Yoshida6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4.69 rtol = 1e-1
sim = test_convergence(dts,prob,KahanLi6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,McAte8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,KahanLi8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sol = solve(prob,Nystrom4(),dt=1/1000)

# Nyström method
dts = 1 .//2 .^(9:-1:6)
sim = test_convergence(dts,prob,Nystrom4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,Nystrom4VelocityIndependent(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,IRKN3(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
sim = test_convergence(dts,prob,IRKN4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
dts = 1.0 ./2.0 .^(5:-1:0)
sim = test_convergence(dts,prob,Nystrom5VelocityIndependent(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

dts = 1.0 ./2.0 .^(2:-1:-2)
sim = test_convergence(dts,prob,SofSpa10(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 10 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

# Methods need BigFloat to test convergence rate
dts = big"1.0"./big"2.0".^(5:-1:1)
prob_big = DynamicalODEProblem(ff_harmonic,[big"1.0",big"1.0"],
                               [big"0.0", big"0.0"],(big"0.",big"70."))
sim = test_convergence(dts,prob_big,DPRKN6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 6 rtol = 1e-1
sim = test_convergence(dts,prob_big,DPRKN8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,DPRKN12(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 12 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,ERKN4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,ERKN5(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

# Adaptive methods regression test
sol = solve(prob, DPRKN6())
@test length(sol.u) < 20
sol = solve(prob, DPRKN8())
@test length(sol.u) < 13
sol = solve(prob, DPRKN12())
@test length(sol.u) < 10
sol = solve(prob, ERKN4(),reltol=1e-8)
@test length(sol.u) < 38
sol = solve(prob, ERKN5(),reltol=1e-8)
@test length(sol.u) < 34

# Test array partition outside of symplectic

function f_ap(du,u,p,t)
  du.x[1] .= -2u.x[2]
  du.x[2] .= u.x[1]
end

u = ArrayPartition((u0,v0))

prob = ODEProblem(f_ap,u,(0.0,5.0))
sol = solve(prob,Euler(),dt=1/100)

################# Out of place symplectic

println("Out of Place")

u0 = 0.0
v0 = 1.0
function f1_harmonic_nip(v,u,p,t)
  -u
end
function f2_harmonic_nip(v,u,p,t)
  v
end

ff_harmonic_nip = DynamicalODEFunction(f1_harmonic_nip,f2_harmonic_nip;
                                       analytic=harmonic_analytic)
prob = DynamicalODEProblem(ff_harmonic_nip,v0,u0,(0.0,5.0))

sol = solve(prob,SymplecticEuler(),dt=1/10)

dts = 1 .//2 .^(6:-1:3)
# Symplectic Euler
sim = test_convergence(dts,prob,SymplecticEuler(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 1 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 1 rtol = 1e-1
# Verlet
sim = test_convergence(dts,prob,VelocityVerlet(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
# Test that position converges faster for Verlet
position_error = :final => [mean(sim[i].u[2].x[1] - sim[i].u_analytic[2].x[1]) for i in 1:length(sim)]
@test first(DiffEqDevTools.calc𝒪estimates(position_error).second) ≈ 4.0 rtol=1e-1

# 2nd Order Tableaus
sim = test_convergence(dts,prob,VerletLeapfrog(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
sim = test_convergence(dts,prob,PseudoVerletLeapfrog(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1
sim = test_convergence(dts,prob,McAte2(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 2 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 2 rtol = 1e-1

# Ruth
sim = test_convergence(dts,prob,Ruth3(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
sim = test_convergence(dts,prob,McAte3(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
sim = test_convergence(dts,prob,CandyRoz4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,McAte4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,McAte42(),dense_errors=true)
@test_broken sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test_broken sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,CalvoSanz4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

dts = 1 .//2 .^(4:-1:0)
sim = test_convergence(dts,prob,McAte5(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,Yoshida6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4.69 rtol = 1e-1
sim = test_convergence(dts,prob,KahanLi6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,McAte8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

sim = test_convergence(dts,prob,KahanLi8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

# Nyström method
dts = 1 .//2 .^(9:-1:6)
sim = test_convergence(dts,prob,Nystrom4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob,Nystrom4VelocityIndependent(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
@test_broken sim = test_convergence(dts,prob,IRKN3(),dense_errors=true)
@test_broken sim.𝒪est[:l2] ≈ 3 rtol = 1e-1
@test_broken sim.𝒪est[:L2] ≈ 3 rtol = 1e-1
@test_broken sim = test_convergence(dts,prob,IRKN4(),dense_errors=true)
#@test_broken sim.𝒪est[:l2] ≈ 4 rtol = 1e-1
#@test_broken sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
dts = 1.0 ./2.0 .^(5:-1:0)
sim = test_convergence(dts,prob,Nystrom5VelocityIndependent(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

dts = 1.0 ./2.0 .^(2:-1:-2)
sim = test_convergence(dts,prob,SofSpa10(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 10 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

# Methods need BigFloat to test convergence rate
dts = big"1.0"./big"2.0".^(5:-1:1)
prob_big = DynamicalODEProblem(ff_harmonic_nip,big"1.0",big"0.0",
                                 (big"0.",big"70."))
sim = test_convergence(dts,prob_big,DPRKN6(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 6 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 6 rtol = 1e-1
sim = test_convergence(dts,prob_big,DPRKN8(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 8 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,DPRKN12(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 12 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,ERKN4(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1
sim = test_convergence(dts,prob_big,ERKN5(),dense_errors=true)
@test sim.𝒪est[:l2] ≈ 5 rtol = 1e-1
@test sim.𝒪est[:L2] ≈ 4 rtol = 1e-1

# Adaptive methods regression test
sol = solve(prob, DPRKN6())
@test length(sol.u) < 20
sol = solve(prob, DPRKN8())
@test length(sol.u) < 13
sol = solve(prob, DPRKN12())
@test length(sol.u) < 10
sol = solve(prob, ERKN4(),reltol=1e-8)
@test length(sol.u) < 38
sol = solve(prob, ERKN5(),reltol=1e-8)
@test length(sol.u) < 34


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

using Base.Test

include("../perfutil.jl")

## recursive fib ##

fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)

@test fib(20) == 6765
@timeit fib(20) "fib" "Recursive fibonacci"

## parse integer ##

function parseintperf(t)
    local n, m
    for i=1:t
        n = rand(UInt32)
        s = hex(n)
        m = UInt32(parse(Int64,s,16))
    end
    @test m == n
    return n
end

@timeit parseintperf(1000) "parse_int" "Integer parsing"

## array constructors ##

@test all(ones(200,200) .== 1)
# @timeit ones(200,200) "ones" "description"

## matmul and transpose ##

A = ones(200,200)
@test all(A*A' .== 200)
# @timeit A*A' "AtA" "description"

## mandelbrot set: complex arithmetic and comprehensions ##

function mandel(z)
    c = z
    maxiter = 80
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

mandelperf() = [ mandel(complex(r,i)) for i=-1.:.1:1., r=-2.0:.1:0.5 ]
@test sum(mandelperf()) == 14791
@timeit mandelperf() "mandel" "Calculation of mandelbrot set"

## numeric vector sort ##

function qsort!(a,lo,hi)
    i, j = lo, hi
    while i < hi
        pivot = a[(lo+hi)>>>1]
        while i <= j
            while a[i] < pivot; i += 1; end
            while a[j] > pivot; j -= 1; end
            if i <= j
                a[i], a[j] = a[j], a[i]
                i, j = i+1, j-1
            end
        end
        if lo < j; qsort!(a,lo,j); end
        lo, j = i, hi
    end
    return a
end

sortperf(n) = qsort!(rand(n), 1, n)
@test issorted(sortperf(5000))
@timeit sortperf(5000) "quicksort" "Sorting of random numbers using quicksort"

## slow pi series ##

function pisum()
    sum = 0.0
    for j = 1:500
        sum = 0.0
        for k = 1:10000
            sum += 1.0/(k*k)
        end
    end
    sum
end

@test abs(pisum()-1.644834071848065) < 1e-12
@timeit pisum() "pi_sum" "Summation of a power series"

## slow pi series, vectorized ##

function pisumvec()
    s = 0.0
    a = [1:10000]
    for j = 1:500
        s = sum(1./(a.^2))
    end
    s
end

#@test abs(pisumvec()-1.644834071848065) < 1e-12
#@timeit pisumvec() "pi_sum_vec"

## random matrix statistics ##

function randmatstat(t)
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i=1:t
        a = randn(n,n)
        b = randn(n,n)
        c = randn(n,n)
        d = randn(n,n)
        P = [a b c d]
        Q = [a b; c d]
        v[i] = trace((P.'*P)^4)
        w[i] = trace((Q.'*Q)^4)
    end
    return (std(v)/mean(v), std(w)/mean(w))
end

(s1, s2) = randmatstat(1000)
@test 0.5 < s1 < 1.0 && 0.5 < s2 < 1.0
@timeit randmatstat(1000) "rand_mat_stat" "Statistics on a random matrix"

## largish random number gen & matmul ##

@timeit rand(1000,1000)*rand(1000,1000) "rand_mat_mul" "Multiplication of random matrices"

## printfd ##

if is_unix()
    function printfd(n)
        open("/dev/null", "w") do io
            for i = 1:n
                @printf(io, "%d %d\n", i, i + 1)
            end
        end
    end

    printfd(1)
    @timeit printfd(100000) "printfd" "Printing to a file descriptor"
end

maxrss("micro")


text:
using Documenter, CollectiveSpins

pages = [
    "index.md",
    "installation.md",
    "geometry.md",
    "interaction.md",
    "effective_interaction.md",
    "system.md",
    "collective_modes.md",
    "descriptions.md",
    "api.md"
]

makedocs(
    modules = [CollectiveSpins],
    checkdocs = :exports,
    format = Documenter.HTML(
        edit_link = nothing,
        canonical = "https://github.io/qojulia/CollectiveSpins.jl/",
        assets = [ asset("assets/favicon.png", class=:ico, islocal = true) ]),
    sitename = "CollectiveSpins.jl",
    pages = pages
    )

deploydocs(
    repo = "github.com/qojulia/CollectiveSpins.jl.git",
    )


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataFrames
using CSV

function plants_report(solution)::DataFrame
    df = DataFrame()
    df."plant type" = String[]
    df."location name" = String[]
    df."year" = Int[]
    df."latitude (deg)" = Float64[]
    df."longitude (deg)" = Float64[]
    df."capacity (tonne)" = Float64[]
    df."amount processed (tonne)" = Float64[]
    df."amount received (tonne)" = Float64[]
    df."amount in storage (tonne)" = Float64[]
    df."utilization factor (%)" = Float64[]
    df."energy (GJ)" = Float64[]
    df."opening cost (\$)" = Float64[]
    df."expansion cost (\$)" = Float64[]
    df."fixed operating cost (\$)" = Float64[]
    df."variable operating cost (\$)" = Float64[]
    df."storage cost (\$)" = Float64[]
    df."total cost (\$)" = Float64[]
    T = length(solution["Energy"]["Plants (GJ)"])
    for (plant_name, plant_dict) in solution["Plants"]
        for (location_name, location_dict) in plant_dict
            for year = 1:T
                capacity = round(location_dict["Capacity (tonne)"][year], digits = 2)
                received = round(location_dict["Total input (tonne)"][year], digits = 2)
                processed = round(location_dict["Process (tonne)"][year], digits = 2)
                in_storage = round(location_dict["Storage (tonne)"][year], digits = 2)
                utilization_factor = round(processed / capacity * 100.0, digits = 2)
                energy = round(location_dict["Energy (GJ)"][year], digits = 2)
                latitude = round(location_dict["Latitude (deg)"], digits = 6)
                longitude = round(location_dict["Longitude (deg)"], digits = 6)
                opening_cost = round(location_dict["Opening cost (\$)"][year], digits = 2)
                expansion_cost =
                    round(location_dict["Expansion cost (\$)"][year], digits = 2)
                fixed_cost =
                    round(location_dict["Fixed operating cost (\$)"][year], digits = 2)
                var_cost =
                    round(location_dict["Variable operating cost (\$)"][year], digits = 2)
                storage_cost = round(location_dict["Storage cost (\$)"][year], digits = 2)
                total_cost = round(
                    opening_cost + expansion_cost + fixed_cost + var_cost + storage_cost,
                    digits = 2,
                )
                push!(
                    df,
                    [
                        plant_name,
                        location_name,
                        year,
                        latitude,
                        longitude,
                        capacity,
                        processed,
                        received,
                        in_storage,
                        utilization_factor,
                        energy,
                        opening_cost,
                        expansion_cost,
                        fixed_cost,
                        var_cost,
                        storage_cost,
                        total_cost,
                    ],
                )
            end
        end
    end
    return df
end

write_plants_report(solution, filename) = CSV.write(filename, plants_report(solution))

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
@testset "Torsion" begin
  # Trivial torsion module

  A = diagonal_matrix(fmpq[1, 1])
  T = Hecke.TorQuadMod(A)
  @test order(T) == 1

  D4_gram = matrix(ZZ, [[2, 0, 0, -1],
                        [0, 2, 0, -1],
                        [0, 0, 2, -1],
                        [-1 ,-1 ,-1 ,2]])

  L = Zlattice(gram = D4_gram)
  T = @inferred discriminant_group(L)
  @test order(T) == 4
  @test elementary_divisors(T) == fmpz[2, 2]

  Ld = dual(L)

  T = @inferred torsion_quadratic_module(Ld, L, snf = false)
  @test ngens(T) == 4
  G = gens(T)
  for i in 1:4
    @test T(lift(G[i])) == G[i]
  end

  TT, mTT = @inferred sub(T, [T([1, 1//2, 1//2, 1])])
  @test order(TT) == 2

  L = Zlattice(ZZ[1 0; 0 1])
  M = lattice(ambient_space(L), ZZ[2 0; 0 3])
  T = @inferred torsion_quadratic_module(L, M, gens = [[1, 1]])
  @test order(T) == 6

  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[0, 0]])
  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[1//2, 0]])
  @test_throws ArgumentError torsion_quadratic_module(L, M, gens = [[1, 1, 2]])
  @test_throws ArgumentError torsion_quadratic_module(L, lattice(ambient_space(L), QQ[1//2 0; 0 0]))

  #primary part of a TorQuadMod
  L = Zlattice(matrix(ZZ, [[2,0,0],[0,2,0],[0,0,2]]))
  T = Hecke.discriminant_group(L)
  @test basis_matrix(Hecke.cover(Hecke.primary_part(T,fmpz(2))[1])) == matrix(QQ, 3, 3, [1//2, 0, 0, 0, 1//2, 0, 0, 0, 1//2])

  #orthogonal submodule to a TorQuadMod
  L = Zlattice(matrix(ZZ, [[2,0,0],[0,2,0],[0,0,2]]))
  T = Hecke.discriminant_group(L)
  S = sub(T, gens(T))[1]
  @test basis_matrix(Hecke.cover(Hecke.orthogonal_submodule_to(T, S)[1])) == matrix(QQ, 3, 3, [1//2,0,0,0,1//2,0,0,0,1//2])

  #checks if a TorQuadMod is degenerate
  @test Hecke.isdegenerate(T) == false

  #test for rescaled torsion quadratic module
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(T, -1)) == matrix(QQ, 3, 3, [7//4,0,0,0,7//4,0,0,0,7//4])
  t = Hecke.TorQuadMod(QQ[1//3 0; 0 1//9])
  @test Hecke.gram_matrix_quadratic(Hecke.rescale(t, 2)) == matrix(QQ, 2, 2, [2//3,0,0,2//9])

  #test for normal form
  L1 = Zlattice(gram=matrix(ZZ, [[-2,0,0],[0,1,0],[0,0,4]]))
  T1 = Hecke.discriminant_group(L1)
  @test Hecke.gram_matrix_quadratic(Hecke.normal_form(T1)[1]) == matrix(QQ, 2, 2, [1//2,0,0,1//4])

  L1 = Zlattice(gram=ZZ[-4 0 0; 0 4 0; 0 0 -2])
  AL1 = discriminant_group(L1)
  L2 = Zlattice(gram=ZZ[-4 0 0; 0 -4 0; 0 0 2])
  AL2 = discriminant_group(L2)
  n1 = normal_form(AL1)[1]
  g1 = QQ[1//2   0   0;
            0 1//4   0;
            0   0 5//4]
  @test Hecke.gram_matrix_quadratic(n1) == g1
  n2 = normal_form(AL2)[1]
  g2 = QQ[1//2   0   0;
            0 1//4   0;
            0   0 5//4]
  @test Hecke.gram_matrix_quadratic(n2) == g2

  #test for brown invariant
  L = Zlattice(gram=matrix(ZZ, [[2,-1,0,0],[-1,2,-1,-1],[0,-1,2,0],[0,-1,0,2]]))
  T = discriminant_group(L)
  @test Hecke.brown_invariant(T) == 4

  #test for genus
  L = Zlattice(gram=diagonal_matrix(fmpz[1,2,3,4]))
  D = discriminant_group(L)
  @test genus(D, (4,0)) == genus(L)

  #test for isgenus
  L = Zlattice(gram=diagonal_matrix(fmpz[1,2,3,4]))
  D = discriminant_group(L)
  @test_throws ErrorException isgenus(D, (4,0))
  L1 = Zlattice(gram=matrix(ZZ, [[2, -1, 0, 0, 0, 0],[-1, 2, -1, -1, 0, 0],[0, -1, 2, 0, 0, 0],[0, -1, 0, 2, 0, 0],[0, 0, 0, 0, 6, 3],[0, 0, 0, 0, 3, 6]]))
  T1 = discriminant_group(L1)
  @test isgenus(T1, (6,0))
  G = genus(diagonal_matrix(fmpz[2, 6, 6]))
  D = discriminant_group(G)
  @test isgenus(D, (3,0))

end


text:
# A custom geometry type to test that we can render arbitrary primitives
# by decomposing them into simple meshes. This replaces the previous test
# which did the same thing using a Polyhedron from Polyhedra.jl. The Polyhedra
# test was removed because it required too many external dependencies just to
# verify a simple interface.
struct CustomGeometry <: GeometryPrimitive{3, Float64}
end

GeometryTypes.isdecomposable(::Type{<:Face}, ::CustomGeometry) = true
function GeometryTypes.decompose(::Type{F}, c::CustomGeometry) where {F <: Face}
    [convert(F, Face(1,2,3))]
end
GeometryTypes.isdecomposable(::Type{<:Point}, ::CustomGeometry) = true
function GeometryTypes.decompose(::Type{P}, c::CustomGeometry) where {P <: Point}
    convert(Vector{P}, [Point(0., 0, 0), Point(0., 1, 0), Point(0., 0, 1)])
end

@testset "Custom decomposable geometry" begin
    Visualizer(CustomGeometry())
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ad471560-28e9-11eb-16b5-f7025684bc5f
using SymEngine

# ╔═╡ 3f007f22-4148-11eb-1c8b-9b6e06ccc541
using OrdinaryDiffEq

# ╔═╡ 42931bc0-4148-11eb-26b9-e5fee42c7580
using Plots

# ╔═╡ 470b0f50-4148-11eb-3e5a-810c88e0c029
using LaTeXStrings

# ╔═╡ 30db9be0-5fda-11eb-3899-278be0c6f536
using PlutoUI

# ╔═╡ 80213b30-28e7-11eb-388d-5501de9711e7
md"""# Tema 6 Backstepping

## Ejemplo final de la teoría, estabilizar:

$$\dot x_1=x_1^2-x_1^3+x_2$$
$$\dot x_2=x_3$$
$$\dot x_3 =u$$

# Plan de acción
* Con $u$ controlo $x_3$
* Con $x_3$ controlo $x_2$
* Con $x_2$ controlo $x_1$

Pero el plan se ejecuta al revés, por eso se llama back-stepping ☺

"""

# ╔═╡ b5d53800-28ea-11eb-1858-4360a7d9fafc
md"De nuevo usamos la función de tema 4"

# ╔═╡ bd1a0730-28ea-11eb-2355-c1d9285733c5
begin
	function derivada(funcion,variables,derivadas)
		d=0
		for i=1:length(variables)
			d=d+diff(funcion,variables[i])*derivadas[i]
		end
		return d
	end

	function derivada(f,x,dx,n) #multiple distpach
		d=f
		for i=1:n
			d=derivada(d,x,dx)
		end
		return d
	end
end

# ╔═╡ 2890d210-28ea-11eb-2f16-23af078db24a
x1,x2,x3,u=symbols("x1,x2,x3,u");

# ╔═╡ e0c749a0-28e9-11eb-2a9d-3b5f6566605a
begin
	dx1=x1^2-x1^3+x2
    dx2=x3
	dx3=u
end;

# ╔═╡ 4cce2a0e-28ea-11eb-26c6-672aeb7bfdf9
md"## Paso 1, controlar $x_1$"

# ╔═╡ 646d6000-28ea-11eb-342d-b1c902896547
V1=x1^2/2

# ╔═╡ 224b2df0-28eb-11eb-2a45-c9eb740f149c
dV1=expand(derivada(V1,[x1],[dx1]))

# ╔═╡ 2d6d8fe0-c925-11eb-3a32-2fd273a45a10
md"Vamos a *tantear* un poco para entender cómo funciona esto, para ello vamos a jugar con $x_2$ a ver si conseguimos que V disminuya como nosotros queremos"

# ╔═╡ 4684dc70-28eb-11eb-1fde-2575c48354e8
md"""De forma directa $x_2$ como control poríamos eliminar el término **malo** y dejar uno cuadrático bueno bastaría con hacer que $x_2x_1 + x_1^3=-k_1x_1^2 \to x_2=-k_1x_1 - x_1^2$"""

# ╔═╡ d589c21e-5fd5-11eb-3c70-535e951c59cf
k1=symbols("k1");

# ╔═╡ 2c671c60-c925-11eb-303d-b398113bdb9b
expand(subs(dV1,x2=>-k1*x1-x1^2))

# ╔═╡ 01b23ad0-c926-11eb-2dc3-f143c1d2b6f4
md"Vamos a guardar lo que hemos encontrado en una funciónde control $\phi_1$"

# ╔═╡ 5159e1d0-28ec-11eb-14f4-99de193e35a1
ϕ1=-k1*x1-x1^2  #el  valor deseado de x2

# ╔═╡ bbd3f930-4181-11eb-1007-7fecf56f39aa
md"""Veamos que no nos hemos equivocado"""

# ╔═╡ 220f19e0-28ec-11eb-057e-0b369e37930c
expand(subs(dV1,x2=>ϕ1)) #si pudiera hacer esto realmente ya habría terminado

# ╔═╡ 883fd650-28ec-11eb-0ba3-014589882c54
md"Esto estaría guay, pero no puedo hacer que $x_2$ valga lo que yo quiera (todavía) así que defino un *error* $z_1=x_2-\phi_1$ y continuo, lo que tengo en realidad es que $x_2=\phi_1+z_1$"

# ╔═╡ c5680400-418b-11eb-32bb-ad1a2f2077ed
z1=symbols("z1");

# ╔═╡ 95969500-418d-11eb-0132-fd7f0ed24652
md"De modo que al derivar lo que realmente sale es"

# ╔═╡ 5f805e82-418b-11eb-05a4-dba1827925ec
expand(subs(dV1,x2=>ϕ1+z1))

# ╔═╡ 5eadc9c0-c926-11eb-2e4b-cb01e83f5e4f
md"""Resulta que el término $x_1z_1$ es **malo**. Para *cargármelo* necesito añadir a $V$ algo que al derivar produzca $z_1\cdot$"algo que pueda contolar" ... eso nos lleva al paso 2"""

# ╔═╡ 42851430-28ed-11eb-2823-a5014cb3eda0
md"""## Paso 2, controlar $x_1$ y $z_1$
Hemos reducido el problema de controlar $x_1$ a hacer $z_1=0$
¿Qué hemos ganado?
Paciencia, la idea es que nos estamos acercando a $u$
"¿Cómo me *desago* del término con $z_1$?
## Ampliemos nuestro V con $\frac{z_1^2}{2}$"""

# ╔═╡ 984ebfb0-28ed-11eb-2cee-fdbf3d493de8
V2=V1+z1^2/2

# ╔═╡ a5582280-28f4-11eb-2c05-51acb08f478b
md"Y derivemos a ver que sale, para ello primero calculo $\dot z=\dot x_2 - \dot \phi_1$"

# ╔═╡ 559eb340-28f3-11eb-3720-059fec2196ee
begin
	dϕ1=symbols("dϕ1")
	dz1=dx2-dϕ1
end;

# ╔═╡ 03ef90e0-4190-11eb-1f86-9f2bf2582d96
md"**Truco**: para no arrastrar *chorizos* creamos una variale que sea la derivada como haríamos **a mano**, luego ya calcularemos la verdadera derivada al final"

# ╔═╡ 7d9e7930-28ed-11eb-20fd-970954cf18b7
let
	dV2=derivada(V2,[x1 z1],[dx1 dz1])
	dV2=subs(dV2,x2=>z1+ϕ1) #La sustiucion es para hacer aparecer z1
	dV2=expand(dV2)
end

# ╔═╡ fb177d00-5fd7-11eb-0636-2f7f9dbd5598
k2=symbols("k2");

# ╔═╡ 2e63aed0-28f7-11eb-27a6-cbf8cb83c09c
ϕ2=-x1 -k2*z1 +dϕ1

# ╔═╡ da519e60-28f6-11eb-3598-c12f14f24180
md"""queremos hacer que $x_1z_1 + x_3z_1 - z_1 \dot ϕ_1 = -k_2z_1^2$, despejamos

$x_1 + x_3 - \dot ϕ_1 = -k_2z_1$

$x_3 = -k_2z_1 -x_1 + \dot ϕ_1$

Y hacemos $ϕ_2=-k_2z_1 -x_1 + \dot ϕ_1$

**Nota:** esto se podría despejar usando un paquete CAS más potente (Symbolics por ejemplo) pero no merece la pena, hacerlo a mano ayuda a *entender* lo que estamos haciendo. El CAS lo usaremos solamente para verificar los cálculos (que per se ya es muy útil)
"""

# ╔═╡ 0dc9ed20-4192-11eb-1039-93d7bfd5d09f
md"veamos que funciona"

# ╔═╡ ebea590e-4182-11eb-2a0b-a9161adc7f97
let
	dV2=derivada(V2,[x1 z1],[dx1 dz1])
	dV2=subs(dV2,x2=>z1+ϕ1)
	dV2=subs(dV2,x3=>ϕ2) #lo nuevo (si pudiese hacerlo claro)
	dV2=expand(dV2)
end

# ╔═╡ a58593c0-28f7-11eb-1674-51101547d596
md"""de nuevo lo mismo **no puedo** hacer eso así que
ahora:
$z_2=x_3-\phi_2$"""

# ╔═╡ 3094ca2e-28f8-11eb-2cc0-a1d2cebc1da6
z2=symbols("z2")

# ╔═╡ 3dbf9b10-4192-11eb-04e5-1f29ccd9da2a
let
	dV2=derivada(V2,[x1 z1],[dx1 dz1])
	dV2=subs(dV2,x2=>z1+ϕ1)
	dV2=subs(dV2,x3=>z2+ϕ2) #esto es lo que SI puedo hacer
	dV2=expand(dV2)
end

# ╔═╡ 8a020db0-c927-11eb-09b9-45afce72c1c8
md"""Al igual que antes queremos *cargarnos* $z_2z_1$ para ello hecesitamos un término adicional que de $z_2\cdot$"cosas" y a estas alturas os imaginaréis ya lo que vamos a poner ;)"""

# ╔═╡ fb7102a0-28f8-11eb-131d-c1792c1aaa27
md"""
## Paso 3, controlar $x_1$, $z_1$ y $z_2$
Hemos llegado a dónde queríamos por fin va a salir $u$
## Ampliemos nuestro V otra vez"""

# ╔═╡ 28437790-28f9-11eb-18f1-11364b9029ee
V3=V2+z2^2/2

# ╔═╡ 9f114912-28f9-11eb-17f1-7fa461ff3371
begin
	dϕ2=symbols("dϕ2")
	dz2=dx3-dϕ2
end;

# ╔═╡ 24332a00-28fa-11eb-2a88-8b5e6f124e6d
md"Aparece $u$ ¡Yupi!"

# ╔═╡ 85a6c900-28f9-11eb-229f-f78f9cd042e4
let
	dV3=derivada(V3,[x1 z1 z2],[dx1 dz1 dz2])
	dV3=subs(dV3,x2=>z1+ϕ1)
	dV3=subs(dV3,x3=>z2+ϕ2)
	dV3=expand(dV3)
end

# ╔═╡ c9b94770-28fa-11eb-25ba-2384e1b7672a
md"""queremos hacer que $uz_2 - z_2 \dot ϕ_2 + z_2z_1 = -k_3z_2^2$

$u - \dot ϕ_2 + z_1 = -k_3z_2$

$u =  \dot ϕ_2 - z_1 -k_3z_2$

Ahora vamos a ver que funciona:"""

# ╔═╡ db3aa380-5fd8-11eb-16fc-6b7ef56e8daa
k3=symbols("k3");

# ╔═╡ f8482150-4183-11eb-05f2-599741c800de
let
	dV3=derivada(V3,[x1 z1 z2],[dx1 dz1 dz2])
	dV3=subs(dV3,x2=>z1+ϕ1)
	dV3=subs(dV3,x3=>z2+ϕ2)
	dV3=subs(dV3,u=>dϕ2-z1-k3*z2)
	dV3=expand(dV3)
end

# ╔═╡ f16d5e90-4183-11eb-16dc-8191ac698ce2
md"""

## Ya tenemos ley de control, ahora hay que implementarla

"""

# ╔═╡ d4991de0-4147-11eb-067b-4f9c143fb673
function control(x,p,t)
	#completemos esto interactivamente
    x1,x2,x3=x
	k1,k2,k3=p

	ϕ1=-k1*x1 - x1^2

	dϕ1=(-k1 - 2*x1)*(x2 + x1^2 - x1^3)

	z1=x2-ϕ1

	ϕ2= dϕ1 - x1 - k2*z1

	z2=x3-ϕ2

	dϕ2= -k2*(-dϕ1 + x3) + (-k1 - 2*x1)*x3 + (x2 + x1^2 - x1^3)*((-k1 - 2*x1)*(2*x1 - 3*x1^2) - 2*(x2 + x1^2 - x1^3)) - (x2 + x1^2 - x1^3)
	#si se hace a mano se puede sustituir todavía menos por ejemplo podemos darnos cuenta de que el término repetido x2 + x1^2 - x1^3 es dx2 y calcularlo una sola vez...

	u=dϕ2-z1-k3*z2

	#bonus para ver que funciona
	V=(x1^2 + z1^2 + z2^2)/2

	return (u,V)
end

# ╔═╡ c259a7d0-28fb-11eb-3cc6-ab3e44ea196c
function derivadas(x,p,t)
	x1,x2,x3=x
	u=control(x,p,t)[1]
	dx1=x1^2-x1^3+x2
    dx2=x3
	dx3=u
	dx=[dx1 dx2 dx3]
end

# ╔═╡ 4704bc10-c963-11eb-24d9-2dd79b9b5f33
#=prueba para ver que devuele números y por lo tanto hemos definido todo
si hemos dejado algo sin definir apareceran las variables simbólicas=#
control([1 2 3],[1 1 1],0)

# ╔═╡ f78bff70-4192-11eb-2f14-05fbee335f4f
md"""ha llegado el momento de hacer los cálculos"""

# ╔═╡ 19ce9560-c939-11eb-0678-6956501cddf9
ϕ1

# ╔═╡ 11015170-c939-11eb-3145-1ff24e53a7f8
dϕ1_b=derivada(ϕ1,[x1],[dx1])

# ╔═╡ b941ace2-c939-11eb-1414-a3befce0a650
ddϕ1=derivada(dϕ1_b,[x1 x2],[dx1 dx2])

# ╔═╡ 83fcb3a0-c938-11eb-1051-a93ef03e0639
ϕ2

# ╔═╡ 70c05310-c96e-11eb-26c8-334d2b13e99b

# ╔═╡ e2c93200-c938-11eb-2d02-752060d9d7d9
dϕ2_b=derivada(ϕ2,[x1 z1 dϕ1],[dx1 dz1 ddϕ1])

# ╔═╡ 4031d6c0-4196-11eb-275f-cda9a32cce36
md"## Vamos a ver cómo queda el *chorizo*"

# ╔═╡ 6ff23950-4195-11eb-250b-898515a5c031
chorizo_u=expand(control([x1 x2 x3],[k1,k2,k3],0)[1])

# ╔═╡ 9e20ff40-c963-11eb-16e2-cd955911eb41
md"## Bonus, depuración:
Para comprobar los cálculos los haremos de otra forma y veremos si concuerdan, para ello trabajaremos directamente con las x, primero ponemos $V$ en términos de $x_1$, $x_2$ y $x_3$"

# ╔═╡ ae7ec160-c963-11eb-07ae-2b1bb1cd78d4
begin
	V_x=V3
	V_x=subs(V_x,z1=>x2-ϕ1)
	V_x=subs(V_x,z2=>x3-ϕ2)
	V_x=subs(V_x,dϕ1=>dϕ1_b)
	V_x=subs(V_x,z1=>x2-ϕ1)
	#El último parece redundante pero hay que insistir hasta que no queden zs
	#Algunas sustituciones hacen aparecer z1 de nuevo
end

# ╔═╡ 242cce30-c968-11eb-2c4b-4f2162e95f54
md"ahora derivamos y substituimos el conrtol quedando todo en términos de las x"

# ╔═╡ ae663050-c963-11eb-24d7-d3dcd3055b4b
begin
	dV_x=derivada(V_x,[x1 x2 x3],[dx1 dx2 dx3])
	dV_x=subs(dV_x,u=>chorizo_u)
end

# ╔═╡ 4ad29ab2-c968-11eb-26c6-f90433e890b1
md"la idea es comparar con lo que esperamos que salga"

# ╔═╡ ae4fe930-c963-11eb-30eb-f16f9e561190
begin
	dV_teórica=-k1*x1^2 - k2*z1^2 - k3*z2^2 - x1^4  #esto dice la teoría
	#para comparar lo ponemos en términos de las x solamente
	dV_teórica=subs(dV_teórica,z1=>x2-ϕ1)
	dV_teórica=subs(dV_teórica,z2=>x3-ϕ2)
	dV_teórica=subs(dV_teórica,dϕ1=>dϕ1_b)
	dV_teórica=subs(dV_teórica,z1=>x2-ϕ1)  #es necesario ya que vuelve a aparecer z1
end

# ╔═╡ cc2f0f30-c968-11eb-0814-bbe5107bdd3b
md"Si sale 0 es que lo hemos hecho bien"

# ╔═╡ ee6c43b0-c963-11eb-13f9-c51ef0a7f9a2
expand(dV_teórica-dV_x)

# ╔═╡ 9e7a974e-c96b-11eb-3a2b-759bf3f51058
md"# Y ahora a simular"

# ╔═╡ 49b79790-5fda-11eb-2917-4b055386f9db
k_1= @bind k_1 Slider(0:0.01:5,default=1,show_value=true)

# ╔═╡ 8ac8b250-5fda-11eb-2acd-4f317055f0cf
k_2= @bind k_2 Slider(0:0.01:5,default=1,show_value=true)

# ╔═╡ 8fbb0600-5fda-11eb-2f43-196140b76f69
k_3= @bind k_3 Slider(0:0.01:5,default=1,show_value=true)

# ╔═╡ fc5e8b30-4147-11eb-30d3-4bef0bac9ca6
begin
	x0 = [1 0 0]
	parametros=[k_1 k_2 k_3]
	tspan = (0.0,10.0)
	prob = ODEProblem(derivadas,x0,tspan,parametros);
	sol = solve(prob,Tsit5());

	Xs=sol.u
	ts=sol.t

	Vs=[]
	Us=[]

	for i=1:length(ts)
		u,V=control(Xs[i],parametros,ts[i])
		push!(Vs,V)
		push!(Us,u)
	end

	# Salida
	fig1=plot(sol,vars=(0,1), xaxis=L"t",yaxis=L"y(t)",label=L"y(t)")

	#Estados
	fig2=plot(sol, xaxis="t",yaxis=L"x(t)", label=:none)

	#Señal de control
	fig3=plot(ts,Us,xaxis="t",yaxis=L"u(t)", label=:none)
	xlims!(tspan)

	#V
	fig4=plot(ts,Vs,xaxis="t",yaxis=L"V(t)", label=:none)
	xlims!(tspan)

	l = @layout [a ; b ; c; d]
	plot(fig1,fig2,fig3,fig4, layout=l)
end

# ╔═╡ Cell order:
# ╟─80213b30-28e7-11eb-388d-5501de9711e7
# ╠═ad471560-28e9-11eb-16b5-f7025684bc5f
# ╟─b5d53800-28ea-11eb-1858-4360a7d9fafc
# ╟─bd1a0730-28ea-11eb-2355-c1d9285733c5
# ╠═2890d210-28ea-11eb-2f16-23af078db24a
# ╠═e0c749a0-28e9-11eb-2a9d-3b5f6566605a
# ╟─4cce2a0e-28ea-11eb-26c6-672aeb7bfdf9
# ╠═646d6000-28ea-11eb-342d-b1c902896547
# ╠═224b2df0-28eb-11eb-2a45-c9eb740f149c
# ╟─2d6d8fe0-c925-11eb-3a32-2fd273a45a10
# ╠═2c671c60-c925-11eb-303d-b398113bdb9b
# ╟─4684dc70-28eb-11eb-1fde-2575c48354e8
# ╠═d589c21e-5fd5-11eb-3c70-535e951c59cf
# ╠═01b23ad0-c926-11eb-2dc3-f143c1d2b6f4
# ╠═5159e1d0-28ec-11eb-14f4-99de193e35a1
# ╟─bbd3f930-4181-11eb-1007-7fecf56f39aa
# ╠═220f19e0-28ec-11eb-057e-0b369e37930c
# ╟─883fd650-28ec-11eb-0ba3-014589882c54
# ╠═c5680400-418b-11eb-32bb-ad1a2f2077ed
# ╟─95969500-418d-11eb-0132-fd7f0ed24652
# ╠═5f805e82-418b-11eb-05a4-dba1827925ec
# ╟─5eadc9c0-c926-11eb-2e4b-cb01e83f5e4f
# ╟─42851430-28ed-11eb-2823-a5014cb3eda0
# ╠═984ebfb0-28ed-11eb-2cee-fdbf3d493de8
# ╟─a5582280-28f4-11eb-2c05-51acb08f478b
# ╠═559eb340-28f3-11eb-3720-059fec2196ee
# ╟─03ef90e0-4190-11eb-1f86-9f2bf2582d96
# ╠═7d9e7930-28ed-11eb-20fd-970954cf18b7
# ╠═fb177d00-5fd7-11eb-0636-2f7f9dbd5598
# ╠═2e63aed0-28f7-11eb-27a6-cbf8cb83c09c
# ╟─da519e60-28f6-11eb-3598-c12f14f24180
# ╟─0dc9ed20-4192-11eb-1039-93d7bfd5d09f
# ╠═ebea590e-4182-11eb-2a0b-a9161adc7f97
# ╟─a58593c0-28f7-11eb-1674-51101547d596
# ╠═3094ca2e-28f8-11eb-2cc0-a1d2cebc1da6
# ╠═3dbf9b10-4192-11eb-04e5-1f29ccd9da2a
# ╟─8a020db0-c927-11eb-09b9-45afce72c1c8
# ╟─fb7102a0-28f8-11eb-131d-c1792c1aaa27
# ╠═28437790-28f9-11eb-18f1-11364b9029ee
# ╠═9f114912-28f9-11eb-17f1-7fa461ff3371
# ╟─24332a00-28fa-11eb-2a88-8b5e6f124e6d
# ╠═85a6c900-28f9-11eb-229f-f78f9cd042e4
# ╟─c9b94770-28fa-11eb-25ba-2384e1b7672a
# ╠═db3aa380-5fd8-11eb-16fc-6b7ef56e8daa
# ╠═f8482150-4183-11eb-05f2-599741c800de
# ╟─f16d5e90-4183-11eb-16dc-8191ac698ce2
# ╠═c259a7d0-28fb-11eb-3cc6-ab3e44ea196c
# ╠═d4991de0-4147-11eb-067b-4f9c143fb673
# ╠═4704bc10-c963-11eb-24d9-2dd79b9b5f33
# ╟─f78bff70-4192-11eb-2f14-05fbee335f4f
# ╠═19ce9560-c939-11eb-0678-6956501cddf9
# ╠═11015170-c939-11eb-3145-1ff24e53a7f8
# ╠═b941ace2-c939-11eb-1414-a3befce0a650
# ╠═83fcb3a0-c938-11eb-1051-a93ef03e0639
# ╠═70c05310-c96e-11eb-26c8-334d2b13e99b
# ╠═e2c93200-c938-11eb-2d02-752060d9d7d9
# ╟─4031d6c0-4196-11eb-275f-cda9a32cce36
# ╠═6ff23950-4195-11eb-250b-898515a5c031
# ╟─9e20ff40-c963-11eb-16e2-cd955911eb41
# ╠═ae7ec160-c963-11eb-07ae-2b1bb1cd78d4
# ╟─242cce30-c968-11eb-2c4b-4f2162e95f54
# ╠═ae663050-c963-11eb-24d7-d3dcd3055b4b
# ╟─4ad29ab2-c968-11eb-26c6-f90433e890b1
# ╠═ae4fe930-c963-11eb-30eb-f16f9e561190
# ╟─cc2f0f30-c968-11eb-0814-bbe5107bdd3b
# ╠═ee6c43b0-c963-11eb-13f9-c51ef0a7f9a2
# ╟─9e7a974e-c96b-11eb-3a2b-759bf3f51058
# ╠═3f007f22-4148-11eb-1c8b-9b6e06ccc541
# ╠═42931bc0-4148-11eb-26b9-e5fee42c7580
# ╠═470b0f50-4148-11eb-3e5a-810c88e0c029
# ╠═30db9be0-5fda-11eb-3899-278be0c6f536
# ╟─49b79790-5fda-11eb-2917-4b055386f9db
# ╟─8ac8b250-5fda-11eb-2acd-4f317055f0cf
# ╟─8fbb0600-5fda-11eb-2f43-196140b76f69
# ╟─fc5e8b30-4147-11eb-30d3-4bef0bac9ca6

max_stars_repo_path: Tema 6.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 6.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 6.jl

max_forks_repo_name: Dictino/CNL


text:
module TestIOHelpers
using Test
using PhaseSpaceIO: setbit, getbit

@testset "getbit, setbit" begin
    for _ in 1:1000
        T = UInt32
        x = rand(T)
        i,j = rand(1:sizeof(T), 2)
        val = rand(Bool)
        xi = setbit(x,val, i)
        @test getbit(xi, i) == val
        if i == j
            @test getbit(xi, j) == val
        else
            @test getbit(xi, j) == getbit(x, j)
        end
    end
end

@testset "getbit" begin
    T = UInt32
    @test getbit(T(1) << 23, 23) == true
    @test getbit(T(1) << 23, 2) == false

    for i in 1:sizeof(T)
        x = one(T) << i
        @test getbit(x, i) == true
        for j in 1:sizeof(T)
            i == j && continue
            @test getbit(x, j) == false
        end
    end

    for _ in 1:100
        x = rand(T)
        y = rand(T)
        i = rand(1:sizeof(T))
        @test getbit(~x, i) === ~getbit(x, i)
        @test getbit(x,i) & getbit(y, i) === getbit(x&y, i)
        @test getbit(x,i) | getbit(y, i) === getbit(x|y, i)
    end
end

end #module


text:
mutable struct Simulation
    status::BUILD_STATUS
    problem::Union{Nothing, DiffEqBase.DAEProblem}
    perturbations::Vector{<:Perturbation}
    x0_init::Vector{Float64}
    initialized::Bool
    tstops::Vector{Float64}
    callbacks::DiffEqBase.CallbackSet
    solution::Union{Nothing, DiffEqBase.DAESolution}
    simulation_folder::String
    simulation_inputs::SimulationInputs
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    multimachine::Bool
end

function Simulation(;
    simulation_inputs,
    perturbations = Vector{Perturbation}(),
    simulation_folder::String = "",
    console_level = Logging.Warn,
    file_level = Logging.Debug,
)
    return Simulation(
        BUILD_INCOMPLETE,
        nothing,
        perturbations,
        Vector{Float64}(),
        false,
        Vector{Float64}(),
        DiffEqBase.CallbackSet(),
        nothing,
        simulation_folder,
        simulation_inputs,
        console_level,
        file_level,
        false,
    )
end

function Simulation!(
    simulation_folder::String,
    system::PSY.System,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
)
    return Simulation!(simulation_folder, system, tspan, [perturbation]; kwargs...)
end

function Simulation(
    simulation_folder::String,
    system::PSY.System,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
)
    return Simulation(simulation_folder, system, tspan, [perturbation]; kwargs...)
end

"""
Builds the simulation object and conducts the indexing process. The initial conditions are stored in the system.

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `system_to_file::Bool`: Serializes the initialized system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation!(
    simulation_folder::String,
    system::PSY.System,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
)
    check_folder(simulation_folder)
    # Instantiates the Simulation object
    sim = Simulation(
        simulation_inputs = SimulationInputs(sys = system, tspan = tspan),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "initialized_system.json"))
    end
    return sim
end

"""
Initializes the simulations and builds the indexing. The input system is not modified during the initialization

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `system_to_file::Bool`: Serializes the original input system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debugg
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation(
    simulation_folder::String,
    system::PSY.System,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
)
    check_folder(simulation_folder)
    simulation_system = deepcopy(system)
    # Instantiates the Simulation object
    sim = Simulation(
        simulation_inputs = SimulationInputs(sys = simulation_system, tspan = tspan),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "input_system.json"))
    end
    return sim
end

function reset!(sim::Simulation)
    @info "Rebuilding the simulation after reset"
    build!(sim; file_mode = "a")
    @info "Simulation reset to status $(sim.status)"
    return
end

function configure_logging(sim::Simulation, file_mode)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = sim.console_level,
        file = true,
        filename = joinpath(sim.simulation_folder, SIMULATION_LOG_FILENAME),
        file_level = sim.file_level,
        file_mode = file_mode,
        tracker = nothing,
        set_global = false,
    )
end

function build!(sim; file_mode = "w", kwargs...)
    logger = configure_logging(sim, file_mode)
    try
        Logging.with_logger(logger) do
            _build!(sim; kwargs...)
        end
    finally
        close(logger)
    end
    return
end

function _build!(sim::Simulation; kwargs...)
    simulation_system = get_system(sim.simulation_inputs)
    sim.status = BUILD_INCOMPLETE
    PSY.set_units_base_system!(simulation_system, "DEVICE_BASE")
    check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
    initialized = false
    build!(sim.simulation_inputs)
    simulation_inputs = sim.simulation_inputs
    @debug "Simulation Inputs Created"
    var_count = get_variable_count(simulation_inputs)

    flat_start = zeros(var_count)
    bus_count = length(PSY.get_components(PSY.Bus, simulation_system))
    flat_start[1:bus_count] .= 1.0
    sim.x0_init = get(kwargs, :initial_guess, flat_start)

    initialize_simulation = get(kwargs, :initialize_simulation, true)
    if initialize_simulation
        @info("Initializing Simulation States")
        _add_aux_arrays!(simulation_inputs, Float64)
        sim.initialized = calculate_initial_conditions!(sim, simulation_inputs)
    end

    dx0 = zeros(var_count)
    _build_perturbations!(sim)
    _add_aux_arrays!(simulation_inputs, Float64)
    sim.problem = DiffEqBase.DAEProblem(
        system!,
        dx0,
        sim.x0_init,
        get_tspan(sim.simulation_inputs),
        simulation_inputs,
        differential_vars = get_DAE_vector(simulation_inputs);
        kwargs...,
    )
    sim.multimachine = (get_global_vars(simulation_inputs)[:ω_sys_index] != 0)
    sim.status = BUILT
    @info "Completed Build Successfully. Simulations status = $(sim.status)"
    return nothing
end

function _add_aux_arrays!(inputs::SimulationInputs, ::Type{T}) where {T <: Number}
    @debug "Auxiliary Arrays created with Type $(T)"
    bus_count = get_bus_count(inputs)
    get_aux_arrays(inputs)[1] = collect(zeros(T, bus_count))                       #I_injections_r
    get_aux_arrays(inputs)[2] = collect(zeros(T, bus_count))                       #I_injections_i
    get_aux_arrays(inputs)[3] = collect(zeros(T, get_n_injection_states(inputs)))  #injection_ode
    get_aux_arrays(inputs)[4] = collect(zeros(T, get_n_branches_states(inputs)))   #branches_ode
    get_aux_arrays(inputs)[5] = collect(zeros(Complex{T}, bus_count))              #I_bus
    get_aux_arrays(inputs)[6] = collect(zeros(T, 2 * bus_count))                   #I_balance
    return
end

function _build_perturbations!(sim::Simulation)
    @info "Attaching Perturbations"
    if isempty(sim.perturbations)
        @debug "The simulation has no perturbations"
        return DiffEqBase.CallbackSet(), [0.0]
    end
    system = get_system(sim.simulation_inputs)
    perturbations = sim.perturbations
    perturbations_count = length(perturbations)
    callback_vector = Vector{DiffEqBase.DiscreteCallback}(undef, perturbations_count)
    tstops = Vector{Float64}(undef, perturbations_count)
    for (ix, pert) in enumerate(perturbations)
        @debug pert
        condition = (x, t, integrator) -> t in [pert.time]
        affect = get_affect(system, pert)
        callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
        tstops[ix] = pert.time
    end
    sim.tstops = tstops
    sim.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
    return
end

function _index_local_states!(
    component_state_index::Vector{Int64},
    local_states::Vector{Symbol},
    component::PSY.DynamicComponent,
)
    for (ix, s) in enumerate(PSY.get_states(component))
        component_state_index[ix] = findfirst(x -> x == s, local_states)
    end
    return
end

function _attach_ports!(component::PSY.DynamicComponent)
    PSY.get_ext(component)[PORTS] = Ports(component)
    return
end

function _attach_inner_vars!(
    dynamic_device::PSY.DynamicGenerator,
    ::Type{T} = Real,
) where {T <: Real}
    dynamic_device.ext[INNER_VARS] = zeros(T, 9)
    return
end

function _attach_inner_vars!(
    dynamic_device::PSY.DynamicInverter,
    ::Type{T} = Real,
) where {T <: Real}
    dynamic_device.ext[INNER_VARS] = zeros(T, 14)
    return
end

function _attach_control_refs!(device::PSY.StaticInjection)
    dynamic_device = PSY.get_dynamic_injector(device)
    dynamic_device.ext[CONTROL_REFS] = [
        PSY.get_V_ref(dynamic_device),
        PSY.get_ω_ref(dynamic_device),
        PSY.get_P_ref(dynamic_device),
        PSY.get_reactive_power(device),
    ]
    return
end

function _index_port_mapping!(
    index_component_inputs::Vector{Int64},
    local_states::Vector{Symbol},
    component::PSY.DynamicComponent,
)
    _attach_ports!(component)
    for i in PSY.get_ext(component)[PORTS][:states]
        tmp = [(ix, var) for (ix, var) in enumerate(local_states) if var == i]
        isempty(tmp) && continue
        push!(index_component_inputs, tmp[1][1])
    end

    return
end

function _get_Ybus(sys::PSY.System)
    n_buses = length(PSY.get_components(PSY.Bus, sys))
    dyn_lines = PSY.get_components(PSY.DynamicBranch, sys)
    if !isempty(PSY.get_components(PSY.ACBranch, sys))
        Ybus_ = PSY.Ybus(sys)
        Ybus = Ybus_[:, :]
        lookup = Ybus_.lookup[1]
        for br in dyn_lines
            ybus!(Ybus, br, lookup, -1.0)
        end
    else
        Ybus = SparseMatrixCSC{Complex{Float64}, Int64}(zeros(n_buses, n_buses))
        lookup = Dict{Int.Int}()
    end
    return Ybus, lookup
end

function _make_device_index!(device::PSY.StaticInjection)
    dynamic_device = PSY.get_dynamic_injector(device)
    states = PSY.get_states(dynamic_device)
    device_state_mapping = Dict{Type{<:PSY.DynamicComponent}, Vector{Int64}}()
    input_port_mapping = Dict{Type{<:PSY.DynamicComponent}, Vector{Int64}}()
    _attach_inner_vars!(dynamic_device)
    _attach_control_refs!(device)

    for c in PSY.get_dynamic_components(dynamic_device)
        device_state_mapping[typeof(c)] = Vector{Int64}(undef, length(PSY.get_states(c)))
        input_port_mapping[typeof(c)] = Vector{Int64}()
        _index_local_states!(device_state_mapping[typeof(c)], states, c)
        _index_port_mapping!(input_port_mapping[typeof(c)], states, c)
        dynamic_device.ext[LOCAL_STATE_MAPPING] = device_state_mapping
        dynamic_device.ext[INPUT_PORT_MAPPING] = input_port_mapping
    end
    return
end

function _add_states_to_global!(
    global_state_index::MAPPING_DICT,
    state_space_ix::Vector{Int64},
    device::PSY.Device,
)
    global_state_index[PSY.get_name(device)] = Dict{Symbol, Int}()
    for s in PSY.get_states(device)
        state_space_ix[1] += 1
        global_state_index[PSY.get_name(device)][s] = state_space_ix[1]
    end

    return
end

function _get_internal_mapping(
    dynamic_device::PSY.DynamicInjection,
    key::AbstractString,
    ty::Type{T},
) where {T <: PSY.DynamicComponent}
    device_index = PSY.get_ext(dynamic_device)[key]
    val = get(device_index, ty, nothing)
    @assert !isnothing(val)
    return val
end

function get_local_state_ix(
    dynamic_device::PSY.DynamicInjection,
    ty::Type{T},
) where {T <: PSY.DynamicComponent}
    return _get_internal_mapping(dynamic_device, LOCAL_STATE_MAPPING, ty)
end

function get_input_port_ix(
    dynamic_device::PSY.DynamicInjection,
    ty::Type{T},
) where {T <: PSY.DynamicComponent}
    return _get_internal_mapping(dynamic_device, INPUT_PORT_MAPPING, ty)
end

function _simulation_pre_step(sim::Simulation, reset_simulation::Bool)
    if sim.status != BUILT && !reset_simulation
        error("The Simulation status is $(sim.status). Use keyword argument reset_simulation = true")
    end

    if reset_simulation
        reset!(sim)
    end
    return
end

function execute!(sim::Simulation, solver; kwargs...)
    reset_simulation = get(kwargs, :reset_simulation, false)
    _simulation_pre_step(sim, reset_simulation)
    sim.status = SIMULATION_STARTED
    sim.solution = DiffEqBase.solve(
        sim.problem,
        solver;
        callback = sim.callbacks,
        tstops = sim.tstops,
        kwargs...,
    )
    if sim.solution.retcode == :Success
        sim.status = SIMULATION_FINALIZED
    else
        sim.status = SIMULATION_FAILED
    end
    return
end

function _change_vector_type!(inputs::SimulationInputs, ::Type{T}) where {T <: Number}
    sys = get_system(inputs)
    for d in PSY.get_components(PSY.DynamicInjection, sys)
        _attach_inner_vars!(d, T)
    end
    _add_aux_arrays!(inputs, Real)
    return
end


text:
export User

"""
A Discord user.
More details [here](https://discordapp.com/developers/docs/resources/user#user-object).
"""
struct User
    id::Snowflake
    # The User inside of a Presence only needs its ID set.
    username::Union{String, Missing}
    discriminator::Union{String, Missing}
    avatar::Union{String, Nothing, Missing}
    bot::Union{Bool, Missing}
    mfa_enabled::Union{Bool, Missing}
    locale::Union{String, Missing}
    verified::Union{Bool, Missing}
    email::Union{String, Nothing, Missing}  # Not supposed to be nullable.
    # TODO: There's a member field here in one case.
end
@boilerplate User :dict :docs :lower :merge


text:
include("../src/hypergraph-affinity-functions.jl")
using StatsBase
using MAT
D = matread("TripAdvisor_NAM_Europe_H.mat")
H = D["H"]
classes = D["classes"]

class1 = "N. America"
class2 = "Europe"
n = size(H,2)
alpha = sum(classes)/n
# B1 = arbitrary_baselines(k,alpha)
# B0 = arbitrary_baselines(k,1-alpha)
## Bootstrap
alpha, B1, B2, H1, H2, tR1, tR2, N = Hypergraph_to_Scores(H,classes,9)

for k =  3:3:9
B = 100
TR1 = tR1[k,1:k]  # True affinity scores
TR0 = tR2[k,1:k]

TA1 = H1[k,1:k]  # True affinity scores
TA0 = H2[k,1:k]

# Bootstraps
R0, R1, MR1, MR0, SR1, SR0, A0, A1, MA1, MA0, SA1, SA0 = Bootstrap_Affinities(H,k,B,classes)

## Plot Ratio Scores
lw = 2.5
ms = 8
gfs = 15
tfs = 15
titlesize = 14
s1 = 400
s2 = 300
sc = 2
color1 = RGB(27/255,158/255,119/255)
color2 = RGB(217/255,95/255,2/255)
title = ""
leg = false
xlab = "Affinity type t"
ylab = "Affinity"
p = plot()
p = Plots.plot(fontfamily = "helvetica",linewidth = 2,yscale = :identity,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/3:k/3:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

# [M0-2*SE0 M0+2*SE]
# [E0m E0]
# [E1m E1]
p = plot!(1:k,[TA0 TA0], fillrange=[MA0-sc*SA0 MA0+sc*SA0], fillalpha=0.3, color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
p = plot!(1:k,[TA1 TA1], fillrange=[MA1-sc*SA1 MA1+sc*SA1],fillalpha=0.3,color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)

savefig("Figures/Affinities_$(k)_TA_Bootstrap.pdf")

## Ratio scores

xlab = "Affinity type t"
ylab = "Affinity / Baseline"
p = plot()
p = Plots.plot([1; k],[1; 1],fontfamily = "helvetica",linewidth = 2,yscale = :log10,legend = false,grid = false, size = (s1,s2),color = :gray)
title!(p,title,titlefontsize = titlesize)
xaxis!(p,xlab, xticks = k/3:k/3:k,tickfontsize=tfs,guidefontsize = gfs)
yaxis!(p,ylab,tickfontsize=tfs,guidefontsize = gfs)

# [E0m E0]
# [E1m E1]
p = plot!(1:k,[TR0 TR0], fillrange=[MR0-sc*SR0 MR0+sc*SR0], fillalpha=0.3, color = color2, linewidth = lw, markersize = ms,
    markershape = :circle, markerstrokecolor = color2, label = class2)
p = plot!(1:k,[TR1 TR1], fillrange=[MR1-sc*SR1 MR1+sc*SR1],fillalpha=0.3,color = color1, linewidth = lw, markersize = ms,
    markershape = :circle,markerstrokecolor = color1, label = class1)
savefig("Figures/Ratios_$(k)_TA_Bootstrap.pdf")
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
settings = class_settings();

context("nlp") do
    # - x - y

    a = sparse([-2.0, -2.0, 0.0])
    b = 0.0
    P = spdiagm([0.0, 0.0, 0.0])
    obj = nl_quad(a,P,b)

    # x^2 + y^2 - s - 4.0
    # x^2 + y^2 >= 4.0
    a1 = sparse([0.0, 0.0, -1.0])
    P1 = spdiagm([2.0, 2.0, 0.0])
    b1 = -4.0
    constr1 = nl_quad(a1,P1,b1)

    a1 = sparse([0.0, 0.0, -1.0])
    P1 = spdiagm([2.0, 2.5, 0.0])
    b1 = -4.0
    constr1_b = nl_quad(a1,P1,b1)

    # x + y + s - 1.0 = 0.0
    # x + y = 1.0
    a2 = sparse([1.0, 1.0, 0.0])
    P2 = spdiagm([0.0, 0.0, 0.0])
    b2 = -1.0 #-4.0
    constr2 = nl_quad(a2,P2,b2)

    # x + x^2 - y = 0
    # x + x^2 = y
    a3 = sparse([1.0, -1.0, 0.0])
    P3 = spdiagm([1.0, 0.0, 0.0])
    b3 = 0.0 #-4.0
    constr3 = nl_quad(a3,P3,b3)

    # x^2 + y^2 - 6x - 6y - s + 17 = 0
    # (x-3)^2 + (y-3)^2 - s - 1 = 0
    # (x-3)^2 + (y-3)^2 <= 1
    # ball function?
    a3 = sparse([-6.0, -6.0, -1.0])
    P3 = spdiagm([1.0, 1.0, 0.0])
    b3 = -17.0 #-4.0
    constr4 = nl_quad(a3,P3,b3)

    # convex problem
    println("CONVEX INFEASIBLE PROBLEM")
    nlp = class_nlp(3, obj, [constr2, constr4])
    vars = class_variables(3,2);
    status, vars = homogeneous_algorithm(nlp, vars, settings)
    @show x_scaled(vars)
    @show y_scaled(vars)

    println("NONCONVEX INFEASIBLE PROBLEM 1")
    nlp = class_nlp(3, obj, [constr1, constr2])
    vars = class_variables(3,2);
    status, vars = homogeneous_algorithm(nlp, vars, settings)
    @show x_scaled(vars)
    @show y_scaled(vars)

    println("NONCONVEX INFEASIBLE PROBLEM 2")
    nlp = class_nlp(3, obj, [constr1_b, constr2])
    vars = class_variables(3,2);
    status, vars = homogeneous_algorithm(nlp, vars, settings)
    @show x_scaled(vars)
    @show y_scaled(vars)

    nlp = class_nlp(3, obj, [constr1, constr3])
    vars = class_variables(3,2);
    status, vars = homogeneous_algorithm(nlp, vars, settings)
    @show x_scaled(vars)
    @show y_scaled(vars)
    #@show internal_eval_hesslag_prod(nlp,x_scaled(vars),y_scaled(vars))

    nlp = class_nlp(3, obj, [constr2, constr3])
    vars = class_variables(3,2);
    status, vars = homogeneous_algorithm(nlp, vars, settings)
    @show x_scaled(vars)
    @show y_scaled(vars)

    #@fact status --> :locally_optimal
    #@fact norm(x_scaled(vars) - [2.0,2.0,0.0]) --> less_than(1e-6)
    #@fact tau(vars) --> greater_than(0.1)
    #@fact 0.0 --> less_than(kappa(vars))
    #@fact kappa(vars) --> less_than(1e-6)
end

max_forks_repo_name: ohinder/Everest


text:
function check_for_numerical_none_value(value::String)::Union{Missing,Float64}

    return_value = missing
    if (value != "None")
        return_value = parse(Float64, value)
    end

    # return -
    return return_value
end