text:
#####################################################################################
# Model
#####################################################################################

const ħoec = ustrip(u"T*nm^2",ħ/e)
const μB = ustrip(u"T^-1*meV",μ_B)
const σ0´ = @SMatrix[1 0; 0 1]
const σx´ = @SMatrix[0 1; 1 0]
const σy´ = @SMatrix[0 -im; im 0]
const σz´ = @SMatrix[1 0; 0 -1]
const σ0τ0´ = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
const σ0τz´ = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
const σxτ0´ = @SMatrix[0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]
const σyτ0´ = @SMatrix[0 -im 0 0; im 0 0 0; 0 0 0 -im; 0 0 im 0]
const σzτ0´ = @SMatrix[1 0 0 0; 0 -1 0 0; 0 0 1 0; 0 0 0 -1]
const σxτz´ = @SMatrix[0 1 0 0; 1 0 0 0; 0 0 0 -1; 0 0 -1 0]
const σyτz´ = @SMatrix[0 -im 0 0; im 0 0 0; 0 0 0 im; 0 0 -im 0]
const σzτz´ = @SMatrix[1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]
const σyτy´ = @SMatrix[0 0 0 -1; 0 0 1 0; 0 1 0 0; -1 0 0 0]

"""
Note if localselfenergy == True || analyticalselfenergy == True, LS = 0 as a preset
Note as it is written you should pass missing args corresponding to the fields: a0, dinter, dintra t0

To do:
  1) fix kfs
"""

# Scaling model: t1 belongs to the low energy sector (high energy sector commented _model1) info:
# https://iopscience.iop.org/article/10.1088/0034-4885/76/5/056503/pdf

@with_kw struct Params @deftype Float64
                                        # Params AB bernal bilayer obtained from DFT calculations:
                                        # params ref -> https://arxiv.org/pdf/1511.06706.pdf (flag <-:dft)
    normal::Bool = false                # determines whether the system is normal or Nambu
    smooth::Bool = false                # U smooth profile if true otherwise constant
    smoothμ::Bool = false               # Chemical potential profile if true otherwise steplike
    mask::Bool = false                  # Determines whether a rectangular hole is made on the N region
    dir::Symbol = :y                    # Determines the direction in which U is smooth. Possible values: :x, :y or :xy
    edges::Symbol = :armchair_zigzag    # edges kwarg sets the nature of the SN and Nvacuum interfaces:
                                        # edges =:typeofxinterface_typeofyinterfaces
                                        # Default :armchair_zigzag.
                                        # Other options :zigzag_armchair, :zigzag_zigzag
    car::Bool = false                   # Interlayer hoppings in the SC leads are by default suppressed
    localselfenergy::Bool = false       # Determines whether we implement a localselfenergy ∝Δσyτy or not
    analyticalselfenergy::Bool = false  # Determines whether the self-energy due to a 3d unbounded SC
                                        # (local + non local terms) is implented or not. See notes
    scale = 1.                          # scaling parameter
    a0 = 0.246
    dinter = 1.36 * a0                  # Note that we leave dinter unaltered by scale
    dintra = a0 /sqrt(3)
    Ls = 3.                             # Length of the SC leads
    Ln = 3.                             # Length of the N region
    W = 3.                              # Width of the JJ junction
    t0 =  2.6                           # <-:dft
    t1 = .36                            # <-:dft
    τL = 1.                             # Transparency left
    τR = 1.                             # Transparency right
    t3 = .28                            # <-:dft #WARN number 1) I NEED TO RETHINK THE SCALING
    μN = 0. * t0
    μSL = t0 / 3.25                     # t0/3 > μS > 2 * t1
    μSR = t0 / 3.25                     # t0/3 > μS > 2 * t1
    U_dimer = 0.0* 0.15                 # Energy difference between dimer and non dimer sites: https://arxiv.org/pdf/1901.01332.pdf
    g1_on = 1.3e-3                      # Comparable to the gap
    g1_hop = 1.3e-3                     # Comparable to the gap
    kfs = 1.                            # Fermi wavelength of the superconductor
    U = 0.0*1e-3                        # Interlayer bias. If smooth = true is the maximal value that is reached at the edges,
                                        # else is the homogeneous value
    Δ = 1.3e-3                          # That of SC MoRe: ref ->  https://arxiv.org/abs/2006.05522
    B_par = 0*0.5                       # Applied on the outof plane
    B_perp = 0.                         # In-plane field and perpendicular to the system.
    λ = 0*0.5e-3                        # Between [-2.5,2.5] meV: ref -> https://arxiv.org/abs/1905.08688
    α = 0.0                             # Rashba coef
    g = 2.                              # g-factor
    V_par = 0.0                         # Zeeman field in y direction
    V_perp = 0.0                        # Zeeman fiel in z direction. You might want to uncomment them and compute it
                                        # with B playing attention to the scaling.
    Vu = 0.                             # modulus of the effective onsite energy Vuvector due to the e-e interactions (spin canting)
    θ = 0.                              # AF canting with respect to the y direction.
    Ex = 0                              # Electric field
    Ey = 0
    Ez = 0
    Δxμ = 10                            # Center of the sigmoidal μSL/R-μN
    λμ = 5                              # Smooth parameter
    λμbarrier = 1
    Δμbarrier = -10                     # Height of the leads barrier (tunnel regime) -> increasing normal reflections
    Δx_mask = 2                         # width of the mask on x
    Δy_mask = 2                         # width of the mask on y
end                                     # Units eV,  T, nm

