text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

using LinearAlgebra
using DoubleFloats
# using SemiLagMPI

using SemiLagrangian

# function landau_old(
#     dt::T,
#     epsilon::T,
#     nbdt,
#     tabcoef::Vector{T},
#     mesh_x::UniformMesh{T},
#     mesh_v::UniformMesh{T},
#     interp_x::AbstractInterpolation{T},
#     interp_v::AbstractInterpolation{T},
# ) where {T,ndims}
#     adv_x = Advection1d(mesh_x, interp_x)
#     adv_v = Advection1d(mesh_v, interp_v)

#     nx = mesh_x.length
#     nv = mesh_v.length

#     t_coef = vcat(tabcoef, tabcoef[end-1:-1:1])
#     nbcoef = size(t_coef, 1)

#     fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))
#     fct_x(x) = epsilon * cos(x / 2) + 1
#     lgn_x = fct_x.(mesh_x.points)
#     lgn_v = fct_v.(mesh_v.points)

#     v = mesh_v.points

#     if ndims == 1
#         nxtp = (nx)
#         nvtp = (nv)
#         fxv = zeros(T, (nx, nv))
#         fvx = zeros(T, (nv, nx))
#         fxv .= lgn_x .* reshape(lgn_v, (1, nv))
#         #        fxv .= fct_x.(mesh_x.points)  .* transpose(fct_v.(mesh_v.points))
#         perm = [2, 1]
#     elseif ndims == 2
#         nxtp = (nx, nx)
#         nvtp = (nv, nv)
#         fxv = zeros(nx, nx, nx, nv)
#         fvx = zeros(nv, nv, nx, nx)
#         fxv =
#             lgn_x .* reshape(lgn_x, (1, nx)) .* reshape(lgn_v, (1, 1, nv)) .*
#             reshape(lgn_v, (1, 1, 1, nv))
#         perm = [3, 4, 1, 2]
#     else
#         println("not yet implemented !!!")
#     end
#     elf = Array{T,ndims}(undef, nxtp)
#     rho = Array{T,ndims}(undef, nxtp)
#     printout(advd, "# dt=$(Float64(dt)) eps=$(Float64(epsilon)) size_x=$nx size_v=$nv")
#     printout(advd, "# x : from $(Float64(mesh_x.start)) to $(Float64(mesh_x.stop))")
#     printout(advd, "# v : from $(Float64(mesh_v.start)) to $(Float64(mesh_v.stop))")
#     printout(advd, "# interpolation : $(get_type(interp_x)) order=$(get_order(interp_x))")
#     printout(advd, "# type=$T precision = $(precision(T))")
#     printout(advd, "#time\tel-energy\tkinetic-energy\tglobal-energy")

#     #   transpose!(fvx, fxv)
#     permutedims!(fvx, fxv, perm)
#     compute_charge!(rho, mesh_v, fvx)
#     compute_elfield!(elf, mesh_x, rho)
#     elenergy = Float64(compute_ee(mesh_x, elf))
#     kinenergy = Float64(compute_ke(mesh_v, mesh_x, fvx))
#     energyall = elenergy + kinenergy
#     printout(advd, "$(Float64(0.))\t$elenergy\t$kinenergy\t$energyall")

#     minall = 10000000
#     maxall = 0
#     t_coef *= dt
#     fltraceend = true
#     for i = 1:nbdt
#         # advection!(adv_x, fxv, v, dt/2)
#         # transpose!(fvx, fxv)
#         # compute_charge!(rho, mesh_v, fvx)
#         # compute_elfield!(elf, mesh_x, rho)
#         # advection!(adv_v, fvx, elf, dt)
#         # transpose!(fxv, fvx)
#         # advection!(adv_x, fxv, v, dt/2)
#         # transpose!(fvx, fxv)
#         # compute_charge!(rho, mesh_v, fvx)
#         # compute_elfield!(elf, mesh_x, rho)
#         # elenergy = Float64(compute_ee(mesh_x, elf))
#         # kinenergy = Float64(compute_ke(mesh_v, mesh_x, fvx))
#         # energyall = elenergy + kinenergy
#         # if (i%modulo == 0)
#         #     println("$(Float64(i*dt))\t$elenergy\t$kinenergy\t$energyall")
#         # end
#         # minall=min(energyall,minall)
#         # maxall=max(energyall,maxall)

#         for k = 1:size(t_coef, 1)
#             if k % 2 == 1
#                 advection!(adv_x, fxv, v, t_coef[k])
#                 permutedims!(fvx, fxv, perm)
#                 if fltraceend || k != size(t_coef, 1)
#                     compute_charge!(rho, mesh_v, fvx)
#                     compute_elfield!(elf, mesh_x, rho)
#                 end
#             else
#                 advection!(adv_v, fvx, elf, t_coef[k])
#                 permutedims!(fxv, fvx, perm)
#             end
#         end
#         if fltraceend
#             elenergy = Float64(compute_ee(mesh_x, elf))
#             kinenergy = Float64(compute_ke(mesh_v, mesh_x, fvx))
#             energyall = elenergy + kinenergy
#             printout(advd, "$(Float64(i*dt))\t$elenergy\t$kinenergy\t$energyall")
#         end
#         minall = min(energyall, minall)
#         maxall = max(energyall, maxall)
#     end
#     printout(advd, "diff=$(maxall-minall)")
# end
function printout(
    advd::AdvectionData{T,N,timeopt},
    str,
) where {T,N,timeopt}
    if timeopt != MPIOpt || advd.adv.mpid.ind == 1
        println(str)
    end
end
printout(str) = println(str)

function trace_energy(
    advd::AdvectionData{T,N,timeopt},
    t,
) where {T,N,timeopt}

    if t == 0
        printout(advd, "#time\tel-energy\tkinetic-energy\tglobal-energy")
    end
    # global cl_obs
    # clockbegin(cl_obs,6)
    compute_charge!(advd)
    compute_elfield!(advd)
    # clockend(cl_obs,6)
    # clockbegin(cl_obs,7)
    elenergy = compute_ee(advd)
    # clockend(cl_obs,7)
    # clockbegin(cl_obs,8)
    kinenergy = compute_ke(advd)
    # clockend(cl_obs,8)
    energyall = elenergy + kinenergy
    printout(
        advd,
        "$(Float32(t))\t$(Float64(elenergy))\t$(Float64(kinenergy))\t$(Float64(energyall))",
    )
    return energyall
end

function landau(advd::AdvectionData, nbdt)

    # global cl_obs
    # clockreset(cl_obs)
    maxdiff = 0
    dt = advd.adv.dt_base
    refel = getenergyall(advd)
    maxel =minel = refel
    #    printall(cl_obs)
    #    clockreset(cl_obs)
    for i = 1:nbdt
        while advection!(advd)
        end
        el = getenergyall(advd)
        maxel = max(maxel, el)
        minel = min(minel,el)
        println("$(maxel-minel)")
        # printall(cl_obs)
        # clockreset(cl_obs)
    end
    println("#  end")
    # printall(cl_obs)
    return maxel-minel
end
function landau1_1(
    T::DataType,
    nbdt,
    timeopt;
    sz = (64, 64),
    dt = big"0.1",
    interp = Lagrange(9, T),
    tab_coef=[1//2, 1, 1//2]
)
    epsilon = T(0.5)
    dt = T(dt)

    spmin, spmax, nsp = T(0), T(4big(pi)), sz[1]
    vmin, vmax, nv = -T(10), T(10), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    tabst = [( [1,2], 1, 1, true),( [2,1], 1, 2, true) ]

    adv = Advection(
        (mesh_sp,mesh_v,),
        [interp,interp],
        dt,
        tabst,
        tab_coef=tab_coef,
        timeopt = timeopt)

    fct_sp(x) = epsilon * cos(x / 2) + 1
    fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))

    lgn_sp = fct_sp.(mesh_sp.points)
    lgn_v = fct_v.(mesh_v.points)

    data = dotprod((lgn_sp, lgn_v))

    pvar = getpoissonvar(adv)

    advd = AdvectionData(adv, data, pvar)

    printout(advd, "# dt=$(Float64(dt)) eps=$(Float64(epsilon)) size_x=$nsp size_v=$nv")
    printout(advd, "# sp : from $(Float64(start(mesh_sp))) to $(Float64(stop(mesh_sp)))")
    printout(advd, "# v : from $(Float64(start(mesh_v))) to $(Float64(stop(mesh_v)))")
    printout(advd, "# interpolation : $interp order=$(get_order(interp))")
    printout(advd, "# tab_coef : $tab_coef")
    printout(advd, "# type=$T precision = $(precision(T))")
    printout(advd, "# timeopt=$timeopt")
    if timeopt == SimpleThreadsOpt || timeopt == SplitThreadsOpt
        printout(advd, "# nb threads : $(Threads.nthreads())")
    elseif timeopt == MPIOpt
        printout(advd, "# nb process : $(adv.mpid.nb)")
    else
        printout(advd, "# monothread version")
    end
    printout(advd, "typeof(data)=$(typeof(data)) size(data)=$(size(data))")

    # advdata = Advection1dData(adv, data, pvar)

    landau(advd, nbdt)
end
function landau2_2(
    T::DataType,
    nbdt,
    timeopt;
    sz = (32, 32, 32, 32),
    dt = big"0.1",
    interpall = ntuple(x -> Lagrange(19, T), 4),
    split=strangsplit,
    tabst = [([3,4,1,2], 1,1, true), ([4,3,1,2], 1, 1, true),  ([1,2,4,3], 1, 2, true),([2,1,3,4], 1, 2, true)]
)
    epsilon = T(0.5)
    dt = T(dt)

    sp1min, sp1max, nsp1 = T(0), T(4big(pi)), sz[1]
    v1min, v1max, nv1 = -T(6.0), T(6.0), sz[3]

    mesh1_sp = UniformMesh(sp1min, sp1max, nsp1)
    mesh1_v = UniformMesh(v1min, v1max, nv1)

    sp2min, sp2max, nsp2 = T(0), T(4big(pi)), sz[2]
    v2min, v2max, nv2 = -T(6.0), T(6.0), sz[4]

    mesh2_sp = UniformMesh(sp2min, sp2max, nsp2)
    mesh2_v = UniformMesh(v2min, v2max, nv2)

    adv = Advection(
        (mesh1_sp, mesh2_sp, mesh1_v, mesh2_v),
        [interpall...],
        dt,
        tabst,
        tab_coef=split(dt),
        timeopt = timeopt
    )

    fct_sp(x) = epsilon * cos(x / 2) + 1
    fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))
    lgn1_sp = fct_sp.(mesh1_sp.points)
    lgn1_v = fct_v.(mesh1_v.points)
    lgn2_sp = fct_sp.(mesh2_sp.points)
    lgn2_v = fct_v.(mesh2_v.points)

    data = dotprod((lgn1_sp, lgn2_sp, lgn1_v, lgn2_v))

    pvar = getpoissonvar(adv)

    advd = AdvectionData(adv, data, pvar)
    # advdata = Advection1dData(adv, data, pvar)
    printout(
        advd,
        "# dt=$(Float32(dt)) eps=$(Float64(epsilon)) size1_sp=$nsp1 size2_sp=$nsp2 size_v1=$nv1 size_v2=$nv2",
    )
    printout(
        advd,
        "# sp1 : from $(Float64(start(mesh1_sp))) to $(Float64(stop(mesh1_sp)))",
    )
    printout(
        advd,
        "# sp2 : from $(Float64(start(mesh2_sp))) to $(Float64(stop(mesh2_sp)))",
    )
    printout(
        advd,
        "# v1 : from $(Float64(start(mesh1_v))) to $(Float64(stop(mesh1_v)))",
    )
    printout(
        advd,
        "# v2 : from $(Float64(start(mesh2_v))) to $(Float64(stop(mesh2_v)))",
    )
    printout(advd, "# interpolation : $interpall")
    printout(advd, "# tab_coef : $split")
    printout(advd, "# tabst : $tabst")
    printout(advd, "# type=$T precision = $(precision(T))")
    printout(advd, "# timeopt=$timeopt")
    if timeopt == SimpleThreadsOpt || timeopt == SplitThreadsOpt
        printout(advd, "# nb threads : $(Threads.nthreads())")
    elseif timeopt == MPIOpt
        printout(advd, "# nb process : $(adv.mpid.nb)")
    else
        printout(advd, "# monothread version")
    end
    printout(advd, "typeof(data)=$(typeof(data)) size(data)=$(size(data))")

    landau(advd, nbdt)
