text:
# control projection problem
mf = 3
nz = mf + 1 + 1 + 1 + 1 + mf
nθ = mf + 1

idx = [3; 1; 2]
function residual(z, θ, κ)
  u = z[1:mf]
  s = z[mf .+ (1:1)]
  y = z[mf + 1 .+ (1:1)]
  w = z[mf + 1 + 1 .+ (1:1)]
  p = z[mf + 1 + 1 + 1 .+ (1:1)]
  v = z[mf + 1 + 1 + 1 + 1 .+ (1:mf)]

  ū = θ[1:mf]
  T = θ[mf .+ (1:1)]

  [
   u - ū - v - [0.0; 0.0; y[1] + p[1]];
   -y[1] - w[1];
   T[1] - u[3] - s[1];
   w .* s .- κ
   p .* u[3] .- κ
   second_order_cone_product(v[idx], u[idx]) .- κ .* [1.0; 0.0; 0.0]
  ]
end

@variables z[1:nz], θ[1:nθ], κ[1:1]

r = residual(z, θ, κ)
r .= simplify(r)
r_func = eval(Symbolics.build_function(r, z, θ, κ)[2])

rz = Symbolics.jacobian(r, z)
rz = simplify.(rz)
rz_func = eval(Symbolics.build_function(rz, z, θ)[2])

rθ = Symbolics.jacobian(r, θ)
rθ = simplify.(rθ)
rθ_func = eval(Symbolics.build_function(rθ, z, θ)[2])

idx_ineq = collect([3, 4, 6, 7])
idx_soc = [collect([3, 1, 2]), collect([10, 8, 9])]#collect([3, 1, 2]), collect([10, 8, 9])]

# ul = -1.0 * ones(m)
# uu = 1.0 * ones(m)
ū = [0.0, 0.0, 0.0]

z0 = 0.1 * ones(nz)
z0[1:mf] = copy(ū)
z0[3] += 1.0
z0[10] += 1.0

θ0 = [copy(ū); 1.0]

# solver
opts_con = InteriorPointOptions(
    κ_init = 1.0,
    κ_tol = 1.0e-4,
    r_tol = 1.0e-8,
    diff_sol = false)

ip_con = interior_point(z0, θ0,
    r! = r_func, rz! = rz_func,
    rz = similar(rz, Float64),
    rθ! = rθ_func,
    rθ = similar(rθ, Float64),
    idx_ineq = idx_ineq,
    idx_soc = idx_soc,
    opts = opts_con)

opts_jac = InteriorPointOptions(
	κ_init = 1.0,
	κ_tol = 1.0e-2,
	r_tol = 1.0e-8,
	diff_sol = true)

ip_jac = interior_point(z0, θ0,
	r! = r_func, rz! = rz_func,
	rz = similar(rz, Float64),
	rθ! = rθ_func,
	rθ = similar(rθ, Float64),
	idx_ineq = idx_ineq,
  idx_soc = idx_soc,
	opts = opts_jac)

interior_point_solve!(ip_con)
interior_point_solve!(ip_jac)

function soc_projection(x, x_min, x_max, scaling)
  proj = second_order_cone_projection(x[[3;1;2]])[1]
	ip_con.z .= [proj[2:3]; proj[1]; 0.1 * ones(7)]
  ip_con.z[3] += max(1.0, norm(proj[2:3])) * 2.0
  ip_con.z[10] += 1.0
	ip_con.θ .= [x; scaling * x_max]

	status = interior_point_solve!(ip_con)

  !status && (@warn "projection failure (res norm: $(norm(ip_con.r, Inf))) \n
		               z = $(ip_con.z), \n
					   θ = $(ip_con.θ)")

	return ip_con.z[1:mf]
end

soc_projection([100.0, 0.0, -1.0], 0.0, 10.0, 0.5)

function soc_projection_jacobian(x, x_min, x_max, scaling)
  proj = second_order_cone_projection(x[[3;1;2]])[1]
	ip_con.z .= [proj[2:3]; proj[1]; 0.1 * ones(7)]
  ip_con.z[3] += max(1.0, norm(proj[2:3])) * 2.0
  ip_con.z[10] += 1.0
	ip_con.θ .= [x; scaling * x_max]

	interior_point_solve!(ip_jac)

	return ip_jac.δz[1:mf, 1:mf]
end

soc_projection_jacobian([10.0, 0.0, 1.0], 0.0, 10.0, 0.5)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
__precompile__(true)

# TODO remove preprocessing stuff, it doesn't belong here
# TODO implement padding for TimeSeriesHandler

module DataHandlers

using Reexport

@reexport using DataTables
using DataUtils

include("DataHandler.jl")
include("TimeSeriesHandler.jl")
include("GroupedDataHandler.jl")

end

ext: jl


text:
using Test, LinearAlgebra, Random, SparseArrays, ExponentialUtilities
using ExponentialUtilities: getH, getV, _exp!

@testset "Exp" begin
    n = 100
    A = randn(n, n)
    expA = exp(A)
    _exp!(A)
    @test A ≈ expA
    A2 = randn(n, n)
    A2 ./= opnorm(A2, 1) # test for small opnorm
    expA2 = exp(A2)
    _exp!(A2)
    @test A2 ≈ expA2
end

@testset "Phi" begin
    # Scalar phi
    K = 4
    z = 0.1
    P = fill(0., K+1); P[1] = exp(z)
    for i = 1:K
        P[i+1] = (P[i] - 1/factorial(i-1))/z
    end
    @test phi(z, K) ≈ P

    # Matrix phi (dense)
    A = [0.1 0.2; 0.3 0.4]
    P = Vector{Matrix{Float64}}(undef, K+1); P[1] = exp(A)
    for i = 1:K
        P[i+1] = (P[i] - 1/factorial(i-1)*I) / A
    end
    @test phi(A, K) ≈ P

    # Matrix phi (Diagonal)
    A = Diagonal([0.1, 0.2, 0.3, 0.4])
    Afull = Matrix(A)
    P = phi(A, K)
    Pfull = phi(Afull, K)
    for i = 1:K+1
        @test Matrix(P[i]) ≈ Pfull[i]
    end
end

@testset "Arnoldi & Krylov" begin
    Random.seed!(0)
    n = 20; m = 5; K = 4
    A = randn(n, n)
    t = 1e-2
    b = randn(n)
    direct = exp(t * A) * b
    @test direct ≈ expv(t, A, b; m=m)
    @test direct ≈ kiops(t, A, b)[1]
    P = phi(t * A, K)
    W = fill(0., n, K+1)
    for i = 1:K+1
        W[:,i] = P[i] * b
    end
    Ks = arnoldi(A, b; m=m)
    W_approx = phiv(t, Ks, K)
    @test W ≈ W_approx
    W_approx_kiops3 = kiops(t, A, hcat([b*inv(t)^i for i in 0:K-1]...))
    @test sum(W[:, 1:K], dims=2) ≈ W_approx_kiops3[1]
    @test_skip begin
        W_approx_kiops4 = kiops(t, A, hcat([b*inv(t)^i for i in 0:K]...))
        @test_broken sum(W[:, 1:K+1], dims=2) ≈ W_approx_kiops4[1]
        @test sum(W[:, 1:K+1], dims=2) ≈ W_approx_kiops4[1] atol=1e-2
    end

    # Happy-breakdown in Krylov
    v = normalize(randn(n))
    A = v * v' # A is Idempotent
    Ks = arnoldi(A, b)
    @test Ks.m == 2

    # Arnoldi vs Lanczos
    A = Hermitian(randn(n, n))
    Aperm = A + 1e-10 * randn(n, n) # no longer Hermitian
    w = expv(t, A, b; m=m)
    wperm = expv(t, Aperm, b; m=m, opnorm=opnorm)
    wkiops = kiops(t, A, b; m=m)[1]
    @test w ≈ wperm
    @test w ≈ wkiops
end

@testset "Complex Value" begin
    n = 20; m = 10;
    for A in [Hermitian(rand(ComplexF64, n, n)), Hermitian(rand(n, n)), rand(ComplexF64, n, n), rand(n, n)]
        for b in [rand(ComplexF64, n), rand(n)], t in [1e-2, 1e-2im, 1e-2 + 1e-2im]
            @test exp(t * A) * b ≈ expv(t, A, b; m=m)
        end
    end
end

@testset "Adaptive Krylov" begin
    # Internal time-stepping for Krylov (with adaptation)
    n = 100
    K = 4
    t = 5.0
    tol = 1e-7
    A = spdiagm(-1=>ones(n-1), 0=>-2*ones(n), 1=>ones(n-1))
    B = randn(n, K+1)
    Phi_half = phi(t/2 * A, K)
    Phi = phi(t * A, K)
    uhalf_exact = sum((t/2)^i * Phi_half[i+1] * B[:,i+1] for i = 0:K)
    u_exact = sum(t^i * Phi[i+1] * B[:,i+1] for i = 0:K)
    U = phiv_timestep([t/2, t], A, B; adaptive=true, tol=tol)
    @test norm(U[:,1] - uhalf_exact) / norm(uhalf_exact) < tol
    @test norm(U[:,2] - u_exact) / norm(u_exact) < tol
    # p = 0 special case (expv_timestep)
    u_exact = Phi[1] * B[:, 1]
    u = expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol, opnorm=opnorm)
    @test_nowarn expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol, opnorm=opnorm(A, Inf))
    @test norm(u - u_exact) / norm(u_exact) < tol
end

@testset "Krylov for Hermitian matrices" begin
    # Hermitian matrices have real spectra. Ensure that the subspace
    # matrix is representable as a real matrix.

    n = 100
    m = 15
    tol = 1e-14

    e = ones(n)
    p = -im*Tridiagonal(-e[2:end], 0e, e[2:end])

    KsA = KrylovSubspace{ComplexF64}(n, m)
    KsL = KrylovSubspace{ComplexF64, Float64}(n, m)

    v = rand(ComplexF64, n)

    arnoldi!(KsA, p, v)
    lanczos!(KsL, p, v)

    AH = view(KsA.H,1:KsA.m,1:KsA.m)
    LH = view(KsL.H,1:KsL.m,1:KsL.m)

    @test norm(AH-LH)/norm(AH) < tol
end

@testset "Alternative Lanczos expv Interface" begin
    n = 300
    m = 30

    A = Hermitian(rand(n,n))
    b = rand(ComplexF64, n)
    dt = 0.1

    atol=1e-10
    rtol=1e-10
    w = expv(-im, dt*A, b, m=m, tol=atol, rtol=rtol, mode=:error_estimate)

    function fullexp(A, v)
        w = similar(v)
        eA = exp(A)
        mul!(w, eA, v)
        w
    end

    w′ = fullexp(-im*dt*A, b)

    δw = norm(w-w′)
    @test δw < atol
    @test δw/abs(1e-16+norm(w)) < rtol
end

struct MatrixFreeOperator{T} <: AbstractMatrix{T}
    A::Matrix{T}
end
Base.eltype(A::MatrixFreeOperator{T}) where T = T
LinearAlgebra.mul!(y::AbstractVector, A::MatrixFreeOperator, x::AbstractVector) = mul!(y, A.A, x)
Base.size(A::MatrixFreeOperator, dim) = size(A.A, dim)
struct OpnormFunctor end
(::OpnormFunctor)(A::MatrixFreeOperator, p::Real) = opnorm(A.A, p)
@testset "Matrix-free Operator" begin
    Random.seed!(123)
    n = 20
    for ishermitian in (false, true)
        A = rand(ComplexF64, n, n)
        M = ishermitian ? A'A : A
        Op = MatrixFreeOperator(M)
        b = rand(ComplexF64, n)
        Ks = arnoldi(Op, b; ishermitian=ishermitian, opnorm=OpnormFunctor(), tol=1e-12)
        pv = phiv(0.01, Ks, 2)
        pv′ = hcat(map(A->A*b, phi(0.01Op.A, 2))...)

        @test pv ≈ pv′ atol=1e-12
    end
end


text:
using ThreeBodyTB
using Plots

#setup chosen Plots backend, i recommend pyplot
#pyplot()
#gr()

#make the crystal object
#here we choose Sc P rocksalt

types=["Sc", "P"];

#positions, crystal units
pos = [0 0 0 ; 0.5000000000  0.5000000000  0.5000000000]
#lattice vectors, in Angstrom units currently
A=[ [4.7 4.7 0]; [4.7 0 4.7 ]; [ 0 4.7 4.7]] * 1.05 * 0.529177;

#makes the crystal
c=makecrys(A, pos, types)

println("Starting crystal is")
println(c)

println("SCF ENERGY ---------------------------------------------------------------")

energy, tbc, flag = scf_energy(c)

println("BAND STRUCTURE ---------------------------------------------------------------")

#plot band structure
kpath=[0.0 0.0 0.0;0.5 0 0; 0.5 0.25 0.75; 0.0 0.0 0.0]
knames=["Γ", "X", "W", "Γ"]

plot_bandstr(tbc, kpath=kpath, names=knames, proj_types=["P"])

savefig("ScP_rocksalt.pdf")

println()
println("DOS ---------------------------------------------------------------")
#Density of states. Default is atomic projection.
dos(tbc)

savefig("ScP_rocksalt_DOS.pdf")

#dos(tbc, proj_type=:orbs) #orbital projection. :none also allowed.
#savefig("ScP_rocksalt_DOS_orbital.pdf")


text:
using Clang

DAILOC="/home/bana/GSP/research/samc/samcnet/deps/libdai"

context = wrap_c.init(output_file="dai_wrap.jl",
    common_file="dai_common.jl",
    header_library=x->"\"src/libdaiwrap.so\"",
    clang_includes=["$DAILOC/include"]
    )

headers = ["src/dai.h"]
wrap_c.wrap_c_headers(context, headers)

# build shared library
cd("src") do
    run(`clang++ -shared -fPIC -o libdaiwrap.so -ldai -lgmp -lgmpxx -L$DAILOC/lib -I$DAILOC/include dai.cpp`)
end


text:
using LatBo

mmHg = 133.3223684211
Npoints = 30, 30
δt, δx = 5e-4, 1e-3
μ₀ = 0.01sqrt(1./3.) * δt / δx
sim = lbgk(:D2Q9, Npoints, δt, δx; viscosity=4e-1, p₀=80.0 * mmHg, μ₀ = ones(2)μ₀)
Playground.initialize(sim.playground) do i, j
    coords = Float64[i, j]
    if j == 1 || j == Npoints[2] || i == 1 || i == Npoints[1]
        return Playground.SOLID
    else
        return Playground.FLUID
    end