"hamiltonian models given `p::Params`"
function models_device(p = Params())
    @unpack a0, dinter, dintra, Ls, Ln, W, μN, μSL, μSR, t0, t1, τL, τR, normal,
            t3, Δ, B_par, B_perp, λ, g, g1_on, g1_hop, kfs, U, U_dimer, α, V_par, V_perp,
            Vu, θ, smooth, smoothμ, scale, edges, dir, car, localselfenergy, analyticalselfenergy, Ex, Ey, Ez,
            Δxμ, λμ, λμbarrier, Δμbarrier = p

    (σ0τ0, σ0τz, σxτz, σyτz, σzτz, σyτy), σvect = sigmas(Val(normal))#normalorNambu(p)

    a0, Ln, Ls, W, regs = modelregs(p)
    edgereg, Redgereg, normalreg, scRreg, scLreg, screg, screg2 = [regs[i] for i in 1:length(regs)]
    t0 = t0 / scale
    t3 = t3 / scale
    dintra = dintra * scale
    dinter_warp = sqrt(dintra^2 + dinter^2)

    # t0, t1, and t3 hoppings
    modelintra = hopping(-t0 * σ0τz, range = dintra,
            sublats = (:At=>:Bt,:Ab=>:Bb,:Bb=>:Ab,:Bt=>:At))
    modelinter = hopping(0*t1 * σ0τz, range = dinter, sublats = (:At=>:Bb,:Bb=>:At)) +
                 hopping(0*t3 * σ0τz, range = dinter_warp, sublats = (:Bt=>:Ab, :Ab=>:Bt))

    ## SOC
    # ISING SOC (different sign of λ in different layers) - "Valley Zeeman"
    # modelKM = λ * (hopKM(:At)-hopKM(:Bt)-hopKM(:Ab)+hopKM(:Bb))
    # we don't have a KM but an Ising with different lambda on each layer
    ϕ0 = π/5
    hopKM(s) = hopping((r,dr) -> im *
        sign(cos(3*atan(dr[2],dr[1]) - ϕ0)) * ifelse(normal === true, σzτz, σzτ0´), sublats = (s => s),
         range = a0)

    # alternatively cos(3*atan(dr[2],dr[1])). Note that in order this to be independent on the
    # lattice orientation phi0 should not coincide with the direction of any NNN in the rotated lattice.
    # The following renormalization of λ comes from the relation between the SOC
    # gap in the k-Hamiltonian and the TB hopping parameter, t2: λ_KM = 3 √3 t2.
    # See: http://www.phys.ufl.edu/~pjh/teaching/phz7427/7427notes/ch6.pdf
    # see also https://journals.aps.org/prb/pdf/10.1103/PhysRevB.82.245412
    # https://journals.aps.org/prb/pdf/10.1103/PhysRevB.100.085412 This is Andor's et al paper

    λ /= (2*3*sqrt(3))
    modelIsing =  hopping((r, dr) ->λ * im * sign(r[3]) *
        sign(cos(3*atan(dr[2],dr[1]) - ϕ0)) * ifelse(normal === true, σzτz, σzτ0´),
        sublats = (:At => :At, :Ab => :Ab, :Bt => :Bt, :Bb => :Bb), range = a0)

    ising! = @hopping!((t, r, dr; λ = λ) -> t + λ * im * sign(r[3]) *
    sign(cos(3*atan(dr[2],dr[1]) - ϕ0)) * ifelse(normal === true, σzτz, σzτ0´);
    sublats = (:At => :At, :Ab => :Ab, :Bt => :Bt, :Bb => :Bb), range = a0)

    # RASHBA SOC (differnt sign of α in different layers)
    hopRashba() = hopping((r,dr) -> sign(r[3])*(dr[2]/dintra * ifelse(normal === true, σxτz, σxτ0´) - dr[1]/dintra * σyτz),
         range = dintra, sublats = (:At=>:Bt, :Ab=>:Bb, :Bt=>:At, :Bb=>:Ab))
    modelRashba = 2*im/3 * α * hopRashba()

    model00 = onsite(0I) + modelintra + modelinter + modelIsing
    model0 = onsite(0I) + modelintra + modelinter + modelIsing + modelRashba

    ## ZEEMAN and canting due to e-e interactions
        # Vu is spatially local (Hartree level, no Fock) and enters as an onsite
        # then it does not contribute to the flux
        # V_par = -1/2*g*μB*B_par
        # V_perp = -1/2*g*μB*B_perp # Uncomment this in params and think about the correct scaling

    Vzeem = @SVector[0, V_par, V_perp]
    Vuvect = Rz([0,Vu,0], θ)
    Vvect = Vzeem + Vuvect

    ##  UNBOUNDED TERMS

    modelNtop = onsite((- μN + U/2 ) * σ0τz + sum(σvect .* Vzeem), sublats= (:At, :Bt)) +
                onsite(sum(σvect .* Vuvect), sublats= (:At)) -
                onsite(sum(σvect .* Vuvect), sublats= (:Bt))

    modelNbot = onsite((- μN - U/2) * σ0τz + sum(σvect .* Vzeem), sublats= (:Ab, :Bb)) +
                onsite(sum(σvect .* Vuvect), sublats= (:Bb)) -
                onsite(sum(σvect .* Vuvect), sublats= (:At))

    modelS = onsite(- μSL * σ0τz - Δ * σyτy)

    ## Peierls phases

    piecewise(x) = ifelse(edges == :armchair_zigzag, clamp(x, 0*-Ln/2, Ln),
        clamp(x, -Ln/2 - a0/(2*√3) , Ln/2 + a0/(2*√3)))

    A(r, B_par, B_perp) = SA[0, B_perp, B_par]/ħoec * piecewise(r[1] + Ln/2)
    eφ(r, dr, B_par, B_perp) = diagphi(dot(A(r, B_par, B_perp), (dr + Δz(r))))
    diagphi(φ) = normal ? Diagonal(SA[cis(φ), cis(φ)]) : Diagonal(SA[cis(φ), cis(φ), cis(-φ), cis(-φ)])
    eφ_hoppings(t, r, dr, B_par, B_perp) = t * eφ(r, dr, B_par, B_perp)

    ## Pairing phases
    #We put φSC/2 to the left and -φSC/2 to the right. This is why we have a 4 factor
    eφ_onsite(o, r, φSC, Bpar, Bperp) =
        eφ(r, r + [0, -W/2, 0], Bpar, Bperp) * diagphi(sign(r[1])*φSC/4) * o * diagphi(sign(r[1])*φSC/4)' *
        eφ(r, r + [0,  -W/2, 0], Bpar, Bperp)'
    deφ_onsite(o, r, φSC, Bpar, Bperp) =
        sign(r[1]) * im/4 * eφ(r,  r + [0, -W/2, 0], Bpar, Bperp) * diagphi(sign(r[1])*φSC/4) *
        (σ0τz * o - o * σ0τz) * diagphi(sign(r[1])*φSC/4)' * eφ(r, r + [0,  -W/2, 0], Bpar, Bperp)'

    # ALTERNATIVE: Now we put φSC on the left and 0 on the right
    # eφ_onsite(o, r, φSC, Bpar, Bperp) =
    #     eφ(r, r, Bpar, Bperp) * diagphi(ifelse(sign(r[1])<0, φSC/2, 0)) * o *
    #     diagphi(ifelse(sign(r[1])<0, φSC/2, 0))' * eφ(r, r, Bpar, Bperp)'
    # deφ_onsite(o, r, φSC, Bpar, Bperp) =
    #     sign(r[1]) * im/4 * eφ(r, r, Bpar, Bperp) * diagphi(ifelse(sign(r[1])<0, φSC/2, 0)) *
    #     (σ0τz * o - o * σ0τz) * diagphi(ifelse(sign(r[1])<0, φSC/2, 0))' * eφ(r, r, Bpar, Bperp)'

    ## Ripples: length scale of 10*a0
    Random.seed!(5);
    kx, ky = rand(0:1e-5:2*2π/10,10,2)
    Δz(r) = [0,0,0*0.5*a0/length(kx) * sum(sin.((kx.*r[1]+ky.*r[2])/a0))] #disabled

    # smooth function describing the band bending due to the SC. λN is the width of the tanh,
    # xμ is the center of the sigmoidal. We want xμ to be inside the normal region
     smoothmu(r) = μN/2*(tanh((r[1]+(Ln/2 - Δxμ))/λμ)-tanh((r[1]-(Ln/2 - Δxμ))/λμ)) +
         ifelse(r[1] < 0, μSL, μSR)/2*(2-tanh((r[1]+(Ln/2 - Δxμ))/λμ)+tanh((r[1]-(Ln/2 - Δxμ))/λμ)) +
         Δμbarrier * ((sigmoidal(r[1] - Ln/2, λμbarrier, λμbarrier/2) * sigmoidal(-r[1] + Ln/2, λμbarrier, λμbarrier/2)) +
                     (sigmoidal(r[1] + Ln/2, λμbarrier, λμbarrier/2) * sigmoidal(-r[1] - Ln/2, λμbarrier, λμbarrier/2)))

    ## CAR processes on or not. Interlayer hoppings suppressed in the leads if CAR == false
    hopSCcondition(r,dr) = ifelse(screg(r - dr/2) == true || screg(r + dr/2) == true ,
        ifelse(car == false, 0, 1), 1)
    ## Transparencies between hopping sites at the SC-N left and right interfaces
    function transparency(r, dr, τL, τR)
        if xor(scRreg(r - dr/2), scRreg(r +  dr/2))
            τR
        elseif xor(scLreg(r - dr/2), scLreg(r +  dr/2))
            τL
        else 1.0 end
    end

    ## Different Δdimer energy (top and bottom):
    Δdimer_top! = @onsite!((o,r) -> o + U_dimer * σ0τz; region = normalreg, sublats = (:At))
    Δdimer_bot! = @onsite!((o,r) -> o + U_dimer * σ0τz; region = normalreg, sublats = (:Bb))

    ###  BOUNDED MODIFIERS
    onsiteNtop! = @onsite!((o, r; μN = μN, U = U) -> (-ifelse(smoothμ == false, μN, smoothmu(r)) + U/2 * ifelse(smooth == true,
        smoothU(r, W, Ln, dir),1)) * σ0τz + sum(σvect .* Vzeem);
        region = normalreg, sublats = (:At, :Bt))
    onsiteNbot! = @onsite!((o, r; μN = μN, U = U) -> (- ifelse(smoothμ == false, μN, smoothmu(r)) - U/2 *
        ifelse(smooth == true, smoothU(r, W, Ln, dir),1)) * σ0τz +
        sum(σvect .* Vzeem); region = normalreg, sublats = (:Ab, :Bb))

    cantingNAtBb! = @onsite!((o) -> o + sum(σvect .* Vuvect);
        region = normalreg, sublats= (:At,:Bb))
    cantingNAbBt! = @onsite!((o) -> o - sum(σvect .* Vuvect);
        region = normalreg, sublats= (:Ab,:Bt))
    field! = @onsite!((o, r) -> o + [Ex, Ey, Ez]'r * σ0τz)
    onsiteS! = @onsite!((o,r; Δ=Δ) -> - ifelse(smoothμ == false, ifelse(scLreg(r), μSL, μSR), smoothmu(r)) * σ0τz - Δ  * σyτy; region = screg)
    peierls! =@hopping!((t, r, dr; B_par=B_par, B_perp=B_perp) ->
        eφ_hoppings(t,r,dr,B_par,B_perp))

    scphases! = @onsite!((o, r; phi=0, B_par=B_par, B_perp=B_perp) ->
         eφ_onsite(o,r,phi,B_par,B_perp);  region = ifelse(localselfenergy == false &&
         analyticalselfenergy == false, screg, edgereg))

    dpeierls! = @onsite!((o, r; phi=0, B_par=B_par, B_perp=B_perp) ->
        deφ_onsite(o,r,phi,B_par,B_perp); region = screg)

    transparency! = @hopping!((t,r,dr) -> t * transparency(r, dr, τL, τR); range = 4a0)

    pairing! = @onsite!((o, r) -> -Δ * σyτy; region = screg)

    localselfenergy! = @onsite!((o,r) -> ifelse(localselfenergy == false, o, o - Δ * σyτy); region = edgereg)
    # Self energy in the static and Andreev limit:
    Σon() = g1_on * σyτy # The Fermi level is renormalized by E -> E - δμ_N term
    Σhop(dr) = 0 * σ0τz + g1_hop * sin(dr * kfs)/dr * σyτy
    # I've neglected the self energy corrections to the normal hopping above
    analyticalselfenergyon! = @onsite!((o,r) ->
        ifelse(analyticalselfenergy == false, o, o + Σon()); region = edgereg)

    analyticalselfenergyhop! = @hopping!((t, r, dr) ->
        ifelse(analyticalselfenergy == false, t, t + Σhop(norm(dr)));
         sublats = (:At =>:At, :Ab=>:Ab, :Bt=>:Bt, :Bb=>:Bb, :At=>:Bb,:Bb=>:At), range = dinter)

    ## Bounded modifiers required for the unbounded phams (note that for infinite y Bperp has to be 0)
    # in order to have a periodic sys.
    onsiteNtop_ribbon! = @onsite!((o, r; μN=μN,U=U, V_par = V_par, V_perp = V_perp) ->
        (-ifelse(smoothμ == false, μN, smoothmu(r))  + U/2 * ifelse(smooth == true, smoothU(r, W, Ln, dir),1) ) * σ0τz
        + V_par * ifelse(normal === true, σyτz, σyτ0´) + V_perp * σzτz; sublats = (:At, :Bt))
        #+ sum(σvect .* Vvect); sublats = (:At, :Bt))
    onsiteNbot_ribbon! = @onsite!((o, r; μN=μN,U=U , V_par = V_par, V_perp = V_perp)  ->
        (-ifelse(smoothμ == false, μN, smoothmu(r)) - U/2 * ifelse(smooth == true, smoothU(r, W, Ln, dir),1) ) * σ0τz
        + V_par * ifelse(normal === true, σyτz, σyτ0´) + V_perp * σzτz; sublats = (:Ab, :Bb))
        #+ sum(σvect .* Vvect); sublats = (:Ab,:Bb))

    onsiteNtop_inf! = @onsite!((o; μN=μN,U=U) -> o + (- μN + U/2 ) * σ0τz +
        sum(σvect .* Vvect); sublats = (:At, :Bt))
    onsiteNbot_inf! = @onsite!((o; μN=μN,U=U)  -> o + (- μN - U/2 ) * σ0τz +
        sum(σvect .* Vvect); sublats = (:Ab,:Bb))

    rashba! = @hopping!((t, r, dr; α = α) -> t + 2*im/3 * α * ifelse(screg(r),0.,1.)*
        sign(r[3])*( dr[2]/dintra * ifelse(normal === true, σxτz, σxτ0´) - dr[1]/dintra * σyτz);
        range = dintra, sublats = (:At=>:Bt, :Ab=>:Bb, :Bt=>:At, :Bb=>:Ab))

    scphasesnoB! = @onsite!((o, r; phi=0) -> eφ_onsite(o,r,phi,0,0);
        region = ifelse(localselfenergy == false && analyticalselfenergy == false, screg, edgereg))
    scphasesnoBperp! = @onsite!((o, r; phi=0, B_par=B_par) ->
        eφ_onsite(o,r,phi,B_par,0); region = ifelse(localselfenergy == false &&
        analyticalselfenergy == false, screg, edgereg))

    interhopSCt1! = @hopping!((t, r, dr) -> t + hopSCcondition(r,dr) * t3 * σ0τz;
        range = dinter_warp, sublats = (:Bt=>:Ab, :Ab=>:Bt))
    interhopSCt3! = @hopping!((t, r, dr) -> t + hopSCcondition(r,dr) * t1 * σ0τz;
        range = dinter, sublats = (:At=>:Bb,:Bb=>:At))
    peierlsnoBperp! = @hopping!((t, r, dr; B_par=B_par) -> eφ_hoppings(t,r,dr,B_par,0))

    ### Build models N-Tuple
    models = (; model0, model00, onsiteNtop!, onsiteNbot!, onsiteNtop_inf!, onsiteNbot_inf!,
              onsiteS!, scphases!, scphasesnoB!, peierls!, dpeierls!,
              scLreg, scRreg, modelNtop, modelNbot, modelS, transparency!, Δdimer_top!, pairing!,
              Δdimer_bot!,scphasesnoBperp!, peierlsnoBperp!, cantingNAtBb!, cantingNAbBt!,
              onsiteNtop_ribbon!, localselfenergy!, analyticalselfenergyon!,
              analyticalselfenergyhop!, rashba!, ising!, onsiteNbot_ribbon!,
              interhopSCt1!,  interhopSCt3!, field!)
    return models
end

"shape of the smooth function f(x,y) s.t U(x, y) = U * f(r) if smooth == true"
function smoothU(r, W, Ln, dir)
    ξ = 10
    if dir == :x
        1 - sigmoidal(r[1], Ln/4, ξ) * sigmoidal(-r[1], Ln/4, ξ)
    elseif dir == :y
        1 - sigmoidal(-r[2], -W/4, ξ) * sigmoidal(r[2], 3W/4, ξ)
    elseif dir == :xy
        1 - sigmoidal(r[1], Ln/4, ξ) * sigmoidal(-r[1], Ln/4, ξ) *
        sigmoidal(-r[2], -W/4, ξ) * sigmoidal(r[2], 3W/4, ξ)
    end
end

sigmoidal(r, length, ξ) = 1/(1 + exp((r-length)/ξ))

"Rotation matrix of an angle θ in around x axis"
Rx(Rvect, θ) = @SMatrix[1 0 0; 0 cos(θ) -sin(θ); 0 sin(θ) cos(θ)] * Rvect
"Rotation matrix of an angle θ in around z axis"
Rz(Rvect, θ) = @SMatrix[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * Rvect

"Unbounded lattice takes `p:Params`"
function lat_infinite(p = Params())
    @unpack a0, dinter, scale, edges = p
    a0 = a0 * scale
    if edges == :armchair_zigzag
        sAbot = sublat((0.0,-1.0a0/sqrt(3.0),  - dinter/2); name = :Ab)
        sBbot = sublat((0.0, 0.0a0/sqrt(3.0),  - dinter/2); name = :Bb)
        sAtop = sublat((0.0, 0.0a0/sqrt(3.0),  + dinter/2); name = :At)
        sBtop = sublat((0.0, 1.0a0/sqrt(3.0),  + dinter/2); name = :Bt)
        br = a0 * SA[cos(pi/3) sin(pi/3) 0; -cos(pi/3) sin(pi/3) 0]'
    elseif edges == :zigzag_armchair
        sAbot = sublat((-1.0a0/sqrt(3.0),0.0, - dinter/2); name = :Ab)
        sBbot = sublat((0.0a0/sqrt(3.0),0.0,  - dinter/2); name = :Bb)
        sAtop = sublat((0.0a0/sqrt(3.0),0.0,  + dinter/2); name = :At)
        sBtop = sublat((1.0a0/sqrt(3.0),0.0,  + dinter/2); name = :Bt)
        br = a0 * Rz(SA[cos(pi/3) sin(pi/3) 0; -cos(pi/3) sin(pi/3) 0]',π/2)
    end
    return lattice(sAtop, sBtop, sAbot, sBbot; bravais = br)
end

# Be careful below I have changed the names to force BdG symmetry but it is not
# intuitive. Simplify this
sigmas(::Val{true}) = [σ0´, σ0´, σx´, σy´, σz´, 0*σ0´], @SVector[σx´, σy´, σz´]
sigmas(::Val{false}) = [σ0τ0´, σ0τz´, σxτz´, σyτz´, σzτz´, σyτy´],
    @SVector[σxτz´, σyτ0´, σzτz´]


text:
@testset "binning.jl" begin
    @testset "EqualSize" begin
        # default constructor
        @test EqualSize().n == 10

        # exception
        @test_throws ErrorException EqualSize(; n=0)
        @test_throws ErrorException EqualSize(; n=-1)

        # binning
        nbins = 11
        b = ReliabilityDiagrams.bins(rand(100), rand(100), EqualSize(; n=nbins))
        @test b isa Histogram{Int}
        @test length(b.weights) == nbins
        @test length(b.edges) == 1
        @test first(b.edges) == range(0; stop=1, length=nbins + 1)

        @test ReliabilityDiagrams.binindex(b, prevfloat(0.0)) == 0
        @test ReliabilityDiagrams.binindex(b, 0) == 1
        @test ReliabilityDiagrams.binindex(b, 0.05) == 1
        @test ReliabilityDiagrams.binindex(b, 0.1) == 2
        @test ReliabilityDiagrams.binindex(b, 1) == nbins
        @test ReliabilityDiagrams.binindex(b, nextfloat(1.0)) == nbins + 1

        @test ReliabilityDiagrams.max_nonempty_bins(b, 1) == 1
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins - 1) == nbins - 1
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins) == nbins
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins + 1) == nbins
    end

    @testset "EqualMass" begin
        # default constructor
        @test EqualMass().n == 10

        # exception
        @test_throws ErrorException EqualMass(; n=0)
        @test_throws ErrorException EqualMass(; n=-1)

        # binning
        x = rand(100)
        nbins = 17
        b = ReliabilityDiagrams.bins(x, rand(100), EqualMass(; n=nbins))
        @test b isa Histogram{Int}
        @test length(b.weights) == nbins
        @test length(b.edges) == 1
        @test first(b.edges) ≈ quantile(x, range(0; stop=1, length=nbins + 1))

        @test ReliabilityDiagrams.binindex(b, prevfloat(minimum(x))) == 0
        @test ReliabilityDiagrams.binindex(b, minimum(x)) == 1
        @test ReliabilityDiagrams.binindex(b, quantile(x, 0.05)) == 1
        @test ReliabilityDiagrams.binindex(b, quantile(x, 0.1)) == 2
        @test ReliabilityDiagrams.binindex(b, maximum(x)) == nbins
        @test ReliabilityDiagrams.binindex(b, nextfloat(maximum(x))) == nbins + 1

        @test ReliabilityDiagrams.max_nonempty_bins(b, 1) == 1
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins - 1) == nbins - 1
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins) == nbins
        @test ReliabilityDiagrams.max_nonempty_bins(b, nbins + 1) == nbins
    end
end

ext: jl


text:
using Colors
using Images