end
# landau2_2(Float64, 50, NoTimeOpt)
# landau2_2(Float64, 50, SimpleThreadsOpt)
# landau2_2(Float64, 50, SplitThreadsOpt)
# landau2_2(Float64, 50, MPIOpt, sz=(32,32,32,32))
# landau2_2(BigFloat, 10000, MPIOpt, sz=(64,64,64,64), dt=big"0.01")
# landau2_2(BigFloat, 10000, MPIOpt, sz=(32,32,32,32), dt=big"0.01")
T = Float64
# landau2_2(T, 10000, NoTimeOpt, sz=(32,32,32,32), dt=big"0.01", interp=B_SplineLU(27,32,T))
# @time landau2_2(T, 1000, NoTimeOpt, sz=(32,32,32,32), dt=big"0.1", interp=Lagrange(5, T))
# @time landau2_2(T, 30, NoTimeOpt, sz=(32,64,36,40), dt=big"0.1")
# sz = (32, 32, 20, 22)
# @time landau2_2(
#     T,
#     10,
#     MPIOpt,
#     sz = sz,
#     dt = big"0.1",
#     interpall = ntuple(x -> B_SplineLU(13, sz[x], T), 4),
# )
## @time landau2_2(T, 640, NoTimeOpt, sz=(32,32,128,128), dt=big"0.125", interp=Lagrange(5, T))

c=BigFloat(2)^(1//3)
c1 = 1/(2(2-c))
c2 = (1-c)/(2(2-c))
d1 = 1/(2-c)
d2 = -c/(2-c)
tc = [c1, d1, c2, d2, c2, d1, c1]
tc = [1, 1]
T=Float64
@time landau2_2(T, 5, NoTimeOpt, sz=(32,32,32,32), dt=T(0.01), interpall=ntuple(x->Lagrange(9,T),4),split=strangsplit)
@time landau2_2(T, 10, NoTimeOpt, sz=(32,32,32,32), dt=T(0.005), interpall=ntuple(x->Lagrange(9,T),4),split=strangsplit)
#  tabst = map(
#     x -> if x%2 != 1
#             ([1,2,3,4], 2, 1, true)
#         else # x%2 == 0
#             ([3,4,1,2], 2, 2, true)
#         end,
#     1:2
# ),
# tab_coef=tc
# )
# @time landau1_1(T, 10000, MPIOpt, sz=(256,256), dt=big"0.01",tab_coef=order6split())
#  tabst = map(
#     x -> if x%2 == 1
#             ([1,2,3,4], 2, x, true)
#         else # x%2 == 0
#             ([3,4,1,2], 2, x, true)
#         end,
#     1:3
# )

#@time landau2_2(T, 30, NoTimeOpt, sz=(32,64,36,40), dt=big"0.1", interpall=ntuple(x->Lagrange(7,T),4), tabst=tabst)
#@time landau2_2(T, 30, NoTimeOpt, sz=(32,64,36,40), dt=big"0.1", interpall=ntuple(x->LagrangeInt(7,T),4))
# landau1_1(T, 50, NoTimeOpt, sz=(64,128))
# landau2_2(T, 10000, MPIOpt, sz=(64,64,64,64), dt=big"0.01", interp=Lagrange(27, T))


text:
using Test, PseudoArcLengthContinuation, LinearAlgebra, SparseArrays, Setfield, Parameters
const PALC = PseudoArcLengthContinuation

f1(u, v) = u^2 * v
norminf = x -> norm(x, Inf)

function Fbru!(f, x, p)
	@unpack α, β, D1, D2, l = p
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h
	c1 = D1 / l^2 / h2
	c2 = D2 / l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	# Dirichlet boundary conditions
	f[1]   = c1 * (α	  - 2u[1] + u[2] ) + α - (β + 1) * u[1] + f1(u[1], v[1])
	f[end] = c2 * (v[n-1] - 2v[n] + β / α)			 + β * u[n] - f1(u[n], v[n])

	f[n]   = c1 * (u[n-1] - 2u[n] +  α  )  + α - (β + 1) * u[n] + f1(u[n], v[n])
	f[n+1] = c2 * (β / α  - 2v[1] + v[2])			 + β * u[1] - f1(u[1], v[1])

	for i=2:n-1
		  f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + α - (β + 1) * u[i] + f1(u[i], v[i])
		f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + β * u[i] - f1(u[i], v[i])
	end
	return f
end

function Fbru(x, p)
	f = similar(x)
	Fbru!(f, x, p)
end

function Jbru_sp(x, p)
	@unpack α, β, D1, D2, l = p
	# compute the Jacobian using a sparse representation
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h

	c1 = D1 / p.l^2 / h2
	c2 = D2 / p.l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	diag   = zeros(eltype(x), 2n)
	diagp1 = zeros(eltype(x), 2n-1)
	diagm1 = zeros(eltype(x), 2n-1)

	diagpn = zeros(eltype(x), n)
	diagmn = zeros(eltype(x), n)

	@. diagmn = β - 2 * u * v
	@. diagm1[1:n-1] = c1
	@. diagm1[n+1:end] = c2

	@. diag[1:n]    = -2c1 - (β + 1) + 2 * u * v
	@. diag[n+1:2n] = -2c2 - u * u

	@. diagp1[1:n-1] = c1
	@. diagp1[n+1:end] = c2

	@. diagpn = u * u
	J = spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
	return J
end

n = 100
par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
	sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))

opt_newton = PALC.NewtonPar(tol = 1e-11, verbose = true)
	out, hist, flag = @time PALC.newton(
		x -> Fbru(x, par_bru),
		x -> Jbru_sp(x, par_bru),
		sol0 .* (1 .+ 0.01rand(2n)),
		opt_newton)

eigls = EigArpack(1.1, :LM)
	opt_newton = PALC.NewtonPar(tol = 1e-11, verbose = false, linsolver = GMRESIterativeSolvers(tol=1e-4, N = 2n), eigsolver = eigls)
	out, hist, flag = @time PALC.newton(
		x -> Fbru(x, par_bru),
		x -> Jbru_sp(x, par_bru),
		sol0 .* (1 .+ 0.01rand(2n)),
		opt_newton)

opts_br0 = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds= 0.0051, pMax = 1.8, theta = 0.01, detectBifurcation = 2, nev = 16)
	br, u1 = @time PALC.continuation(
		(x, p) ->  Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		out,
		0.3,
		opts_br0,
		plot = false,
		printSolution = (x, p) -> norm(x, Inf64), verbosity = 0)
#################################################################################################### Continuation of the Hopf Point using Dense method
ind_hopf = 1
av = randn(Complex{Float64},2n); av = av./norm(av)
bv = randn(Complex{Float64},2n); bv = bv./norm(bv)
hopfpt = PALC.HopfPoint(br, ind_hopf)
# hopfpt[1:2n] .= rand(2n)
hopfvariable = β -> HopfProblemMinimallyAugmented(
					(x, p) ->  Fbru(x, setproperties(par_bru, (l=p, β=β))),
					(x, p) -> Jbru_sp(x, setproperties(par_bru, (l=p, β=β))),
					(x, p) -> transpose(Jbru_sp(x, setproperties(par_bru, (l=p, β=β)))),
					br.eig[br.bifpoint[ind_hopf].idx].eigenvec[:, br.bifpoint[ind_hopf].ind_bif],
					conj.(br.eig[br.bifpoint[ind_hopf].idx].eigenvec[:, br.bifpoint[ind_hopf].ind_bif]),
					# av,
					# bv,
					opts_br0.newtonOptions.linsolver)
	hopfPb = (u, p) -> hopfvariable(p)(u)

hopfvariable(par_bru.β)(hopfpt) |> norm

Bd2Vec(x) = vcat(x.u, x.p)
Vec2Bd(x) = BorderedArray(x[1:end-2], x[end-1:end])
hopfpbVec(x, p) = Bd2Vec(hopfvariable(p)(Vec2Bd(x)))

# finite differences Jacobian
Jac_hopf_fdMA(u0, p) = PALC.finiteDifferences( u-> hopfpbVec(u, p), u0)
# ``analytical'' jacobian
Jac_hopf_MA(u0, pb::HopfProblemMinimallyAugmented) = (return (u0, pb, x -> x))

# # on construit la matrice a inverser pour calculer sigma1 et sigma2
# ω = hopfpt[end]
# # av = hopfvariable(b).a
# # bv = hopfvariable(b).b
# Jh = hopfvariable(b).J(hopfpt[1:end-2], hopfpt[end-1])
# A = Jh + Complex(0, ω ) * I
# A = Array(A)
# Atot = [A av ; (bv') 0.]
# res1 = Atot \ [zeros(2n)' Complex(1,0)]'
# res2 = Atot' \ [zeros(2n)' Complex(1,0)]'
# println("--> σ1 = ", res1[end], ",\n--> σ2 = ", res2[end])
#
# dot(res2, Atot*res1)
#
# v, σ1, _ = PALC.linearBorderedSolver(Jh + Complex(0, ω) * I, av, bv, 0., zeros(2n), 1.0, hopfvariable(b).linsolve)
#
# w, σ2, _ = PALC.linearBorderedSolver(Jh' - Complex(0, ω) * I, bv, av, 0., zeros(2n), 1.0, hopfvariable(b).linsolve)
#
# res1[1:end-1]-v |> norm
# res2[1:end-1]-w |> norm
#
# L1 = hcat(Jh, -ω*I, real.(av), -imag.(av))
# L2 = hcat(ω*I, Jh, imag.(av), real.(av))
# L3 = hcat(real.(bv)', imag.(bv)', 0, 0)
# L4 = hcat(-imag.(bv)', real.(bv)', 0, 0)
# Atot2 = vcat(L1, L2, L3, L4)
# res12 = Atot2 \ hcat(zeros(4n)',[1],[0])'
# res22 = Atot2' \ hcat(zeros(4n)',[1],[0])'
#
# @test res12[1:2n] - real(res1[1:end-1]) |> norm < 1e-10
# @test res12[2n+1:4n] - imag(res1[1:end-1]) |> norm < 1e-10
# @test res22[1:2n] - real(res2[1:end-1]) |> norm < 1e-10
# @test res22[2n+1:4n] - imag(res2[1:end-1]) |> norm < 1e-10
#
# @test v - res1[1:end-1] |> norm < 1e-10
# @test w - res2[1:end-1] |> norm < 1e-10
#
# println("--> version reelle σ1 = ", res12[end-1:end])
# println("--> version reelle σ2 = ", res22[end-1:end])
#
# # on calcule les derivees de sigma1 et sigma2
# jac_hopf_fd = Jac_hopf_fdMA(hopfpt, b)
# println("-->\n FD sigma_omega = ", jac_hopf_fd[end-1:end,end-1:end])
# dot(res22[1:2n], res12[1:2n])
# dot(res22[2n+1:4n], res12[2n+1:4n])
#
# newton convergence toward

outhopf, _, flag, _ = @time PALC.newton(u -> hopfpbVec(u, par_bru.β),
							# u -> Jac_hopf_fdMA(u, par_bru.β),
							Bd2Vec(hopfpt), NewtonPar(verbose = true))
	flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf[end-1], ", ω = ", outhopf[end], " from ", hopfpt.p, "\n")

rhs = rand(length(hopfpt))
jac_hopf_fd = Jac_hopf_fdMA(Bd2Vec(hopfpt), par_bru.β)
sol_fd = jac_hopf_fd \ rhs
# create a linear solver
hopfls = HopfLinearSolveMinAug()
sol_ma, _, _, sigomMA  = hopfls(Jac_hopf_MA(hopfpt, hopfvariable(par_bru.β)), BorderedArray(rhs[1:end-2],rhs[end-1:end]), debug_ = true)

# TODO TODO fix these two lines
println("--> test jacobian expression for Hopf Minimally Augmented")
@test Bd2Vec(sol_ma) - sol_fd |> x-> norm(x, Inf64) < 1e3

@test (Bd2Vec(sol_ma) - sol_fd)[1:end-2] |> x->norm(x, Inf64) < 1e-1