end

run!(sim, nsteps=2max(Npoints...))

Profile.clear_malloc_data()
run!(sim, nsteps=500, doinit=false)


text:
#these four macros convert the Vnum value into the appropriate type of rounding
#sigmoid by calling the call constructors which are designed to handle bumping
#up or down the ubit.

typeas(value::Valid{N,ES}, stype::Symbol) where {N,ES} = Sigmoid{N,ES,stype}
Base.:~(value::Valid{N,ES}, stype::Symbol) where {N,ES} = typeas(value, stype)

macro upper(value)
  esc(:(isulp(($value).upper) ? reinterpret(($value ~ :upper), ($value).upper |> lub) : reinterpret(($value ~ :exact), ($value).upper)))
end

macro lower(value)
  esc(:(isulp(($value).lower) ? reinterpret(($value ~ :lower), ($value).lower |> glb) : reinterpret(($value ~ :exact), ($value).lower)))
end

macro exact(value)
  esc(:(reinterpret(($value ~ :exact), ($value).lower)))
end

"""
    SigmoidNumbers.@d_lower(value)

retrieves a directionally typecast lower value.  This could be typecast as one
of: :inward_exact, :inward_ulp, :outward_exact, :outward_ulp, depending on
the nature of the values.

#Examples

"""
macro d_lower(value)
    esc(quote
        if Base.iszero($value.lower)
            Base.iszero($value.upper) ? zero($value ~ :inward_exact)                            : zero($value ~ :outward_exact)
        elseif !isfinite($value.lower)
            isfinite($value.upper)    ? reinterpret(($value ~ :inward_exact), $value.lower)     : reinterpret(($value ~ :outward_exact), $value.lower)
        elseif ($value.lower < zero($value ~ :ubit))
            isulp($value.lower)       ? reinterpret(($value ~ :inward_ulp),  glb($value.lower)) : reinterpret(($value ~ :inward_exact),  $value.lower)
        else
            isulp($value.lower)       ? reinterpret(($value ~ :outward_ulp), glb($value.lower)) : reinterpret(($value ~ :outward_exact), $value.lower)
        end
    end)
end

macro d_upper(value)
    esc(quote
        if Base.iszero($value.upper)
            Base.iszero($value.lower) ? zero($value ~ :inward_exact)                            : zero($value ~ :outward_exact)
        elseif !isfinite($value.upper)
            isfinite($value.lower)    ? reinterpret(($value ~ :inward_exact), $value.upper)     : reinterpret(($value ~ :outward_exact), $value.upper)
        elseif ($value.upper < zero($value ~ :ubit))
            isulp($value.upper)       ? reinterpret(($value ~ :outward_ulp), lub($value.upper)) : reinterpret(($value ~ :outward_exact), $value.upper)
        else
            isulp($value.upper)       ? reinterpret(($value ~ :inward_ulp),  lub($value.upper)) : reinterpret(($value ~ :inward_exact),  $value.upper)
        end
    end)
end

recast_as_lower(x::Sigmoid{N,ES,:exact}) where {N,ES} = x
recast_as_lower(x::Sigmoid{N,ES,:lower}) where {N,ES} = x
recast_as_lower(x::Sigmoid{N,ES,:cross}) where {N,ES} = reinterpret(Sigmoid{N,ES,:lower}, x)
recast_as_lower(x::Sigmoid{N,ES,:upper}) where {N,ES} = reinterpret(Sigmoid{N,ES,:lower}, x)
macro rl(value)
  esc(:(recast_as_lower($value)))
end

recast_as_upper(x::Sigmoid{N,ES,:exact}) where {N,ES} = x
recast_as_upper(x::Sigmoid{N,ES,:upper}) where {N,ES} = x
recast_as_upper(x::Sigmoid{N,ES,:cross}) where {N,ES} = reinterpret(Sigmoid{N,ES,:upper}, x)
recast_as_upper(x::Sigmoid{N,ES,:lower}) where {N,ES} = reinterpret(Sigmoid{N,ES,:upper}, x)
macro ru(value)
  esc(:(recast_as_upper($value)))
end

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
using FastAI
using Test

@testset "FastAI.jl" begin
    include("utils.jl")
    include("test_learner.jl")
    include("test_recorder.jl")
end


text:
# run from REPL

using Makie

φ = 0:π/100:2π

θ = 0:π/200:π

x = [cos(θ) * sin(φ) for θ in θ, φ in φ]
y = [sin(θ)*sin(φ) for θ in θ, φ in φ]
z = [cos(φ) for θ in θ, φ in φ]

surface(x, y, z, backgroundcolor = :black, show_axis = false)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Printf
"""
Funções para Análise de Pórtico Espacial

Sistema de Coordenadas

* X = Horizontal, positivo para direita
* Y = Vertical, positivo para cima
* Z = Horizontal, positivo para direita

"""
module PorticoEspacial

"""Constrói matriz de transformação do elemento"""
function matriztransformacao(noinicial::Vector{Float64}, nofinal::Vector{Float64}, gama::Float64 = 1.0)::Matrix{Float64}
    Tr::Matrix{Float64} = zeros(3, 3)
    TrAux::Matrix{Float64} = zeros(3, 3)
    Cx = (nofinal[1] - noinicial[1])
    Cy = (nofinal[2] - noinicial[2])
    Cz = (nofinal[3] - noinicial[3])
    L = sqrt(sqrt(Cx) + sqrt(Cy) + sqrt(Cz))
    Cx /= L
    Cy /= L
    Cz /= L
    cosBeta = sqrt(sqrt(Cx) + sqrt(Cz))
    sinBeta = Cy

    if cosBeta > 0.0
        cosAlfa = Cx / cosBeta
        sinAlfa = Cz / cosBeta
    else
        cosAlfa = 1.0
        sinAlfa = 0.0
    end

    cosGama = cos(gama)
    sinGama = sin(gama)

    Tr[1, 1] = cosAlfa
    Tr[1, 3] = sinBeta
    Tr[2, 2] = 1.0
    Tr[3, 1] = -sinAlfa
    Tr[3, 3] = cosAlfa

    TrAux[1, 1] = cosBeta
    TrAux[1, 2] = sinBeta
    TrAux[2, 1] = -sinBeta
    TrAux[2, 2] = cosBeta
    TrAux[3, 3] = 1.0
    Tr = TrAux * Tr

    TrAux[1, 1] = 1.0
    TrAux[2, 2] = cosGama
    TrAux[2, 3] = sinGama
    TrAux[3, 2] = -sinGama
    TrAux[3, 3] = cosGama
    Tr = TrAux * Tr

    return Tr
end

"""Constrói matriz de rotação da barra"""
function matrizrotacaoespacial(Tr::Matrix{Float64})::Matrix{Float64}
    rotBar::Matrix{Float64} = zeros(12, 12)

    m = 0
    for _ in 1:4
        for i in 1:3
            for j in 1:3
                rotBar[m+i, m+j] = Tr[i, j]
            end
        end
        m += 3
    end

    return rotBar
end

"""Inclui articulação na matriz de rigidez do elemento"""
function addarticulacao(matrizrigidez::Matrix{Float64})::Matrix{Float64}
    return nothing
end

"""Constrói matriz de rigidez local"""
function matrizrigidezlocal(noinicial::Vector{Float64}, nofinal::Vector{Float64})::Matrix{Float64}
    return nothing
end

"""Constrói matriz de rigidez global"""
function matrizrigidezglobal(matrizrigidezlocal::Matrix{Float64})::Matrix{Float64}
    return nothing
end

"""Constrói vetor apontador"""
function vetorapontador()::Vector{Float64}
    return nothing
end

"""Constrói matriz de rigidez da estrutura"""
function matrizrigidezestrutura()::Matrix{Float64}
    return nothing
end

"""Verifica condições de contorno"""
function condicoesgeometricascontorno() end

"""Entry point para executar a análise"""
function execute()
    # Carregar Materiais
    # Configurar propriedades das seções
    # Carregar coordenadas nodais
    # Carregar Elementos
    # Carregar apoios
    # Carregar Cargas

    vetorapontador()
    matrizrigidezestrutura()
    condicoesgeometricascontorno()
    # Resolver sistema de equações
end
end

println("========================================== ANÁLISE DE ESTRUTURAS ==========================================\n")
println("1. Cargas\n")
println("2. Ações\n")
println("3. Deslocamentos Iniciais\n")
println("4. Matrizes de Rigidez Locais\n")
println("5. Matriz de Rigidez Global\n")

println()
println("=============================================== RESULTADO =================================================\n")
println("1. Materiais\n")
println("ID       Módulo E   Módulo G   Densidade")
println("-------- ---------- ---------- ----------")

@printf "%s % 10.6f % 10.2f % 10.2f\n" "CONCR20M" 123E-6 3.4 2345.45
@printf "%s % 10.6f % 10.2f % 10.2f\n" "CONCR30M" 123E-6 3.4 2345.45
@printf "%s % 10.6f % 10.2f % 10.2f\n" "CONCR40M" 123E-6 3.4 2345.45
println()
println()
println("2. Elementos\n")

println("ID        Material Comprimento  Área Seção   Nó Inicial Nó Final   Anotação")
println("--------- -------- ------------ ------------ ---------- ---------- ----------------------------------------")
@printf "%s %s % 12.2f % 12.2f %s %s\n" "V101     " "CONCR40M" 12.34 12.34 "NO00000001" "NO00000002"
@printf "%s %s % 12.2f % 12.2f %s %s\n" "V101     " "CONCR40M" 12.34 12.34 "NO00000001" "NO00000002"
@printf "%s %s % 12.2f % 12.2f %s %s\n" "V101     " "CONCR40M" 12.34 12.34 "NO00000001" "NO00000002"
@printf "%s %s % 12.2f % 12.2f %s %s\n" "V101     " "CONCR40M" 12.34 12.34 "NO00000001" "NO00000002"

println()
println("NO00000001                                             NO00000002")
println("----------------------------------------------------   ----------------------------------------------------")
@printf "X1: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m   X2: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Y1: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m   Y2: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Z1: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m   Z2: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
println()

println("NO00000001                                             NO00000002")
println("----------------------------------------------------   ----------------------------------------------------")
@printf "X1: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m   X2: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Y1: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m   Y2: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Z1: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m   Z2: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
println()

println("NO00000001                                             NO00000002")
println("----------------------------------------------------   ----------------------------------------------------")
@printf "X1: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m   X2: % 8.2f m  Fx: % 8.2f kgf  Mx: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Y1: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m   Y2: % 8.2f m  Fy: % 8.2f kgf  My: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56
@printf "Z1: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m   Z2: % 8.2f m  Fz: % 8.2f kgf  Mz: % 8.2f kgf.m\n" 12.34 1234.56 1234.56 12.34 1234.56 1234.56

println()

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
"""
This is the [`DocumenterShowcase`](@ref), which contains

* [`DocumenterShowcase.foo`](@ref)
* [`DocumenterShowcase.bar`](@ref)

!!! compat "Documenter 0.24"

    This showcase page is included in Documenter version 0.24.0 and above.

# Contents

Docstrings can contain admonitions and other block-level nodes.

!!! info

    This admonition is in a docstring. It itself can contain block levels nodes such
    as code blocks:

    ```julia
    println("Hello World")
    ```

    ... or block quotes:

    > Lorem ipsum.

In fact, while not recommended, you can actually have a matryoshka of admonitions:

!!! danger
    !!! warning
        !!! tip
            !!! note
                Stack overflow.
"""
module DocumenterShowcase

"""
    foo(::Integer)

Docstring for `foo(::Integer)`.
"""
foo(::Integer) = nothing

"""
    foo(::AbstractString)

Docstring for `foo(::AbstractString)`.
"""
foo(::AbstractString) = nothing

"""
    bar(::Integer)

Docstring for `bar(::Integer)`.
"""
bar(::Integer) = nothing

"""
    bar(::AbstractString)

Docstring for `bar(::AbstractString)`.
"""
bar(::AbstractString) = nothing

function hello(who)
    println("Hello, $(who)!")
end

struct SVGCircle
    stroke :: String
    fill :: String
end
function Base.show(io, ::MIME"image/svg+xml", c::SVGCircle)
    write(io, """
    <svg width="50" height="50">
      <g style="stroke-width: 3">
        <circle cx="25" cy="25" r="24" stroke-width="2" style="stroke: #$(c.stroke); fill: #$(c.fill)" />
      </g>
    </svg>
    """)
end

"The type definition."
struct Foo{T,S} end

"Constructor `Foo()` with no arguments."
Foo() = Foo{Nothing,Nothing}()

"Constructor `Foo{T}()` with one parametric argument."
Foo{T}() where T = Foo{T,Nothing}()

end # module


text:
function get_inl_sections()
    return ["GOVERNOR PERMANANT DROOP"]
end

function get_inl_section_info()
    governor_droop = [["I", 1, Int], ["ID", 2, String], ["R", 6, Float64]]

    section_info = Dict{String,Array}(
        "GOVERNOR PERMANANT DROOP" => governor_droop)

    return section_info
end

function idx_inl()
    INL_BUS = 1
    INL_ID = 2
    INL_R = 3

    return INL_BUS, INL_ID, INL_R
end

function parse_inl(filename::AbstractString; delim=",")
    f = open(filename, "r")
    lines = readlines(f)
    close(f)

    sections = get_inl_sections()
    section_info = get_inl_section_info()
    pos = 1

    inldata = Dict{String,Array}()

    # There is only one section in case.inl file.
    s = sections[1]
    pos, data = parse_data(lines, pos, section_info[s])
    inldata[s] = data

    return inldata
end


text:
""" Test the Syntax module.

The unit tests are sparse because many of the Doctrine tests are really just
tests of the Syntax module.
"""
module TestSyntax

using Test
using Catlab

# Monoid
########

""" Theory of monoids.
"""
@signature Monoid(Elem) begin
  Elem::TYPE
  munit()::Elem
  mtimes(x::Elem,y::Elem)::Elem
end

""" Syntax for the theory of monoids.
"""
@syntax FreeMonoid Monoid

Elem(mod::Module, args...) = Elem(mod.Elem, args...)

@test isa(FreeMonoid, Module)
@test occursin("theory of monoids", string(Docs.doc(FreeMonoid)))
@test sort(names(FreeMonoid)) == sort([:FreeMonoid, :Elem])