function calc_image(width::UInt8, height::UInt8, bitwise_function)
    println("Starting calculation")
    image = zeros(RGB{Float64}, width, height)

    r::UInt8 = 0
    g::UInt8 = 126
    b::UInt8 = 255

    for x = 1:width, y = 1:height
        r = 0 + bitwise_function(x, y)
        g = 255 - bitwise_function(x, y)
        b = 64 + div(bitwise_function(x, y), 2)

        image[x, y] = RGB{N0f8}(reinterpret(N0f8, UInt8(r)), reinterpret(N0f8, UInt8(g)), reinterpret(N0f8, UInt8(b)))
    end

    println("Calculation complete")

    image
end

save("xor.png", calc_image(UInt8(255), UInt8(255), (x, y) -> xor(x, y)))
save("and.png", calc_image(UInt8(255), UInt8(255), (x, y) -> x & y))
save("or.png", calc_image(UInt8(255), UInt8(255), (x, y) -> x | y))

max_stars_repo_name: James-P-D/BitwiseImages

max_issues_repo_name: James-P-D/BitwiseImages

max_forks_repo_name: James-P-D/BitwiseImages


text:
@testset "pointer" begin

# constructors

voidptr_a = CuPtr{Cvoid,AS.Host}(Int(0xDEADBEEF))
@test voidptr_a == CuHostPtr{Cvoid}(Int(0xDEADBEEF))

voidptr_b = CuPtr{Cvoid,AS.Host}(Int(0xCAFEBABE))
@test voidptr_b == CuHostPtr{Cvoid}(Int(0xCAFEBABE))

# getters
@test eltype(voidptr_a) == Cvoid

# comparisons
@test voidptr_a != voidptr_b

@testset "conversions" begin

# between regular and CUDA pointers

@test_throws ArgumentError convert(Ptr{Cvoid}, voidptr_device_a)

# between CUDA pointers

intptr_a = CuPtr{Int,AS.Host}(Int(0xDEADBEEF))
@test convert(typeof(intptr_a), voidptr_a) == intptr_a

end

end


text:
module CoinGeckoAPI

# dependacies

using HTTP
using JSON3

# exports
export ping
export get_price
export get_token_price
export get_supported_vs_currencies
export get_coins_list
export get_coins_markets
export get_coin_by_id
export get_coin_ticker_by_id
export get_coin_history_by_id
export get_coin_market_chart_by_id
export get_coin_market_chart_range_by_id
export get_coin_status_updates_by_id

# globals constants
const global url_base = "https://api.coingecko.com/api/v3/"

"""
    ping()
Function to check if the API is up and running

# Arguments
    None

# Returns
"""
function ping()
    apiurl ="ping"
    r =  HTTP.request("GET", url_base * apiurl)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    _api_url_params(api_url, params, api_url_has_params = false)

Function to build the API url with the parameters

# Arguments
    `api_url: string` : the API url
    `params: dict` : the parameters to be added to the API url
    `api_url_has_params: bool` : if the API url has parameters

# Returns
    `api_url: String` : the API url with the parameters

"""
function _api_url_params(api_url, params, api_url_has_params = false)

    if api_url_has_params == true
        api_url = api_url * "&"
    else
        api_url = api_url * "?"
    end

    for a_pair in params
        api_url *=  a_pair.first * "=" * a_pair.second * "&"
    end
    @show api_url[1:end-1]
    return api_url[1:end-1]
end

"""
    get_price(ids, vs_currencies, kargs...)

Get the price of a cryptocurrency in any other supported currencies needed

# Arguments
    `ids: list` : the cryptocurrency ids
    `vs_currencies: list` : the currencies to get the price in
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `prices: dict` : the prices of the cryptocurrencies in the currencies

# Example:

```julia
r = get_price("bitcoin", "usd", "include_market_cap" => "true", "include_24hr_vol" => "true", "include_24hr_change" => "true")
```

"""
function get_price(ids, vs_currencies, kargs...)
    apiurl = "simple/price"
    kwards= Dict(kargs)
    kwards["ids"] = ids
    kwards["vs_currencies"] = vs_currencies

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

function get_token_price(id, contract_addresses, vs_currencies, kargs...)
    apiurl = "simple/token_price"
    kwards= Dict(kargs)
    kwards["id"] = id
    kwards["contract_addresses"] = contract_addresses
    kwards["vs_currencies"] = vs_currencies

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_supported_vs_currencies()

Get list of supported_vs_currencies

# Arguments
    None

# Returns
    `supported_vs_currencies: Dict` : the list of supported_vs_currencies
"""
function get_supported_vs_currencies()
    apiurl = "simple/supported_vs_currencies"
    r =  HTTP.request("GET", url_base * apiurl)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coins(kargs...)

List all coins with data (name, price, market, developer, community, etc)

# Arguments
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coins: Dict` : the list of coins with id, name and symbol

"""
function get_coins_list(kargs...)
    apiurl = "coins/list"
    kwards= Dict(kargs)

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coins_markets(vs_currencies, kargs...)

Use this to obtain all the coins market data (price, market cap, volume)

# Arguments
    `vs_currencies: list` : the currencies to get the price in
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coins_markets: Dict` : the list of coins with market data

# Example:

```julia
r = get_coins_markets("usd")

````
"""
function get_coins_markets(vs_currency, kargs...)
    apiurl = "coins/markets"
    kwards= Dict(kargs)
    kwards["vs_currency"] = vs_currency

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_by_id(id, kargs...)

Get current data (name, price, market, including exchange tickers for a coin)

# Arguments
    `id: string` : the coin id
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin: Dict` : the current data for a coin

"""
function get_coin_by_id(id, kargs...)
    apiurl = "coins"
    kwards= Dict(kargs)
    kwards["id"] = id

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_ticker_by_id(id, kargs...)

Get coin tickers (paginated to 100 items)

# Arguments
    `id: string` : the coin id
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin_ticker: Dict` : the coin tickers

"""
function get_coin_ticker_by_id(id, kargs...)
    apiurl = "coins/$id/tickers"
    kwards= Dict(kargs)
    kwards["id"] = id

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_history_by_id(id, date, kargs...)

Get Historical data (name, price, market, stats ) at a given date for a coin.

# Arguments
    `id: string` : the coin id
    `date: string` : the date to get the data for
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin_history: Dict` : the historical data for a coin

# Example:
```julia
r = get_coin_history_by_id("bitcoin", "30-12-2017")

````
"""
function get_coin_history_by_id(id, date, kargs...)
    apiurl = "coins/$id/history"
    kwards= Dict(kargs)
    kwards["id"] = id
    kwards["date"] = date

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_market_chart_by_id(id, vs_currency, days, kargs...)

Get historical market data include price, market cap, and 24h volume (granularity auto)
Minutely data will be used for duration within 1 day, Hourly data will be used for duration between 1 day and 90 days,
    Daily data will be used for duration above 90 days.

# Arguments
    `id: string` : the coin id
    `vs_currency: string` : the currency to get the data in
    `days: String` : the number of days to get the data for
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin_market_chart: Dict` : the historical market data for a coin

# Example:
```julia
r = get_coin_market_chart_by_id("bitcoin", "usd", "30")
```

"""
function get_coin_market_chart_by_id(id, vs_currency, days, kargs...)
    apiurl = "coins/$id/market_chart"
    kwards= Dict(kargs)
    kwards["id"] = id
    kwards["vs_currency"] = vs_currency
    kwards["days"] = days

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_market_chart_range_by_id(id, vs_currency, from_timestamp, to_timestamp, kargs...)

Get historical market data include price, market cap, and 24h volume within a range of timestamp (granularity auto)
Data granularity is automatic (cannot be adjusted)
1 day from query time = 5 minute interval data
1 - 90 days from query time = hourly data
above 90 days from query time = daily data (00:00 UTC)