# dF = jac_hopf_fd[:,end-1]
# sig_vec_re = jac_hopf_fd[end-1,1:end-2]
# sig_vec_im = jac_hopf_fd[end,1:end-2]
#
# println("-->\n FD sigma_omega = ", jac_hopf_fd[end-1:end,end-1:end]')
#
# norm(sig_vec_re - real.(sigomMA[1]), Inf64)
# norm(sig_vec_im - imag.(sigomMA[1]), Inf64)
# norm(dF[1:end-2] - sigomMA[end], Inf64)
#
# # on essaie de resoudre jac_hopf_fd \ rhs
# Jh = jac_hopf_fd[1:end-2,1:end-2]
# sig1 = jac_hopf_fd[end-1,1:end-2]
# sig2 = jac_hopf_fd[end,1:end-2]
# sig1p = jac_hopf_fd[end-1,end-1]
# sig2p = jac_hopf_fd[end,end-1]
# sig1o = jac_hopf_fd[end-1,end]
# sig2o = jac_hopf_fd[end,end]
#
# sigma = hcat(sig1,sig2)
#
# X1 = Jh \ rhs[1:end-2]
# X2 = Jh \ dF[1:end-2]#jac_hopf_fd[1:end-2,end-2]
# X2m = hcat(X2, 0*X2)
# C = jac_hopf_fd[end-1:end,end-1:end]
# println("--> dp, dom = ",(C - sigma' * X2m) \ (rhs[end-1:end] - sigma' * X1))
# println("--> dp, dom FD = ",sol_fd[end-1:end])
outhopf, hist, flag = @time PALC.newtonHopf(
		(x, p) ->  Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		(x, p) -> transpose(Jbru_sp(x, @set par_bru.l = p)),
		br, 1,
		NewtonPar(verbose = true))
		flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[end], ", from l = ",hopfpt.p[1],"\n")

outhopf, _, flag, _ = @time PALC.newton(u -> hopfvariable(par_bru.β)(u),
							x -> Jac_hopf_MA(x, hopfvariable(par_bru.β)),
							hopfpt, NewtonPar(verbose = true, linsolver = HopfLinearSolveMinAug()))
	flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[2], " from ", hopfpt.p, "\n")

# version with analytical Hessian = 2 P(du2) P(du1) QU + 2 PU P(du1) Q(du2) + 2PU P(du2) Q(du1)
function d2F(x, p1::T, p2::T, du1, du2) where T
	n = div(length(x),2)
	out = similar(du1)
	@views out[1:n] .= 2 .* x[n+1:end] .* du1[1:n] .* du2[1:n] .+
				2 .* x[1:n] .* du1[1:n] .* du2[n+1:end] .+
				2 .* x[1:n] .* du2[1:n] .* du1[1:n]
	@inbounds for ii=1:n
		out[ii+n] = -out[ii]
	end
	return out
end

outhopf, hist, flag = @time PALC.newtonHopf(
		(x, p) ->  Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		(x, p) -> transpose(Jbru_sp(x, @set par_bru.l = p)),
		(x, p1, v1, v2) -> d2F(x, 0., 0., v1, v2),
		br, 1,
		NewtonPar(verbose = true))
		flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[end], ", from l = ",hopfpt.p[1],"\n")

br_hopf, u1_hopf = @time PALC.continuationHopf(
			(x, p, β) ->  Fbru(x, setproperties(par_bru, (l=p, β=β))),
			(x, p, β) ->  Jbru_sp(x, setproperties(par_bru, (l=p, β=β))),
			br, ind_hopf,
			par_bru.β,
			ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, a = 2., theta = 0.4, maxSteps = 3, newtonOptions = NewtonPar(verbose = false)), verbosity = 1, plot = false)

br_hopf, u1_hopf = @time PALC.continuationHopf(
			(x, p, β) ->  Fbru(x, setproperties(par_bru, (l=p, β=β))),
			(x, p, β) ->  Jbru_sp(x, setproperties(par_bru, (l=p, β=β))),
			(x, p, β) ->  transpose(Jbru_sp(x, setproperties(par_bru, (l=p, β=β)))),
			p2 -> (x, p1, v1, v2) -> d2F(x, 0., 0., v1, v2),
			br, ind_hopf,
			par_bru.β,
			ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, a = 2., theta = 0.4, maxSteps = 3, newtonOptions = NewtonPar(verbose = false)), verbosity = 0, plot = false)
#################################################################################################### Continuation of Periodic Orbit
ind_hopf = 1
hopfpt = PALC.HopfPoint(br, ind_hopf)

l_hopf  = hopfpt.p[1]
ωH		= hopfpt.p[2] |> abs
M = 20

orbitguess = zeros(2n, M)
phase = []; scalphase = []
vec_hopf = geteigenvector(opt_newton.eigsolver, br.eig[br.bifpoint[ind_hopf].idx][2], br.bifpoint[ind_hopf].ind_bif-1)
for ii=1:M
	t = (ii-1)/(M-1)
	orbitguess[:, ii] .= real.(hopfpt.u + 26*0.1 * vec_hopf * exp(-2pi * complex(0, 1) * (t - .252)))
end

orbitguess_f = vcat(vec(orbitguess), 2pi/ωH) |> vec

# test guess using function
l_hopf, Th, orbitguess2, hopfpt, vec_hopf = PALC.guessFromHopf(br, ind_hopf, opt_newton.eigsolver, M, 2.6; phase = 0.252)

poTrap = p-> PeriodicOrbitTrapProblem(
			x ->	Fbru(x, @set par_bru.l = p),
			x -> Jbru_sp(x, @set par_bru.l = p),
			real.(vec_hopf),
			hopfpt.u,
			M)

jac_PO_fd = PALC.finiteDifferences(x -> poTrap(l_hopf + 0.01)(x), orbitguess_f)
jac_PO_sp =  poTrap(l_hopf + 0.01)(Val(:JacFullSparse), orbitguess_f)

# test of the Jacobian for PeriodicOrbit via Finite differences VS the FD associated jacobian
println("--> test jacobian expression for Periodic Orbit solve problem")
@test norm(jac_PO_fd - jac_PO_sp, Inf64) < 1e-4

# newton to find Periodic orbit
opt_po = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 150)
	outpo_f, _, flag = @time PALC.newton(
		x ->  poTrap(l_hopf + 0.01)(x),
		x ->  poTrap(l_hopf + 0.01)(Val(:JacFullSparse),x),
		orbitguess_f,
		opt_po)
	println("--> T = ", outpo_f[end])
flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", PALC.amplitude(outpo_f, n, M; ratio = 2),"\n")

# jacobian of the functional
Jpo2 = poTrap(l_hopf + 0.01)(Val(:JacCyclicSparse), orbitguess_f)

# calcul des exposants de Floquet
floquetES = PALC.FloquetQaDTrap(DefaultEig())

# continuation of periodic orbits using :BorderedLU linear algorithm
opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.001, pMax = 2.3, maxSteps = 3, theta = 0.1, newtonOptions = NewtonPar(verbose = false), detectBifurcation = 1)
	br_pok2, upo , _= @time PALC.continuationPOTrap(
		p ->  poTrap(p),
		outpo_f, l_hopf + 0.01,
		opts_po_cont, :BorderedLU;
		plot = false,
		verbosity = 0)

# test of simple calls to newton / continuation
deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zero(orbitguess_f)])
opt_po = PALC.NewtonPar(tol = 1e-8, verbose = true, maxIter = 10)
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01, pMax = 3.0, maxSteps = 3, newtonOptions = (@set opt_po.verbose = false), computeEigenValues = true, nev = 2, precisionStability = 1e-8, detectBifurcation = 1)
for linalgo in [:FullLU, :BorderedLU]
	@show linalgo
	outpo_f, hist, flag = @time PALC.newton(poTrap(l_hopf + 0.01),
			orbitguess_f, opt_po, deflationOp, linalgo; normN = norminf)
	outpo_f, hist, flag = @time PALC.newton(poTrap(l_hopf + 0.01),
			orbitguess_f, opt_po, linalgo; normN = norminf)
	br_pok2, upo , _= @time PALC.continuationPOTrap(poTrap,
			outpo_f, l_hopf + 0.01,
			opts_po_cont; verbosity = 0,
			plot = false, normC = norminf)
end

ext: jl


text:
using pylib::FileWriteString
using std::fs::File
using std::fs::OpenOptions
# Code to construct puncture initial data for single black hole.

using numpy: zeros, size, sqrt, linspace
import scipy.linalg
struct EllipticSolver
n_grid::
delta::
rhs_1d::
A::
sol::
rad::
end

function __init__{T0, T1, T2}(self::EllipticSolver, x::T0, y::T1, z::T2)
println(join([" Setting up Poisson solver..."], " "));
self.n_grid = size(x)
self.delta = (x[1] - x[0])
nnn = pow(self.n_grid, 3)
self.rhs_1d = zeros(nnn)
self.A = zeros((nnn, nnn))
self.sol = zeros((self.n_grid, self.n_grid, self.n_grid))
self.rad = zeros((self.n_grid, self.n_grid, self.n_grid))
for i in (0:self.n_grid - 1)
for j in (0:self.n_grid - 1)
for k in (0:self.n_grid - 1)
rad2 = ((pow(x[i], 2) + pow(y[j], 2)) + pow(z[k], 2))
self.rad[(i, j, k)] = sqrt(rad2)
end
end
end
end

function setup_matrix{T0}(self::EllipticSolver, fct::T0)
n_grid = self.n_grid
i = 0
for j in (0:n_grid - 1)
for k in (0:n_grid - 1)
index = super_index(self, convert(, i), j, k)
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index + 1))] = -(self.rad[((i + 1), j, k)])
end
end
i = (n_grid - 1);
for j in (0:n_grid - 1)
for k in (0:n_grid - 1)
index = super_index(self, convert(, i), j, k)
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index - 1))] = -(self.rad[((i - 1), j, k)])
end
end
j = 0
for i in (1:(n_grid - 1) - 1)
for k in (0:n_grid - 1)
index = super_index(self, i, convert(, j), k)
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index + n_grid))] = -(self.rad[(i, (j + 1), k)])
end
end
j = (n_grid - 1);
for i in (1:(n_grid - 1) - 1)
for k in (0:n_grid - 1)
index = super_index(self, i, convert(, j), k)
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index - n_grid))] = -(self.rad[(i, (j - 1), k)])
end
end
k = 0
for i in (1:(n_grid - 1) - 1)
for j in (1:(n_grid - 1) - 1)
index = super_index(self, i, j, convert(, k))
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index + (n_grid*n_grid)))] = -(self.rad[(i, j, (k + 1))])
end
end
k = (n_grid - 1);
for i in (1:(n_grid - 1) - 1)
for j in (1:(n_grid - 1) - 1)
index = super_index(self, i, j, convert(, k))
self.A[(index, index)] = self.rad[(i, j, k)]
self.A[(index, (index - (n_grid*n_grid)))] = -(self.rad[(i, j, (k - 1))])
end
end
for i in (1:(n_grid - 1) - 1)
for j in (1:(n_grid - 1) - 1)
for k in (1:(n_grid - 1) - 1)
index = super_index(self, i, j, k)
self.A[(index, index)] = (-6.0 + (pow(self.delta, 2)*fct[(i, j, k)]))
self.A[(index, (index - 1))] = 1.0
self.A[(index, (index + 1))] = 1.0
self.A[(index, (index - n_grid))] = 1.0
self.A[(index, (index + n_grid))] = 1.0
self.A[(index, (index - (n_grid*n_grid)))] = 1.0
self.A[(index, (index + (n_grid*n_grid)))] = 1.0
end
end
end
end

function setup_rhs{T0}(self::EllipticSolver, rhs::T0)
n_grid = self.n_grid
for i in (1:(n_grid - 1) - 1)
for j in (1:(n_grid - 1) - 1)
for k in (1:(n_grid - 1) - 1)
index = super_index(self, i, j, k)
self.rhs_1d[index] = (pow(self.delta, 2)*rhs[(i, j, k)])
end
end
end
end

function solve{RT}(self::EllipticSolver)::RT
sol_1d = solve(la)
for i in (0:self.n_grid - 1)
for j in (0:self.n_grid - 1)
for k in (0:self.n_grid - 1)
index = super_index(self, i, j, k)
self.sol[(i, j, k)] = sol_1d[index]
end
end
end
return self.sol
end

function super_index{T0, T1, T2, RT}(self::EllipticSolver, i::T0, j::T1, k::T2)::RT
return (i + (self.n_grid*(j + (self.n_grid*k))))
end

struct Puncture
bh_loc::
lin_mom::
n_grid::
x_out::
delta::Float64
x::
y::
z::
solver::EllipticSolver
alpha::
beta::
u::
res::
end