x, y, z = Elem(FreeMonoid,:x), Elem(FreeMonoid,:y), Elem(FreeMonoid,:z)
@test isa(mtimes(x,y), FreeMonoid.Elem)
@test isa(munit(FreeMonoid.Elem), FreeMonoid.Elem)
@test mtimes(mtimes(x,y),z) != mtimes(x,mtimes(y,z))

# Test equality
@test x == Elem(FreeMonoid,:x)
@test x != y
@test Elem(FreeMonoid,"X") == Elem(FreeMonoid,"X")
@test Elem(FreeMonoid,"X") != Elem(FreeMonoid,"Y")

# Test hash
@test hash(x) == hash(x)
@test hash(x) != hash(y)
@test hash(mtimes(x,y)) == hash(mtimes(x,y))
@test hash(mtimes(x,y)) != hash(mtimes(x,z))

@syntax FreeMonoidAssoc Monoid begin
  mtimes(x::Elem, y::Elem) = associate(new(x,y))
end

x, y, z = [ Elem(FreeMonoidAssoc,sym) for sym in [:x,:y,:z] ]
e = munit(FreeMonoidAssoc.Elem)
@test mtimes(mtimes(x,y),z) == mtimes(x,mtimes(y,z))
@test mtimes(e,x) != x && mtimes(x,e) != x

@syntax FreeMonoidAssocUnit Monoid begin
  mtimes(x::Elem, y::Elem) = associate_unit(new(x,y), munit)
end

x, y, z = [ Elem(FreeMonoidAssocUnit,sym) for sym in [:x,:y,:z] ]
e = munit(FreeMonoidAssocUnit.Elem)
@test mtimes(mtimes(x,y),z) == mtimes(x,mtimes(y,z))
@test mtimes(e,x) == x && mtimes(x,e) == x

abstract type MonoidExpr{T} <: GATExpr{T} end
@syntax FreeMonoidTyped(MonoidExpr) Monoid

x = Elem(FreeMonoidTyped.Elem, :x)
@test FreeMonoidTyped.Elem <: MonoidExpr
@test isa(x, FreeMonoidTyped.Elem) && isa(x, MonoidExpr)

@signature Monoid(Elem) => MonoidNumeric(Elem) begin
  elem_int(x::Int)::Elem
end
@syntax FreeMonoidNumeric MonoidNumeric

x = elem_int(FreeMonoidNumeric.Elem, 1)
@test isa(x, FreeMonoidNumeric.Elem)
@test first(x) == 1

""" A monoid with two distinguished elements.
"""
@signature Monoid(Elem) => MonoidTwo(Elem) begin
  one()::Elem
  two()::Elem
end
""" The free monoid on two generators.
"""
@syntax FreeMonoidTwo MonoidTwo begin
  Elem(::Type{Elem}, value) = error("No extra generators allowed!")
end

x, y = one(FreeMonoidTwo.Elem), two(FreeMonoidTwo.Elem)
@test all(isa(expr, FreeMonoidTwo.Elem) for expr in [x, y, mtimes(x,y)])
@test_throws ErrorException Elem(FreeMonoidTwo, :x)

# Category
##########

@signature Category(Ob,Hom) begin
  Ob::TYPE
  Hom(dom::Ob, codom::Ob)::TYPE

  id(X::Ob)::Hom(X,X)
  compose(f::Hom(X,Y), g::Hom(Y,Z))::Hom(X,Z) ⊣ (X::Ob, Y::Ob, Z::Ob)
end

@syntax FreeCategory Category begin
  compose(f::Hom, g::Hom) = associate(new(f,g))
end

@test isa(FreeCategory, Module)
@test sort(names(FreeCategory)) == sort([:FreeCategory, :Ob, :Hom])

X, Y, Z, W = [ Ob(FreeCategory.Ob, sym) for sym in [:X, :Y, :Z, :W] ]
f, g, h = Hom(:f, X, Y), Hom(:g, Y, Z), Hom(:h, Z, W)
@test isa(X, FreeCategory.Ob) && isa(f, FreeCategory.Hom)
@test_throws MethodError FreeCategory.Hom(:f)
@test dom(f) == X
@test codom(f) == Y

@test isa(id(X), FreeCategory.Hom)
@test dom(id(X)) == X
@test codom(id(X)) == X

@test isa(compose(f,g), FreeCategory.Hom)
@test dom(compose(f,g)) == X
@test codom(compose(f,g)) == Z
@test isa(compose(f,f), FreeCategory.Hom) # Doesn't check domains.
@test compose(compose(f,g),h) == compose(f,compose(g,h))

@syntax FreeCategoryStrict Category begin
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

X, Y = Ob(FreeCategoryStrict.Ob, :X), Ob(FreeCategoryStrict.Ob, :Y)
f, g = Hom(:f, X, Y), Hom(:g, Y, X)

@test isa(compose(f,g), FreeCategoryStrict.Hom)
@test_throws SyntaxDomainError compose(f,f)

# Functor
#########

@instance Monoid(String) begin
  munit(::Type{String}) = ""
  mtimes(x::String, y::String) = string(x,y)
end

F(expr; kw...) = functor((String,), expr; kw...)

x, y, z = Elem(FreeMonoid,:x), Elem(FreeMonoid,:y), Elem(FreeMonoid,:z)
gens = Dict(x => "x", y => "y", z => "z")
@test F(mtimes(x,mtimes(y,z)); generators=gens) == "xyz"
@test F(mtimes(x,munit(FreeMonoid.Elem)); generators=gens) == "x"

terms = Dict(:Elem => (x) -> string(first(x)))
@test F(mtimes(x,mtimes(y,z)); terms=terms) == "xyz"
@test F(mtimes(x,munit(FreeMonoid.Elem)); terms=terms) == "x"

# Serialization
###############

# To JSON
X, Y, Z = [ Ob(FreeCategory.Ob, sym) for sym in [:X, :Y, :Z] ]
f = Hom(:f, X, Y)
g = Hom(:g, Y, Z)
@test to_json_sexpr(X) == ["Ob", "X"]
@test to_json_sexpr(f) == ["Hom", "f", ["Ob", "X"], ["Ob", "Y"]]
@test to_json_sexpr(compose(f,g)) == [
  "compose",
  ["Hom", "f", ["Ob", "X"], ["Ob", "Y"]],
  ["Hom", "g", ["Ob", "Y"], ["Ob", "Z"]],
]
@test to_json_sexpr(f, by_reference=x->true) == "f"
@test to_json_sexpr(compose(f,g), by_reference=x->true) == ["compose","f","g"]

# From JSON
@test parse_json_sexpr(FreeMonoid, [:Elem, "x"]) == Elem(FreeMonoid, :x)
@test parse_json_sexpr(FreeMonoid, [:munit]) == munit(FreeMonoid.Elem)
@test parse_json_sexpr(FreeCategory, ["Ob", "X"]) == X
@test parse_json_sexpr(FreeCategory, ["Ob", "X"]; symbols=false) ==
  Ob(FreeCategory.Ob, "X")
@test parse_json_sexpr(FreeCategory, ["Hom", "f", ["Ob", "X"], ["Ob", "Y"]]) == f
@test parse_json_sexpr(FreeCategory, ["Hom", "f", ["Ob", "X"], ["Ob", "Y"]]; symbols=false) ==
  Hom("f", Ob(FreeCategory.Ob, "X"), Ob(FreeCategory.Ob, "Y"))

# Round trip
@test parse_json_sexpr(FreeCategory, to_json_sexpr(compose(f,g))) == compose(f,g)
@test parse_json_sexpr(FreeCategory, to_json_sexpr(id(X))) == id(X)

end


text:
@safetestset facecontinuous_smoke = "FaceContinuous basic test" begin
# Test using the vector form of the Hilbert index so that we don't
# introduce another layer of error. When this works, delete it.
using BijectiveHilbert: FaceContinuous, H_encode!, H_decode!, msb
# Works for UInt32 or UInt16.
A = UInt32
H = UInt32
fc1 = FaceContinuous(H, 9, 3)
fc2 = FaceContinuous(H, 8, 3)
X1 = zeros(A, 3)
X1[1] = A(0x8f)
X1[2] = A(0x8f)
X1[3] = A(0x8f)
hvec1 = zeros(A, 3)
hvec2 = zeros(A, 3)
H_encode!(fc1, X1, hvec1)
H_encode!(fc2, X1, hvec2)
@assert hvec1 == hvec2
X2 = zeros(A, 3)
H_decode!(fc1, hvec1, X1)
H_decode!(fc2, hvec2, X2)
@test X1 == X2
end

@safetestset facecontinuous_api = "FaceContinuous api test" begin
# Test using the vector form of the Hilbert index so that we don't
# introduce another layer of error. When this works, delete it.
using BijectiveHilbert: FaceContinuous, encode_hilbert_zero, decode_hilbert_zero!
# Works for UInt32 or UInt16.
A = UInt32
H = UInt32
fc = FaceContinuous(H, 9, 3)
X = zeros(A, 3)
X[1] = A(0x8f)
X[2] = A(0x85)
X[3] = A(0x43)
h = encode_hilbert_zero(fc, X)
X2 = zeros(A, 3)
decode_hilbert_zero!(fc, X2, h)
@test X == X2
X[1] = A(0x5)
X[2] = A(0x8)
X[3] = A(0x3)
h = encode_hilbert_zero(fc, X)
X3 = zeros(A, 3)
decode_hilbert_zero!(fc, X3, h)
@test X == X3
end

@safetestset facecontinuous_own_inverse = "FaceContinuous is its own inverse" begin
using BijectiveHilbert: FaceContinuous
using ..HilbertTestSuite: check_own_inverse
for n in [2, 3, 5]
    for b in [2, 3]
        fc = FaceContinuous(b, n)
        @test check_own_inverse(fc, b, n)
    end
end
end

@safetestset facecontinuous_complete_set = "FaceContinuous is complete set" begin
using BijectiveHilbert: FaceContinuous
using ..HilbertTestSuite: check_complete_set
for n in [2, 3, 5]
    for b in [2, 3]
        fc = FaceContinuous(b, n)
        @test check_complete_set(fc, b, n)
    end
end
end


text:
export run_scunittfopf, run_ac_scunittfopf, run_dc_scunittfopf

""
function run_ac_scunittfopf(data::Dict{String,Any}, solver; kwargs...)
    return run_scunittfopf(data, ACPPowerModel, solver; kwargs...)
end

""
function run_dc_scunittfopf(data::Dict{String,Any}, solver; kwargs...)
    return run_scunittfopf(data, DCPPowerModel, solver; kwargs...)
end

""
function run_scunittfopf(data::Dict{String,Any}, model_constructor, solver; kwargs...)
    process_additional_data!(data)
    pm = PowerModels.build_generic_model(data, model_constructor, post_scunittfopf; kwargs...)
    return PowerModels.solve_generic_model(pm, solver; solution_builder = get_solution_tf)
end

""
function post_scunittfopf(pm::GenericPowerModel)
    variable_risk(pm)
    variable_dispatch_cost(pm)
    variable_redispatch_cost(pm)
    variable_loadshedding_cost(pm)
    first_stage_network_id = 1;
    second_stage_network_ids = Dict();
    for (n, network) in pm.ref[:nw]
        if n != 1
            push!(second_stage_network_ids, n => n);
        end
    end
    first_stage_model(pm, first_stage_network_id)
    second_stage_model(pm, first_stage_network_id, second_stage_network_ids)
    objective_total_risk(pm, first_stage_network_id, second_stage_network_ids)
end

function first_stage_model(pm::GenericPowerModel, first_stage_network_id)
    n = first_stage_network_id;
    add_load_model!(pm, n) # To add load data
    add_power_factor!(pm, n) # To add load data
    PowerModels.variable_voltage(pm; nw = n)
    PowerModels.variable_generation(pm;  nw = n)
    PowerModels.variable_branch_flow(pm;  nw = n)
    PowerModels.variable_dcline_flow(pm;  nw = n)
    variable_transformation(pm;  nw = n)
    #variable_dispatch_cost(pm, n)
    variable_node_aggregation(pm;  nw = n)
    variable_load(pm;  nw = n)
    variable_action_indicator(pm;  nw = n)
    variable_auxiliary_power(pm;  nw = n)

    PowerModels.constraint_voltage(pm;  nw = n)
    for i in PowerModels.ids(pm, n, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i;  nw = n)
    end

    for i in PowerModels.ids(pm, n, :bus)
        constraint_kcl_shunt_aggregated(pm, i;  nw = n)
        constraint_load_gen_aggregation_sheddable(pm, i; nw = n)
    end

    for i in PowerModels.ids(pm, n, :load)
        constraint_fixed_load(pm, i; nw = n)
    end

    for i in PowerModels.ids(pm, n, :gen)
        constraint_flexible_gen(pm, i; nw = n)
        constraint_redispatch_power_gen(pm, i; nw = n)
    end

    for i in PowerModels.ids(pm, n, :branch)
        branch = PowerModels.ref(pm, n, :branch, i)

        if branch["shiftable"] == false && branch["tappable"] == false
            PowerModels.constraint_ohms_yt_from(pm, i; nw = n)
            PowerModels.constraint_ohms_yt_to(pm, i; nw = n)
            constraint_link_voltage_magnitudes(pm, i; nw = n )
        else
            constraint_variable_transformer_y_from(pm, i; nw = n)
            constraint_variable_transformer_y_to(pm, i; nw = n)
        end
        PowerModels.constraint_voltage_angle_difference(pm, i; nw = n)

        PowerModels.constraint_thermal_limit_from(pm, i; nw = n)
        PowerModels.constraint_thermal_limit_to(pm, i; nw = n)
    end
    for i in PowerModels.ids(pm, n, :dcline)
        PowerModels.constraint_dcline(pm, i; nw = n)
    end
end