# Arguments
    `id: string` : the coin id
    `vs_currency: string` : the currency to get the data in
    `from_timestamp: String` : the timestamp to start the data from (UNIX Timestamp (eg 1392577232))
    `to_timestamp: String` : the timestamp to end the data from (UNIX Timestamp (eg 1422577232))
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin_market_chart: Dict` : the historical market data for a coin

# Example:
```julia
r = get_coin_market_chart_range_by_id("bitcoin", "usd", "1392577232", "1422577232")
```
"""
function get_coin_market_chart_range_by_id(id, vs_currency, from_timestamp, to_timestamp, kargs...)
    apiurl = "coins/$id/market_chart/range"
    kwards= Dict(kargs)
    # kwards["id"] = id
    kwards["vs_currency"] = vs_currency
    kwards["from"] = from_timestamp
    kwards["to"] = to_timestamp

    api_url = _api_url_params(apiurl, kwards)

    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

"""
    get_coin_status_updates_by_id(id, kargs...)

Get status updates for a given coin.

# Arguments
    `id: string` : the coin id  (eg bitcoin)
    `kargs: dict` : the parameters to be added to the API url

# Returns
    `coin_status_updates: Dict` : the status updates for a coin

"""
function get_coin_status_updates_by_id(id, kargs...)
    apiurl = "coins/$id/status_updates"
    kwards= Dict(kargs)
    kwards["id"] = id

    api_url = _api_url_params(apiurl, kwards)
    r =  HTTP.request("GET", url_base * api_url)
    r = String(r.body)
    return JSON3.read(r)

end

end


text:
import FractionalDiffEq.solve, FractionalDiffEq.FractionalDiffEqAlgorithm

using LinearAlgebra, InvertedIndices

"""
Using [triangular strip matrices](https://en.wikipedia.org/wiki/Triangular_matrix) to discrete fractional differential equations to simple algebra system and solve the system.

```tex
@inproceedings{Podlubny2000MATRIXAT,
  title={MATRIX APPROACH TO DISCRETE FRACTIONAL CALCULUS},
  author={Igor Podlubny},
  year={2000}
}
```
"""
struct MatrixDiscrete <: FractionalDiffEqAlgorithm end

"""
@inproceedings{Podlubny1998FractionalDE,
  title={Fractional differential equations},
  author={Igor Podlubny},
  year={1998}
}
"""

"""
@inproceedings{Podlubny2000MATRIXAT,
  title={MATRIX APPROACH TO DISCRETE FRACTIONAL CALCULUS},
  author={Igor Podlubny},
  year={2000}
}
"""

"""
    solve(p1, α, p2, c, h, T, MatrixDiscrete())

Using the **Matrix Discretization algorithm** proposed by [Prof Igor Podlubny](http://people.tuke.sk/igor.podlubny/index.html) to approximate the numerical solution.
"""
function solve(p1, α, p2, c, h, T, ::MatrixDiscrete)
    n=Int64(floor(α))
    rows=collect(1:n)
    N=Int64(T/h+1)
    equation = p1*B(N, α, h)+p2*(zeros(N, N)+I)
    equation = eliminator(N, rows)*equation*eliminator(N, rows)'
    righthand = c*eliminator(N, rows)*ones(N)
    result = equation\righthand
    return vcat(zeros(n), result)
end

function solve(equation, right, h, T)
    N=Int64(T/h+1)
    equation = eliminator(N, rows)*equation*eliminator(N, rows)'
    right = right*eliminator(N, rows)*ones(N)
    result = equation\right
    return vcat(zeros(n), result)
end

"""
    eliminator(n, row)

Compute the eliminator matrix Sₖ by omiting n-th row
"""
function eliminator(n, row)
    temp = zeros(n, n)+I
    return temp[Not(row), :]
end

"""
Generating elements in Matrix.
"""
function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    for i in range(1, n, step=1)
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega

end

function B(N, p, h)
    result=zeros(N, N)
    temp=omega(N, p)

    for i in range(1, N, step=1)
        result[i, 1:i]=reverse(temp[1:i])
    end

    return h^(-p)*result
end

function F(N, p, h)
    result=zeros(N, N)
    temp = omega(N, p)

    for i in range(1, N, step=1)
        result[i, 1:i]=reverse(temp[1:i])
    end

    result=reverse(reverse(result, dims=1), dims=2)

    return (-1)^(ceil(p))*h^(-p)*result
end

"""
    bagleytorvik(p1, p2, p3, T, h)

By specifying the parameters of Bagley Torvik Equation, we can use **bagleytorvik** to directly obtain the numerical approximation.

!!! info p2 ≠ 0
    Please note that the parameter of fractional derivative item must not be 0
"""
function bagleytorvik(p1, p2, p3, T, h)
    N=Int64(T/h+1)
    equation = p1*B(N, 2, h)+p2*B(N, 1.5, h)+p3*(zeros(N, N)+I)
    equation = eliminator(N, [1,2])*equation*eliminator(N, [1,2])'
    right = eliminator(N, [1,2])*ones(N)
    result = equation\right

    return vcat(zeros(2), result)
end


text:
@testset "AVLTree" begin
    @testset "inserting values" begin
        t = AVLTree{Int}()
        for i in 100:-1:1
            insert!(t, i)
        end

        @test length(t) == 100

        for i in 1:100
            @test haskey(t, i)
        end

        for i = 101:150
            @test !haskey(t, i)
        end
    end

    @testset "deleting values" begin
        t = AVLTree{Int}()
        for i in 1:100
            insert!(t, i)
        end
        for i in 1:2:100
            delete!(t, i)
        end

        @test length(t) == 50

        for i in 1:100
            if iseven(i)
                @test haskey(t, i)
            else
                @test !haskey(t, i)
            end
        end

        for i in 1:2:100
            insert!(t, i)
        end

        @test length(t) == 100
    end

    @testset "handling different cases of delete!" begin
        t2 = AVLTree{Int}()
        for i in 1:100000
            insert!(t2, i)
        end

        @test length(t2) == 100000

        nums = rand(1:100000, 8599)
        visited = Set()
        for num in nums
            if !(num in visited)
                delete!(t2, num)
                push!(visited, num)
            end
        end

        for i in visited
            @test !haskey(t2, i)
        end
        @test (length(t2) + length(visited)) == 100000
    end

    @testset "handling different cases of insert!" begin
        nums = rand(1:100000, 1000)
        t3 = AVLTree{Int}()
        uniq_nums = Set(nums)
        for num in uniq_nums
            insert!(t3, num)
        end
        @test length(t3) == length(uniq_nums)
    end

    @testset "in" begin
        t4 = AVLTree{Char}()
        push!(t4, 'a')
        push!(t4, 'b')
        @test length(t4) == 2
        @test in('a', t4)
        @test !in('c', t4)
    end

    @testset "search_node" begin
        t5 = AVLTree{Int}()
        for i in 1:32
            push!(t5, i)
        end
        n1 = search_node(t5, 21)
        @test n1.data == 21
        n2 = search_node(t5, 35)
        @test n2.data == 32
        n3 = search_node(t5, 0)
        @test n3.data == 1
    end

    @testset "getindex" begin
        t6 = AVLTree{Int}()
        for i in 1:10
            push!(t6, i)
        end
        for i in 1:10
            @test getindex(t6, i) == i
        end
        @test_throws BoundsError getindex(t6, 0)
        @test_throws BoundsError getindex(t6, 11)
    end

    @testset "key conversion in push!" begin
        t7 = AVLTree{Int}()
        push!(t7, Int8(1))
        @test length(t7) == 1
        @test haskey(t7, 1)
    end

    @testset "minimum_node" begin
        t8 = AVLTree{Int}()
        @test minimum_node(t8.root) == nothing
        for i in 1:32
            push!(t8, i)
        end
        m1 = minimum_node(t8.root)
        @test m1.data == 1
        node = t8.root
        while node.leftChild != nothing
            m = minimum_node(node.leftChild)
            @test m == m1
            node = node.leftChild
        end
    end

    @testset "rank" begin
        t9 = AVLTree{Int}()
        for i in 1:20
            push!(t9, i)
        end
        for i in 1:20
            @test sorted_rank(t9, i) == i
        end
        @test_throws KeyError sorted_rank(t9, 21)
    end
end


text:
using JuLIP, JuLIPMaterials

sym = :Cu
r0 = rnn(sym)
V = lennardjones(r0 = r0)

CC1 = []
for R in [3, 4, 6, 9, 13, 20]
   athom = cluster(sym, R*r0+2*cutoff(V))              # homogeneous reference configuration
   at = deleteat!(deepcopy(athom), 1)   # create a vacancy
   x̄ = athom[1]
   r = [norm(x-x̄) for x in positions(at)]
   set_constraint!(at, FixedCell(at, free = find(r .< R *r0)))
   set_calculator!(at, V)

   C0, C1 = JuLIPMaterials.PoleExpansion.dipole_tensor(at, athom;
            method = :lbfgs, precond = FF(at, V), verbose=0, gtol = 1e-4)
   println("R = $R, C0 = $C0")
   println("        C1 = $(round.(C1,2))")
   push!(CC1, C1)
end

CC1

err = [ norm(C - CC1[end]) for C in CC1[1:end-1] ]
RR = [3, 4, 6, 9, 13] * r0

using Plots
Plots.gr()
P = Plots.plot(RR, err, lw=2, ms=12, m=:o, label = "err C1",
         xaxis = (:log, [0.8*RR[1],1.2*RR[end]]), yaxis = (:log, [1e-3, 4]))
Plots.plot!(RR, 10*RR.^(-2.5), lw=1, ls = :dash, ms = 0, label = "~ R^{-2.5}")
# Plots.plot!(RR, RR.^(-2.5), lw=1, linestyle = :dash, ms = 0, label = "~ R^{-2.5}")

;

# # TEST DECAY
#
# R = 10
# athom = cluster(sym, R*r0 + 2*cutoff(V))
# Xhom = positions(athom)
# x̄ = Xhom[1]
# at = deleteat!(deepcopy(athom), 1)
# @show length(at)
# X0 = positions(at)
# @assert norm(x̄ - mean(X0)) < 1e-10
# r = [ norm(x - x̄) for x in X0 ]
# set_constraint!( at, FixedCell(at, free = find(r .< R*r0)) )
# set_calculator!( at, V )
# minimise!(at, method = :lbfgs, precond = FF(at, V), verbose=2, gtol = 1e-6)
# X = positions(at)
# U = X - X0
#
# e = strains(U, set_positions!(at, X0); rcut = 1.5 * r0)
#
# h = 1e-5
# F1 = (forces(set_positions!(at, X0 + h*U)) -
#       forces(set_positions!(at, X0 - h*U)) ) / (2*h)
# Iin = find(r .< 6*r0)
# sum(F1[Iin])
# sum( f * (x - x̄)'   for (f, x) in zip(F1[Iin], X0[Iin]) )
#
#
# Uhom = [ [zero(JVecF)]; U ]
# X0hom = positions(athom)
# F2 = (forces(V, set_positions!(athom, Xhom + h*Uhom)) -
#       forces(V, set_positions!(athom, Xhom - h*Uhom)) ) / (2*h)
# sum(F2[Iin])
# sum( f * (x - x̄)'   for (f, x) in zip(F2[Iin], X0hom[Iin]) )
#
#
#
# using Plots
# gr()
# t = [extrema(r)...]
# P = plot(r, 1e-15+e, lw = 0, m=:o, ms=2,
#          xaxis = (:log, [0.8,1.2].*t),
#          yaxis = (:log, [1e-5, 1e-2]))
# plot!(P, t, t.^(-3), lw=1)
#
#
#
#
#
#
# R = 10
# athom = cluster(sym, R*r0 + 2*cutoff(V))
# x̄ = at[1]
# at = deleteat!(deepcopy(at), 1)
# X0 = positions(at)
# set_calculator!(at, V)
# C0, C1 = JuLIPMaterials.PoleExpansion.dipole_tensor(at, athom)
# X = positions(at)
# maximum(norm.(X - X0))
# U = X - X0
# e = strains(U, set_positions!(at, X0); rcut = 1.5 * r0)
# r = [ norm(x - x̄) for x in X0 ]
# t = [extrema(r)...]
# P = plot(r, 1e-15+e, lw = 0, m=:o, ms=2,
#          xaxis = (:log, [0.8,1.2].*t),
#          yaxis = (:log, [1e-5, 1e-2]))
# plot!(P, t, t.^(-3), lw=1)


text:
module SparsityDetection

using SpecialFunctions
using Cassette, LinearAlgebra, SparseArrays
using Cassette: tag, untag, Tagged, metadata, hasmetadata, istagged, canrecurse
using Cassette: tagged_new_tuple, ContextTagged, BindingMeta, DisableHooks, nametype
using Core: SSAValue

export Sparsity, hsparsity, sparsity!

include("program_sparsity.jl")
include("sparsity_tracker.jl")
include("path.jl")
include("take_all_branches.jl")
include("terms.jl")
include("linearity.jl")
include("hessian.jl")
include("blas.jl")
include("linearity_special.jl")

end


text:
using MetaGraphsNext
import Documenter: doctest

doctest(MetaGraphsNext)


text:
using MPI
using PETSc_jll
using Libdl

export Scalar
export comm_type
export MPI_Comm
export PetscErrorCode
export PetscBool
export PetscInt
export have_petsc, petsc_libs, petsc_type
const depfile = joinpath(dirname(@__FILE__), "..", "..", "deps", "deps.jl")
isfile(depfile) || error("GridapDistributedPETScWrappers not properly installed. Please run Pkg.build(\"GridapDistributedPETScWrappers\")")
include(depfile)

if (deps_file_petscRealDouble == "")
  const petscRealDouble=PETSc_jll.libpetsc_path
else
  const petscRealDouble=deps_file_petscRealDouble
end
const petscRealSingle=""
const petscComplexDouble=""

const libs=[petscRealDouble]
const expected_petsc_scalar_types=[Float64]
const expected_petsc_real_types=[Float64]

function initialize(libhdl::Ptr{Cvoid})
  PetscInitializeNoArguments_ptr = dlsym(libhdl, :PetscInitializeNoArguments)
  result=ccall(PetscInitializeNoArguments_ptr, PetscErrorCode, ())
  @assert result==0
  return nothing
end

function DataTypeFromString(libhdl::Ptr{Cvoid}, name::AbstractString)
    PetscDataTypeFromString_ptr = dlsym(libhdl, :PetscDataTypeFromString)
    dtype_ref = Ref{UInt32}()
    found_ref = Ref{UInt32}()
    result=ccall(PetscDataTypeFromString_ptr, PetscErrorCode,
             (Cstring, Ptr{UInt32}, Ptr{UInt32}),
             name, dtype_ref, found_ref)
    @assert result==0
    @assert found_ref[] == UInt32(1)
    return dtype_ref[]
end

function PetscDataTypeGetSize(libhdl::Ptr{Cvoid}, dtype::UInt32)
    PetscDataTypeGetSize_ptr = dlsym(libhdl, :PetscDataTypeGetSize)
    datasize_ref = Ref{Csize_t}()
    result=ccall(PetscDataTypeGetSize_ptr, PetscErrorCode,
             (PetscDataType, Ptr{Csize_t}),
             dtype, datasize_ref)
    @assert result==0
    return datasize_ref[]
end

const libtypes = map(libs) do lib
    libhdl = dlopen(lib)
    initialize(libhdl)
    PETSC_REAL = DataTypeFromString(libhdl, "Real")
    PETSC_SCALAR = DataTypeFromString(libhdl, "Scalar")
    PETSC_INT_SIZE = PetscDataTypeGetSize(libhdl, PETSC_INT)

    PetscReal =
        PETSC_REAL == PETSC_DOUBLE ? Cdouble :
        PETSC_REAL == PETSC_FLOAT ? Cfloat :
        error("PETSC_REAL = $PETSC_REAL not supported.")

    PetscScalar =
        PETSC_SCALAR == PETSC_REAL ? PetscReal :
        PETSC_SCALAR == PETSC_COMPLEX ? Comlex{PetscReal} :
        error("PETSC_SCALAR = $PETSC_SCALAR not supported.")

    PetscInt =
        PETSC_INT_SIZE == 4 ? Int32 :
        PETSC_INT_SIZE == 8 ? Int64 :
        error("PETSC_INT_SIZE = $PETSC_INT_SIZE not supported.")

    # TODO: PetscBLASInt, PetscMPIInt ?
    return (lib, PetscScalar, PetscReal, PetscInt)
end

for (i,lib) in enumerate(libs)
  if (expected_petsc_scalar_types[i] != libtypes[i][2])
    error("""Expected PetscScalar ($(expected_petsc_scalar_types[i]))
          does NOT match lib $(lib) PetscScalar ($(libtypes[i][2]))
          """)
  end
  if (expected_petsc_real_types[i] != libtypes[i][3])
    error("""Expected PetscReal ($(expected_petsc_real_types[i]))
          does NOT match lib $(lib) PetscReal ($(libtypes[i][3]))
          """)
  end
end

const PetscInt=libtypes[1][4]

const petsc_libs = [:petscRealDouble, :petscRealSingle, :petscComplexDouble]
const petsc_type = [Float64, Float32, ComplexF64]

const Scalar=Union{Float32, Float64, ComplexF64}

# some auxiliary functions used by ccall wrappers
# get an array of pointers to UInt8s that is the same shape as
# the Symbol array
# does *not* allocate the pointers
function symbol_get_before(sym_arr)
  ptr_arr = Array{Ptr{UInt8}}(undef,length(sym_arr))
#  println("ptr_arr = ", ptr_arr)
#  for i=1:length(sym_arr)
#    println("ptr_arr[$i] = ", ptr_arr[i])
#  end

  return pointer(ptr_arr), ptr_arr
end

# turn array of strings (UInt8 *) and  puts them into Symbol array
function symbol_get_after(ptr, sym_arr)
  ptr_arr = unsafe_wrap(Array, ptr, length(sym_arr))

  for i=1:length(sym_arr)
    sym_arr[i] = Symbol(unsafe_string(ptr_arr[i]))
  end

end

function symbol_set_before(sym_arr)
  str_arr = similar(sym_arr, String)

  for i=1:length(str_arr)
    str_arr[i] = string(sym_arr[i])
  end

end

# TODO: auto-generate these
function PetscObjectComm(::Type{Float64},arg1::Ptr)
   ccall((:PetscObjectComm,petscRealDouble),comm_type,(Ptr{Cvoid},),arg1)
end
function PetscObjectComm(::Type{Float32},arg1::Ptr)
   ccall((:PetscObjectComm,petscRealSingle),comm_type,(Ptr{Cvoid},),arg1)
end
function PetscObjectComm(::Type{ComplexF64},arg1::Ptr)
   ccall((:PetscObjectComm,petscComplexDouble),comm_type,(Ptr{Cvoid},),arg1)
end


text:
@everywhere function fake_slow_condition(x, delay = 0.10, expected = 30)
  sleep(delay)
  sum(x) == expected ? x : false # Return x iff it fulfills the poperty
end

# Similar to Base.findfirst(predicate, A) but process the array "A" columnwise in parallel
# and opt out of further processing if found. Returns both the index at which found
# as well as the result returned.
@everywhere function pfind_colwise(predicate, A)
  np = nprocs()  # determine the number of processes available
  n = size(A, 2)
  i = 1
  results = Any[false for i in 1:n]
  found = false
  wasfound(idx) = found = true
  isfound() = found
  # function to produce the next work item from the queue.
  # in this case it's just an index to a column.
  nextidx() = (idx=i; i+=1; idx)
  @sync begin
    for p=1:np
      if p != myid() || np == 1
        @async begin
          while true
            idx = nextidx()
            if idx > n
              break
            end
            results[idx] = remotecall_fetch(p, predicate, A[:,idx])
            if isfound()
              break # Someone else already found it
            end
            if results[idx] != false
              wasfound(idx)
              break # We found it
            end
          end
        end
      end
    end
  end
  fidx = findfirst((r)->(r != false), results)
  (fidx, (fidx != 0 ? results[fidx] : nothing))
end


text:
module YorkFit

export fit

"""Calculate the weight of the error `σX`."""
ω(σX::Real)::Real = σX^-2
α{T<:Real}(ωX::T, ωY::T)::T = sqrt(ωX * ωY)
W{T<:Real}(ωX::T, ωY::T, r::T, b::T, α::T)::T =
    (ωX * ωY) / (ωX + b^2 * ωY - 2 * b * r * α)

"""Compute Σ(Ws * Xs) / Σ(Xs)"""
barvar{T<:Real}(Ws::Vector{T}, Xs::Vector{T})::T = sum(Ws .* Xs) / sum(Ws)

function β{T<:Real}(W::T, U::T, V::T, ωX::T, ωY::T, b::T, r::T, α)::T
    W * ((U / ωY) + (b * V / ωX) - (r / α) * (b * U + V))
end

"""Calculate `a` and `b` using the method of least squares. Returns the tuple
`(a, b)` where `a` is the intercept and `b` is the gradient."""
function lsq{T<:Real}(Xs::Vector{T}, Ys::Vector{T})::Tuple{T,T}
    A = hcat(ones(length(Xs)), Xs)
    m = (A' * A) \ A' * Ys

    (m[1], m[2])
end

"""
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σYs::Vector{T}, rs::Vector{T}; kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σYs::Vector{T}, r::T=T(0); kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σY::T, rs::Vector{T}; kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σY::T, r::T=T(0); kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σY::T, rs::Vector{T}; kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σY::T, r::T=T(0); kwargs...)
    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σYs::Vector{T}, r::T=T(0); kwargs...)

    fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σYs::Vector{T}, rs::Vector{T}; kwargs...)::Tuple{T,T,T,T}

Fits a slope and intercept to a set of observations using the method of York et
al. (2004).

Arguments
=========

- `Xs` -- A vector of observations.
- `Ys` -- A vector of observations related to `Xs` by `Ys = a + b * Xs`.
- `σXs` -- Errors in the observations in `Xs`. A vector of equal length to `Xs` or
  a scalar value if all the errors are equal.
- `σYs` -- Errors in the observations in `Ys`. A vector of equal length to `Ys` or
  a scalar value if all the errors are equal.
- `rs` -- Correlation coefficient between errors. A vector of equal length to
  `Xs` and `Ys` or a scalar value if the coefficient is constant. Defaults to 0.

Keyword Arguments
==================

- `niter=400` -- Number of iterations to calculate the slope `b` with.
- `tol=1E-15` -- Stopping criteria. When the difference between two successive
  iterations of `b` falls below this, the iterations stop.

Returns
=======

A 4-tuple containing `(a, b, σa, σb)` where `a` is the intercept, `b` is the
slope and `σa` and `σb` are the errors in the intercept and slope.

References
==========

York, D., Evensen, N.M., Martınez, M.L. and Delgado, J.D.B., 2004. Unified
equations for the slope, intercept, and standard errors of the best straight
line. American Journal of Physics, 72(3), pp.367-375.
"""
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σYs::Vector{T}, rs::Vector{T}; kwargs...) =
    fit(Xs, Ys, fill(σX, length(Xs)), σYs, rs; kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σYs::Vector{T}, r::T=T(0); kwargs...) =
    fit(Xs, Ys, fill(σX, length(Xs)), σYs, fill(r, length(Xs)); kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σY::T, rs::Vector{T}; kwargs...) =
    fit(Xs, Ys, σXs, fill(σY, length(Ys)), rs; kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σY::T, r::T=T(0); kwargs...) =
    fit(Xs, Ys, σXs, fill(σY, length(Ys)), fill(r, length(Xs)); kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σY::T, rs::Vector{T}; kwargs...) =
    fit(Xs, Ys, fill(σX, length(Xs)), fill(σY, length(Ys)), rs; kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σX::T, σY::T, r::T=T(0); kwargs...) =
    fit(Xs, Ys, fill(σX, length(Xs)), fill(σY, length(Ys)), fill(r, length(Xs)); kwargs...)
fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σYs::Vector{T}, r::T=T(0); kwargs...) =
    fit(Xs, Ys, σXs, σYs, fill(r, length(Xs)); kwargs...)
function fit{T<:Real}(Xs::Vector{T}, Ys::Vector{T}, σXs::Vector{T}, σYs::Vector{T}, rs::Vector{T}; niter=400, tol=1E-15)::Tuple{T,T,T,T}
    # TODO: Check for dimension mismatch of vectors
    (a, b) = lsq(Xs, Ys)

    ωXs = ω.(σXs)
    ωYs = ω.(σYs)
    αs = α.(ωXs, ωYs)

    Ws = zeros(length(Xs))
    βs = zeros(length(Xs))
    barX = 0
    barY = 0

    for i in 1:niter
        Ws = W.(ωXs, ωYs, rs, b, αs)
        barX = barvar(Ws, Xs)
        barY = barvar(Ws, Ys)
        Us = Xs - barX
        Vs = Ys - barY
        βs = β.(Ws, Us, Vs, ωXs, ωYs, b, rs, αs)

        b_old = b
        b = sum(Ws .* βs .* Vs) / sum(Ws .* βs .* Us)

        if abs((b - b_old) / b) < tol
            break
        end
    end

    a = barY - b * barX
    # Adjusted x values
    xs = barX + βs
    barx = sum(Ws .* xs) / sum(Ws)
    us = xs - barx

    σb = sqrt(sum(Ws .* us.^2)^-1)
    σa = sqrt(sum(Ws)^-1 + barx^2 * σb^2)

    return (a, b, σa, σb)
end

end # module


text:
module TestModels

    using LinearAlgebra

    export p_fun,
        model,
        set_locations,
        set_indices

    export HyperCube,
        Polytope,
        HyperSphere,
        OddShape

    """
        model(parms, args...; kwargs...)

    Model function for test models. It is an identity function.

    # Arguments

    - `parms`: a proposed location

    """
    model(parms, args...; kwargs...) = parms

    """
        HyperCube

    A hypercube oject divided into partitions.

    # Field

    - `p_bounds`: a vector of boundaries of partitions
    """
    struct HyperCube
        p_bounds::Vector{Float64}
    end

    """
        p_fun(location, hypercube::HyperCube, args...; kwargs...)

    Returns the qualitative pattern associated with a hypercube.

    # Arguments

    - 'location`: location of chain in the hypercube
    - `p_bounds`: boundaries of partitions for each dimension

    """
    function p_fun(location, hypercube::HyperCube, args...; kwargs...)
        p_bounds = hypercube.p_bounds
        nb = length(p_bounds)
        nd = length(location)
        vals = fill(-100, nd)
        for j in 1:nd
            for i in 1:(nb-1)
                if (location[j] ≥ p_bounds[i]) && (location[j] ≤ p_bounds[i+1])
                    vals[j] = i
                    continue
                end
            end
        end
        return vals
    end

    """
        Polytope

    A polytope object containing centroid locations for each partition

    # Fields

    - `location`: location of centroid for each partition
    """
    struct Polytope
        location::Vector{Float64}
    end

    """
        p_fun(location, points::Vector{Polytope}, args...; kwargs...)

    Returns the qualitative pattern associated with a polytope.

    # Arguments

    - 'location`: location of chain in the hypercube
    - `points`: a vector of polytopes containing locations

    """
    function p_fun(location, points::Vector{Polytope}, args...; kwargs...)
        distances = map(p -> norm(location .- p.location), points)
        _,idx = findmin(distances)
        return idx
    end

    """
        HyperSphere

    A hypersphere object with arbitrary dimensionality.

    # Fields

    - `location`: a vector representing the centroid of the hypersphere
    - `radius`: the radius of the hypersphere
    """
    mutable struct HyperSphere
        location::Vector{Float64}
        radius::Float64
    end

    function HyperSphere(;location, radius)
        return HyperSphere(location, radius)
    end

    function in_contact(v1::HyperSphere, v2::HyperSphere)
        dist = norm(v1.location - v2.location)
        return dist ≤ max(v1.radius, v2.radius)
    end

    function in_contact(v1::HyperSphere, proposal::Vector{<:Number})
        dist = norm(v1.location - proposal)
        return dist ≤ v1.radius
    end

    not_valid(v, vs)  = any(x -> in_contact(v, x), vs) || out_of_bounds(v)

    out_of_bounds(v) = any((v.location .- v.radius) .< 0) || any((v.location .+ v.radius) .> 1)

    """
        set_locations(r_dist, n_dims, n_obj)

    Sets the location of hypspheres within a unit hypercube

    # Arguments

    - `r_dist`: a function that samples a value for the radius
    - `n_dims`: the number of dimensions of the hypersphere
    - `n_obj`: the number of hyperspheres

    """
    function set_locations(r_dist, n_dims, n_obj)
        hs = [HyperSphere(;location = rand(n_dims), radius = r_dist())]
        for i in 1:(n_obj-1)
            h = HyperSphere(;location = rand(n_dims), radius = r_dist())
            while not_valid(h, hs)
                h.location = rand(n_dims)
            end
            push!(hs, h)
        end
        return hs
    end

    """
        p_fun(location, hyperspheres::Vector{HyperSphere}, args...; kwargs...)

    Returns the qualitative data pattern for hypspheres embedded in a hypercube. -100 is
    returned if the proposed location is not in a hypsphere. Otherwise, an index of the hypersphere
    in which the proposed location is contained.

    # Arguments

    - `location`: a proposed location in a unit hypercube
    - `hyperspheres`: a vector of `HyperSphere` objects
    """
    function p_fun(location, hyperspheres::Vector{HyperSphere}, args...; kwargs...)
        for (i,h) in enumerate(hyperspheres)
            in_contact(h, location) ? (return i) : nothing
        end
        # invalid area (not in a hypersphere)
        return -100
    end

"""
        OddShape

    An odd shape the is found by (1) initializing a shape at random partition and adding it to set of indices, (2)
    selecting random indices from the set, (3) adding an ajecent index to the set of indices, and repeating 2-3.

    # Fields

    - `indices`: a vector of vectors that represent the location of the odd shape
    - `p_bounds`: the partition boundaries for each dimension
    """
    mutable struct OddShape{T1,T2}
        indices::T1
        p_bounds::T2
        shape_ids::Vector{Int}
        n_shapes::Int
    end

    function OddShape(indices, p_bounds, n_cells::Vector{Int})
        shape_ids = vcat(fill.(1:length(n_cells), n_cells)...)
        n_shapes = length(n_cells)
        return OddShape(indices, p_bounds, shape_ids, n_shapes)
    end

    """
        set_locations(r_dist, n_dims, n_obj)

    Sets the location of `OddShape` within a unit hypercube

    # Arguments

    - `n_dims`: the number of dimensions of the hypersphere
    - `n_part`: number of partitions along a dimension
    - `n_cells`: the number of cells per shape
    """
    function set_indices(n_dims, n_part, n_cells)
        n_shapes = length(n_cells)
        indices = Vector{Vector{Int}}()
        for s in 1:n_shapes
            index = rand(1:n_part, n_dims)
            push!(indices, index)
            for i in 1:(n_cells[s]-1)
                while not_valid(indices, index, n_part)
                    index = propose(indices)
                end
                push!(indices, index)
            end
        end
        return indices
    end

    function not_valid(xs, x, n_part)
        return (x ∈ xs) || any(x .< 1) || any(x .> n_part)
    end

    function propose(indices)
        index = deepcopy(rand(indices))
        i = rand(1:length(index))
        index[i] += rand() ≤ .5 ? 1 : -1
        return index
    end

    """
        p_fun(location, shape::OddShape, args...; kwargs...)

    Returns 1 if location is in OddShape and 0 otherwise.

    # Arguments

    - `location`: a proposed location in a unit hypercube
    - `Oddshape`: an odd contiguous shape
    """
    function p_fun(location, shape::OddShape, args...; kwargs...)
        bounds = shape.p_bounds
        indices = shape.indices
        shape_ids = shape.shape_ids
        n_dims = length(location)
        n_b = length(bounds)
        index = fill(0, n_dims)

        for j in 1:n_dims
            for i in 1:(n_b-1)
                if (location[j] ≥ bounds[i]) && (location[j] ≤ bounds[i+1])
                    index[j] = i
                    continue
                end
            end
        end
        return get_id(index, indices, shape_ids)
    end

    function get_id(index, indices, shape_ids)
        for i in 1:length(indices)
            if indices[i] == index
                return shape_ids[i]
            end
        end
        return 0
    end

end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
export
    propagate,
    get_actions!,
    tick!,
    reset_hidden_states!,
    simulate!

include("driver_models.jl")
include("simulation.jl")
include("callbacks.jl")
include("actions.jl")
include("features.jl")


text:
@testset "BoundConstraints" begin

    @testset "Basic x F64, bc $T " for T in (Float64, Int)
        n=5
        bc = BoundConstraints(T,n)
        x = rand(n)

        @test lower_bound(bc)==zeros(T,n)
        @test upper_bound(bc)==ones(T,n)
        @test x ∈ bc
    end

    @testset "BC operations type=$T " for T in (Float64, Int)
        n=5
        bc = BoundConstraints(T,n)

        bc_copy = deepcopy(bc)

        bc = bc - ones(Int,n)

        @test lower_bound(bc)==-ones(T,n)
        @test upper_bound(bc)==zeros(T,n)
        @test eltype(bc)==T

        bc = bc + ones(T,n)

        @test lower_bound(bc)≈zeros(T,n)
        @test upper_bound(bc)≈ones(T,n)
    end

end


text:
include("./params.jl")
include("./fitgamma.jl")

function intensity_gamma(data)
    # Example usage:
    # x = 0.5-0.5*log(rand(1,1e3));
    # rs = intensity_gamma(x);
    #
    # Input argument
    # data: Spike-rate.
    #
    # Output argument
    # rs:   Estimated intensity function

    isi = data[2:end] - data[1:end-1]
    lambda, k = fitgamma(isi)

    Nspike = length(data)

    P = 0;
    prev_flg = 0;
    rs = zeros(L)
    for i=1:L

        if prev_flg == 0
            tau = i*dt
        else
            tau = i*dt - data[prev_flg]
        end

        p = 1.0 / gamma(k) * (k * lambda).^k * tau.^(k-1) .* exp(-k*lambda*tau); # ISI
        P = P + p*dt;

        rs[i] = p / (1.0 - P); # Intensity

        if prev_flg+1 < Nspike && data[prev_flg+1] <= i*dt

            P = 0;
            prev_flg += 1;

        end

    end

    rs
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
# nums = [[parse(Int64, l) for l in line] for line in readlines("inputs/04.txt")]
lines = readlines("inputs/04.txt")

nums = lines[1]
nums = [parse(Int64, x ) for x in split(nums, ',')]

lines = lines[2:end]

grids = []
for i in 2:6:length(lines)-4
  println(lines[i:i+5])
  grid = zeros(Int64, 5,5)
  gridstrs = lines[i:i+4]
  for (r_i, r) in enumerate(gridstrs)
    row_vec = [parse(Int64, x) for x in filter(q -> length(q) > 0 ,split(strip(r), ' '))]
    grid[r_i,:] = row_vec
  end
  push!(grids, grid)
end

function day3(nums, grids_)
  grids = copy(grids_)
  for num in nums
    grids = map(g -> replace(g, num => 100), grids)
    for g in grids
      if is_bingo(g)
        print(g)
        println(sum(g) - (100 * count(i -> i==100, g1)))
        println("last num: ", num)
        return g
      end
    end
  end
end
gx = day3(nums, grids)

function is_bingo(grid)
  if 500 in sum(grid, dims = 1) || 500 in sum(grid, dims = 2)
    return true
  end
  # if sum([grid[i,i] for i in 1:5]) == 500 || sum([grid[5,1], grid[4,2], grid[3,3], grid[2,4], grid[1,5]]) == 500
  #   return true
  # end
  return false
end


text:
using Test
using AtiyahBott

# first test
println( "Testing O1" )
P = (g,c,w,s,m) -> O1(g,c,w,s,m)^2
res = AtiyahBottFormula(2,2,6-1,P,true,true,false)
if res[1] == 1//1
    @test true
else
    @test false
end

#second test
println( "Testing Incidency and Contact" )
d = 4
P = Vector{Any}(undef, d+1)
for a in 0:d
    P[a+1] = (g,c,w,s,m) -> Incidency(g,c,w,s,3)^a*
        Incidency(g,c,w,s,2)^(2*(d-a)+1)*
        Contact(g,c,w,s)
end
res = AtiyahBottFormula(3,d,0,P,true,true,false);

if res == [ 1089024, 96512, 9408, 1024, 128 ]
    @test true
else
    @test false
end

#third test
println( "Testing Hypersurface" )
d = 4 #for other values of d, change this line
P = (g,c,w,s,m) -> Hypersurface(g,c,w,s,5)
res = AtiyahBottFormula(4,d,0,P,true,true,false);

if res == [ 15517926796875//64 ]
    @test true
else
    @test false
end

#fourth test
println( "Testing R1" )
d = 4 #for other values of d, change this line
P = (g,c,w,s,m) -> R1(g,c,w,s,1)^2
res = AtiyahBottFormula(1,d,0,P,true,true,false)

if res == [ 1//64 ]
    @test true
else
    @test false
end


text:
using InteractiveUtils

mutable struct Point
    x::Float64
    y::Float64
    z::Float64
end

typeof(Point) #> DataType

fieldnames(Point) #> (:x, :y, :z)

orig = Point(0, 0, 0) #> Point(0.0,0.0,0.0)
p1 = Point(2, 4, 1.3) #> Point(2.0,4.0,1.3)
fieldnames(typeof(p1)) #> (:x, :y, :z)

methods(Point)
# 3 methods for generic function "(::Type)":
# [1] Point(x::Float64, y::Float64, z::Float64) in Main at REPL[12]:2
# [2] Point(x, y, z) in Main at REPL[12]:2
# [3] (::Type{T})(arg) where T in Base at deprecated.jl:461

typeof(p1) #> Point (constructor with 1 method)
subtypes(Point) #> 0-element Array{Type,1}

p1.y #> 4.0
p1.z = 3.14 #> 3.14
p1 #> Point(2.0, 4.0, 3.14)
# p1.z = "A" #> ERROR: MethodError: Cannot `convert` an object of type String
# to an object of type Float64

struct Vector3D
    x::Float64
    y::Float64
    z::Float64
end
typeof(Vector3D) #> DataType
p = Vector3D(1, 2, 3) #> Vector3D(1.0,2.0,3.0)
# p.y = 5 #> ERROR: type Vector3D is immutable

Point3D = Point
p31 = Point3D(1, 2, 3) #> Point(1.0,2.0,3.0)

5 == 5 #> true
5 == 5.0 #> true
isequal(5, 5) #> true
isequal(5, 5.0) #> true
5 === 5   #> true
5 === 5.0 #> false

q = Vector3D(4.0, 3.14, 2.71)
r = Vector3D(4.0, 3.14, 2.71)
isequal(q, r) #> true
q === r #> true

mutable struct MVector3D
    x::Float64
    y::Float64
    z::Float64
end

q = MVector3D(4.0, 3.14, 2.71)
r = MVector3D(4.0, 3.14, 2.71)
isequal(q, r) #> false
q === r #> false

# constructors, subtyping and multiple dispatch
abstract type Employee
end

# Employee() #> ERROR: MethodError: no constructors have been defined for Employee

mutable struct Developer <: Employee
    name::String
    iq
    favorite_lang
end

# mutable struct MobileDeveloper <: Developer
#    platform
# end
# ERROR: invalid subtyping in definition of MobileDeveloper

# outer constructor:
Developer(name, iq) = Developer(name, iq, "Java")
#> Developer (constructor with 3 methods)

mutable struct Manager
    name::String
    iq
    department::String
end

devel1 = Developer("Bob", 110) #> Developer("Bob",110,"Java")
devel2 = Developer("William", 145, "Julia")
#> Developer("William",145,"Julia")
man1 = Manager("Julia", 120, "ICT") #> Manager("Julia",120,"ICT")

cleverness(emp::Employee) = emp.iq
cleverness(devel1) #> 110
# cleverness(man1) #> ERROR: MethodError: `cleverness` has no method
# matching cleverness(::Manager)

function cleverer(m::Manager, e::Employee)
    println("The manager $(m.name) is cleverer!")
end

cleverer(man1, devel1) #> The manager Julia is cleverer!
cleverer(man1, devel2) #> The manager Julia is cleverer!

function cleverer(d::Developer, e::Employee)
    println("The developer $(d.name) is cleverer!")
end

cleverer(devel1, devel2) #> The developer Bob is cleverer!
# cleverer(devel1, man1) #> ERROR: `cleverer` has no method matching cleverer(::Developer,::Manager)

function cleverer(e::Employee, d::Developer)
    if e.iq <= d.iq
        println("The developer $(d.name) is cleverer!")
    else
        println("The employee $(e.name) is cleverer!")
    end
end

# cleverer(devel1, devel2)
#> ERROR: MethodError: cleverer(::Developer, ::Developer) is ambiguous. Candidates:
#    cleverer(e::Employee, d::Developer) in Main at REPL[32]:2
#    cleverer(d::Developer, e::Employee) in Main at REPL[29]:2
#  Possible fix, define
#    cleverer(::Developer, ::Developer)

function cleverer(d1::Developer, d2::Developer)
    if d1.iq <= d2.iq
        println("The developer $(d2.name) is cleverer!")
    else
        println("The developer $(d1.name) is cleverer!")
    end
end

cleverer(devel1, devel2) #> The developer William is cleverer!
cleverer(devel2, devel1) #> The developer William is cleverer!


text:
# # Eady model of Baroclinic Instability
#
# A simulation of the growth of barolinic instability in the Eady n-layer model
# when we impose a vertical mean flow shear as a difference $\Delta U$ in the
# imposed, domain-averaged, zonal flow at each layer.

using FourierFlows, Plots, Printf

using FFTW: rfft, irfft
using Statistics: mean

import GeophysicalFlows.MultiLayerQG
import GeophysicalFlows.MultiLayerQG: energies

# ## Numerical parameters and time-stepping parameters
nx = 256                  # 2D resolution = nx^2
ny = nx

stepper = "FilteredRK4"   # timestepper
dt = 2e-3      # timestep
nsteps = 20000            # total number of time-steps
nsubs  = 100              # number of time-steps for plotting (nsteps must be multiple of nsubs)

# ## Physical parameters
Lx = 2π         # domain size
 μ = 1e-2       # bottom drag
 β = 0.1          # the y-gradient of planetary PV
 g = 1.0
f₀ = 2.0

# Vertical grid
nlayers = 5     # number of layers

total_depth = 1
Δh = total_depth / nlayers
z = [ (i - 1/2) * Δh for i = nlayers:-1:1 ] .- total_depth
H = [ Δh for i = 1:nlayers ]

# Density stratification
ρ₀ = 1
N² = 5

ρ_resting(z) = ρ₀ * (1 - N² * z / g)
ρ = ρ_resting.(z)

@info @sprintf("The largest Rossby radius of deformation is %.2f",
               sqrt(N²) * total_depth / f₀)

@info @sprintf("The smallest Rossby radius of deformation is %.2f",
               sqrt(N²) * Δh / f₀)

# Background shear
ΔU = 1 / (nlayers - 1)
U = [i * ΔU for i=nlayers-1:-1:0]

twod_grid = TwoDGrid(nx, Lx)

X, Y = gridpoints(twod_grid)
σx, σy = 0.7, 0.5
eta = f₀ / sum(H) * Δh * @. exp(-X^2/σx^2 -Y^2/σy^2)

# ## Problem setup
#
# We initialize a `Problem` by providing a set of keyword arguments,

prob = MultiLayerQG.Problem(nlayers,
                             nx = nx,
                             Lx = Lx,
                             f₀ = f₀,
                              g = g,
                              H = H,
                              ρ = ρ,
                              U = U,
                             dt = dt,
                        stepper = stepper,
                            eta = eta,
                              μ = μ,
                              β = β
                           )

# and define some shortcuts.
sol, clock, params, vars, grid = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
x, y = grid.x, grid.y

# ## Setting initial conditions

# Our initial condition is some small amplitude random noise. We smooth our initial
# condidtion using the `timestepper`'s high-wavenumber `filter`.

qᵢ  = randn((nx, ny, nlayers))
for j in 1:nlayers
  qᵢ[:, :, j] = qᵢ[:, :, j] .- mean(qᵢ[:, :, j])
end
qhᵢ = prob.timestepper.filter .* rfft(qᵢ, (1, 2)) # only apply rfft in dims=1, 2
qᵢ  = irfft(qhᵢ, grid.nx, (1, 2)) # only apply irfft in dims=1, 2

MultiLayerQG.set_q!(prob, qᵢ)

KEᵢ, PEᵢ = MultiLayerQG.energies(prob)
KEᵢ

ψhᵢ = similar(qhᵢ)

MultiLayerQG.streamfunctionfrompv!(ψhᵢ, qhᵢ, prob.params, prob.grid)

for j in 1:nlayers
  ψhᵢ[:, :, j] = sqrt(5e-8 / KEᵢ[j]) * ψhᵢ[:, :, j]
end
ψᵢ  = irfft(ψhᵢ, grid.nx, (1, 2))

MultiLayerQG.set_ψ!(prob, ψᵢ)
nothing # hide

# ## Diagnostics

# Create Diagnostics -- `energies` function is imported at the top.
E = Diagnostic(energies, prob; nsteps=nsteps)
diags = [E] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.

nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_eady"
plotname = "snapshots"
filename = joinpath(filepath, "eady.jld2")
nothing # hide

# Do some basic file management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# And then create Output
get_sol(prob) = prob.sol # extracts the Fourier-transformed solution

out = Output(prob, filename, (:sol, get_sol))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the potential vorticity field and the evolution
# of energy and enstrophy.

function plot_output(prob)
  q = prob.vars.q
  x, y = prob.grid.x, prob.grid.y
  Lx, Ly = prob.grid.Lx, prob.grid.Ly

  layout = @layout Plots.grid(2, Int(round(nlayers/2)+1)+1)

  p = plot(layout=layout, size=(1000, 600), dpi=180)

  for m in 1:nlayers

    heatmap!(p[m], x, y, q[:, :, m],
         aspectratio = 1,
              legend = :none,
                   c = :deep,
               xlims = (-Lx/2, Lx/2),
               ylims = (-Ly/2, Ly/2),
               title = "q layer "*string(m),
              xticks = :none,
              yticks = :none,
            colorbar = false,
          framestyle = :box

         )
  end

  plot!(p[nlayers+1], nlayers,
             # label = ["KE1" "KE2"],
            legend = false,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.05*μ*nsteps*dt, 1.01*μ*nsteps*dt),
             ylims = (1e-8, 1e2),
            yscale = :log10,
            yticks = 10.0.^(-8:2:2),
            xlabel = "μt",
             title = "KE at layers")

  plot!(p[nlayers+2], nlayers-1,
             # label = "PE",
            legend = false,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.05*μ*nsteps*dt, 1.01*μ*nsteps*dt),
             ylims = (1e-9, 1e1),
            yscale = :log10,
            yticks = 10.0.^(-8:2:2),
            xlabel = "μt",
             title = "PE at interfaces")

  return p
end

nothing # hide

# ## Time-stepping the `Problem` forward

# Finally, we time-step the `Problem` forward in time.

p = plot_output(prob)

startwalltime = time()

anim = @animate for j = 0:Int(nsteps / nsubs)

    cfl = clock.dt * maximum([maximum(abs, vars.u) / grid.dx, maximum(abs, vars.v) / grid.dy])

    log = @sprintf("step: %04d, dt: %.3f, t: %0.2f, cfl: %.2f, walltime: %.2f min",
                    clock.step, clock.dt, clock.t, cfl, (time()-startwalltime)/60)

    log_KE = @sprintf("KE₁: %.3e     KE₂: %.3e     KE₃: %.3e     KE₄: %.3e     KE₅: %.3e",
                  E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][1][3], E.data[E.i][1][4], E.data[E.i][1][5])

    log_PE = @sprintf("      PE₁.₅: %.3e   PE₂.₅: %.3e   PE₃.₅: %.3e   PE₄.₅: %.3e",
                  E.data[E.i][2][1], E.data[E.i][2][2], E.data[E.i][2][3], E.data[E.i][2][4])

    # log_KE = @sprintf("KE₁: %.3e     KE₂: %.3e     KE₃: %.3e",
    #               E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][1][3])
    #
    # log_PE = @sprintf("      PE₁.₅: %.3e   PE₂.₅: %.3e",
    #               E.data[E.i][2][1], E.data[E.i][2][2])
    #

    if j % (500 / nsubs) == 0
      println(log)
      println(log_KE)
      println(log_PE)
    end

    for m in 1:nlayers
      p[m][1][:z] = @. vars.q[:, :, m] + params.β * grid.y'
      if m==nlayers
        p[m][1][:z] = @. vars.q[:, :, m] + params.β * grid.y' + params.eta
      end
      push!(p[nlayers+1][m], μ*E.t[E.i], E.data[E.i][1][m])
    end

    for m in 1:nlayers-1
      push!(p[nlayers+2][m], μ*E.t[E.i], E.data[E.i][2][m])
    end

    stepforward!(prob, diags, nsubs)

    MultiLayerQG.updatevars!(prob)

end

mp4(anim, "multilayerqg_eady.mp4", fps=18)

# ## Save

#=
# Finally save the last snapshot.
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
savefig(savename)
=#

max_stars_repo_name: navidcy/MultilayerQG-example

max_issues_repo_name: navidcy/MultilayerQG-example

max_forks_repo_name: navidcy/MultilayerQG-example


text:
@testset "kinetic energy zero at zero velocity" begin
    θ_t = deg2rad(30.0)
    θ′_t = 0.0
    θ_k = deg2rad(22.0)
    θ′_k = 0.0
    γ_OA = deg2rad(60.0)
    γ′_OA = 0.0
    T = @test_nowarn kinetic_energy(0.0, 0.0, 0.0, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA)
    @test T == 0.0
end

@testset "kinetic energy non-negative on various inputs" begin
    θ_t = deg2rad(30.0)
    θ′_t = 0.0
    θ_k = deg2rad(22.0)
    θ′_k = 0.1
    γ_OA = deg2rad(60.0)
    γ′_OA = -0.2
    @test kinetic_energy(-0.1, 0.0, 0.0, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(0.1, 0.0, 0.0, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(0.0, -0.2, 0.0, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(0.1, 0.12, 0.0, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(0.0, 0.0, -0.22, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(0.0, 0.4, 0.22, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(-1000.0, 100.4, 10.22, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    @test kinetic_energy(-1000.0, -100.4, -1000.22, θ_t, θ′_t, θ_k, θ′_k, γ_OA, γ′_OA) ≥ 0
    x = randn(9, 12) .* [ones(3) .* 1e3; ones(6) .* 0.1]
    for p in eachcol(x)
        @test kinetic_energy(p...) ≥ 0
    end
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Basics (Modia3D/Basics/_module.jl)
#

"""
    path = PTP_path(names;
                    positions = [zeros(size(names,1))'; ones(size(names,1))'],
                    startTime=0.0,
                    v_max = ones(size(names,1)),
                    a_max = ones(size(names,2))))