function __init__{T0, T1, T2, T3}(self::Puncture, bh_loc::T0, lin_mom::T1, n_grid::T2, x_out::T3)
self.bh_loc = bh_loc
self.lin_mom = lin_mom
println(join([" Constructing class Puncture for single black hole"], " "));
println(join(["    at bh_loc = (", bh_loc[0], ",", bh_loc[1], ",", bh_loc[2], ")"], " "));
println(join(["    with momentum p = (", lin_mom[0], ",", lin_mom[1], ",", lin_mom[2], ")"], " "));
println(join([" Using", n_grid, "^3 gridpoints with outer boundary at", x_out], " "));
self.n_grid = n_grid
self.x_out = x_out
self.delta = ((2.0*x_out)/n_grid)
half_delta = (self.delta/2.0)
self.x = linspace((half_delta - x_out), (x_out - half_delta), n_grid)
self.y = linspace((half_delta - x_out), (x_out - half_delta), n_grid)
self.z = linspace((half_delta - x_out), (x_out - half_delta), n_grid)
self.solver = EllipticSolver(self.x, self.y, self.z)
self.alpha = zeros((n_grid, n_grid, n_grid))
self.beta = zeros((n_grid, n_grid, n_grid))
self.u = zeros((n_grid, n_grid, n_grid))
self.res = zeros((n_grid, n_grid, n_grid))
end

function construct_solution{T0, T1}(self::Puncture, tol::T0, it_max::T1)
setup_alpha_beta(self);
residual_norm = residual(self)
println(join([" Initial Residual = ", residual_norm], " "));
println(join([" Using up to", it_max, "iteration steps to reach tolerance of", tol], " "));
it_step = 0
while residual_norm > tol&&it_step < it_max
it_step += 1
update_u(self);
residual_norm = residual(self);
println(join([" Residual after", it_step, "iterations :", residual_norm], " "));
end
if residual_norm < tol
println(join([" Done!"], " "));
else

println(join([" Giving up..."], " "));
end
end

function update_u(self::Puncture)
n_grid = self.n_grid
fct = zeros((n_grid, n_grid, n_grid))
rhs = zeros((n_grid, n_grid, n_grid))
for i in (1:(n_grid - 1) - 1)
for j in (1:(n_grid - 1) - 1)
for k in (1:(n_grid - 1) - 1)
temp = ((self.alpha[(i, j, k)]*(1.0 + self.u[(i, j, k)])) + 1.0)
fct[(i, j, k)] = (((-7.0*self.beta[(i, j, k)])*self.alpha[(i, j, k)])/pow(temp, 8))
rhs[(i, j, k)] = -(self.res[(i, j, k)])
end
end
end
setup_matrix(self.solver, fct);
setup_rhs(self.solver, rhs);
delta_u = solve(self.solver)
self.u += delta_u
end

function residual(self::Puncture)::Float64
residual_norm = 0.0
for i in (1:(self.n_grid - 1) - 1)
for j in (1:(self.n_grid - 1) - 1)
for k in (1:(self.n_grid - 1) - 1)
ddx = ((self.u[((i + 1), j, k)] - (2.0*self.u[(i, j, k)])) + self.u[((i - 1), j, k)])
ddy = ((self.u[(i, (j + 1), k)] - (2.0*self.u[(i, j, k)])) + self.u[(i, (j - 1), k)])
ddz = ((self.u[(i, j, (k + 1))] - (2.0*self.u[(i, j, k)])) + self.u[(i, j, (k - 1))])
lhs = (((ddx + ddy) + ddz)/pow(self.delta, 2))
temp = ((self.alpha[(i, j, k)]*(1.0 + self.u[(i, j, k)])) + 1.0)
rhs = (-(self.beta[(i, j, k)])/pow(temp, 7))
self.res[(i, j, k)] = (lhs - rhs)
residual_norm += pow(self.res[(i, j, k)], 2)
end
end
end
residual_norm = (sqrt(residual_norm)*pow(self.delta, 3));
return residual_norm
end

function setup_alpha_beta(self::Puncture)
n_grid = self.n_grid
p_x = self.lin_mom[0]
p_y = self.lin_mom[1]
p_z = self.lin_mom[2]
for i in (0:n_grid - 1)
for j in (0:n_grid - 1)
for k in (0:n_grid - 1)
s_x = (self.x[i] - self.bh_loc[0])
s_y = (self.y[j] - self.bh_loc[1])
s_z = (self.z[k] - self.bh_loc[2])
s2 = ((pow(s_x, 2) + pow(s_y, 2)) + pow(s_z, 2))
s_bh = sqrt(s2)
l_x = (s_x/s_bh)
l_y = (s_y/s_bh)
l_z = (s_z/s_bh)
lP = (((l_x*p_x) + (l_y*p_y)) + (l_z*p_z))
fac = (3.0/(2.0*s2))
A_xx = (fac*(((2.0*p_x)*l_x) - ((1.0 - (l_x*l_x))*lP)))
A_yy = (fac*(((2.0*p_y)*l_y) - ((1.0 - (l_y*l_y))*lP)))
A_zz = (fac*(((2.0*p_z)*l_z) - ((1.0 - (l_z*l_z))*lP)))
A_xy = (fac*(((p_x*l_y) + (p_y*l_x)) + ((l_x*l_y)*lP)))
A_xz = (fac*(((p_x*l_z) + (p_z*l_x)) + ((l_x*l_z)*lP)))
A_yz = (fac*(((p_y*l_z) + (p_z*l_y)) + ((l_y*l_z)*lP)))
A2 = (((pow(A_xx, 2) + pow(A_yy, 2)) + pow(A_zz, 2)) + (2.0*((pow(A_xy, 2) + pow(A_xz, 2)) + pow(A_yz, 2))))
self.alpha[(i, j, k)] = (2.0*s_bh)
self.beta[(i, j, k)] = ((pow(self.alpha[(i, j, k)], 7)*A2)/8.0)
end
end
end
end