function second_stage_model(pm::GenericPowerModel, first_stage_network_id, second_stage_network_ids)
    for (n, contingency) in second_stage_network_ids
        add_load_model!(pm, n) # To add load data
        add_power_factor!(pm, n) # To add load data
        PowerModels.variable_voltage(pm; nw = n)
        PowerModels.variable_generation(pm; nw = n)
        PowerModels.variable_branch_flow(pm; nw = n)
        PowerModels.variable_dcline_flow(pm; nw = n)
        variable_transformation(pm; nw = n)
        #variable_dispatch_cost(pm, n)
        variable_node_aggregation(pm; nw = n)
        variable_load(pm; nw = n)
        variable_action_indicator(pm; nw = n)
        variable_auxiliary_power(pm; nw = n)

        PowerModels.constraint_voltage(pm; nw = n)
        for i in PowerModels.ids(pm, n, :ref_buses)
            PowerModels.constraint_theta_ref(pm, i; nw = n)
        end

        for i in PowerModels.ids(pm, n, :bus)
            constraint_kcl_shunt_aggregated(pm, i;  nw = n)
            constraint_load_gen_aggregation_sheddable(pm, i;  nw = n)
        end

        for i in PowerModels.ids(pm, n, :load)
            constraint_flexible_load(pm, i;  nw = n)
            constraint_second_stage_redispatch_power_load(pm, n, i, first_stage_network_id)
        end

        contingencies = PowerModels.ref(pm, n, :contingencies)
        contingency_id = second_stage_network_ids[n]
        for i in PowerModels.ids(pm, n, :gen)
            if contingencies[contingency_id]["gen_id1"] == i || contingencies[contingency_id]["gen_id2"] == i || contingencies[contingency_id]["gen_id3"] == i
                constraint_gen_contingency(pm, n, i)
            else
                constraint_flexible_gen(pm, i;  nw = n)
                constraint_second_stage_redispatch_power_gen(pm, n, i, first_stage_network_id)
            end
        end

        for i in PowerModels.ids(pm, n, :branch)
            branch = PowerModels.ref(pm, n, :branch, i)
            if branch["shiftable"] == false && branch["tappable"] == false
                if contingencies[contingency_id]["branch_id1"] == i || contingencies[contingency_id]["branch_id2"] == i || contingencies[contingency_id]["branch_id3"] == i
                    constraint_branch_contingency(pm, n, i)
                else
                    PowerModels.constraint_ohms_yt_from(pm, i; nw = n)
                    PowerModels.constraint_ohms_yt_to(pm, i; nw = n)
                    constraint_link_voltage_magnitudes(pm, i; nw = n)
                end
            else
                if contingencies[contingency_id]["branch_id1"] == i || contingencies[contingency_id]["branch_id2"] == i || contingencies[contingency_id]["branch_id3"] == i
                    constraint_branch_contingency(pm, n, i)
                else
                    constraint_variable_transformer_y_from(pm, i, nw = n)
                    constraint_variable_transformer_y_to(pm, i, nw = n)
                end
            end
            PowerModels.constraint_voltage_angle_difference(pm, i; nw = n)

            PowerModels.constraint_thermal_limit_from(pm, i; nw = n)
            PowerModels.constraint_thermal_limit_to(pm, i; nw = n)
        end
        for i in PowerModels.ids(pm, n, :dcline)
            PowerModels.constraint_dcline(pm, i; nw = n)
        end
    end
end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
module EMpht

using Dates
using LinearAlgebra
using Random
using Statistics

export empht
# export loglikelihoodcensored
# export parse_settings

include("sample.jl")
include("phasetype.jl")
include("utils.jl")
include("ode.jl")
include("uniformization.jl")

function em_iterate(name, s, fit, ph_structure, method, max_iter, timeout,
        verbose)
    p = fit.p

    # Count the total of all weight.
    sumOfWeights = sum(s.obsweight) + sum(s.censweight) + sum(s.intweight)

    # Sort the observations.
    perm = sortperm(s.obs)
    obsSort = s.obs[perm]
    obsweightSort = s.obsweight[perm]
    s = Sample(obsSort, obsweightSort, s.cens, s.censweight, s.int, s.intweight)

    start = now()
    endTime = start+Dates.Minute(round(timeout))

    save_progress(name, s, fit, start)

    ll = 0

    numPlots = 0
    iter = 1
    while iter < max_iter
        ##  The expectation step!
        Bs = zeros(p); Zs = zeros(p); Ns = zeros(p, p+1)

        if length(s.obs) > 0
            if method == :unif
                e_step_observed_uniform!(s, fit, Bs, Zs, Ns)
            elseif method == :ode
                e_step_observed_ode!(s, fit, Bs, Zs, Ns)
            else
                error("Method should be :unif or :ode")
            end
        end

        if length(s.cens) > 0 || length(s.int) > 0
            if method == :unif
                e_step_censored_uniform!(s, fit, Bs, Zs, Ns)
            elseif method == :ode
                e_step_censored_ode!(s, fit, Bs, Zs, Ns)
            else
                error("Method should be :unif or :ode")
            end
        end

        ## The maximisation step!
        π_next = max.(Bs ./ sumOfWeights, 0)
        t_next = max.(Ns[:,end] ./ Zs, 0)
        t_next[isnan.(t_next)] .= 0

        T_next = zeros(p,p)
        for i=1:p
            T_next[i,:] = max.(Ns[i,1:end-1] ./ Zs[i], 0)
            T_next[i,isnan.(T_next[i,:])] .= 0
            T_next[i,i] = -(t_next[i] + sum(T_next[i,:]))
        end

        # Remove any numerical instabilities.
        π_next = max.(π_next, 0)
        π_next /= sum(π_next)

        fit = PhaseType(π_next, T_next, t_next)

        # The CanonicalForm1 phase-type is a bit special. It needs to have an
        # ordering to its diagonal of T, and the EM steps will probably break
        # this, so it needs to be sorted on every iteration.
        if ph_structure == "CanonicalForm1"
            fit = sort_into_canonical_form(fit)
        end

        remaining = endTime - now()
        if remaining.value < 0 || iter >= max_iter
            ll = save_progress(name, s, fit, start)
            seconds = round(now() - start, Dates.Second)
            duration = Dates.canonicalize(Dates.CompoundPeriod(seconds))
            if verbose
                println("Quitting after $duration/$iter iterations with ll = ",
                    round(Float64(ll), sigdigits=5))
            end
            break
        else
            if verbose
                ll = loglikelihoodcensored(s, fit)
                seconds = round(remaining, Dates.Second)
                duration = Dates.canonicalize(Dates.CompoundPeriod(seconds))
                println("Iteration: $iter/$max_iter\t",
                    "Timeout in $duration\t",
                    "Loglikelihood: ", round(Float64(ll), sigdigits=5))
            end
        end

        iter += 1
    end

    fit
end

function empht(s; p=1, ph_structure="Coxian", name="", continueFit=false,
        method=:unif, max_iter=100, timeout=1, verbose=false)

    if verbose
        println("(p,ph_structure,name,method,continueFit,max_iter,timeout) = ",
            ((p, ph_structure, name, method, continueFit, max_iter, timeout)))
    end

    # Check we don't just have right-censored obs, since this blows things up.
    if length(s.obs) == 0 && length(s.cens) > 0 && length(s.int) == 0
        error("Can't just have right-censored observations!")
    end

    # If not continuing previous fit, remove any left-over output files.
    if ~continueFit
        rm(string(name, "_loglikelihood.csv"), force=true)
        rm(string(name, "_fit.csv"), force=true)
    end

    fit = initial_phasetype(name, p, ph_structure, continueFit, s, verbose)

    if verbose && p <= 10
        println("first pi is $(fit.π), first T is $(fit.T)\n")
    end

    em_iterate(name, s, fit, ph_structure, method, max_iter, timeout, verbose)
end

function empht(settings_filename::String; method=:unif, verbose=false)
    # Read in details for the fit from the settings file.
    s, p, ph_structure, name, continueFit, max_iter, timeout =
        parse_settings(settings_filename, verbose)
    empht(s, p=p, ph_structure=ph_structure, name=name, continueFit=continueFit,
            method=method, max_iter=max_iter, timeout=timeout, verbose=verbose)
end

end # module


text:
using ModelsGenerators
srand(1)
y,x,z=randiv(1000)
xb = randn(1000,2)
x = [x xb]
z = [z xb]
n, m = size(z)

f(theta) = z.*(y-x*theta);

function f3(theta)
  g = Array(Float64, n, m)
  broadcast!(*, g, z, y-x*theta)
end

@time f3([.1, 0, .1])
@time f([.1, 0, .1])
@time f_p([.1, 0, .1])

function f_p(theta)
    n, m = size(z)
    v = Array(Float64, n, m)
    u = (y-x*theta)
    @inbounds for k = 1:m
        for j = 1:n
            v[j,k] = z[j,k]*u[j]
        end
    end
    v
end

function f_p!(v, theta)
    n, m = size(z)
    u = (y-x*theta)
    @inbounds for k = 1:m
        for j = 1:n
            v[j,k] = z[j,k]*u[j]
        end
    end
    v
end

n, m = size(z)
v = Array(Float64, n, m)

f([.1, 0, 0])
f_p([.1, 0, 0])
f_p!(v, [.1, 0, 0])
@time f([.1, 0, 0])
@time f_p([.1, 0, 0])
@time f_p!(v, [.1, 0, 0])

p = ones(10000);
l = ones(7,1);

ws() = p'v*l

@time ws()

using Base.LinAlg.BLAS

using FastAnonymous

af = @anon theta->z.*(y-x*theta)

af([.1,0,0])

type dummy
    a::DataType
end

uu = dummy(af)

type dummy_
    a::Function
end

uuu = dummy_(f)

fdummy() = for j = 1:1000 f([.1,0,0]) end
fdummy_2() = for j = 1:1000 uu.a([.1,0,0]) end
fdummy_() = for j = 1:1000 uuu.a([.1,0,0]) end
@time fdummy()
@time fdummy_2()
@time fdummy_()


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

"""
Monte Carlo search by drawing uniform samples from the root of the grammar.
Returns the sample with the best result.
"""
module MC

export MCESParams, MCESResult, mc_search, exprsearch, SearchParams, SearchResult, get_derivtree

using Reexport
using ExprSearch
using RLESUtils, GitUtils, CPUTimeUtils, Observers, LogSystems
@reexport using LinearDerivTrees  #for pretty strings
using GrammaticalEvolution
using Iterators
using RLESTypes.SymbolTable

import ..ExprSearch: exprsearch, get_fitness, get_derivtree, get_expr
import Base: isless, copy!

const DEFAULT_EXPR = :()

include("logdefs.jl")

type MCESParams <: SearchParams
    #tree params
    maxsteps::Int64

    #MC
    n_samples::Int64 #samples
    default_expr        #if derivation is incompelete, use default expr
    logsys::LogSystem
    userargs::SymbolTable
end
MCESParams(maxsteps::Int64, n_samples::Int64, logsys::LogSystem=logsystem();
    userargs::SymbolTable=SymbolTable()) =
    MCESParams(maxsteps, n_samples, DEFAULT_EXPR, logsys, userargs)

type MCState
    tree::LinearDerivTree
    fitness::Float64
    expr
end

MCState(tree::LinearDerivTree) = MCState(tree, realmax(Float64), 0)

type MCESResult <: SearchResult
    tree::LinearDerivTree
    actions::Vector{Int64}
    fitness::Float64
    expr
    best_at_eval::Int64
    totalevals::Int64

    function MCESResult()
        result = new()
        result.actions = Int64[]
        result.fitness = realmax(Float64)
        result.expr = 0
        result.best_at_eval = 0
        result.totalevals = 0
        result
    end

    function MCESResult(tree::LinearDerivTree)
        result = MCESResult()
        result.tree = tree
        result
    end
end

exprsearch(p::MCESParams, problem::ExprProblem) = mc_search(p, problem)

get_derivtree(result::MCESResult) = get_derivtree(result.tree)
get_expr(result::MCESResult) = result.expr
get_fitness(result::MCESResult) = result.fitness