Generate a new path object to move as fast as possible from
positions[i,:] to positions[i+1,:]. The `positions[i,:]` can be a set of translational
positions in [m] (= absolute distances) and/or rotational positions in [rad]
(= angles). In robotics such a movement is called PTP (Point-To-Point).
The signals are constructed in such a way that it is not possible
to move faster, given the maximally allowed velocity `v_max[j]` and
the maximally allowed acceleration `a_max[j]` for signal `names[j]`
and have a velocity of zero at the given `positions`.

If there are two or more signals (that is length(names) > 1) the
the path is constructed such that all signals
are in the same periods in the acceleration, constant velocity
and deceleration phase. This means that only one of the signals
is at its limits whereas the others are synchronized in such a way
that the end point is reached at the same time instant.

For example, this means that the signals have a velocity of zero at positions[1,:],
one of the signals is accelerated with its maximally allowed acceleration until
one of the signals reaches its maximally allowed velocity. At a proper time instant,
one of the signals is decelerated with the negative value of its maximally allowed
acceleration, so that all signals reach positions[2,:] with velocity zero.

This element is useful to generate a reference signal for a controller
which controls, e.g., a drive train, or to drive
a flange according to a given acceleration.

# Example
```julia
using Modia3D
using ModiaLang
@usingModiaPlot

const ptp_path = PTP_path(["angle1", "angle2", "angle3"],
                          positions = [0.0 2.0 3.0;  # angle1=0.0, angle2=2.0, angle3=3.0
                                       0.5 3.0 4.0;
                                       0.8 1.5 0.3;
                                       0.2 1.5 0.8],
                          startTime = 0.1,
                          v_max = 2*ones(3),
                          a_max = 3*ones(3))
angles = zeros(3)
getPosition!(ptp_path, 0.5, angles)   # angles = [0.12, 2.24, 3.24]
plotPath(ptp_path, plot)   # used plot(..) defined with @usingModiaPlot
```
"""
mutable struct PTP_path
   names::Vector{AbstractString}
   startTime::Float64
   v_max::Vector{Float64}
   a_max::Vector{Float64}
   positions::Matrix{Float64}

   delta::Matrix{Float64}
   hasPath::Vector{Bool}
   sd_max::Vector{Float64}
   sdd_max::Vector{Float64}
   Ta1::Vector{Float64}
   Ta2::Vector{Float64}
   noWphase::Vector{Bool}
   Tv::Vector{Float64}
   Te::Vector{Float64}
   Ta1s::Vector{Float64}
   Ta2s::Vector{Float64}
   Tvs::Vector{Float64}
   Tes::Vector{Float64}
   sd_max2::Vector{Float64}
   s1::Vector{Float64}
   s2::Vector{Float64}
   s3::Vector{Float64}
   Tend::Float64
   posTemp::Vector{Float64}   # Temporary storage that can be used by the functions operation on PTP_path

   function PTP_path(names::AbstractVector;
                     positions::Matrix{Float64} = [zeros(size(names,1))'; ones(size(names,1))'],
                     startTime=0.0,
                     v_max = ones(size(names,1)),
                     a_max = ones(size(names,1)))

        #@assert(size(positions,1) > 1)
        #@assert(size(positions,2) == Base.length(names))
        #@assert(Base.length(v_max) == Base.length(names))
        #@assert(Base.length(a_max) == Base.length(names))
        np    = Base.length(names)
        npath = size(positions,1) - 1   # number of path points

        for i in eachindex(v_max)
            @assert(v_max[i] > 0.0)
            @assert(a_max[i] > 0.0)
        end

        delta    = zeros(npath, np)
        hasPath  = fill(false,npath)
        sd_max   = zeros(npath)
        sdd_max  = zeros(npath)
        Ta1      = zeros(npath)
        Ta2      = zeros(npath)
        noWphase = fill(false,npath)
        Tv       = zeros(npath)
        Te       = zeros(npath)
        Ta1s     = zeros(npath)
        Ta2s     = zeros(npath)
        Tvs      = zeros(npath)
        Tes      = zeros(npath)
        sd_max2  = zeros(npath)
        s1       = zeros(npath)
        s2       = zeros(npath)
        s3       = zeros(npath)
        aux1     = zeros(np)
        aux2     = zeros(np)
        small    = 1000*eps()

        for i in 1:npath
            delta[i,:] = positions[i+1,:] - positions[i,:]

            for j in 1:np
                aux1[j] = delta[i,j]/v_max[j]
                aux2[j] = delta[i,j]/a_max[j]
            end

            sd_max_inv  = maximum(abs.(aux1))
            sdd_max_inv = maximum(abs.(aux2))

            if sd_max_inv <= small || sdd_max_inv <= small
                hasPath[i] = false
            else
                hasPath[i] = true
                sd_max[i]  = 1/sd_max_inv
                sdd_max[i] = 1/sdd_max_inv

                Ta1[i] = sqrt(1/sdd_max[i])
                Ta2[i] = sd_max[i]/sdd_max[i]
                noWphase[i] = Ta2[i] >= Ta1[i]

                Tv[i]   = noWphase[i] ? Ta1[i] : 1/sd_max[i]
                Te[i]   = noWphase[i] ? Ta1[i] + Ta1[i] : Tv[i] + Ta2[i]

                Tbegin  = i==1 ? startTime : Tes[i-1]
                Ta1s[i] = Ta1[i] + Tbegin
                Ta2s[i] = Ta2[i] + Tbegin
                Tvs[i]  = Tv[i]  + Tbegin
                Tes[i]  = Te[i]  + Tbegin
                sd_max2[i] = sdd_max[i]*Ta1[i]
                s1[i] = sdd_max[i]*(noWphase[i] ? Ta1[i]*Ta1[i] : Ta2[i]*Ta2[i])/2
                s2[i] = s1[i] + (noWphase[i] ? sd_max2[i]*(Te[i] - Ta1[i]) - (sdd_max[i]/2)*(Te[i] - Ta1[i])^2 : sd_max[i]*(Tv[i] - Ta2[i]))
                s3[i] = s2[i] + sd_max[i]*(Te[i] - Tv[i]) - (sdd_max[i]/2)*(Te[i] - Tv[i])^2
            end
        end

        Tend = Tes[end]

        new(names, startTime, v_max, a_max, positions, delta,
            hasPath, sd_max, sdd_max, Ta1, Ta2, noWphase, Tv, Te, Ta1s, Ta2s, Tvs, Tes, sd_max2, s1, s2, s3, Tend, zeros(np))
    end
end

"""
    Tend = pathEndTime(path)