function write_to_file(self::Puncture)
filename = ((("Puncture_" + string(self.n_grid)) + "_") + string(self.x_out))
filename = (filename + ".data");
out = OpenOptions::new().write(true).open(filename)
if out
k = (self.n_grid / 2)
write(out, ("# Data for black hole at x = (%f,%f,%f)
" % (self.bh_loc[0], self.bh_loc[1], self.bh_loc[2])));
write(out, ("# with linear momentum P = (%f, %f, %f)
" % self.lin_mom));
write(out, ("# in plane for z = %e
" % self.z[k]));
write(out, "# x            y              u
");
write(out, "#============================================
");
for i in (0:self.n_grid - 1)
for j in (0:self.n_grid - 1)
write(out, ("%e  %e  %e
" % (self.x[i], self.y[j], self.u[(i, j, k)])));
end
write(out, "
");
end
close(out);
else

println(join([" Could not open file", filename, "in write_to_file()"], " "));
println(join([" Check permissions?"], " "));
end
end

function main_func()
println(join([" -------------------------------------------------------"], " "));
println(join([" --- puncture.py --- use flag -h for list of options ---"], " "));
println(join([" -------------------------------------------------------"], " "));
loc_x = 0.0
loc_y = 0.0
loc_z = 0.0
p_x = 1.0
p_y = 0.0
p_z = 0.0
n_grid = 16
x_out = 4.0
tol = 1e-12
it_max = 50
for i in (0:length(append!([PROGRAM_FILE], ARGS)) - 1)
if append!([PROGRAM_FILE], ARGS)[i] == "-h"
usage();
return
end
if append!([PROGRAM_FILE], ARGS)[i] == "-n_grid"
n_grid = Int64(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-x_out"
x_out = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-loc_x"
loc_x = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-loc_y"
loc_y = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-loc_z"
loc_z = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-p_x"
p_x = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-p_y"
p_y = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-p_z"
p_z = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-tol"
tol = float(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
if append!([PROGRAM_FILE], ARGS)[i] == "-it_max"
it_max = Int64(append!([PROGRAM_FILE], ARGS)[(i + 1)]);
end
end
bh_loc = (loc_x, loc_y, loc_z)
lin_mom = (p_x, p_y, p_z)
black_hole = Puncture(bh_loc, lin_mom, n_grid, x_out)
construct_solution(black_hole, tol, it_max);
write_to_file(black_hole);
end

function usage()
println(join(["Constructs puncture initial data for single black hole."], " "));
println(join([""], " "));
println(join(["The following options can be used to over-write default parameters"], " "));
println(join(["	-n_grid: number of grid points [default: 16]"], " "));
println(join(["	-x_out: location of outer boundary [4.0]"], " "));
println(join(["	-loc_x, -loc_y, -loc_z: location of black hole [(0.0, 0.0, 0.0)]"], " "));
println(join(["	-p_x, -p_y, -p_z: lin. momentum of black hole [(1.0, 0.0, 0.0)]"], " "));
println(join(["	-tol: tolerance for elliptic solver [1.e-12]"], " "));
println(join(["	-it_max: maximum number of iterations [50]"], " "));
println(join(["For example, to construct data with x_out = 6.0, call"], " "));
println(join(["	python3 puncture.py -x_out 6.0"], " "));
end

function main()
main_func();
end

main()


text:
@everywhere function mark_graph!(graph, mark, k, islice, jslice)
    for i in islice, j in jslice
        gs = min(graph[i,k],graph[k,j])
        gij = graph[i,j]
        if  gs != 0 && gij != 0 && gij <= gs
            graph[i,j] = gs
            mark[i,j] = true
        end
    end
end

@everywhere function change_graph!(graph, mark, islice, jslice)
    for i in islice, j in jslice
        if mark[i,j]
            graph[i,j] = 0
        end
    end
end

@everywhere function start_kernel_mark!(graph, mark, k)
    mark_graph!(graph, mark, k, 1:size(graph,1), 1:size(graph,2))
end

@everywhere function start_kernel_change!(graph, mark, interval_split)
    idx = indexpids(graph)
    change_graph!(graph, mark, interval_split[idx]+1:interval_split[idx+1], 1:size(graph,2))
end

function transitive_reduction_parallel!(graph::SharedArray{Int8,2}, markedgraph::SharedArray{Bool})
    graphsize = size(graph)

    # processes
    procs_list = procs(graph)
    # graph split
    interval_split = [round(Int, s) for s in linspace(0,size(graph,1),length(procs_list)+1)]

    @sync begin
        for k in 1:graphsize[1]
            @async remotecall_wait(start_kernel_mark!, procs_list[(k%length(procs_list))+1], graph, markedgraph, k)
        end
    end

    @sync begin
        for p in procs_list
            @async remotecall_wait(start_kernel_change!, p, graph, markedgraph, interval_split)
        end
    end
end

function transitive_reduction_parallel(graphin::Array{Int8,2})
    graph = SharedArray{Int8,2}(copy(graphin))
    markedgraph = SharedArray{Bool}(size(graph), init = S -> S[Base.localindexes(S)] = false)

    # Time without alloc here
    transitive_reduction_parallel!(graph, markedgraph)

    return graph
end

function transitive_reduction_serial(graphin::Array{Int8,2})
    graph = Array{Int8,2}(copy(graphin))
    n = size(graphin,1)
    for j in 1:n
        for i in 1:n
            if graph[i, j] > 0
                for k in 1:n
                    if graph[j, k] > 0
                        graph[i, k] = 0
                    end
                end
            end
        end
    end
end

containment_graph_little = Array{Int8,2}([
      [0 1 1 1 1]
      [0 0 1 0 0]
      [0 0 0 0 0]
      [0 0 0 0 1]
      [0 0 0 0 0]])

containment_graph = Array{Int8,2}([
      [0 0 0 0 0 0 0 0]
      [0 0 0 0 0 0 1 1]
      [1 1 0 1 1 1 1 1]
      [0 0 0 0 0 0 0 0]
      [0 0 0 0 0 1 0 0]
      [0 0 0 0 0 0 0 0]
      [0 0 0 0 0 0 0 0]
      [0 0 0 0 0 0 1 0]])

# Run once to warm up the GC
transitive_reduction_serial(ones(Int8,10,10))
transitive_reduction_parallel(ones(Int8,10,10))

# Test
@time res1 = transitive_reduction_parallel(containment_graph_little)
@time res2 = transitive_reduction_parallel(containment_graph)
@time res3 = transitive_reduction_parallel(ones(Int8,10,10))
@time res4 = transitive_reduction_parallel(ones(Int8,100,100))
@time res5 = transitive_reduction_parallel(ones(Int8,1000,1000))
# @time res5 = transitive_reduction_parallel(ones(Int8,10000,10000))

@time sres1 = transitive_reduction_serial(containment_graph_little)
@time sres2 = transitive_reduction_serial(containment_graph)
@time sres3 = transitive_reduction_serial(ones(Int8,10,10))
@time sres4 = transitive_reduction_serial(ones(Int8,100,100))
@time sres5 = transitive_reduction_serial(ones(Int8,1000,1000))

# Verify results
println(res1)
println(res2)
# println(res3)
# println(res4)

max_stars_repo_name: furio/lar-parallel-tests

max_issues_repo_name: furio/lar-parallel-tests

max_forks_repo_name: furio/lar-parallel-tests


text:
module Visualize_Solids

using Modia3D
using Modia3D.ModiaMath
using Modia3D.Unitful

filename = joinpath(Modia3D.path, "objects", "pascal", "pascal.obj")

# Properties
vmat1 = Modia3D.Material(color="LightBlue", transparency=0.5, wireframe=false)
vmat2 = deepcopy(vmat1)
vmat2.transparency = 0.7
smat1 = Modia3D.MassProperties(m=1.0)
font  = Modia3D.Font(fontFamily="TimesNewRoman", color="Black", bold=false, charSize=0.2)
rtext = [0.0, 0.0, 1.0]
ltext = [0.0, 0.0, -0.5]
Rmesh = ModiaMath.rot123(90u"°", 180u"°", 0.0)

# Solid objects
@assembly Solids begin
   world     = Modia3D.Object3D(visualizeFrame=false)
   box       = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(0.9,0.5,0.3)      , "Aluminium", vmat1); r=[ 4.5,0.0,0.0])
   sphere    = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(0.7)           , "Aluminium", vmat1); r=[ 3.0,0.0,0.0])
   cylinder  = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCylinder(0.5,0.8)     , smat1      , vmat1); r=[ 1.5,0.0,0.0])
   cone      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCone(0.3,0.7)         , smat1      , vmat1); r=[ 0.0,0.0,0.0])
   capsule   = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCapsule(0.4,0.45)     , smat1      , vmat1); r=[-1.5,0.0,0.0])
   pipe      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidPipe(0.5,0.8)         , smat1      , vmat1); r=[ 4.5,0.0,-2.5])
   beam      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBeam(0.4,0.5,0.3)     , smat1      , vmat1); r=[ 3.0,0.0,-2.5])
   fileMesh1 = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidFileMesh(filename,0.2), smat1      , vmat1);  R=Rmesh, r=[ 1.5,0.0,-2.5])
   fileMesh2 = Modia3D.Object3D(world, Modia3D.SolidWithConvexDecomposition(Modia3D.SolidFileMesh(filename,0.2), smat1, vmat1, vmat2); r=[-1.5,0.0,-2.5], R=Rmesh)

   # coord Systems
   coordBox     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 4.5,0.0,0.0])
   coordSphere  = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 3.0,0.0,0.0])
   coordCyl     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 1.5,0.0,0.0])
   coordCone    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 0.0,0.0,0.0])
   coordCap     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ -1.5,0.0,0.0])
   coordpipe    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 4.5,0.0,-2.5])
   coordBeam    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 3.0,0.0,-2.5])
   coordMesh1   = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 1.5,0.0,-2.5])
   coordMesh2   = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ -1.5,0.0,-2.5])

   # Place text above the shapes
   boxText       = Modia3D.Object3D(box      , Modia3D.TextShape("SolidBox"     ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   sphereText    = Modia3D.Object3D(sphere   , Modia3D.TextShape("SolidSphere"  ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   cylinderText  = Modia3D.Object3D(cylinder , Modia3D.TextShape("SolidCylinder"; font=font); r=rtext, visualizeFrame=Modia3D.False)
   coneText      = Modia3D.Object3D(cone     , Modia3D.TextShape("SolidCone"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   capsuleText   = Modia3D.Object3D(capsule  , Modia3D.TextShape("SolidCapsule" ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   pipeText      = Modia3D.Object3D(pipe     , Modia3D.TextShape("SolidPipe"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   beamText      = Modia3D.Object3D(beam     , Modia3D.TextShape("SolidBeam"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   meshText1     = Modia3D.Object3D(coordMesh1, Modia3D.TextShape("SolidFileMesh"; font=font); r=rtext, visualizeFrame=Modia3D.False)
   meshText2     = Modia3D.Object3D(coordMesh2, Modia3D.TextShape("SolidFileMesh with convex decomposition"; font=font); r=rtext, visualizeFrame=Modia3D.False)
   meshText3     = Modia3D.Object3D(coordMesh2, Modia3D.TextShape("(.obj)"        ; font=font); r=ltext, visualizeFrame=Modia3D.False)
   meshText4     = Modia3D.Object3D(coordMesh1, Modia3D.TextShape("(.obj)"        ; font=font); r=ltext, visualizeFrame=Modia3D.False)
end

Modia3D.visualizeAssembly!( Solids() )

println("... success of Visualize_Solids.jl!")
end


text:
using LinearAlgebra

mutable struct materialsData{T}
    youngModulus::T
    poissonRatio::T
end

struct elementsData
    type::String
    crossSectionArea::Float64
end

mutable struct boundaryCondsData
    DirichletNodalDOFs::Vector{Int}
    DirichletNodalVals::Vector{Float64}
    NeumannNodalDOFs::Vector{Int}
    NeumannNodalVals::Vector{Float64}
end

struct BCsIndexes
    diriDofs::Vector{Int}
    neumDofs::Vector{Int}
end

struct systemMatrix{T}
    matrix::T
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function calc_dt_for_cleaning_speed(cfl::Real, mesh,
                                    equations::Union{AbstractIdealGlmMhdEquations, AbstractIdealGlmMhdMulticomponentEquations}, dg::DG, cache)
# compute time step for GLM linear advection equation with c_h=1 for the DG discretization on
# Cartesian meshes
  max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(equations)
  # OBS! This depends on the implementation details of the StepsizeCallback and needs to be adapted
  #      as well if that callback changes.
  return cfl * 2 / (nnodes(dg) * max_scaled_speed_for_c_h)
end

function calc_dt_for_cleaning_speed(cfl::Real, mesh,
                                    equations::Union{AbstractIdealGlmMhdEquations, AbstractIdealGlmMhdMulticomponentEquations},
                                    dg::DGMulti, cache)
  rd = dg.basis
  md = mesh.md

  # Compute time step for GLM linear advection equation with c_h=1 for a DGMulti discretization.
  # Copies implementation behavior of `calc_dt_for_cleaning_speed` for DGSEM discretizations.
  max_scaled_speed_for_c_h = (1 / minimum(md.J)) * ndims(equations)

  # This mimics `max_dt` for `TreeMesh`, except that `nnodes(dg)` is replaced by
  # `polydeg+1`. This is because `nnodes(dg)` returns the total number of
  # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns
  # the number of 1D nodes for `DGSEM` solvers.
  polydeg = rd.N
  return cfl * 2 / ((polydeg + 1) * max_scaled_speed_for_c_h)
end

end # @muladd


text:
include("src/CellularAutomata.jl")
using Random

init_config = bitrand(10)

automaton = CellularAutomata.CellularAutomaton(init_config)

for i in 1:50
    CellularAutomata.step!(automaton, 30)
end

CellularAutomata.plot(automaton)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Unicorn"
version = v"1.0.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/oblivia-simplex/unicorn.git", "cda971566980d3a051a730f63019e74eeeb44703")
]

# Bash recipe for building across all platforms
script = raw"""
ln -s /usr/bin/make /usr/bin/gmake
cd ${WORKSPACE}/srcdir/unicorn
make -j $(nproc)
make PREFIX=${prefix} install
install_license COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libunicorn", :libunicorn)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
"Your optional docstring here"
function distance(a, b)

end


text:
""" Functions that write and read 2OKP instances on files """

# All functions are public

"""
Returns the name of the file associated with the given problem `prob` and `id`

`id`, unique number to avoid writing on existing files (generally between 1 and 5)
"""
function FILES_getProbName(prob::BiOKP,id::Int;isInteger=true)
    if isInteger
        return string("Inst_",string(prob.nbVar),"_",string(id))
    else
        return string("Inst_",string(prob.nbVar),"_",string(id),"_float")
    end
end

"""
Write the given problem `prob` attributes on a file

`id`, unique number to avoid writing on existing files (generally between 1 and 5)

# Example
```jldoctest
julia> FILES_writeProblem(myProblem,10)
```
"""
function FILES_writeProblem(prob::BiOKP,id::Int; isInteger = true)
    fileName = FILES_getProbName(prob,id,isInteger=isInteger)

    f = open(string("instances/",fileName,".dat"),"w")

    println(f,string(prob.nbVar," ",prob.nbObj))
    for i in 1:2
        for profit in prob.profits[i,:]
            print(f,string(profit," "))
        end
        println(f,"")
    end
    for weight in prob.weights
        print(f,string(weight," "))
    end
    println(f,"")
    println(f,string(prob.maxWeight))
    println(f,string(prob.isInteger))

    close(f)
end

"""
Returns the problem stored in the given file

The `fileName` is the name of the given file (no extension, no path, just the name)

# Example
```jldoctest
julia> p = FILES_readProblem("Inst_5_1")
BiOKP(5, 2, [29.0 28.0 … 18.0 31.0; 11.0 36.0 … 27.0 37.0], [16.0, 50.0, 32.0, 13.0, 50.0], 65.0, true)
```
"""
function FILES_readProblem(fileName::String)

    f = open(string("instances/",fileName,".dat"),"r")

    first_line::Array{SubString{String},1} = split(readline(f))
    nbVars = parse(Int64, first_line[1])
    nbObj = parse(Int64, first_line[2])
    profits_1 = split(readline(f))
    profits_2 = split(readline(f))
    profits = zeros(Float64, 2, nbVars)
    weights_str = split(readline(f))
    weights = zeros(Float64,nbVars)
    for i in 1:nbVars
        profits[1,i] = parse(Float64,profits_1[i])
        profits[2,i] = parse(Float64,profits_2[i])
        weights[i] = parse(Float64,weights_str[i])
    end
    maxWeight = parse(Float64,readline(f))
    isInteger = parse(Bool,readline(f))

    return BiOKP(nbVars,nbObj,profits,weights,maxWeight,isInteger)
end


text:
using PyCall
@pyimport numpy
X = numpy.load("test_x.npy");
Y = numpy.load("test_y.npy");

# Get the predicted value
round_output(x) = x[1]>=0.5?1:0

using Keras

#Change this value to test results on smaller datasets
num_tests = 25000

println("Loading model...")
model = Keras.load("model_structure.json", "model_weight.h5")
println("Model loaded successfully.")
# Test
println("Testing model, Please wait. This may take a few minutes.")
count  = 0
for i=1:num_tests
    if round_output(model(X[i,:])) == Y[i]
        count += 1
    end
end

println("Test complete: $(count*100/num_tests) % accuracy achieved.")


text:
@testset "construction and promotion" begin
    @test basetype(ZZ{Int}) == Int
    @test depth(ZZ{Int8}) == 1
    @test iscoprime(ZZ(18), ZZ(35))
    @test isirreducible(ZZ(17))
    @test !isirreducible(ZZ(18))
    @test typeof(copy(ZZ(big"2"))) == ZZ{BigInt}
    a = ZZ(16)
    @test ZZ{Int}(a) === a
    @test ZZ{Int8}(100) + ZZ(100) == ZZ(200)
    @test ZZ(10) + 10 == 20
    @test ZZ(1) + QQ(1, 2) == QQ(3, 2)
    @test ZZ(Int8(1)) + QQ(1, 2) == QQ(3, 2)
    @test typeof(ZZ(1) + QQ(1, 2)) == QQ{Int}
    @test typeof(ZZ(1) + 1//2) == QQ{Int}
    @test hash(ZZ(big"123")) == hash(123)
end

@testset "ZZ{$T}" for T in (Int32, Int64, BigInt)

    Z = ZZ{T}
    z = zero(Z)
    @test z == zero(z)
    @test iszero(z)
    @test !isone(z)
    o = one(z)
    @test isunit(o)
    @test isone(o)
    @test !iszero(o)
    n1 = T(19)
    n2 = Base.hastypemax(T) ? typemax(T) : T(big"987654321987654321987654321")
    z1 = Z(n1)
    z2 = ZZ(n2)
    @test !isunit(z1)
    @test z1 + z1 == ZZ(T(2n1))
    T != BigInt && @test_throws OverflowError z2 + z2

    @test z1 - z1 == z
    @test z1 - z2 == ZZ(n1 - n2)
    @test -z2 == ZZ(-n2)
    @test isunit(-o)
    T != BigInt && @test_throws OverflowError -z1 - z2

    @test z1 * z1 == ZZ(T(n1 * n1))
    @test z1 * n1 == ZZ(T(n1 * n1))
    @test n1 * z1 == ZZ(T(n1 * n1))
    T != BigInt && @test_throws OverflowError z1 * z2

    @test z1^2 == ZZ(n1^2)
    T != BigInt && @test_throws OverflowError z1^1000

    z3 = 2z1
    z4 = z3 + o
    @test z3 ÷ z1 == ZZ(T(2))
    @test z1 ÷ z3 == z
    @test z3 / z1 == ZZ(T(2))
    @test z1 \ z3 == ZZ(T(2))
    @test_throws DomainError z4 / z2

    @test inv(-o) == -o
    @test_throws DomainError inv(z1)

    @test gcd(z3, z1) == z1
    @test gcd(z3, z1, z3) == z1
    @test gcdx(z3, z1) == (z1, zero(z3), one(z1))
    @test gcdx(ZZ(60), ZZ(28), ZZ(6)) == (2, [-1, 2, 1])
    @test lcm(z3, z1) == z3
    @test rem(ZZ(12), ZZ(5)) == ZZ(2)

    @test sprint(show, z3) == "$(Int(2n1))"

    @test value(z2) === n2

    @test factor(z3) == [Z(2) => 1, n1 => 1]
    @test eltype(factor(z3)) == Pair{Z,Int}
end


text:
Msun = 1.99e33
km = 1e5
kg = 1e3

rad = 10 * km
mass = 1.4 * Msun
vol = (4*pi/3)*rad^3

dens = mass/vol

max_stars_repo_name: natj/thesis

max_issues_repo_name: natj/thesis

max_forks_repo_name: natj/thesis


text:
module DefferedAcceptance

function converter(a::Vector{Vector{Int}}, b::Vector{Vector{Int}})
    m = length(a)
    n = length(b)

    a_2d = Matrix{Int64}(n+1, m)
    b_2d = Matrix{Int64}(m+1, n)

    for (t,i) in enumerate(a)
        if length(i) != n
            a_2d[1:length(i), t] = i
            a_2d[length(i)+1, t] = 0
            a_2d[(length(i)+2):end, t] = Array([j for j in 1:n if !(j in i)])
        else
            a_2d[1:length(i), t] = i
            a_2d[length(i)+1, t] = 0
        end
    end

    for (t,i) in enumerate(b)
        if length(i) != m
            b_2d[1:length(i), t] = i
            b_2d[length(i)+1, t] = 0
            b_2d[length(i)+2:end, t] = Array([j for j in 1:m if !(j in i)])
        else
            b_2d[1:length(i), t] = i
            b_2d[length(i)+1, t] = 0
        end
    end
    return a_2d, b_2d
end

# many-to-manyの関数を入れる
function ikegami(a::Vector{Vector{Int}}, b::Vector{Vector{Int}}, caps::Vector{Int})
    w = converter(a, b)
    a_2d = w[1]
    b_2d = w[2]
    return ikegamida_mm(a_2d, b_2d, caps)
end

function ikegami(a::Vector{Vector{Int}}, b::Vector{Vector{Int}})
    caps = ones(Int, length(b))
    w = converter(a, b)
    a_2d = w[1]
    b_2d = w[2]
    a, b, c =  ikegami(a_2d, b_2d, caps)
    return a, b
end

# one-to-one
function ikegamida(m_prefs, f_prefs)

    m_num, f_num = size(m_prefs)[2], size(f_prefs)[2]
    m_pool, f_pool = collect(1:m_num), collect(1:f_num)
    match_or_unmatch = Array(Bool, m_num)
    for i in 1:m_num
        match_or_unmatch[i] = true
    end
    m_matched, f_matched = ones(Int64, m_num) + f_num, ones(Int64, f_num) + m_num

    m_rank = Array(Int64, (f_num + 1, m_num))
    f_rank = Array(Int64, (m_num + 1, f_num))

    sorting = Array(Int64, (f_num + 1, 2))
    sorting[:, 1] = collect(1:(f_num + 1))

    for i in 1:m_num
        sorting[:, 2] = m_prefs[:, i]
        m_rank[:,i] = sortrows(sorting, by = x->(x[2]))[:,1]
        sorting = sortrows(sorting, by = x->(x[1]))
    end

    sorting2 = Array(Int64, (m_num + 1, 2))
    sorting2[:, 1] = collect(1:(m_num + 1))

    for i in 1:f_num
        sorting2[:, 2] = f_prefs[:, i]
        f_rank[:,i] = sortrows(sorting2, by = x->(x[2]))[:,1]
        sorting2 = sortrows(sorting2, by = x->(x[1]))
    end

    while sum(match_or_unmatch) > 0
        for i in 1:m_num
            if match_or_unmatch[i] == true

                for j in m_prefs[:, i]

                    if j == 0
                        m_matched[i] = 0
                        match_or_unmatch[i] = false
                        break

                    else
                        no_marriage = f_rank[1, j]

                        if f_rank[i+1, j] < no_marriage

                            if f_matched[j] == m_num + 1
                                f_matched[j] = i
                                m_matched[i] = j
                                match_or_unmatch[i] = false
                                break

                            else

                                if f_rank[f_matched[j]+1, j] > f_rank[i+1, j]
                                    match_or_unmatch[f_matched[j]] = true
                                    f_matched[j] = i
                                    m_matched[i] = j
                                    match_or_unmatch[i] = false
                                    break
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    for (i,t) in enumerate(f_matched)
        if t == m_num + 1
            f_matched[i] = 0
        end
    end

     for (i,t) in enumerate(m_matched)
        if t == f_num + 1
            m_matched[i] = 0
        end
    end

    return m_matched, f_matched
end

# many-to-many
# capsがrespの時とpropの時とを関数内で場合分けする。

function ikegamida_mm(prop_prefs::Matrix{Int}, resp_prefs::Matrix{Int}, caps::Vector{Int})

    # set up
    prop_num, resp_num = size(prop_prefs)[2], size(resp_prefs)[2]
    # capsがrespの持つcapsの場合
    if length(caps) == resp_num

        match_or_unmatch = Array(Bool, prop_num)
        for i in 1:prop_num
            match_or_unmatch[i] = true
        end
        prop_matched, resp_matched = ones(Int64, prop_num) + resp_num, ones(Int64, sum(caps)) + prop_num

        prop_rank = Array(Int64, (resp_num + 1, prop_num))
        resp_rank = Array(Int64, (prop_num + 1, resp_num))

        sorting = Array(Int64, (resp_num + 1, 2))
        sorting[:, 1] = collect(1:(resp_num + 1))
        for i in 1:prop_num
            sorting[:, 2] = prop_prefs[:, i]
            prop_rank[:,i] = sortrows(sorting, by = x->(x[2]))[:,1]
            sorting = sortrows(sorting, by = x->(x[1]))
        end

        sorting2 = Array(Int64, (prop_num + 1, 2))
        sorting2[:, 1] = collect(1:(prop_num + 1))
        for i in 1:resp_num
            sorting2[:, 2] = resp_prefs[:, i]
            resp_rank[:,i] = sortrows(sorting2, by = x->(x[2]))[:,1]
            sorting2 = sortrows(sorting2, by = x->(x[1]))
        end

        insertarray = ones(Int64, maximum(caps) + 1) * (prop_num+1)

        #making indptr
        indptr = Array(Int, resp_num+1)
        indptr[1] = 1
        for i in 1:resp_num
            indptr[i+1] = indptr[i] + caps[i]
        end

        #main loop
        while sum(match_or_unmatch) > 0
            for i in 1:prop_num

                if match_or_unmatch[i] == true

                    for j in prop_prefs[:, i]

                        if j == 0
                            prop_matched[i] = 0
                            match_or_unmatch[i] = false
                            break

                        else
                            cutoff = resp_rank[1, j]

                            if resp_rank[i+1, j] < cutoff

                                #capを超えてないなら無条件まっち
                                if prop_num + 1 in resp_matched[indptr[j]:indptr[j+1]-1]

                                    for (m,q) in enumerate(resp_matched[indptr[j]:indptr[j+1]-1])

                                        if q == prop_num + 1
                                            resp_matched[indptr[j]+m-1] = i
                                            prop_matched[i] = j
                                            match_or_unmatch[i] = false
                                            break

                                        else
                                            if resp_rank[q+1, j] > resp_rank[i+1, j]
                                                if m == 1
                                                    insertarray[m+1:length(resp_matched[indptr[j]:indptr[j+1]-1])] = resp_matched[indptr[j]:indptr[j+1]-1][m:end-1]
                                                    insertarray[m] = i
                                                else
                                                    insertarray[m+1:length(resp_matched[indptr[j]:indptr[j+1]-1])] = resp_matched[indptr[j]:indptr[j+1]-1][m:end-1]
                                                    insertarray[1:m-1] = resp_matched[indptr[j]:indptr[j+1]-1][1:m-1]
                                                    insertarray[m] = i
                                                end
                                                resp_matched[indptr[j]:indptr[j+1]-1] = insertarray[1:length(resp_matched[indptr[j]:indptr[j+1]-1])]
                                                prop_matched[i] = j
                                                match_or_unmatch[i] = false
                                                break
                                            end
                                        end
                                    end
                                    break

                                #capを超えてたら競争
                                else
                                    for (n,p) in enumerate(resp_matched[indptr[j]:indptr[j+1]-1])
                                        if resp_rank[p+1, j] > resp_rank[i+1, j]
                                            if n == 1
                                                insertarray[n+1:length(resp_matched[indptr[j]:indptr[j+1]-1])] = resp_matched[indptr[j]:indptr[j+1]-1][n:end-1]
                                                insertarray[n] = i
                                            else
                                                insertarray[n+1:length(resp_matched[indptr[j]:indptr[j+1]-1])] = resp_matched[indptr[j]:indptr[j+1]-1][n:end-1]
                                                insertarray[1:n-1] = resp_matched[indptr[j]:indptr[j+1]-1][1:n-1]
                                                insertarray[n] = i
                                            end
                                            match_or_unmatch[resp_matched[indptr[j]:indptr[j+1]-1][end]] = true
                                            prop_matched[resp_matched[indptr[j]:indptr[j+1]-1][end]] = resp_num + 1
                                            resp_matched[indptr[j]:indptr[j+1]-1] = insertarray[1:length(resp_matched[indptr[j]:indptr[j+1]-1])]
                                            prop_matched[i] = j
                                            match_or_unmatch[i]=false
                                            break
                                        end
                                    end
                                    if match_or_unmatch[i] == false
                                        break
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end

    # propがcapsを持つ時
    else
        match_or_unmatch = Array(Bool, prop_num)
        for i in 1:prop_num
            match_or_unmatch[i] = true
        end
        prop_matched, resp_matched = ones(Int64, sum(caps)) + resp_num, ones(Int64, resp_num) + prop_num

        prop_rank = Array(Int64, (resp_num + 1, prop_num))
        resp_rank = Array(Int64, (prop_num + 1, resp_num))

        sorting = Array(Int64, (resp_num + 1, 2))
        sorting[:, 1] = collect(1:(resp_num + 1))
        for i in 1:prop_num
            sorting[:, 2] = prop_prefs[:, i]
            prop_rank[:,i] = sortrows(sorting, by = x->(x[2]))[:,1]
            sorting = sortrows(sorting, by = x->(x[1]))
        end

        sorting2 = Array(Int64, (prop_num + 1, 2))
        sorting2[:, 1] = collect(1:(prop_num + 1))
        for i in 1:resp_num
            sorting2[:, 2] = resp_prefs[:, i]
            resp_rank[:,i] = sortrows(sorting2, by = x->(x[2]))[:,1]
            sorting2 = sortrows(sorting2, by = x->(x[1]))
        end

        insertarray = ones(Int64, maximum(caps) + 1) * (prop_num+1)

        #making indptr
        indptr = Array(Int, prop_num+1)
        indptr[1] = 1
        for i in 1:prop_num
            indptr[i+1] = indptr[i] + caps[i]
        end

        #main loop
        while sum(match_or_unmatch) > 0
            for i in 1:prop_num

                if match_or_unmatch[i] == false
                    break

                else
                    for j in prop_prefs[:, i]

                        # それ以降は0にして良い
                        if j == 0
                            for (n,t) in enumerate(prop_matched[indptr[i]:indptr[i+1]-1])
                                if t == resp_num + 1
                                    prop_matched[indptr[i]:indptr[i+1]-1][n] = 0
                                end
                            end
                            match_or_unmatch[i] = false
                            break

                        else
                            cutoff = resp_rank[1, j] # jが今proposeしている相手の持つ選考における0のランク

                            # proposeを受け入れてもいいなら…
                            if resp_rank[i+1, j] < cutoff

                                #マッチしてなければ無条件でマッチ
                                if prop_num + 1 == resp_matched[j]
                                    resp_matched[j] = i
                                    # ここでprop_matchedの場所を指定する必要
                                    m = sum(prop_matched[indptr[i]:indptr[i+1]-1] .!= resp_num+1)
                                    prop_matched[indptr[i]+m] = j

                                    # capが埋まったらfalseにしてbreak
                                    if m+1 == caps[i]
                                        match_or_unmatch[i] = false
                                        break
                                    end

                                #マッチしていたら競争
                                else
                                    rival = resp_matched[j]

                                    # rivalがかつ場合はするー
                                    # rivalが負ける場合
                                    if resp_rank[rival+1, j] > resp_rank[i+1, j]
                                        # rivalとか指定しなくても、prop_matchedにはもともとjは一人しかいないから、そこをreplaceすればいい？
                                        for (k,l) in enumerate(prop_matched)
                                            if l == j
                                                prop_matched[k] = resp_num + 1
                                            end
                                        end
                                        prop_matched[indptr[rival]:indptr[rival+1]-1] = sort(prop_matched[indptr[rival]:indptr[rival+1]-1])
                                        match_or_unmatch[rival] = true

                                        m = sum(prop_matched[indptr[i]:indptr[i+1]-1] .!= resp_num+1)
                                        prop_matched[indptr[i]+m] = j
                                        resp_matched[j] = i
                                        if m + 1 == caps[i]
                                            match_or_unmatch[i] = false
                                            break
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    for (i,t) in enumerate(prop_matched)
        if t == resp_num + 1
            prop_matched[i] = 0
        end
    end

    for (i,t) in enumerate(resp_matched)
        if t == prop_num + 1
            resp_matched[i] = 0
        end
    end

    return prop_matched, resp_matched, indptr
end

function ikegamida_mm(prop_prefs::Matrix{Int}, resp_prefs::Matrix{Int})
    caps = ones(Int, size(resp_prefs, 2))
    prop_matches, resp_matches, indptr =
    ikegamida_mm(prop_prefs, resp_prefs, caps)
    return prop_matches, resp_matches
end

end


text:
# decompose formula into normal vs iv part
function decompose_iv(f::FormulaTerm)
	formula = f
	formula_endo = nothing
	formula_iv = nothing
	for term in eachterm(f.rhs)
		if isa(term, FormulaTerm)
			formula_endo = FormulaTerm(ConstantTerm(0), term.lhs)
			formula_iv = FormulaTerm(ConstantTerm(0), term.rhs)
		end
		formula = FormulaTerm(f.lhs, tuple((t for t in eachterm(f.rhs) if !isa(t, FormulaTerm))...))
	end
	return formula, formula_endo, formula_iv
end

##############################################################################
##
## build model
##
##############################################################################
eachterm(x::AbstractTerm) = (x,)
eachterm(x::NTuple{N, AbstractTerm}) where {N} = x

allvars(::Nothing) = Array{Symbol}(undef, 0)
allvars(term::Union{AbstractTerm, NTuple{N, AbstractTerm}}) where {N} = termvars(term)
allvars(f::Union{Expr, Symbol}) = termvars(@eval(@formula(nothing ~$(f))).rhs)


text:
"Autumn Language"
module LambdaC
using Reexport

include("aexpr.jl")
@reexport using .AExpressions

include("util.jl")
@reexport using .Util

include("sexpr.jl")
@reexport using .SExpr

include("subexpr.jl")
@reexport using .SubExpressions

# include("passes.jl")

# include("compileutils.jl")
# @reexport using .CompileUtils

# include("compile.jl")
# @reexport using .Compile

# include("abstractinterpretation.jl")
# @reexport using .AbstractInterpretation

include("interpret/interpret.jl")
@reexport using .InterpretBig

# include("interpret/interpret2.jl")
# @reexport using .Interpret2

# include("scope.jl")
# @reexport using .Scope

# include("transform.jl")
# @reexport using .Transform

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# TODO: fix and test smoothing

module Simulations1D

using ..Curlcurls
using ..Dispersions
using ..Domains
using ..Laplacians
using ..Lattices
using ..Points
using ..SelfEnergies
using ..Shapes
using ..VectorFields
# using LinearAlgebra
using SparseArrays
using Statistics

import ..AbstractDomain
import ..Symmetric, ..Unsymmetric
import ..Simulation
import ..smooth_dielectric!
import ..smooth_pump!
import ..which_domains
import ..simulation_dielectric!
import ..simulation_pump!

function Simulation(
			ω₀::Real,
			lattice_domains::NTuple{K,LatticeDomain{1}},
			nondispersive_domains::NTuple{L,NondispersiveDomain{1}},
			dispersive_domains::NTuple{M,DispersiveDomain{1}};
			k₂₀::Real=0,
			k₃₀::Real=0,
			k₁₀::Real=sqrt(ω₀^2-k₂₀^2-k₃₀^2)
			) where {K,L,M}

	# simple checking for 1D symmetric
	K==1 || throw("in 1D, can only provide one LatticeDomain, but $K provided instead")
	lattice_domain = lattice_domains[1]

	smoothed = Ref(false)
	x = lattice_domain.x

	lattice_domain_indices = fill(1,length(x))
	nondispersive_domain_indices = which_domains(nondispersive_domains,x)
	dispersive_domain_indices = which_domains(dispersive_domains,x)

	ε = Vector{ComplexF64}(undef,length(x))
	if isempty(nondispersive_domains)
		for i ∈ eachindex(ε) ε[i] = lattice_domain.ε end
	else
		dielectrics = map(n -> getfield(n,:dielectric), nondispersive_domains)
		for i ∈ eachindex(ε)
			d = nondispersive_domain_indices[i]
			if iszero(d)
				ε[i] = lattice_domain.ε
			else
				ε[i] = dielectrics[d](x[i])
			end
		end
	end

	F = Vector{Float64}(undef,length(x))
	χ = Vector{AbstractDispersion}(undef,length(x))
	Fs = Vector{Vector{Float64}}(undef,length(dispersive_domains))
	for i ∈ eachindex(Fs) Fs[i] = zeros(Float64,length(x)) end
	if isempty(dispersive_domains)
		for i ∈ eachindex(F) F[i] = 0 end
		for i ∈ eachindex(F) χ[i] = NoDispersion() end
	else
		pumps = map(n -> n.pump, dispersive_domains)
		χs = map(n -> n.χ, dispersive_domains)
		for i ∈ eachindex(F)
			d = dispersive_domain_indices[i]
			if iszero(d)
				F[i] = 0
				χ[i] = NoDispersion()
			else
				F[i] = pumps[d](x[i])
				χ[i] = χs[d]
				Fs[d][i] = pumps[d](x[i])
			end
		end
	end

	indices = lattice_domain.indices
	nnm = lattice_domain.nnm
	nnp = lattice_domain.nnp
	interior = lattice_domain.interior
	surface = lattice_domain.surface
	bulk = lattice_domain.bulk
	corner = lattice_domain.corner

	# checking
	surface[1] || throw("first site should be an endpoint, something is amiss")
	surface[end] || throw("last site should be an endpoint, something is amiss")
	sum(surface)>2 && throw("no point except first and last should be an endpoint, something is amiss")

	x_half = _generate_half_x(lattice_domain.lattice, x, indices)

	# σ = (Vector{ComplexF64}(undef,length(x)),)
	# σ_half = (Vector{ComplexF64}(undef,length(x)+1),)
	# boundary_layer!(σ, lattice_domains, x, lattice_domain_indices)
	# boundary_layer!(σ_half, lattice_domains, x_half[1], vcat(lattice_domain_indices, lattice_domain_indices[end]))
	σ = boundary_layer(lattice_domain,x)
	σ_half = boundary_layer(lattice_domain,x_half[1])
	α = (1 .+ 1im*σ[1]/k₁₀,)
	α_half = (1 .+ 1im*σ_half[1]/k₁₀,)

	laplacian = Laplacian{Symmetric}(lattice_domain.lattice,α[1],α_half[1])
	# curlcurl = Curlcurl{Symmetric}(lattice_domain.lattice,α[1],α_half[1],nnm,nnp,indices,interior,surface)

	Σ = SelfEnergy{Symmetric}(lattice_domain,α_half[1])

	return Simulation{1,Symmetric,ComplexF64,typeof(lattice_domains),typeof(nondispersive_domains),typeof(dispersive_domains),typeof(Σ)}(
		lattice_domains,
		nondispersive_domains,
		dispersive_domains,
		x,
		lattice_domain_indices,
		nondispersive_domain_indices,
		dispersive_domain_indices,
		ε,
		F,
		χ,
		Fs,
		laplacian,
		# curlcurl,
		Σ,
		α,
		σ,
		x_half,
		σ_half,
		ω₀,
		k₁₀,
		k₂₀,
		k₃₀,
		smoothed)
end

################################################################################
# SIMULATION{1} building utilities

# used in building Simulation{1}
function _generate_half_x(lattice::Lattice{1}, x::Vector{Point{1,TP}}, indices) where TP
	half_x = (Vector{Point{1,Cartesian}}(undef,length(x)+1),)
	for i ∈ eachindex(x) half_x[1][i] = lattice[indices[i][1]-.5] end
	half_x[1][end] = lattice[indices[length(x)][1]+.5]
	return half_x
end

# 1-D hook for boundary_layer found in Simulations.jl
@inline function boundary_layer(domain::LatticeDomain{1}, x::Vector{Point{1,C}}) where C
	typeof(domain.boundary.shape)<:Interval || throw("1D only defined for shape=Interval, but shape=", domain.boundary.shape)
	return (domain.boundary.bls[1].(x)+domain.boundary.bls[2].(x),)
end

################################################################################
# 1-d smoothing functions

function smooth_dielectric!(sim::Simulation{1}, num_sub_pixel::Integer = NUM_SUBPIXELS)
    indices = sim.nondispersive_domain_indices
    lattice = sim.lattice
    X = Vector{Point{1}}(undef,num_sub_pixel)
    E = Vector{ComplexF64}(undef,num_sub_pixel)
	r = LinRange(-.5,.5,num_sub_pixel)
    for i ∈ eachindex(indices)
		if 1<i<length(indices)
			idx, = Tuple(sim.lattice_domain.indices[i])
			if !(indices[i] == indices[i-1] == indices[i+1])
				for j ∈ eachindex(r) X[j] = lattice[idx+r[j]] end
		        simulation_dielectric!(E,sim,X)
				sim.ε[i] = mean(E)
			end
		else
			# placeholder for periodic lattice domains
		end
    end
    return nothing
end

function smooth_pump!(sim::Simulation{1}, num_sub_pixel::Integer = NUM_SUBPIXELS)
	indices = sim.dispersive_domain_indices
    lattice = sim.lattice
    X = Vector{Point{1}}(undef,num_sub_pixel)
    F = Vector{ComplexF64}(undef,num_sub_pixel)
	r = LinRange(-.5,.5,num_sub_pixel)
    for i ∈ eachindex(indices)
		if 1<i<length(indices)
			idx, = Tuple(sim.lattice_domain.indices[i])
			if !(indices[i] == indices[i-1] == indices[i+1])
				for j ∈ eachindex(r) X[j] = lattice[idx+r[j]] end
		        simulation_pump!(F,sim,X)
				sim.F[i] = mean(F)
			end
		end
    end
    return nothing
end

################################################################################
# extras

function Base.getproperty(sim::Simulation{1}, sym::Symbol)
	if sym==:lattice_domain
		return getfield(sim,:lattice_domains)[1]
	elseif Base.sym_in(sym,(:k10,:k1₀))
		return getfield(sim,:k₁₀)
	elseif Base.sym_in(sym,(:k20,:k2₀))
		return getfield(sim,:k₂₀)
	elseif Base.sym_in(sym,(:k30,:k3₀))
		return getfield(sim,:k₃₀)
	elseif Base.sym_in(sym,(:shape,:boundary,:lattice))
		return getproperty(getfield(sim,:lattice_domains)[1],sym)
	elseif Base.sym_in(sym,propertynames(getproperty(getfield(sim,:lattice_domains)[1],:lattice)))
		return getproperty(getfield(sim,:lattice_domains)[1],sym)
	elseif Base.sym_in(sym,(:lat,:lattice,:Lat,:Lattice))
		return getfield(sim,:lattice_domains)[1].lattice
	else
		return getfield(sim,sym)
	end
end

function Base.propertynames(sim::Simulation{1}, private=false)
	if private
		return fieldnames(Simulation)
	else
		return (:shape, :boundary, :lattice, :lattice_domain, :lattice_domains, :nondispersive_domains, :dispersive_domains, :ε, :F, :x, :χ, :ω₀, :k1₀, :k2₀, :k3₀, :lattice, propertynames(sim.lattice)...)
	end
end

################################################################################
# Pretty Printing

import ...PRINTED_COLOR_GOOD
import ...PRINTED_COLOR_WARN
import ...PRINTED_COLOR_DARK
import ...PRINTED_COLOR_NUMBER

function Base.show(io::IO,sim::Simulation{1,CLASS}) where CLASS
	if sim.smoothed[]
		printstyled(io,"Smoothed ",color=PRINTED_COLOR_GOOD)
	else
		printstyled(io,"Unsmoothed ",color=PRINTED_COLOR_WARN)
	end
	CLASS<:Symmetric ? print(io,"Symmetric") : nothing
	CLASS<:Unsymmetric ? print(io,"Unsymmetric") : nothing
	printstyled(io," Simulation\n",color=PRINTED_COLOR_DARK)
    print(io,"\tn sites: ")
	printstyled(io,length(sim),"\n",color=PRINTED_COLOR_NUMBER)
	println(io,"\t====================")
	println(IOContext(io,:tabbed2=>true),sim.boundary)
	println(io)
	println(IOContext(io,:tabbed2=>true),sim.lattice)
	println(io)
	domains = (sim.nondispersive_domains...,sim.dispersive_domains...)
	for d ∈ eachindex(domains)
		print(io,"\t\tDomain ", d, " ")
		print(IOContext(io,:tabbed2=>true),domains[d])
		d < length(domains) ? println(io) : nothing
	end
end

end

using .Simulations1D


text:
"""
    Petri

Provides a modeling framework for representing and solving stochastic petri nets

"""
module Petri

export Model, NullPetri, EmptyPetri, vectorfield, Graph

include("types.jl")
include("solvers.jl")
include("visualization.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule mlpack_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("mlpack")
JLLWrappers.@generate_main_file("mlpack", UUID("156a7db3-27ba-586e-a86b-f061da4f95ea"))
end  # module mlpack_jll


text:
using Documenter, OptiMimi

makedocs(
         modules=[OptiMimi],
         pages = Any[
                     "Introduction" => "index.md",
                     "Linear Programming" => "linproghouse.md",
                     ]
         )


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#     Functions to manage the SPI.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export init_spi, spi_transfer, spi_transfer!

################################################################################
#                                Initialization
################################################################################

"""
    init_spi(devices; mode = 0, max_speed_hz = 4_000_000, bits_per_word = 8)

Initialize the SPI devices. `devices` can be a string with the path to `spidev`
or a vector of strings with a set of SPI devices that will be initialized.

# Keywords

* `mode`: Set the mode of the SPI. (**Default** = 0)
* `max_speed_hz`: Maximum allowed speed in SPI communication [Hz].
                  (**Default** = 4_000_000)
* `bits_per_word`: Number of bits per word in SPI communication.
                   (**Default** = 8)

Notice that all keywords can be a `Integer`, when the configuration will be
applied to all SPI devices, or a `Vector` of `Integers`, when different
configurations can be applied to the initialized devices.

"""
@inline init_spi(devices::String; kwargs...) = init_spi([devices]; kwargs...)

function init_spi(devices::AbstractVector{String};
                  mode::Union{Integer,Vector{Integer}} = 0,
                  max_speed_hz::Union{Integer,Vector{Integer}} = 4_000_000,
                  bits_per_word::Union{Integer,Vector{Integer}} = 8)

    # Number of devices that the user wants to initialize.
    num_init_devices = length(devices)

    # If the arguments are not vectors, then transform them to vectors
    # considering the size of `devices`.
    if !(typeof(mode) <: AbstractVector)
        mode = [mode for _ = 1:num_init_devices]
    end

    if !(typeof(max_speed_hz) <: AbstractVector)
        max_speed_hz = [max_speed_hz for _ = 1:num_init_devices]
    end

    if !(typeof(bits_per_word) <: AbstractVector)
        bits_per_word = [bits_per_word for _ = 1:num_init_devices]
    end

    # If we have already initialized, then close all the devices first.
    objects.spi_init && close_spi()

    # Initialize each device.
    spidev = Vector{SPIDEV}(undef, 0)

    @inbounds for i = 1:num_init_devices
        @assert (mode[i] ≥ 0) error("Invalid SPI mode.")
        @assert (max_speed_hz[i] > 0) error("Invalid maximum speed [Hz].")
        @assert (bits_per_word[i] > 0) error("Invalid number of bits per word.")

        try
            # Open the device
            # ==================================================================

            io = open(devices[i], "w+")

            # Configure SPI interface
            # ==================================================================

            _ioctl(io, SPI_IOC_WR_MODE, Ref(mode[i]))
            _ioctl(io, SPI_IOC_RD_MODE, Ref(mode[i]))

            _ioctl(io, SPI_IOC_WR_BITS_PER_WORD, Ref(bits_per_word[i]))
            _ioctl(io, SPI_IOC_RD_BITS_PER_WORD, Ref(bits_per_word[i]))

            _ioctl(io, SPI_IOC_WR_MAX_SPEED_HZ, Ref(max_speed_hz[i]))
            _ioctl(io, SPI_IOC_RD_MAX_SPEED_HZ, Ref(max_speed_hz[i]))

            # Pack values
            # ==================================================================
            spidev_i = SPIDEV(io, max_speed_hz[i], bits_per_word[i])
            push!(spidev, spidev_i)
        catch
            error("Could not open `$(devices[i])`. Make sure you have the required permissions.")
        end
    end

    # Get the size of SPI buffer.
    buffer_size = parse(Int, read("/sys/module/spidev/parameters/bufsiz", String))

    objects.spidev          = spidev
    objects.spi_buffer_size = buffer_size
    objects.spi_init        = true

    return nothing
end

"""
    close_spi()

Close all SPI connections.

"""
function close_spi()
    if objects.spi_init
        @inbounds for d in objects.spidev
            close(d.io)
        end
        objects.spidev = Vector{SPIDEV}(undef, 0)
        objects.spi_init = false
    end

    return nothing
end

################################################################################
#                                   Transfer
################################################################################

"""
    spi_transfer(devid, tx_buf; kwargs...)

Execute a full duplex transfer to SPI device `devid`. `devid` is the ID of the
SPI device considering the initialization order when the function `init_spi` was
called.

`tx_buf` can be a vector of `Integer`, in which only one message will be sent,
or a vector of vectors of `Integer`, in which multiple messages will be sent.

The result is returned in an object with the same type of `tx_buf` together with
the number of words received.

!!! note "Allocations"

    This function will perform allocations because they create the vector that
    will be returned. If this is not desired, then use the in-place version
    `spi_transfer!`.

# Keywords

The same keywords of `spi_transfer!` can be used.

"""
function spi_transfer(devid::Integer, tx_buf::AbstractVector{T}; kwargs...) where
    T<:Integer

    # Allocate the vector that will hold the returned words.
    rx_buf = zeros(T, length(tx_buf))

    # Perform the transfer.
    ret = spi_transfer!(devid, tx_buf, rx_buf; kwargs...)

    return rx_buf, ret
end

function spi_transfer(devid::Integer, tx_buf::AbstractVector{T}; kwargs...) where
    T<:Vector{U} where U<:Integer

    # Allocate the vectors that will hold the returned words.
    rx_buf = zeros.(U, length.(tx_buf))

    # Perform the transfer.
    ret = spi_transfer!(devid, tx_buf, rx_buf; kwargs...)

    return rx_buf, ret
end

"""
    spi_transfer!(devid, tx_buf, rx_buf; kwargs...)

Execute a full duplex transfer to SPI device `devid`. `devid` is the ID of the
SPI device considering the initialization order when the function `init_spi` was
called.

`tx_buf` can be a vector of `Integer`, in which only one message will be sent,
or a vector of vectors of `Integer`, in which multiple messages will be sent.

The received data will be stored in `rx_buf` that must have the same type of
`tx_buf` and enough size.

This function returns the number of bytes received.

!!! note "Allocations"

    This function will not allocate only if the number of messages sent is lower
    than the constant `BaremetalPi._SPI_BUFFER_SIZE`. Otherwise, it will perform
    an allocation because it must allocate a vector of
    `struct_spi_ioc_transfer`.

# Keywords

* `max_speed_hz`: If > 0, then override the default maximum transfer speed with
                  this value [Hz]. (**Default** = 0)
* `delay_usecs`: If ≥ 0, then override the default delay with this value.
                 (**Default** = -1)
* `bits_per_word`: If > 0, then override the number of bits per word with this
                   value. (**Default** = 0)
* `cs_change`: If `false`, the deselect the device at the end of the transfer.

"""
function spi_transfer!(devid::Integer,
                       tx_buf::AbstractVector{T},
                       rx_buf::AbstractVector{T};
                       max_speed_hz::Integer = 0,
                       delay_usecs::Integer = -1,
                       bits_per_word::Integer = 8,
                       cs_change::Bool = false) where T<:Vector{U} where U <:Integer

    @assert objects.spi_init "SPI not initialized. Run init_spi()."
    @assert (0 < devid ≤ length(objects.spidev)) "SPI device ID is out of bounds."

    spidev = objects.spidev[devid]

    # Number of messages to be transmitted.
    num_msgs = length(tx_buf)

    @assert (length(rx_buf) ≥ num_msgs) "The number of buffers in `rx_buf` must be equal or bigger than the number of buffers in `tx_buf`."

    # Check default parameters.
    max_speed_hz  ≤ 0 && (max_speed_hz  = spidev.max_speed_hz)
    delay_usecs   < 0 && (delay_usecs   = 0)
    bits_per_word ≤ 0 && (bits_per_word = spidev.bits_per_word)

    # Check if we can use the allocated buffer to transfer the message.
    if num_msgs > _SPI_BUFFER_SIZE
        descs = Vector{struct_spi_ioc_transfer}(undef, num_msgs)
    else
        descs = spidev.bdescs
    end

    @inbounds for i = 1:num_msgs
        msg_size = length(tx_buf[i])

        @assert (msg_size*bits_per_word ≤ objects.spi_buffer_size*8) "The message to be transmitted is larger than the SPI buffer."
        @assert (length(rx_buf[i]) ≥ msg_size) "The length of `rx_buf[i]` must be equal or bigger than that of `tx_buf[i]`."

        # Create the structure that contains the information of the SPI transfer.
        descs[i] = struct_spi_ioc_transfer(pointer(tx_buf[i]),
                                           pointer(rx_buf[i]),
                                           # In SPI, the number of transmitted
                                           # and received words are always the
                                           # same.
                                           msg_size*sizeof(U),
                                           max_speed_hz,
                                           delay_usecs,
                                           bits_per_word,
                                           cs_change)
    end

    # Execute the transfer.
    return _ioctl(fd(spidev.io), SPI_IOC_MESSAGE(num_msgs), descs)
end

function spi_transfer!(devid::Integer,
                       tx_buf::AbstractVector{T},
                       rx_buf::AbstractVector{T};
                       max_speed_hz::Integer = 0,
                       delay_usecs::Integer = -1,
                       bits_per_word::Integer = 8,
                       cs_change::Bool = false) where T<:Integer

    @assert objects.spi_init "SPI not initialized. Run init_spi()."
    @assert (0 < devid ≤ length(objects.spidev)) "SPI device ID is out of bounds."

    spidev = objects.spidev[devid]

    # Check default parameters.
    max_speed_hz  ≤ 0 && (max_speed_hz  = spidev.max_speed_hz)
    delay_usecs   < 0 && (delay_usecs   = 0)
    bits_per_word ≤ 0 && (bits_per_word = spidev.bits_per_word)

    msg_size = length(tx_buf)

    @assert (msg_size*bits_per_word ≤ objects.spi_buffer_size*8) "The message to be transmitted is larger than the SPI buffer."
    @assert (length(rx_buf) ≥ msg_size) "The length of `rx_buf` must be equal or bigger than that of `tx_buf`."

    # Create the structure that contains the information of the SPI transfer.
    desc = struct_spi_ioc_transfer(pointer(tx_buf),
                                   pointer(rx_buf),
                                   # In SPI, the number of transmitted and
                                   # received words are always the same.
                                   msg_size*sizeof(T),
                                   max_speed_hz,
                                   delay_usecs,
                                   bits_per_word,
                                   cs_change)

    # Execute the transfer.
    return _ioctl(fd(spidev.io), SPI_IOC_MESSAGE(1), Ref(desc))
end