function mc_search(p::MCESParams, problem::ExprProblem)
    @notify_observer(p.logsys.observer, "verbose1", ["Starting MC search"])
    @notify_observer(p.logsys.observer, "computeinfo", ["starttime", string(now())])

    initialize!(problem)
    grammar = get_grammar(problem)
    tree_params = LDTParams(grammar, p.maxsteps)

    s = MCState(LinearDerivTree(tree_params))
    result = MCESResult(LinearDerivTree(tree_params))

    tstart = CPUtime_start()
    for i = 1:p.n_samples
        @notify_observer(p.logsys.observer, "iteration", [i])
        ###############
        #MC algorithm
        sample!(s, p, problem)
        update!(result, s)
        ###############

        @notify_observer(p.logsys.observer, "elapsed_cpu_s", [i, CPUtime_elapsed_s(tstart)])
        @notify_observer(p.logsys.observer, "current_best", [i, result.fitness, string(result.expr)])
    end
    @notify_observer(p.logsys.observer, "result", [result.fitness, string(result.expr),
        result.best_at_eval, result.totalevals])

    #meta info
    @notify_observer(p.logsys.observer, "computeinfo", ["endtime",  string(now())])
    @notify_observer(p.logsys.observer, "computeinfo", ["hostname", gethostname()])
    @notify_observer(p.logsys.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
    @notify_observer(p.logsys.observer, "computeinfo", ["cpu_time", CPUtime_elapsed_s(tstart)])
    @notify_observer(p.logsys.observer, "parameters", ["maxsteps", p.maxsteps])
    @notify_observer(p.logsys.observer, "parameters", ["n_samples", p.n_samples])

    result
end

#initialize to random state
function sample!(s::MCState, p::MCESParams, problem::ExprProblem)
    rand!(s.tree) #sample uniformly
    try
        s.expr = get_expr(s.tree)
        s.fitness = get_fitness(problem, get_derivtree(s.tree), p.userargs)
    catch e
        if !isa(e, IncompleteException)
            rethrow(e)
        end
        s.expr = p.default_expr
        s.fitness = realmax(Float64)
    end
    s
end

#update the global best trackers with the current state
function update!(result::MCESResult, s::MCState)
    result.totalevals += 1 #assumes an eval was called prior

    #update globals
    if s.fitness < result.fitness
        copy!(result.tree, s.tree)
        resize!(result.actions, length(s.tree.actions))
        result.actions = convert(Array, s.tree.actions)
        result.fitness = s.fitness
        result.expr = s.expr
        result.best_at_eval = result.totalevals
    end
end

isless(r1::MCESResult, r2::MCESResult) = r1.fitness < r2.fitness

function copy!(dst::MCState, src::MCState)
    copy!(dst.tree, src.tree)
    dst.fitness = src.fitness
    dst.expr = src.expr
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
function extract_docstrings(fPaths, outFile::ASCIIString)
    lnsToWrite = (ByteString)[]
    for fn=1:length(fPaths)
        fPath = fPaths[fn]
        fHandle = open(fPath, "r")
        lns = readlines(fHandle)
        close(fHandle)

        startPnts = (Int)[]
        stopPnts = (Int)[]

        for i=1:length(lns)
            if strip(lns[i]) == "@doc doc\"\"\""
                push!(startPnts, i)
            elseif strip(lns[i]) == "\"\"\"->"
                push!(stopPnts, i)
            end
        end

        for s=1:length(startPnts)
            for nl=1:stopPnts[s]-startPnts[s]-1
                push!(lnsToWrite, lns[startPnts[s]+nl])
            end
            push!(lnsToWrite, "\n")
        end
    end

    fOut = open(outFile, "w")
    write(fOut, lnsToWrite)
    close(fOut)
end


text:
const XREGISTERS = 2

struct Xoshift{P} <: AbstractVRNG{P}
    ptr::Ptr{UInt64}
end
@inline Base.pointer(rng::Xoshift) = rng.ptr

struct XoshiftState{P,W} <: AbstractState{P,W}
    eins::VecUnroll{P,W,UInt64,Vec{W,UInt64}}
    zwei::VecUnroll{P,W,UInt64,Vec{W,UInt64}}
    drei::VecUnroll{P,W,UInt64,Vec{W,UInt64}}
    vier::VecUnroll{P,W,UInt64,Vec{W,UInt64}}
end

Xoshift(ptr, ::StaticInt{X}) where {X} = Xoshift{X}(ptr)
Xoshift(ptr) = Xoshift(ptr, pick_vector_width(UInt64) * StaticInt(XREGISTERS))
function randnonzero()
    while true
        r = rand(UInt64)
        iszero(r) || return r
    end
end
function initXoshift!(ptr::Ptr{UInt64}, P)
    e = randnonzero(); z = randnonzero();
    d = randnonzero(); v = randnonzero();
    initXoshift!(ptr, P, e, z, d, v)
end
function initXoshift!(ptr::Ptr{UInt64}, P, e::UInt64, z::UInt64, d::UInt64, v::UInt64) # P here means number of streams
    for j ∈ 1:P-1
        i = P - j
        vstoreu!(ptr, e, 8i); vstoreu!(ptr, z, 8*(i + P)); vstoreu!(ptr, d, 8*(i + 2P)); vstoreu!(ptr, v, 8*(i + 3P));
        e, z, d, v = jump(e, z, d, v)
    end
    vstoreu!(ptr, e); vstoreu!(ptr, z, 8P); vstoreu!(ptr, d, 8*(2P)); vstoreu!(ptr, v, 8*(3P));
    vstoreu!(Base.unsafe_convert(Ptr{UInt32}, ptr), 0x00000000, 8*(4P));
end
function jump(eins, zwei, drei, vier)
    e = zero(UInt64); z = zero(UInt64); d = zero(UInt64); v = zero(UInt64)
    for u ∈ (0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c)
        for b ∈ 0:63
            if u % Bool
                e ⊻= eins
                z ⊻= zwei
                d ⊻= drei
                v ⊻= vier
            end
            u >>>= 1
            e, z, d, v = nextstate(e, z, d, v)
        end
    end
    e, z, d, v
end
function seed!(seed::Integer)
    i = seed % UInt64
    e = z = d = v = zero(UInt64)
    increment = 0xa04de531e612e1b9
    while any(iszero, (e, z, d, v))
        e = ((i * 0x90ce6ecbad5e33b5) + increment)
        z = ((e * 0x90ce6ecbad5e33b5) + increment)
        d = ((z * 0x90ce6ecbad5e33b5) + increment)
        v = ((d * 0x90ce6ecbad5e33b5) + increment)
        increment += 0x0000000000000002
    end
    nstreams = XREGISTERS * Base.Threads.nthreads() * pick_vector_width(UInt64)
    initXoshift!(GLOBAL_vRNGs[], nstreams, e, z, d, v)
end

# @inline function getstate(rng::Xoshift{P}, ::Val{N}, ::Val{W}) where {P,N,W}
#     ptr = pointer(rng)
#     XoshiftState(
#         ntuple(n -> vloada(Vec{W,UInt64}, ptr, W64*(n - 1)), Val{N}()),
#         ntuple(n -> vloada(Vec{W,UInt64}, ptr, W64*(n - 1) + W64*P), Val{N}()),
#         ntuple(n -> vloada(Vec{W,UInt64}, ptr, W64*(n - 1) + 2W64*P), Val{N}()),
#         ntuple(n -> vloada(Vec{W,UInt64}, ptr, W64*(n - 1) + 3W64*P), Val{N}())
#     )
# end
@inline function getrandu64counter(rng::Xoshift{P}) where {P}
    vloadu(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), 4simd_integer_register_size()*P)
end
@inline function getrand64counter(rng::Xoshift{P}) where {P}
    vloadu(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), 4simd_integer_register_size()*P + 2)
end
@inline function getrandn64counter(rng::Xoshift{P}) where {P}
    vloadu(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), 4simd_integer_register_size()*P + 3)
end
@inline function setrandu64counter!(rng::Xoshift{P}, v::UInt8) where {P}
    vstoreu!(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), v, 4simd_integer_register_size()*P)
end
@inline function setrand64counter!(rng::Xoshift{P}, v::UInt8) where {P}
    vstoreu!(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), v, 4simd_integer_register_size()*P + 2)
end
@inline function setrandn64counter!(rng::Xoshift{P}, v::UInt8) where {P}
    vstoreu!(Base.unsafe_convert(Ptr{UInt8}, pointer(rng)), v, 4simd_integer_register_size()*P + 3)
end

@inline function getstate(rng::Xoshift{P}, ::Val{1}, ::StaticInt{W}) where {P,W}
    ptr = pointer(rng)
    XoshiftState(
        VecUnroll((vloada(ptr, MM{W,8}(StaticInt{0}())),)),
        VecUnroll((vloada(ptr, MM{W,8}(simd_integer_register_size()*P)),)),
        VecUnroll((vloada(ptr, MM{W,8}(2simd_integer_register_size()*P)),)),
        VecUnroll((vloada(ptr, MM{W,8}(3simd_integer_register_size()*P)),))
    )
end
@inline function getstate(rng::Xoshift{P}, ::Val{2}, ::StaticInt{W}) where {P,W}
    ptr = pointer(rng)
    XoshiftState(
        VecUnroll((vloada(ptr, MM{W,8}(StaticInt{0}()  )), vloada(ptr, MM{W,8}(simd_integer_register_size()         )))),
        VecUnroll((vloada(ptr, MM{W,8}( P*simd_integer_register_size())), vloada(ptr, MM{W,8}(simd_integer_register_size()*(1 +  P))))),
        VecUnroll((vloada(ptr, MM{W,8}(2P*simd_integer_register_size())), vloada(ptr, MM{W,8}(simd_integer_register_size()*(1 + 2P))))),
        VecUnroll((vloada(ptr, MM{W,8}(3P*simd_integer_register_size())), vloada(ptr, MM{W,8}(simd_integer_register_size()*(1 + 3P)))))
    )
end
@inline function getstate(rng::Xoshift{P}, ::Val{4}, ::StaticInt{W}) where {P,W}
    ptr = pointer(rng)
    RS = simd_integer_register_size()
    XoshiftState(
        VecUnroll((vloada(ptr, MM{W,8}(StaticInt{0}())), vloada(ptr, MM{W,8}(RS)         ), vloada(ptr, MM{W,8}(RS* 2)      ), vloada(ptr, MM{W,8}(RS* 3      )))),
        VecUnroll((vloada(ptr, MM{W,8}( P*RS)),          vloada(ptr, MM{W,8}(RS*(1 +  P))), vloada(ptr, MM{W,8}(RS*(2 +  P))), vloada(ptr, MM{W,8}(RS*(3 +  P))))),
        VecUnroll((vloada(ptr, MM{W,8}(2P*RS)),          vloada(ptr, MM{W,8}(RS*(1 + 2P))), vloada(ptr, MM{W,8}(RS*(2 + 2P))), vloada(ptr, MM{W,8}(RS*(3 + 2P))))),
        VecUnroll((vloada(ptr, MM{W,8}(3P*RS)),          vloada(ptr, MM{W,8}(RS*(1 + 3P))), vloada(ptr, MM{W,8}(RS*(2 + 3P))), vloada(ptr, MM{W,8}(RS*(3 + 3P)))))
    )
end
@inline function storestate!(rng::Xoshift{P}, s::XoshiftState{N,W}) where {P,N,W}
    ptr = pointer(rng)
    @unpack eins, zwei, drei, vier = s
    @inbounds for n ∈ 0:N
        vstorea!(rng, data(eins)[n], simd_integer_register_size()*n)
    end
    @inbounds for n ∈ 0:N
        vstorea!(rng, data(zwei)[n], simd_integer_register_size()*(n +  P))
    end
    @inbounds for n ∈ 0:N
        vstorea!(rng, data(drei)[n], simd_integer_register_size()*(n + 2P))
    end
    @inbounds for n ∈ 0:N
        vstorea!(rng, data(vier)[n], simd_integer_register_size()*(n + 3P))
    end
end
@inline function storestate!(rng::Xoshift{P}, s::XoshiftState{0,W}) where {P,W}
    ptr = pointer(rng)
    @unpack eins, zwei, drei, vier = s;
    _eins = data(eins); _zwei = data(zwei); _drei = data(drei); _vier = data(vier);
    @inbounds begin
        vstorea!(ptr, _eins[1],       )
        vstorea!(ptr, _zwei[1],  P*simd_integer_register_size())
        vstorea!(ptr, _drei[1], 2P*simd_integer_register_size())
        vstorea!(ptr, _vier[1], 3P*simd_integer_register_size())
    end
end
@inline function storestate!(rng::Xoshift{P}, s::XoshiftState{1,W}) where {P,W}
    ptr = pointer(rng)
    @unpack eins, zwei, drei, vier = s;
    _eins = data(eins); _zwei = data(zwei); _drei = data(drei); _vier = data(vier);
    @inbounds begin
        vstorea!(ptr, _eins[1],                            )
        vstorea!(ptr, _eins[2],   simd_integer_register_size()          )
        vstorea!(ptr, _zwei[1],   simd_integer_register_size()*       P )
        vstorea!(ptr, _zwei[2],   simd_integer_register_size()*(1 +   P))
        vstorea!(ptr, _drei[1],   simd_integer_register_size()*      2P )
        vstorea!(ptr, _drei[2],   simd_integer_register_size()*(1 +  2P))
        vstorea!(ptr, _vier[1],   simd_integer_register_size()*      3P )
        vstorea!(ptr, _vier[2],   simd_integer_register_size()*(1 +  3P))
    end
end
@inline function storestate!(rng::Xoshift{P}, s::XoshiftState{3,W}) where {P,W}
    ptr = pointer(rng)
    @unpack eins, zwei, drei, vier = s;
    _eins = data(eins); _zwei = data(zwei); _drei = data(drei); _vier = data(vier);
    @inbounds begin
        vstorea!(ptr, _eins[1],      )
        vstorea!(ptr, _eins[2],   simd_integer_register_size())
        vstorea!(ptr, _eins[3],   simd_integer_register_size()*2)
        vstorea!(ptr, _eins[4],   simd_integer_register_size()*3)
        vstorea!(ptr, _zwei[1],   simd_integer_register_size()*       P)
        vstorea!(ptr, _zwei[2],   simd_integer_register_size()*(1 +   P))
        vstorea!(ptr, _zwei[3],   simd_integer_register_size()*(2 +   P))
        vstorea!(ptr, _zwei[4],   simd_integer_register_size()*(3 +   P))
        vstorea!(ptr, _drei[1],   simd_integer_register_size()*      2P)
        vstorea!(ptr, _drei[2],   simd_integer_register_size()*(1 +  2P))
        vstorea!(ptr, _drei[3],   simd_integer_register_size()*(2 +  2P))
        vstorea!(ptr, _drei[4],   simd_integer_register_size()*(3 +  2P))
        vstorea!(ptr, _vier[1],   simd_integer_register_size()*      3P)
        vstorea!(ptr, _vier[2],   simd_integer_register_size()*(1 +  3P))
        vstorea!(ptr, _vier[3],   simd_integer_register_size()*(2 +  3P))
        vstorea!(ptr, _vier[4],   simd_integer_register_size()*(3 +  3P))
    end
end

@inline function nextstate(eins, zwei, drei, vier)
    t = zwei << 0x0000000000000011
    drei = drei ⊻ eins
    vier = vier ⊻ zwei
    zwei = zwei ⊻ drei
    eins = eins ⊻ vier
    drei = drei ⊻ t
    vier = rotate_right(vier, 0x000000000000002d)
    eins, zwei, drei, vier
end

@generated function _unpack(s::XoshiftState{P}, ::Val{U}) where {P,U}
    if U ≤ P
        quote
            $(Expr(:meta,:inline))
            @unpack eins, zwei, drei, vier = s
            _eins = data(eins); _zwei = data(zwei); _drei = data(drei); _vier = data(vier);
            (
                VecUnroll( Base.Cartesian.@ntuple $U u -> _eins[u] ),
                VecUnroll( Base.Cartesian.@ntuple $U u -> _zwei[u] ),
                VecUnroll( Base.Cartesian.@ntuple $U u -> _drei[u] ),
                VecUnroll( Base.Cartesian.@ntuple $U u -> _vier[u] )
            )
        end
    elseif U == P+1
        quote
            $(Expr(:meta,:inline))
            s.eins, s.zwei, s.drei, s.vier
        end
    else
        throw("$U > $(P+1).")
    end
end

@inline XoshiftState(eins::VecUnroll{N}, zwei::VecUnroll{N}, drei::VecUnroll{N}, vier::VecUnroll{N}, s::XoshiftState{N}) where {N} = XoshiftState( eins, zwei, drei, vier )
@generated function XoshiftState(eins::VecUnroll{N}, zwei::VecUnroll{N}, drei::VecUnroll{N}, vier::VecUnroll{N}, s::XoshiftState{P}) where {N,P}
    @assert P > N
    q = quote
        $(Expr(:meta,:inline))
        e = data(eins); z = data(zwei); d = data(drei); v = data(vier)
        _e = data(s.eins); _z = data(s.zwei); _d = data(s.drei); _v = data(s.vier)
    end
    _eins = Expr(:tuple)
    _zwei = Expr(:tuple)
    _drei = Expr(:tuple)
    _vier = Expr(:tuple)
    for n ∈ 1:N+1
        push!(_eins.args, Expr(:ref, :e, n))
        push!(_zwei.args, Expr(:ref, :z, n))
        push!(_drei.args, Expr(:ref, :d, n))
        push!(_vier.args, Expr(:ref, :v, n))
    end
    for n ∈ N+2:P+1
        push!(_eins.args, Expr(:ref, :_e, n))
        push!(_zwei.args, Expr(:ref, :_z, n))
        push!(_drei.args, Expr(:ref, :_d, n))
        push!(_vier.args, Expr(:ref, :_v, n))
    end
    push!(q.args, :(XoshiftState( VecUnroll($_eins), VecUnroll($_zwei), VecUnroll($_drei), VecUnroll($_vier) )))
    q