Given a `path::PTP_path` return the end time `Tend` of the path.
"""
pathEndTime(path::PTP_path) = path.Tend

"""
    getPosition!(path, time, position)

Given a `path::PTP_path` and a time instant `time`, return the actual
position at time `time` in vector `position`.
"""
function getPosition!(path::PTP_path, time::Number, position::Vector{Float64})
    time = ustrip(time)
    @assert(length(position) == size(path.positions,2))
    npath = length(path.hasPath)
    np    = length(position)

    # Search correct time interval
    i = 0
    if time <= path.startTime
        i = 1
        s = 0

    else
        while i < npath
            i = i+1
            if time <= path.Tes[i]
               break
            end
        end
        if time >= path.Tes[end]
            i = npath
            s = path.noWphase[i] ? path.s2[end] : path.s3[end]
            #println("... time=$time i=$i s=$s qbegin=", path.positions[i,1], ", qdelta = ", path.delta[i,1])
        else
            Tbegin = i==1 ? path.startTime : path.Tes[i-1]
            if path.noWphase[i]
                if time < path.Ta1s[i]
                    s = (path.sdd_max[i]/2)*(time - Tbegin)^2
                elseif time < path.Tes[i]
                    s = path.s1[i] + path.sd_max2[i]*(time - path.Ta1s[i]) - (path.sdd_max[i]/2)*(time - path.Ta1s[i])^2
                else
                    s = path.s2[i]
                end
            elseif time < path.Ta2s[i]
                s = (path.sdd_max[i]/2)*(time - Tbegin)^2
            elseif time < path.Tvs[i]
                s = path.s1[i] + path.sd_max[i]*(time - path.Ta2s[i])
            elseif time < path.Tes[i]
                s = path.s2[i] + path.sd_max[i]*(time - path.Tvs[i]) - (path.sdd_max[i]/2)*(time - path.Tvs[i])*(time - path.Tvs[i])
            else
                s = path.s3[i]
            end
        end
    end

    for j in 1:np
        # println("... i = $i, j= $j, np = $np, s=$s, positions[i,1]=",path.positions[i,1])
        position[j]= path.positions[i,j] + path.delta[i,j]*s
    end
end

"""
    getPosition!(path, time, position, velocity, acceleration)