end

@inline function nextstate(s::XoshiftState{P}, ::Val{U}) where {P,U}
    eins, zwei, drei, vier = _unpack(s, Val{U}())
    out = eins + vier
    out = rotate_right(out, 0x0000000000000017)
    eins, zwei, drei, vier = nextstate(eins, zwei, drei, vier)
    # out += eins
    XoshiftState( eins, zwei, drei, vier, s ), out
end

function randbuffer64(r::Xoshift{P}) where {P}
    ptr = pointer(r)
    Buffer256(Base.unsafe_convert(Ptr{Float64}, ptr + P * 5simd_integer_register_size()))
end
function randnbuffer64(r::Xoshift{P}) where {P}
    ptr = pointer(r)
    Buffer256(Base.unsafe_convert(Ptr{Float64}, ptr + P * 5simd_integer_register_size() + 2048))
end
# function randbuffer32(r::Xoshift{P}) where {P}
#     ptr = pointer(r)
#     Buffer256(Base.unsafe_convert(Ptr{Float64}, ptr + P * 5simd_integer_register_size() + 4096))
# end
function randubuffer64(r::Xoshift{P}) where {P}
    ptr = pointer(r)
    Buffer256(Base.unsafe_convert(Ptr{UInt64}, ptr + P * 5simd_integer_register_size() + 4096))
end
# function randnbuffer32(r::Xoshift{P}) where {P}
#     ptr = pointer(r)
#     Buffer256(Base.unsafe_convert(Ptr{Float64}, ptr + P * 5simd_integer_register_size() + 5120))
# end


text:
const UA_TYPES = Ref{Ptr{UA_DataType}}(0) # Initilize with C_NULL and initialize correct address during __init__
const UA_TYPES_PTRS = OffsetVector{Ptr{UA_DataType}}(undef, 0:UA_TYPES_COUNT-1) # Initilize vector of UA_TYPES pointer undefined and write values during __init__
const UA_TYPES_MAP = Vector{DataType}(undef, UA_TYPES_COUNT) # Initilize vector of mapping between UA_TYPES and Julia types as undefined and write values during __init__

function juliadatatype(p::Ptr{UA_DataType})
    global UA_TYPES_PTRS
    global UA_TYPES_MAP
    ind = Int(Int((p - UA_TYPES_PTRS[0])) / sizeof(UA_DataType))
    return UA_TYPES_MAP[ind + 1]
end

function UA_init(p::Ref{T}) where T
    @ccall memset(p::Ptr{Cvoid}, 0::Cint, (sizeof(T))::Csize_t)::Ptr{Cvoid}
    return nothing
end

function Base.unsafe_wrap(s::Ptr{T}, field::Symbol) where T
    size_fieldname = Symbol(field, :Size)
    ptr = unsafe_load(getproperty(s,field))
    datasize = unsafe_load(getproperty(s, size_fieldname))
    return unsafe_wrap(Array, ptr, (datasize,))
end

struct UA_Array{T<:Ptr} <: AbstractArray{T, 1}
    ptr::T
    length::Int64
end

## UA_Array
# Julia wrapper for C array types
function UA_Array(s::T, field::Symbol) where T
    size_fieldname = Symbol(field, :Size)
    ptr = getfield(s,field)
    datasize = getfield(s, size_fieldname)
    return UA_Array(ptr, Int64(datasize))
end

Base.size(a::UA_Array) = (a.length,)
Base.length(a::UA_Array) = a.length
Base.IndexStyle(::Type{<:UA_Array}) = IndexLinear()
function Base.getindex(a::UA_Array{Ptr{T}}, i::Int) where T
    1 <= i <= a.length || throw(BoundsError(a, i))
    return a.ptr + (i-1) * sizeof(T)
end
Base.firstindex(a::UA_Array) = 1
Base.lastindex(a::UA_Array) = a.length
Base.setindex!(a::UA_Array, v, i::Int) = (a[i] = v)
Base.unsafe_wrap(a::UA_Array) = unsafe_wrap(Array, a[begin], size(a))
Base.pointer(a::UA_Array) = a[begin]
Base.convert(::Type{Ptr{T}}, a::UA_Array{Ptr{T}}) where T = a[begin]

function UA_Array_init(p::UA_Array)
    for i in p
        UA_init(i)
    end
end

function UA_Array_new(v::AbstractVector{T}, type_ptr::Ptr{UA_DataType}) where T
    v_typed = convert(Vector{juliadatatype(type_ptr)}, v) # Implicit check if T can be converted to type_ptr
    arr_ptr = convert(Ptr{T}, UA_Array_new(length(v), type_ptr))
    GC.@preserve v_typed unsafe_copyto!(arr_ptr, pointer(v_typed), length(v))
    return arr_ptr
end

function UA_print(p::Ref{T}, type_ptr::Ptr{UA_DataType} = ua_data_type_ptr_default(T)) where T
    buf = UA_String_new()
    UA_print(p, type_ptr, buf)
    s = unsafe_string(buf)
    UA_String_clear(buf)
    UA_String_delete(buf)
    return s
end

UA_print(v::T, type_ptr = ua_data_type_ptr_default(T)) where T = UA_print(Ref(v), type_ptr)

for (i, type_name) in enumerate(type_names)
    type_ind_name = Symbol("UA_TYPES_", uppercase(String(type_name)[4:end]))
    julia_type = julia_types[i]
    val_type = Val{type_name}

    @eval begin
        ua_data_type_ptr(::$(val_type)) = UA_TYPES_PTRS[$(i-1)]

        if !(type_names[$(i)] in types_ambiguous_denylist)
            ua_data_type_ptr_default(::Type{$(julia_type)}) = UA_TYPES_PTRS[$(i-1)]

            Base.show(io::IO, ::MIME"text/plain", v::$(julia_type)) = print(io, UA_print(v))
        end

        function $(Symbol(type_name, "_new"))()
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            return convert(Ptr{$(type_name)}, UA_new(data_type_ptr))
        end

        $(Symbol(type_name, "_init"))(p::Ptr{$(type_name)}) = UA_init(p)

        function $(Symbol(type_name, "_copy"))(src::Ptr{$(type_name)}, dst::Ptr{$(type_name)})
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            return UA_copy(src, dst, data_type_ptr)
        end

        function $(Symbol(type_name, "_clear"))(p::Ptr{$(type_name)})
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            UA_clear(p, data_type_ptr)
        end

        function $(Symbol(type_name, "_delete"))(p::Ptr{$(type_name)})
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            UA_delete(p, data_type_ptr)
        end

        function $(Symbol(type_name, "_Array_new"))(length::Integer)
            length <= 0 && error("Length of new array must be larger than zero.")
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            arr_ptr = convert(Ptr{$(type_name)}, UA_Array_new(length, data_type_ptr))
            return UA_Array(arr_ptr, length)
        end

        function $(Symbol(type_name, "_Array_new"))(v::AbstractVector)
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            v_typed = convert(Vector{$(type_name)}, v)
            arr_ptr = convert(Ptr{$(type_name)}, UA_Array_new(length(v), data_type_ptr))
            GC.@preserve v_typed unsafe_copyto!(arr_ptr, pointer(v_typed), length(v))
            return UA_Array(arr_ptr, length(v))
        end

        function $(Symbol(type_name, "_Array_init"))(p::UA_Array{Ptr{$(type_name)}})
            UA_Array_init(p)
        end

        function $(Symbol(type_name, "_Array_copy"))(
            src::Ptr{$(type_name)},
            dst::Ptr{$(type_name)},
            length::Integer,
        )
            length < 0 && error("Length of copied array must be larger than zero.")
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            return UA_Array_copy(src, length, Ref(dst), data_type_ptr)
        end

        function $(Symbol(type_name, "_Array_delete"))(p::Ptr{$(type_name)}, length::Integer)
            length < 0 && error("Length of deleted array must be larger than zero.")
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            UA_Array_delete(p, length, data_type_ptr)
        end

        function $(Symbol(type_name, "_Array_delete"))(p::UA_Array{Ptr{$(type_name)}})
            data_type_ptr = UA_TYPES_PTRS[$(type_ind_name)]
            UA_Array_delete(p, p.length, data_type_ptr)
        end
    end
end

## StatusCode

UA_StatusCode_name_print(sc::Integer) = return unsafe_string(UA_StatusCode_name(UA_StatusCode(sc)))

## String

# String `s` is copied to newly allocated memory that needs to be freed. Returns a pointer to a new `UA_String`.
function UA_String_set_alloc(data::AbstractString, ua_str::Ptr{UA_String})
    s = UA_String_fromChars(data)
    ua_str.data = s.data
    ua_str.length = s.length
    return ua_str
end

# String `s` must be kept valid using GC.@preserve as long as the return value is used. It is recommended to use UA_STRING_ALLOC with a subsequent call to UA_String_delete.
function UA_STRING_unsafe(s::AbstractString)
    GC.@preserve s begin
        isempty(s) && return UA_String(0, C_NULL)
        return UA_String(length(s), pointer(s))
    end
end

# String `s` is copied to newly allocated memory. The result's field `data` needs to be freed. Returns a `UA_String` struct.
UA_STRING_ALLOC(s::AbstractString) = UA_String_fromChars(s)

UA_String_delete(s::UA_String) = UA_Byte_delete(s.data)

Base.unsafe_string(s::UA_String) = unsafe_string(s.data, s.length)
Base.unsafe_string(s::Ref{UA_String}) = unsafe_string(s[])
Base.unsafe_string(s::Ptr{UA_String}) = unsafe_string(unsafe_load(s))

## DateTime

function UA_DateTime_toUnixTime(date::UA_DateTime)
    return (date - UA_DATETIME_UNIX_EPOCH) / UA_DATETIME_SEC
end

function UA_DateTime_fromUnixTime(unixDate::Integer)
    return UA_DateTime(unixDate * UA_DATETIME_SEC) + UA_DATETIME_UNIX_EPOCH
end

datetime2ua_datetime(dt::DateTime) = UA_DateTime_fromUnixTime(round(Int, datetime2unix(dt)))
ua_datetime2datetime(dt::UA_DateTime) = unix2datetime(UA_DateTime_toUnixTime(dt))

## Guid

function UA_GUID(s::AbstractString)
    guid = Ref{UA_Guid}()
    ua_s = UA_STRING_unsafe(s)
    retval = GC.@preserve s UA_Guid_parse(guid, ua_s)
    retval != UA_STATUSCODE_GOOD && error("Parsing of Guid \"$(s)\" failed with statuscode \"$(UA_StatusCode_name_print(retval))\".")
    return guid[]
end

## NodeId

function UA_NODEID_NUMERIC(nsIndex::Integer, identifier::Integer)
    identifier_tuple = anonymous_struct_tuple(UInt32(identifier), fieldtype(UA_NodeId, :identifier))
    return UA_NodeId(nsIndex, UA_NODEIDTYPE_NUMERIC, identifier_tuple)
end

# String `s` must be kept valid using GC.@preserve as long as the return value is used. It is recommended to use UA_NODEID_STRING_ALLOC with a subsequent call to UA_NodeId_delete.
function UA_NODEID_STRING_unsafe(nsIndex::Integer, s::AbstractString)
    GC.@preserve s identifier_tuple = anonymous_struct_tuple(UA_STRING_unsafe(s), fieldtype(UA_NodeId, :identifier))
    return UA_NodeId(nsIndex, UA_NODEIDTYPE_STRING, identifier_tuple)
end

# String `s` is copied to newly allocated memory that needs to be freed
function UA_NODEID_STRING_ALLOC(nsIndex::Integer, s::AbstractString)
    identifier_tuple = anonymous_struct_tuple(UA_String_fromChars(s), fieldtype(UA_NodeId, :identifier))
    return UA_NodeId(nsIndex, UA_NODEIDTYPE_STRING, identifier_tuple)
end

function UA_NODEID_GUID(nsIndex::Integer, guid::UA_Guid)
    identifier_tuple = anonymous_struct_tuple(guid, fieldtype(UA_NodeId, :identifier))
    return UA_NodeId(nsIndex, UA_NODEIDTYPE_GUID, identifier_tuple)
end

function UA_NodeId_delete(n::UA_NodeId)
    if n.identifier == UA_NODEIDTYPE_STRING
        UA_String_delete(n.identifier.string)
    end
    return nothing
end

UA_NodeId_equal(n1::Ref{UA_NodeId}, n2::Ref{UA_NodeId}) = UA_NodeId_order(n1, n2) == UA_ORDER_EQ

## ExpandedNodeId

# String `ns_uri` is copied to newly allocated memory that needs to be freed.
function UA_EXPANDEDNODEID_NUMERIC_ALLOC(
    nsIndex::Integer,
    identifier::Integer,
    ns_uri::AbstractString,
    server_ind::Integer,
)
    nodeid = UA_NODEID_NUMERIC(nsIndex, identifier)
    ua_ns_uri = UA_STRING_ALLOC(ns_uri)
    return UA_ExpandedNodeId(nodeid, ua_ns_uri, server_ind)
end

# Strings `s` and `ns_uri` are copied to newly allocated memory that needs to be freed.
function UA_EXPANDEDNODEID_STRING_ALLOC(
    nsIndex::Integer,
    s::AbstractString,
    ns_uri::AbstractString,
    server_ind::Integer,
)
    nodeid = UA_NODEID_STRING_ALLOC(nsIndex, s)
    ua_ns_uri = UA_STRING_ALLOC(ns_uri)
    return UA_ExpandedNodeId(nodeid, ua_ns_uri, server_ind)
end

# String `ns_uri` is copied to newly allocated memory that needs to be freed.
function UA_EXPANDEDNODEID_GUID_ALLOC(
    nsIndex::Integer,
    guid::UA_Guid,
    ns_uri::AbstractString,
    server_ind::Integer,
)
    nodeid = UA_NODEID_GUID(nsIndex, guid)
    ua_ns_uri = UA_STRING_ALLOC(ns_uri)
    return UA_ExpandedNodeId(nodeid, ua_ns_uri, server_ind)