Given a `path::PTP_path` and a time instant `time`, return the actual
position, velocity and acceleration at time `time` in vectors
`position, velocity, acceleration`.
"""
function getPosition!(path::PTP_path, time::Number,
                      position::Vector{Float64},
                      velocity::Vector{Float64},
                      acceleration::Vector{Float64})::Nothing
    time = ustrip(time)
    @assert(length(position)     == size(path.positions,2))
    @assert(length(velocity)     == length(position))
    @assert(length(acceleration) == length(acceleration))
    npath = length(path.hasPath)
    np    = length(position)
    s   = 0.0
    sd  = 0.0
    sdd = 0.0

    # Search correct time interval
    i = 0
    if time <= path.startTime
        i   = 1
        s   = 0
        sd  = 0.0
        sdd = 0.0
    else
        while i < npath
            i = i+1
            if time <= path.Tes[i]
               break
            end
        end
        if time >= path.Tes[end]
            i   = npath
            s   = path.noWphase[i] ? path.s2[end] : path.s3[end]
            sd  = 0.0
            sdd = 0.0
            #println("... time=$time i=$i s=$s qbegin=", path.positions[i,1], ", qdelta = ", path.delta[i,1])
        else
            Tbegin = i==1 ? path.startTime : path.Tes[i-1]
            if path.noWphase[i]
                if time < path.Ta1s[i]
                    s   = (path.sdd_max[i]/2)*(time - Tbegin)^2
                    sd  = path.sdd_max[i]*(time - Tbegin)
                    sdd = path.sdd_max[i]
                elseif time < path.Tes[i]
                    s   = path.s1[i] + path.sd_max2[i]*(time - path.Ta1s[i]) - (path.sdd_max[i]/2)*(time - path.Ta1s[i])^2
                    sd  = path.sd_max2[i] - path.sdd_max[i]*(time - path.Ta1s[i])
                    sdd = -path.sdd_max[i]
                else
                    s   = path.s2[i]
                    sd  = 0.0
                    sdd = 0.0
                end
            elseif time < path.Ta2s[i]
                s   = (path.sdd_max[i]/2)*(time - Tbegin)^2
                sd  = path.sdd_max[i]*(time - Tbegin)
                sdd = path.sdd_max[i]
            elseif time < path.Tvs[i]
                s   = path.s1[i] + path.sd_max[i]*(time - path.Ta2s[i])
                sd  = path.sd_max[i]
                sdd = 0.0
            elseif time < path.Tes[i]
                s   = path.s2[i] + path.sd_max[i]*(time - path.Tvs[i]) - (path.sdd_max[i]/2)*(time - path.Tvs[i])^2
                sd  = path.sd_max[i] - path.sdd_max[i]*(time - path.Tvs[i])
                sdd = -path.sdd_max[i]
            else
                s   = path.s3[i]
                sd  = 0.0
                sdd = 0.0
            end
        end
    end

    for j in 1:np
        # println("... i = $i, j= $j, np = $np, s=$s, positions[i,1]=",path.positions[i,1])
        position[j]     = path.positions[i,j] + path.delta[i,j]*s
        velocity[j]     = path.delta[i,j]*sd
        acceleration[j] = path.delta[i,j]*sdd
    end

    return nothing
end

"""
    pos = getPosition(path, index, time)

Given a `path::PTP_path`, the `index` of a signal, and a time instant `time`, return the actual
position at time `time`.
"""
function getPosition(path::PTP_path, index, time::Number)
    getPosition!(path, time, path.posTemp)
    return path.posTemp[index]
end

"""
    index = getIndex(path, name)

Return the index of `name` in `path` or trigger an error, if not present.
"""
function getIndex(path, name)
    index = findfirst(x -> x==name, path.names)
    if isnothing(index)
        error("getIndex(path,name): \"", name, "\" not in path")
    end
    return index
end

"""
    plotPath(path, plot::Function;
             names=path.names, heading="PTP plots",
             tend=1.1*path.Tend, figure=1, ntime=101, onlyPositions=true)

Given a `path::PTP_path`, plot the path over `time` up to `tend` for all points
identified by the vector or tuple `names` to figure `figure`
using `ntime` time points.

# Example

```julia
using Modia3D
using ModiaLang
@usingModiaPlot

const ptp_path = PTP_path(["angle1", "angle2", "angle3"],
                          positions = [0.0 2.0 3.0;  # angle1=0.0, angle2=2.0, angle3=3.0
                                       0.5 3.0 4.0;
                                       0.8 1.5 0.3;
                                       0.2 1.5 0.8],
                          startTime = 0.1,
                          v_max = 2*ones(3),
                          a_max = 3*ones(3))
angles = zeros(3)
getPosition!(ptp_path, 0.5, angles)   # angles = [0.12, 2.24, 3.24]
plotPath(ptp_path, plot)   # used plot(..) defined with @usingModiaPlot
```
"""
function plotPath(path::PTP_path, plot::Function; names=path.names, heading="PTP plots", figure=1,
                  ntime=101, tend = 1.1*path.Tend, onlyPositions=true)::Nothing
    time = range(0u"s",(tend)u"s",length=ntime)
    indices = indexin(names, path.names)
    names2  = deepcopy(names)
    for i in eachindex(indices)
        if isnothing(i)
            @warn "plotPath(path, ...): \""*names[i]*"\" is ignored, because not in path"
            deleteat!(indices,i)
            deleteat!(names2,i)
        end
    end

    np   = length(indices)
    q    = zeros(length(time), np)
    qt   = zeros(length(path.names))

    series = Dict{AbstractString,Any}()
    series["time"] = time

    if onlyPositions
        for i in eachindex(time)
            getPosition!(path, time[i], qt)
            q[i,:] = qt[indices]
        end

        for i in eachindex(names2)
            series[names2[i]] = q[:,i]
        end

        plot(series, Tuple(names2), heading=heading, figure=figure)
    else
        der_names2  = "der(" .* names2 .* ")"
        der2_names2 = "der2(" .* names2 .* ")"
        qd   = zeros(length(time), np)
        qdd  = zeros(length(time), np)
        qtd  = zeros(length(path.names))
        qtdd = zeros(length(path.names))
        for i in eachindex(time)
            getPosition!(path, time[i], qt, qtd, qtdd)
            q[i,:]   = qt[indices]
            qd[i,:]  = qtd[indices]
            qdd[i,:] = qtdd[indices]
        end

        for i in eachindex(names2)
            series[names2[i]]      = q[:,i]
            series[der_names2[i]]  = qd[:,i]
            series[der2_names2[i]] = qdd[:,i]
        end

        plot(series, [Tuple(names2), Tuple(der_names2), Tuple(der2_names2)],
             heading=heading, figure=figure)
    end

    return nothing
end

"""
    component = Modia3D.getPathComponent(rootComponent, path::String)

Return the component defined by rootComponent and path.

"""
function getPathComponent(rootComponent, path::String)
    component = rootComponent
    lenPath = length(path)
    if lenPath == 0
        @goto PathInvalid
    end
    i = 1
    while true
        j = findnext(".", path, i)
        key = isnothing(j) ? Symbol(path[i:end]) : Symbol(path[i:j[1]-1])
        if !isdefined(component, key) # hasfield(typeof(component), key)
            @goto PathInvalid
        end
        component = getfield(component, key)
        if isnothing(j)
            return component
        end
        i=j[1]+1
        if i > lenPath
            @goto PathInvalid
        end
    end

    @label PathInvalid
    T = typeof(rootComponent)
    error("getPathComponent(rootComponent::$T, \"$path\"): path is wrong.\n",
          "Allowed symbols on the last correct level:\n",
          fieldnames(typeof(component)))
end


text:
# -------------------------------------------------------------------
# JULIA TOOL TO CALCULATE MAJORANA DENSITIES OF STATES
# FOR DIFFERENT 2D AND 3D KITAEV MODELS
# (written for Julia 0.6)
# Author: Tim Eschmann
# Date: 2019-11-11
# Modified for new version of LatticePhysics: 2019-01-31
# -------------------------------------------------------------------

using PyCall

# Get unit cell size
function get_uc_length(lat)
    if lat == "6a"
        N = 2
    elseif lat == "8a"
        N = 6
    elseif lat == "8b"
        N = 6
    elseif lat == "8n"
        N = 16
    elseif lat == "10a"
        N = 4
    elseif lat == "10b"
        N = 4
    elseif lat == "10c"
        N = 12
    elseif lat == "10d"
        N = 8
    end
    return N
end

# Generate Bloch matrix for lattice 'lat' with couplings 'J' at k-point 'k'
function get_ham(k,J,lat)
    if lat == "6a"
        f(k,J) = 2*(J[1]*exp(2*im*k[1]*pi) + J[2]*exp(2*im*k[2]*pi) + J[3])
        ham = [0 im*f(k,J); -im*conj(f(k,J)) 0]
    elseif lat == "8a"
        A1(k,J) = exp(2*im*k[1]*pi)*J[2]
        A2(k,J) = exp(-2*im*k[2]*pi)*J[3]
        A3(k,J) = exp(-2*im*k[3]*pi)*J[1]
        B = [im*J[2] im*A2(k,J) 0; im*J[3] 0 im*A1(k,J); im*J[1] im*J[2] im*J[3]]
        C1 = [0 0 0; 0 0 im*A3(k,J); 0 -im*conj(A3(k,J)) 0]
        C2 = [0 -im*A3(k,J) 0; im*conj(A3(k,J)) 0 0; 0 0 0]
        ham = [C2 B; transpose(conj(B)) C1]
    elseif lat == "8b"
        a13(k,J) = exp(-2*im*(k[1]+k[3])*pi)*J[2]
        a2(k,J) = exp(2*im*k[2]*pi)*J[1]
        a3(k,J) = exp(-2*im*k[3]*pi)*J[1]
        B2 = [im*J[3] 0 -im*a13(k,J); -im*J[2] im*J[3] 0; -im*a2(k,J) -im*J[2] im*J[3]]
        C1 = [0 -im*a3(k,J) 0; im*conj(a3(k,J)) 0 0 ; 0 0 0]
        C2 = [0 0 0; 0 0 im*a3(k,J); 0 -im*conj(a3(k,J)) 0]
        ham = [C1 B2; transpose(conj(B2)) C2]
    elseif lat == "8n"
        b1(k,J) = exp(-2*im*k[1]*pi)
        b13(k,J) = exp(-2*im*(k[1]-k[3])*pi)
        b2(k,J) = exp(2*im*k[2]*pi)
        b23(k,J) = exp(2*im*(k[2]-k[3])*pi)
        b3(k,J) = exp(2*im*k[3]*pi)
        D = [im*J[3] im*J[1] 0 0 0 0 -im*J[2] 0; 0 -im*J[2] im*J[3] 0 0 0 0 -im*J[1]; 0 0 -im*J[1] im*J[2] 0 0 -im*b3(k,J)*J[3] 0; 0 0 0 -im*J[3] im*J[1] 0 0 -im*b2(k,J)*J[2]; 0 0 0 0 -im*J[2] im*J[3] im*b2(k,J)*J[1] 0; im*J[2] 0 0 0 0 -im*J[1] 0 im*b13(k,J)*J[3]; 0 im*b23(k,J)*J[3] 0 -im*b1(k,J)*J[1] 0 im*J[2] 0 0; -im*J[1] 0 -im*b1(k,J)*J[2] 0 -im*conj(b3(k,J))*J[3] 0 0 0]
        E = zeros(8,8)
        ham = [E D; transpose(conj(D)) E]
    elseif lat == "10a"
        c1(k,J) = exp(-2*im*k[1]*pi)*J[3]
        c2(k,J) = exp(-2*im*k[2]*pi)*J[1]
        c3(k,J) = exp(-2*im*k[3]*pi)*J[1]
        ham = [0 -im*c2(k,J) -im*J[2] -im*c1(k,J); im*conj(c2(k,J)) 0 -im*J[3] im*J[2]; im*J[2] im*J[3] 0 -im*c3(k,J); im*conj(c1(k,J)) -im*J[2] im*conj(c3(k,J)) 0]
    elseif lat == "10b"
        A13(k,J) = exp(-2*im*k[3]*pi)*(J[1] + exp(2*im*k[1]*pi) * J[2])
        A4(k,J) = J[1] + exp(2*im*k[2]*pi) * J[2]
        B = [0 0; 0 0]
        C = [im*J[3] im*A13(k,J); im*A4(k,J) im*J[3]]
        ham = [B C; transpose(conj(C)) B]
    elseif lat == "10c"
        d1(k,J) = J[1] + exp(-2*im*k[1]*pi)*J[2]
        d2(k,J) = exp(-2*im*k[3]*pi)*J[3]
        d3(k,J) = exp(-2*im*k[2]*pi)*J[2]
        d4(k,J) = exp(2*im*k[1]*pi)*J[2]
        d5(k,J) = -J[1] + exp(-2*im*k[1]*pi)*J[2]
        d6(k,J) = exp(2*im*(k[1]+k[2])*pi)*J[2]
        A = [im*d1(k,J) 0 -im*d2(k,J) 0 0 0; im*J[3] -im*J[1] 0 0 im*d3(k,J) 0; 0 im*J[3] -im*J[1] 0 0 im*d4(k,J); 0 0 0 im*d5(k,J) 0 -im*d2(k,J); 0 im*J[2] 0 im*J[3] im*J[1] 0; 0 0 im*d6(k,J) 0 im*J[3] -im*J[1]]
        ham = [zeros(6,6) A; transpose(conj(A)) zeros(6,6)]
    elseif lat == "10d"
        A_ = im*[0 J[3]*exp(2*im*k[2]*pi) J[1] J[2]*exp(-2*im*k[3]*pi); J[3]*exp(-2*im*k[1]*pi) 0 J[2] J[1]; J[1] J[2] 0 J[3]; J[2]*exp(-2*im*k[3]*pi) J[1] J[3]*exp(2*im*(k[1]-k[2])*pi) 0]
        ham = [zeros(4,4) A_; transpose(conj(A_)) zeros(4,4)]
    end
    return 2*ham
end

function main()
    # Define lattice, couplings and number of k-points
    lats = ["8a"]
    J = [1/3., 1/3., 1/3.]
    ni = 400
    nj = 400
    nk = 400

    a1 = [1,0,0]
    a2 = [0,1,0]
    a3 = [0,0,1]

    # 2D:
    #bzone = [i*a1 + j*a2 for i=0:1/ni:1 for j=0:1/nj:1]

    # 3D:
    bzone = [i*a1 + j*a2 + k*a3 for i=0:1/ni:1 for j=0:1/nj:1 for k=0:1/nk:1]
    #println(bzone)

    for lat in lats
        N = get_uc_length(lat)
        #evals = Array{Float64,2}((ni+1)*(nj+1),N) #2D
        evals = Array{Float64,2}((ni+1)*(nj+1)*(nk+1),N) #3D
        ct = 1
        for k in bzone
            evals[ct,:] = eigvals(get_ham(k,J,lat))
            ct += 1
            if ct % 10000 == 0
                println("$(lat):  $(ct) / $(size(bzone,1))\n")
            end
        end

        #println(minimum(evals), " ", maximum(evals))

        open("dos$(lat).txt", "w") do f
            for i in 1:size(evals,1)
                for j in 1:size(evals,2)
                    write(f, "$(evals[i,j])\n")
                end
            end
        end
    end

end

main()

max_stars_repo_name: timeschmann/KitaevDOS

max_issues_repo_name: timeschmann/KitaevDOS

max_forks_repo_name: timeschmann/KitaevDOS