end

function UA_ExpandedNodeId_equal(n1::Ref{UA_ExpandedNodeId}, n2::Ref{UA_ExpandedNodeId})
    return UA_ExpandedNodeId_order(n1, n2) == UA_ORDER_EQ
end

function UA_ExpandedNodeId_delete(n::UA_ExpandedNodeId)
    UA_NodeId_delete(n.nodeId)
    UA_String_delete(n.namespaceUri)
    return nothing
end

## QualifiedName

UA_QualifiedName_isNull(q::UA_QualifiedName) = (q.namespaceIndex == 0 && q.name.length == 0)
UA_QualifiedName_isNull(q::Ref{UA_QualifiedName}) = UA_QualifiedName_isNull(q[])

# String `s` must be kept valid using GC.@preserve as long as the return value is used. It is recommended to use UA_QUALIFIEDNAME_ALLOC with a subsequent call to UA_QualifiedName_delete.
function UA_QUALIFIEDNAME(nsIndex::Integer, s::AbstractString)
    GC.@preserve s return UA_QualifiedName(nsIndex,  UA_STRING_unsafe(s))
end

# String `s` is copied to newly allocated memory that needs to be freed.
function UA_QUALIFIEDNAME_ALLOC(nsIndex::Integer, s::AbstractString)
    return UA_QualifiedName(nsIndex, UA_String_fromChars(s))
end

UA_QualifiedName_delete(q::UA_QualifiedName) = UA_String_delete(q.name)

## LocalizedText

# Strings `locale` and `text` must be kept valid using GC.@preserve as long as the return value is used. It is recommended to use UA_QUALIFIEDNAME_ALLOC with a subsequent call to UA_LocalizedText_delete.
function UA_LOCALIZEDTEXT_unsafe(locale::AbstractString, text::AbstractString)
    GC.@preserve locale text begin
        return UA_LocalizedText(UA_STRING_unsafe(locale), UA_STRING_unsafe(text))
    end
end

# Strings `locale` and `text` are copied to newly allocated memory that needs to be freed
function UA_LOCALIZEDTEXT_ALLOC(locale::AbstractString, text::AbstractString)
    return UA_LocalizedText(UA_STRING_ALLOC(locale), UA_STRING_ALLOC(text))
end

UA_LocalizedText_delete(l::UA_LocalizedText) = UA_String_delete(l.text)

## NumericRange

function UA_NUMERICRANGE(s::AbstractArray)
    nr = Ref{UA_NumericRange}()
    retval = GC.@preserve s UA_NumericRange_parse(nr, UA_STRING_unsafe(s))
    retval != UA_STATUSCODE_GOOD && error("Parsing of NumericRange \"$(s)\" failed with statuscode \"$(UA_StatusCode_name_print(retval))\".")
    return nr[]
end

## Variant

function unsafe_size(v::UA_Variant)
    UA_Variant_isScalar(v) && return ()
    v.arrayDimensionsSize == 0 && return (Int(v.arrayLength),)
    return Tuple([Int(unsafe_load(v.arrayDimensions, d+1)) for d = 0:unsafe_load(v.arrayDimensionsSize)])
end

unsafe_size(p::Ref{UA_Variant}) = unsafe_size(unsafe_load(p))

Base.length(v::UA_Variant) = Int(v.arrayLength)
Base.length(p::Ref{UA_Variant}) = length(unsafe_load(p))

function UA_Variant_new_copy(value::T, type_ptr::Ptr{UA_DataType}) where T
    var = UA_Variant_new()
    var.type = type_ptr
    var.storageType = UA_VARIANT_DATA
    var.arrayLength = length(value)
    var.arrayDimensionsSize = length(size(value))
    if isempty(size(value)) # Scalar value
        var.arrayDimensions = C_NULL
        value_ptr = convert(Ptr{T}, UA_new(type_ptr))
        unsafe_store!(value_ptr, value)
    else # Array value
        value_row_major = permutedims(value, reverse(1:length(size(value))))
        var.data = UA_Array_new(value, type_ptr)
        var.arrayDimensions = UA_UInt32_Array_new([size(value_row_major)...])
    end
    return var
end

UA_Variant_new_copy(value::T) where T = UA_Variant_new_copy(value, ua_data_type_ptr_default(eltype(T)))
UA_Variant_new_copy(value, type_sym::Symbol) = UA_Variant_new_copy(value, ua_data_type_ptr(Val(type_sym)))

function Base.unsafe_wrap(v::UA_Variant)
    type = juliadatatype(v.type)
    UA_Variant_isScalar(v) && return unsafe_load(reinterpret(Ptr{type}, v.data))

    values = unsafe_wrap(Array, reinterpret(Ptr{type}, v.data), unsafe_size(v))
    values_row_major = reshape(values, unsafe_size(v))
    return permutedims(values_row_major, reverse(1:length(unsafe_size(v)))) # To column major format
end

Base.unsafe_wrap(p::Ref{UA_Variant}) = unsafe_wrap(unsafe_load(p))

UA_Variant_isEmpty(v::UA_Variant) = v.type == C_NULL
UA_Variant_isEmpty(p::Ref{UA_Variant}) = UA_Variant_isEmpty(unsafe_load(p))

UA_Variant_isScalar(v::UA_Variant) = v.arrayLength == 0 && v.data > UA_EMPTY_ARRAY_SENTINEL
UA_Variant_isScalar(p::Ref{UA_Variant}) = UA_Variant_isScalar(unsafe_load(p))

function UA_Variant_hasScalarType(v::UA_Variant, type::Ref{UA_DataType})
    return UA_Variant_isScalar(v) && type == v.type
end

function UA_Variant_hasScalarType(p::Ref{UA_Variant}, type::Ref{UA_DataType})
    return UA_Variant_hasScalarType(unsafe_load(p), type)
end

function UA_Variant_hasArrayType(v::UA_Variant, type::Ref{UA_DataType})
    return !UA_Variant_isScalar(v) && type == v.type
end

function UA_Variant_hasArrayType(p::Ref{UA_Variant}, type::Ref{UA_DataType})
    return UA_Variant_hasArrayType(unsafe_load(p), type)
end


text:
using Base.Test
using Meshes

# run lint if run with --lint
if "--lint" in ARGS
    using Lint
    println("Running Lint...")
    lintpkg("Meshes")
end


text:
abstract type AbstractMatrixOperator{Tconv} end

@params struct MatrixOperator{Tconv} <: AbstractMatrixOperator{Tconv}
    K
    f
    conv::Tconv
end
LinearAlgebra.mul!(c, op::MatrixOperator, b) = mul!(c, op.K, b)
Base.size(op::MatrixOperator, i...) = size(op.K, i...)
Base.eltype(op::MatrixOperator) = eltype(op.K)
LinearAlgebra.:*(op::MatrixOperator, b) = mul!(similar(b), op.K, b)

@params struct MatrixFreeOperator{Tconv, T, dim} <: AbstractMatrixOperator{Tconv}
    f::AbstractVector{T}
    elementinfo::ElementFEAInfo{dim, T}
    meandiag::T
    vars::AbstractVector{T}
    xes
    fixed_dofs
    free_dofs
    xmin::T
    penalty
    conv::Tconv
end
GPUUtils.whichdevice(m::MatrixFreeOperator) = whichdevice(m.vars)
Base.size(op::MatrixFreeOperator) = (size(op, 1), size(op, 2))
Base.size(op::MatrixFreeOperator, i) = 1 <= i <= 2 ? length(op.elementinfo.fixedload) : 1
Base.eltype(op::MatrixFreeOperator{<:Any, T}) where {T} = T

import LinearAlgebra: *, mul!

function *(A::MatrixFreeOperator, x)
    y = similar(x)
    mul!(y, A::MatrixFreeOperator, x)
    y
end

function mul!(y::TV, A::MatrixFreeOperator, x::TV) where {TV <: AbstractVector}
    T = eltype(y)
    nels = length(A.elementinfo.Kes)
    ndofs = length(A.elementinfo.fixedload)
    dofspercell = size(A.elementinfo.Kes[1], 1)
    meandiag = A.meandiag

    @unpack Kes, metadata, black, white, varind = A.elementinfo
    @unpack cell_dofs, dof_cells = metadata
    @unpack penalty, xmin, vars, fixed_dofs, free_dofs, xes = A

    for i in 1:nels
        px = ifelse(black[i], one(T),
                    ifelse(white[i], xmin,
                            density(penalty(vars[varind[i]]), xmin)
                        )
                    )
        xe = xes[i]
        for j in 1:dofspercell
            xe = @set xe[j] = x[cell_dofs[j,i]]
        end
        if eltype(Kes) <: Symmetric
            xes[i] = px * (bcmatrix(Kes[i]).data * xe)
        else
            xes[i] = px * (bcmatrix(Kes[i]) * xe)
        end
    end

    for i in 1:length(fixed_dofs)
        dof = fixed_dofs[i]
        y[dof] = meandiag * x[dof]
    end
    for i in 1:length(free_dofs)
        dof = free_dofs[i]
        yi = zero(T)
        r = dof_cells.offsets[dof] : dof_cells.offsets[dof+1]-1
        for ind in r
            k, m = dof_cells.values[ind]
            yi += xes[k][m]
        end
        y[dof] = yi
    end
    y
end

function mul!(y::TV, A::MatrixFreeOperator, x::TV) where {TV <: CuArrays.CuVector}
    T = eltype(y)
    nels = length(A.elementinfo.Kes)
    ndofs = length(A.elementinfo.fixedload)
    meandiag = A.meandiag

    @unpack Kes, metadata, black, white, varind = A.elementinfo
    @unpack cell_dofs, dof_cells = metadata
    @unpack penalty, xmin, vars, fixed_dofs, free_dofs, xes = A

    args1 = (xes, x, black, white, vars, varind, cell_dofs, Kes, xmin, penalty, nels)
    callkernel(dev, mul_kernel1, args1)
    CUDAdrv.synchronize(ctx)

    args2 = (y, x, dof_cells.offsets, dof_cells.values, xes, fixed_dofs, free_dofs, meandiag)
    callkernel(dev, mul_kernel2, args2)
    CUDAdrv.synchronize(ctx)

    return y
end

# CUDA kernels
function mul_kernel1(xes::AbstractVector{TV}, x, black, white, vars, varind, cell_dofs, Kes, xmin, penalty, nels) where {N, T, TV<:SVector{N, T}}
    i = @thread_global_index()
    offset = @total_threads()
    while i <= nels
        px = ifelse(black[i], one(T),
                    ifelse(white[i], xmin,
                            density(penalty(vars[varind[i]]), xmin)
                        )
                    )
        xe = xes[i]
        for j in 1:N
            xe = @set xe[j] = x[cell_dofs[j, i]]
        end
        if eltype(Kes) <: Symmetric
            xe = SVector{1, T}((px,)) .* (bcmatrix(Kes[i]).data * xe)
        else
            xe = SVector{1, T}((px,)) .* (bcmatrix(Kes[i]) * xe)
        end

        xes[i] = xe
        i += offset
    end

    return
end

function mul_kernel2(y, x, dof_cells_offsets, dof_cells_values, xes, fixed_dofs, free_dofs, meandiag)
    T = eltype(y)
    offset = @total_threads()
    n_fixeddofs = length(fixed_dofs)
    ndofs = length(y)

    i = @thread_global_index()
    while i <= n_fixeddofs
        dof = fixed_dofs[i]
        y[dof] = meandiag * x[dof]
        i += offset
    end

    i = @thread_global_index()
    while n_fixeddofs < i <= ndofs
        dof = free_dofs[i - n_fixeddofs]
        yi = zero(T)
        r = dof_cells_offsets[dof] : dof_cells_offsets[dof+1]-1
        for ind in r
            k, m = dof_cells_values[ind]
            yi += xes[k][m]
        end
        y[dof] = yi
        i += offset
    end
    return
end


text:
"""
splits long vector with all the outputs concatented into array of
vectors of individule outputs
"""
function split_outputs(y::Array{<:AbstractFloat}, t::Array{Array{<:AbstractFloat,1},1})
    s = [size(ti)[1] for ti in t]
    out = fill(Float64[], size(t)[1])
    out[1] = y[1:s[1]]
    sc = s[1]
    for i = 2:size(t)[1]
        out[i] = y[sc + 1:sc + s[i]]
        sc += s[i]
    end
    out
end

sample_mean(s) = sum(s) / size(s)[1]
sample_std(s) = sqrt(sum((s .- sample_mean(s)).^2) / size(s)[1])

function summary_stats(s)
    m = sample_mean(s)
    std = sample_std(s)
    println("$m ± $std ")
end

Base.copy(gp::GaussianProcess) = GaussianProcess(gp.base_kernel, gp.D, gp.C, gp.P, gp.data, gp.dpars)


text:
# some simple tests with periodic NN style functoins for central pattern generators.
# The goal is to effectively learn the parameters of sine and cosine functions via gradient descent on a time series
using Plots

arr = zeros(Float32,1000)
arr2 = zeros(Float32, 1000)
arr3 = zeros(Float32, 1000)
for i in 1:1000
    arr[i] = sin(0.05 * i)
    arr2[i] = cos(0.05 * i)
    arr3[i] = tan(0.05 * i)
end
plot(arr)
plot!(arr2)
plot!(arr3)

using Flux.Tracker: update!
data = [sin(0.05f0* i) for i in 1:1000]

prediction(period,i) = sin(i * period)

MSE(x,pred) = (x .- pred).^2 # i.e. flux mse loss
loss(data, preds) = mean(MSE.(data, preds))
period = param(0.001)
periods = []
ls = []
predslist = []
grads = []
for i in 1:100
    global period
    println("epoch $i, period: $period")
    preds = [prediction(period,j) for j in 1:1000]
    l = loss(data, preds)
    gs = Tracker.gradient(()->loss(data,preds), params(period))
    #period -= 0.01 * grad
    #period = period.data
    #period = param(period)
    grad = min(max(gs[period].data,1f0),-1f0)
    println("grad: $grad")
    update!(period, -0.01 *grad)
    push!(periods, period.data)
    push!(ls, l.data)
    push!(predslist, Tracker.data.(preds))
    push!(grads, grad)

end
plot(ls)
plot(periods)
plot(predslist[10])
plot(data)
plot(grads)

W = rand(2, 5)
b = rand(2)

prediction(x) = W*x .+ b

function loss(x, y)
  ŷ = prediction(x)
  sum((y .- ŷ).^2)
end

x, y = rand(5), rand(2) # Dummy data
loss(x, y) # ~ 3

# for rhythmical NNs
period = param(0.001)
periods = []
ls = []
predslists = []
grads = []
datalists =[]

for i in 1:100

    global period,L
    println("epoch $i, period: $period")
    for j in 1:40
        predslist = []
        datalist = []
        pred = prediction(period, j)
        data = sin(0.1 * j)
        L = MSE(pred, data)
        gs = Tracker.gradient(()->MSE(pred, data), params(period))
        grad = min(max(gs[period].data,1f0),-1f0)
        #grad = gs[period]
        print("grad: $grad")
        #println("grad: $grad")
        update!(period, -0.001 *grad)
        push!(predslist, pred.data)
        #push!(grads, grad)
        push!(datalist, data)
        println(period)
        println("L: $L")
    end
    push!(datalists, datalist)
    push!(predslists, predslist)
    datalist = []
    predslist = []
    println("outside: $period")
    println("outside l: $L")
    push!(periods, period.data)
    push!(ls, L.data)

end
plot(ls)
plot(periods)

d = [sin(0.001 * i) for i in 1:10000]
plot(d)
period = param(0.01)
pred2(period) = [sin(period*i) for i in 1:10000]
loss(x,y) = mean((x .- y).^2)
periods = []
ls = []
grads = []
preds = []
##
for i in 1:1000
    global period
    l = loss(d, pred2(period))
    pred = pred2(period)
    push!(preds, Tracker.data.(pred))
    gs = Tracker.gradient(()-> loss(d, pred2(period)), params(period))
    grad = gs[period]
    update!(period, -0.001 * grad)
    push!(ls, l.data)
    push!(periods, period.data)
    push!(grads, grad.data)
end
plot(ls)
plot(periods)
plot(grads)
plot(d)
plot!(preds[end])
mutable struct CPG
    i::Int
    p::AbstractFloat
end

function forward(c::CPG)
    pred =  sin(c.i * c.p)
    c.i = c.i+1
    return pred
end

a = CPG(0,0.01)
forward(a)
a.i
##
using Flux
params(a)
d = [sin(0.001 * i) for i in 1:10000]
plot(d)
period = param(0.01)
pred2(c)
loss(x,y) = mean((x .- y).^2)
periods = []
ls = []
grads = []
preds = []
##
for i in 1:1000
    global period
    l = loss(d, pred2(period))
    pred = pred2(period)
    push!(preds, Tracker.data.(pred))
    gs = Tracker.gradient(()-> loss(d, pred2(period)), params(period))
    grad = gs[period]
    update!(period, -0.001 * grad)
    push!(ls, l.data)
    push!(periods, period.data)
    push!(grads, grad.data)
end
plot(ls)
plot(periods)
plot(grads)
plot(d)
plot!(preds[end])


text:
include("common.jl")
using LinearAlgebra: rank, eigvals

function balanced_ρ(p, grav, mesh, elm)
  zdg = mesh.zdg
  Nq, K = size(zdg)
  𝟙 = ones(Bool, Nq)

  pdg = reshape(mesh.scatter * p, Nq, K)
  Q = elm.Q

  # Form the RHS vector
  pdg_diff = similar(zdg)
  Fρdg = similar(zdg, Nq * K, Nq * K)
  for e in 1:K
    p̄ = (pdg[:, e] .+ pdg[:, e]') / 2
    pdg_diff[:, e] = ((Q - Q') .* p̄) * 𝟙
  end
  # BCs hack
  pdg_diff[1,1] -= p[1,1]
  pdg_diff[end,end] += p[end,end]
  p_diff = mesh.scatter' * pdg_diff[:]
  @show p_diff

  # Form the ρsystem
  T = eltype(zdg)
  Fρdg = zeros(T, Nq * K, Nq * K)
  for e in 1:K
    rng = (e - 1) * Nq .+ (1: Nq)
    Δϕ = grav * (zdg[:, e]' .- zdg[:, e]) / 2
    A = (Q - Q') .* Δϕ
    Fρdg[rng, rng] .= ((Diagonal(A * 𝟙) + A) / 2 )
  end
  Fρ = mesh.scatter' * Fρdg * mesh.scatter
  ρ_bal = -Fρ \ p_diff
  return ρ_bal
end

let
  T = Float64
  K = 10
  N = 5
  #
  elm = element_operators(N)
  mesh = create_mesh(K, elm, T(0), T(30e4))
  zcg = mesh.zcg
  #
  Nq = N + 1
  #
  params = Params{T}()
  p = similar(zcg)
  ρ = similar(zcg)
  for (i, z) in enumerate(zcg)
    (p[i], ρ[i]) = decaying_temperature_profile(params, z)
  end
  #
  ρdg = reshape(mesh.scatter * ρ, Nq, K)
  pdg = reshape(mesh.scatter * p, Nq, K)
  #
  vdg = fill!(similar(ρdg), 0)
  #
  Q = elm.Q
  zdg = mesh.zdg
  grav = params.grav
  #
  pdg_diff_0 = fill!(similar(pdg), 0)
  ρdg_diff_0 = fill!(similar(ρdg), 0)
  for e in 1:K
    for i in 1:Nq, j in 1:Nq
      p_i, p_j = pdg[i, e], pdg[j, e]
      ρ_i, ρ_j = ρdg[i, e], ρdg[j, e]
      ϕ_i, ϕ_j = grav * zdg[i, e], grav * zdg[j, e]
      # flux terms
      ρ̄ = (ρ_i + ρ_j) / 2
      Δϕ = (ϕ_j - ϕ_i) / 2
      p̄ = (p_i + p_j) / 2
      # flux differences
      pdg_diff_0[i, e] += (Q[i, j] - Q[j, i]) * p̄
      ρdg_diff_0[i, e] += (Q[i, j] - Q[j, i]) * ρ̄ * Δϕ
    end
  end
  # BCs hack
  pdg_diff_0[1,1] -= p[1,1]
  pdg_diff_0[end,end] += p[end,end]
  #
  # @show extrema(mesh.scatter'*(pdg_diff_0 + ρdg_diff_0)[:])
  #
  # Remove the inner for loops
  pdg_diff_1 = fill!(similar(pdg), 0)
  ρdg_diff_1 = fill!(similar(ρdg), 0)
  𝟙 = ones(Bool, Nq)
  for e in 1:K
    p̄ = (pdg[:, e] .+ pdg[:, e]') / 2
    ρ̄ = (ρdg[:, e] .+ ρdg[:, e]') / 2
    Δϕ = grav * (zdg[:, e]' .- zdg[:, e]) / 2
    pdg_diff_1[:, e] = ((Q - Q') .* p̄) * 𝟙
    ρdg_diff_1[:, e] = ((Q - Q') .* ρ̄ .* Δϕ) * 𝟙
  end
  # BCs hack
  pdg_diff_1[1,1] -= p[1,1]
  pdg_diff_1[end,end] += p[end,end]
  @assert all(pdg_diff_0 .≈ pdg_diff_1)
  @assert all(ρdg_diff_0 .≈ ρdg_diff_1)
  #
  # Build the ρ matrix
  Fρdg = zeros(T, Nq * K, Nq * K)
  for e in 1:K
    rng = (e - 1) * Nq .+ (1: Nq)
    Δϕ = grav * (zdg[:, e]' .- zdg[:, e]) / 2
    A = (Q - Q') .* Δϕ
    Fρdg[rng, rng] .= ((Diagonal(A * 𝟙) + A) / 2 )
  end
  ρdg_diff_2 = reshape(Fρdg * ρdg[:], Nq, K)
  @assert all(ρdg_diff_0 .≈ ρdg_diff_2)
  #
  Fρ = mesh.scatter' * Fρdg * mesh.scatter
  ρ_bal = -Fρ \ (mesh.scatter' * pdg_diff_0[:])
  ρdg_bal = reshape(mesh.scatter * ρ_bal, Nq, K)
  display(ρdg)
  display(ρdg_bal)
  display(ρdg_bal - ρdg)
  ρ_bal2 = balanced_ρ(p, grav, mesh, elm)
  @show extrema(ρ_bal2 - ρ_bal)
  nothing
end


text:
macro interface(ex)
    return interfacem(__module__, __source__, ex)
end

function interfacem(__module__::Module, __source__::LineNumberNode, ex::Expr)
    r = handle(ex)

    if r === nothing
        return :(error("expect a function definition or a function call as interface definition"))
    end

    name, args, body = r

    if body === nothing && args === nothing
        return quote
            export $(esc(name))
            Core.@__doc__ $(esc(ex)) = throw(NotImplementedError($(QuoteNode(name))))
        end
    elseif body === nothing
        return quote
            export $(esc(name))
            Core.@__doc__ $(esc(ex)) = throw(NotImplementedError($(QuoteNode(name)), tuple($(esc.(args)...))))

            # push!(METADATA, Interface($(QuoteNode(handle(ex))), $(Base.Docs.signature(ex))))
        end
    else
        return quote
            export $(esc(name))
            Core.@__doc__ $(esc(ex))

            # push!(METADATA, Interface($(QuoteNode(handle(ex))), $(Base.Docs.signature(ex))))
        end
    end
end

struct Interface
    name::Symbol
    signatures::DataType
end

handle(ex::Expr, isvalid=false) = handle(Val(ex.head), ex, isvalid)
handle(ex::Symbol, isvalid) = isvalid ? (ex, nothing, nothing, nothing) : error("expect a function")
handle(::Val{:where}, ex::Expr, isvalid) = handle(ex.args[1], isvalid)

# valid syntax
function handle(::Val{:(=)}, ex::Expr, isvalid)
    name, args, _ = handle(ex.args[1], isvalid)
    return name, args, ex.args[2]
end

function handle(::Val{:function}, ex::Expr, isvalid)
    name, args, _ = handle(ex.args[1], true)
    if isempty(ex.args[2:end])
        return name, args, nothing
    else
        return name, args, ex.args[2:end]
    end
end

function handle(::Val{:call}, ex::Expr, isvalid)
    args = []
    for each in ex.args[2:end]
        if !(each isa Union{Expr, Symbol})
            return nothing
        end
        if each isa Expr && each.head === :parameters
            continue
        else
            push!(args, name_handle(each))
        end
    end
    return name_handle(ex.args[1]), args, nothing
end

name_handle(x::Symbol) = x
name_handle(x::QuoteNode) = x.value
name_handle(x::Expr) = name_handle(Val(x.head), x)
# type annotation
function name_handle(::Val{:(::)}, x::Expr)
    if length(x.args) == 2
        return name_handle(x.args[1])
    else
        temp = x.args[1]
        x.args[1] = gensym(:x)
        push!(x.args, temp)
        return x.args[1]
    end
end

name_handle(::Val{:(...)}, x::Expr) = Expr(:(...), name_handle(x.args[1]))
name_handle(::Val{:(.)}, x::Expr) = name_handle(x.args[2])
name_handle(::Val{:kw}, x::Expr) = name_handle(x.args[1])
name_handle(::Val, x::Expr) = error("expect a valid variable name got $x")


text:
q = Quaternion(rand() * 2pi, normalize(ℝ³(rand(3))))
cp = ComplexPlane(q)
s = SU2(cp)
θ = rand() * pi
u = normalize(ℝ³(rand(3)))
g = Quaternion(θ, u)

@test isapprox(cp, s)
@test isapprox(cp, q)
@test isapprox(g * q, g * Quaternion(s))

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
@testset "matching $G" begin

g =G(2)
add_edge!(g,1,2)
w = EdgeMap(g, Int)
w[1,2]=500
match = minimum_weight_perfect_matching(g, w)
@test match.mate[1] == 2

g = CompleteGraph(4, G)
w = EdgeMap(g, Int)
w[1,2]=500
w[1,3]=600
w[2,3]=700
w[3,4]=100
w[2,4]=1000
match = minimum_weight_perfect_matching(g, w)
@test match.mate[1] == 2
@test match.mate[2] == 1
@test match.mate[3] == 4
@test match.mate[4] == 3
@test match.weight == 600

g = CompleteGraph(4, G)
w = EdgeMap(g, Int)
w[1, 2]=500
w[1, 3]=400
w[2, 3]=300
w[3, 4]=1000
w[2, 4]=1000
match = minimum_weight_perfect_matching(g, w)
@test match.mate[1] == 3
@test match.mate[2] == 4
@test match.mate[3] == 1
@test match.mate[4] == 2
@test match.weight == 1400

g =CompleteBipartiteGraph(2,2, G)
w = EdgeMap(g, Float64)
w[1,3] = 10
w[1,4] = 19.5
w[2,3] = 9
w[2,4] = 19
match = minimum_weight_perfect_matching(g, w)
@test match.mate[1] == 4
@test match.mate[4] == 1
@test match.mate[2] == 3
@test match.mate[3] == 2
@test match.weight == 28.5

g =CompleteGraph(4, G)
w = EdgeMap(g, Float64)
w[1,3] = 10
w[1,4] = 0.5
w[2,3] = 11
w[2,4] = 2
w[1,2] = 100
match = minimum_weight_perfect_matching(g, w, cutoff=50)
@test match.mate[1] == 4
@test match.mate[4] == 1
@test match.mate[2] == 3
@test match.mate[3] == 2
@test match.weight == 11.5

w = (s=sprand(100,100,0.2); s=s+s'; s-Diagonal(s))
g = Graph(w)
match = minimum_weight_perfect_matching(g, EdgeMap(g, w))
for i=1:nv(g)
    j = match.mate[i]
    @test match.mate[j] == i
end

# TEST matching 1D
dist_periodic_bc(a,b,R) = abs(a-b)<0.5*R ? abs(a-b) : R-abs(a-b)

N=100
p=2.
g = CompleteGraph(N)
for r in 1:100
    w = EdgeMap(g, Float64)
    points = sort(rand(N))
    for e in edges(g)
        i, j = src(e), dst(e)
        w[i,j] = dist_periodic_bc(points[i], points[j], 1)^p
    end
    match = minimum_weight_perfect_matching(g, w)
    if p>1
        if match.mate[1]==2
            for i=2:nv(g)
                @test match.mate[i] ==(i - (-1)^i) % (N+1)
            end
        elseif match.mate[1] == N
            for i=2:nv(g)
                @test match.mate[i] == (i + (-1)^i) % N
            end
        end
    elseif p < 0
        for i=1:nv(g)
            @test match.mate[i] == (i + N/2) % N
        end
    end
end

end # testset